lcsim/src/org/lcsim/recon/muon
diff -u -r1.7 -r1.8
--- BarrelCalSegmentFinder.java 31 Oct 2005 22:31:14 -0000 1.7
+++ BarrelCalSegmentFinder.java 6 Dec 2005 16:11:21 -0000 1.8
@@ -1,15 +1,10 @@
package org.lcsim.recon.muon;
-//import hep.analysis.*; // C.M.-27Jan03- Folder
-//import org.lcsim.geometry.*;
-//import org.lcsim.geometry.component.*;
-//import org.lcsim.recon.tracking.*;
import org.lcsim.event.*;
import java.util.*;
import org.lcsim.util.step.TrackStepper; // C.M. Stepper
-import org.lcsim.material.*;
-//import org.lcsim.material.RadLengths;
-import org.lcsim.recon.muon.*;
+import org.lcsim.util.step.DeDx;
+
import org.lcsim.geometry.Detector;
import org.lcsim.geometry.FieldMap;
import org.lcsim.geometry.subdetector.CylindricalBarrelCalorimeter;
@@ -19,9 +14,7 @@
import org.lcsim.geometry.util.IDDescriptor;
import org.lcsim.geometry.layer.LayerSlice;
import org.lcsim.digisim.CalHitMapMgr;
-
-// import hep.lcd.util.material.DeDx;
-import org.lcsim.util.step.DeDx;
+import org.lcsim.recon.cluster.density.MyTools;
// Muon Segment Finder for barrel calorimeter
// Modified to include the EM Calorimeter- C.Milstene
@@ -30,15 +23,16 @@
// elements as well- C. Milstene-January-02-04
public class BarrelCalSegmentFinder extends SegmentFinder
{
- public BarrelCalSegmentFinder(Detector aDet, String name)
+ public BarrelCalSegmentFinder(Detector aDet, String subdetectorName)
{
if(dedx==null) dedx = DeDx.instance();
+ subdetName = subdetectorName;
// Get the parameters that describe this detector component
det = aDet;
- CylindricalBarrelCalorimeter calsub = (CylindricalBarrelCalorimeter)det.getSubdetectors().get(name);
-
+ CylindricalBarrelCalorimeter calsub = (CylindricalBarrelCalorimeter)det.getSubdetectors().get(subdetName);
+
/**
* Get materials (of 1st layer only) from LayeredCalorimeter class.
*
@@ -49,9 +43,18 @@
for ( Object o : calsub.getLayering().getLayers().getLayer(0).getSlices() )
{
LayerSlice ls = (LayerSlice) o;
- materials.put(ls.getMaterial().getName(), ls.getThickness() );
- }
- //calsub.getMaterials();
+ String sublayerMat = ls.getMaterial().getName();
+ if( materials.containsKey(sublayerMat) ) {
+ // if material exists, just add to thickness
+ double thickness = materials.get( sublayerMat );
+ thickness += ls.getThickness();
+ materials.put( sublayerMat, thickness );
+ }
+ else {
+ // a new material in this cal component
+ materials.put(ls.getMaterial().getName(), ls.getThickness() );
+ }
+ }
rmin = calsub.getInnerRadius();
rmax = calsub.getOuterRadius();
@@ -67,25 +70,26 @@
nTheta = projCal.getThetaBins();
phiBinSize = 2*Math.PI/nPhi;
thetaBinSize = Math.PI/nTheta;
- aname = name; // C.M.-June-03
- System.out.println("MY NAME IS "+aname+" rin="+rmin+" rout="+rmax+" zmax="+zmax);
+ System.out.println("MY NAME IS "+subdetName+" rin="+rmin+" rout="+rmax+" zmax="+zmax+", phiBinSize="+phiBinSize+", thetaBinSize="+thetaBinSize);
// Set the nearest neighbor cuts to default values
- if(aname.equals(EcalName))
+ if(subdetName.equals(ecalSubdetName))
{
phiNNCut =2;
thetaNNCut = 2;
}
- else if(aname.equals(HcalName))
+ else if(subdetName.equals(hcalSubdetName))
{
phiNNCut =3;
thetaNNCut = 3;
}
}
- // FIXME: This code is specific to SDJan03 - it would be nice to code
- // this up independent of detector geometry - to be read from
- // compact.xml or .lcdd
- String description = "layer:7,system:3,barrel:3,theta:32:11,phi:11";
+ // FIXME: This code is specific to SDJan03 - it would be nice to code
+ // this up independent of detector geometry - to be read from
+ // compact.xml or .lcdd
+ // Look up on <detector>/compact.xml for descriptions (under EcalBarrHits and similars)
+ // This is sidaug05 description (for all cal components)
+ String description = "layer:7,system:6,barrel:3,theta:32:11,phi:11";
try {
IDDescriptor desc = new IDDescriptor(description);
encoder = new IDEncoder(desc);
@@ -95,8 +99,10 @@
}
// Set system and barrel/endcap values
int sys = 0;
- if( name.contains("Ecal") ) sys = 2;
- if( name.contains("Hcal") ) sys = 3;
+ System.out.println("BCSF constructor: subdetName="+subdetName);
+ if( subdetName.equals(ecalSubdetName) ) sys = 2;
+ if( subdetName.equals(hcalSubdetName) ) sys = 3;
+ if( subdetName.equals(muSubdetName) ) sys = 4;
encoder.setValue("system", sys);
encoder.setValue("barrel", 0);
}
@@ -116,8 +122,8 @@
// C.M. Dec-03- Tkparams to stepper
double tmpxyz ; double tmpp;
r=rmin;
- System.out.println(aname);
- if(aname.equals("EMBarrel"))
+ System.out.println("findSegment: subd="+subdetName);
+ if(subdetName.equals(ecalSubdetName))
{
rpVect = stpr.getNewRp();
// Then do the Tracker up to r=rmin
@@ -133,14 +139,15 @@
System.out.println("Barrel:out of Tracker");
}
- if(aname.equals("HADBarrel"))
+ if(subdetName.equals(hcalSubdetName))
{
+ System.out.println(" !!!=== rmin="+rmin+" r="+r);
// Then air to r=rmin
double[] pos = new double[]{rmin, 0.0, 0.0};
double[] field = det.getFieldMap().getField(pos); // C.M.-June-03
BField = field[2];
rpVect = stpr.getNewRp();
- double rnow= Math.abs(Math.sqrt(rpVect[0]*rpVect[0]+rpVect[1]));
+ double rnow= Math.abs(Math.sqrt(rpVect[0]*rpVect[0]+rpVect[1]*rpVect[1]));
double d_EMHD = Math.abs(r-rnow);
System.out.println(" distance between calorimeters EM-HD "+d_EMHD);
stepConditions = steprConditions((10*d_EMHD),BField,0.0);
@@ -149,8 +156,8 @@
rpVect = stpr.getNewRp();
xNumSteps = 10.;
- stepConditions = steprConditions(xNumSteps,BField,meanDEdxPerLay(aname));
- System.out.println("Barrel:MeanDEdxPerLay="+meanDEdxPerLay(aname));
+ stepConditions = steprConditions(xNumSteps,BField,meanDEdxPerLay(subdetName));
+ System.out.println("Barrel:MeanDEdxPerLay="+meanDEdxPerLay(subdetName));
// Forms a muon segment by looking for calorimeter hits that match
// cells which the track passes through
// Begin by finding intersection of track with innermost layer
@@ -159,7 +166,7 @@
stpr.tkSteps(r, zmax,stepConditions);
rpVect = stpr.getNewRp();
- System.out.println(aname+" !!!!!! rmin="+rmin);
+ System.out.println(subdetName+" !!!!!! rmin="+rmin);
// Get intersections with successive layers as long as track is
// within the detector, and create a list of the cells that
// the track passes through
@@ -170,12 +177,19 @@
int iphi = xyzToPhiBin(rpVect);
int itheta = xyzToThetaBin(rpVect);
+ double xx = rpVect[0];
+ double yy = rpVect[1];
+ double rho = Math.sqrt(xx*xx+yy*yy);
+// System.out.println("matchHits: ("+rpVect[0]+"; "+rpVect[1]+"; "+rpVect[2]+"), rho="+rho);
+
+// System.out.println(" CellID: layer="+nr+", iphi="+iphi+", itheta="+itheta);
encoder.setValue("layer", nr);
encoder.setValue("theta", itheta);
encoder.setValue("phi", iphi);
- Long cellid = encoder.getID();
+ long cellid = encoder.getID();
+ System.out.println("adding to cellList: cellid="+MyTools.printID(cellid));
cellList.add(cellid);
-
+
if(stpr.getStopTkELow()) break;
if(nr!=(nLayers-1)) r = r + dr;
else {
@@ -183,8 +197,7 @@
} // Go on to the next layer
stpr.tkSteps(r, zmax,stepConditions);
rpVect = stpr.getNewRp();
- nr++;
-
+ nr++;
}
// Save number of layers track passed through
numLayersHit = nr;
@@ -195,16 +208,21 @@
// Now match the cells on this list to calorimeter hits
- if(aname.equals(HcalName)) //C.M.-27Jan03
+// System.out.println("Names: subdetName=<"+subdetName+">"
+// +", EcalSubdetName="+ecalSubdetName
+// +", HcalSubdetName="+hcalSubdetName);
+ if(subdetName.equals(hcalSubdetName)) //C.M.-27Jan03
{
// HistogramFolder.setDefaultFolder("/HDCal");
- matchHitsFast(segm, dataMgr.getCollHitMap(HcalName) );
+// System.out.println("BCSF: Doing matchHits in HCal");
+ matchHitsFast(segm, dataMgr.getCollHitMap(hcalHitmapName) );
// HistogramFolder.setDefaultFolder("..");
}
- else if(aname.equals(EcalName)) //
+ else if(subdetName.equals(ecalSubdetName)) //
{
// HistogramFolder.setDefaultFolder("/EMCal");
- matchHitsFast(segm, dataMgr.getCollHitMap(EcalName) );
+// System.out.println("BCSF: Doing matchHits in ECal");
+ matchHitsFast(segm, dataMgr.getCollHitMap(ecalHitmapName) );
// HistogramFolder.setDefaultFolder(".."); // C.M.-27Jan03- Folder
}
@@ -212,8 +230,7 @@
public void reset()
{
- // Erase lists of hits and cells
-
+ // Erase lists of hits and cells
hitList.removeAllElements();
cellList.removeAllElements();
}
@@ -302,7 +319,7 @@
//--------------------------------------------C.M.- Get DE/dx per Calorimeter Layer
protected double meanDEdxPerLay(String name)
{
-// String compoundMaterial = name+"_Material";
+// String compoundMaterial = name+"_Material";
double matterDEdx=0.;
// MaterialList emml = (MaterialList) det.get(name+"_Material");
int i=0;
@@ -311,9 +328,10 @@
// Material emmat = emml.getMaterial(i);
// String mat = emml.getName();
// double tthk = emml.getThickness();
- i++;
+// i++;
double tthk = materials.get(mat).doubleValue();
double e1 = tthk * (dedx.getDeDx(mat)/10); // energy loss this material
+// System.out.println("mat="+mat+", tth="+tthk+", e1="+e1);
//System.out.println("matName="+mat+", thickness="+tthk+", dEdx="+dedx.getDeDx(mat));
matterDEdx += e1; // sum energy loss per material
//System.out.println(" DedxPerMateriel (GeV/cm) ="+dedx.getDeDx(mat)+" thickness(mm)="+tthk+ " element="+i);
@@ -328,79 +346,60 @@
protected void matchHitsFast(Segmentation cell,
Map<Long,SimCalorimeterHit> calHits)
{
+ this.matchHitsFast(cell, calHits, false);
+ }
+
+ protected void matchHitsFast(Segmentation cell,
+ Map<Long,SimCalorimeterHit> calHits,
+ boolean debug)
+ {
+ int nhitsTotal = 0;
+ ProjectiveCylinder projCal = (ProjectiveCylinder)segm;
+ if(debug) System.out.println("matchhits: nhits = "+calHits.size());
for( Long cellid : cellList ) {
if(segm instanceof ProjectiveCylinder) {
- ProjectiveCylinder projCal = (ProjectiveCylinder)segm;
- long[] neighs = projCal.getNeighbourIDs(0, thetaNNCut, phiNNCut );
+ // check cell itself
+ SimCalorimeterHit hit = calHits.get(cellid);
+ if(hit!=null) {
+ if(debug) System.out.println("matchHits: Adding hit: "+Long.toHexString(cellid));
+ hitList.add(hit);
+ nhitsTotal++;
+ }
+ // now check for the neighbors
+ projCal.setID(cellid);
+ long[] neighs = projCal.getNeighbourIDs(0, thetaNNCut, phiNNCut );
for(int i=0; i<neighs.length; ++i) {
- SimCalorimeterHit hit = calHits.get(neighs[i]);
- hitList.add(hit);
+ hit = calHits.get(neighs[i]);
+ if(debug) System.out.println("neighb: cellid="+MyTools.printID(neighs[i])+", hit="+hit);
+ if(hit!=null) {
+ if(debug) System.out.println("matchHits: Adding hit: "+Long.toHexString(cellid));
+ hitList.add(hit);
+ nhitsTotal++;
+ }
}
}
}
+ if(debug) System.out.println("Total # hits: virtual="+cellList.size()
+ +", real="+nhitsTotal);
}
-// protected void matchHits(Segmentation cell, List<SimCalorimeterHit> calHits)
-// {
-
-// // Match cells the track passes through with calorimeter hits
-
-// // Begin by looping over calorimeter hits
-
-// for (Enumeration eCal = calHits.getHits(); eCal.hasMoreElements();)
-// {
-// CalorimeterHit hit = (CalorimeterHit)eCal.nextElement();
-
-// // Use this event's CalorimeterCell object to decode
-// // the TowerIDs
-
-// cell.setTowerID(hit.getTowerID());
-// int cellPhi = cell.getPhiBin();
-// int cellTheta = cell.getThetaBin();
-// int cellLayer = cell.getLayer();
-
-// // Loop over cells the track passed through. If one of these is
-// // close enough to this calorimeter hit, add it to the list
-// int dPhiTkHit;
-// for (Enumeration eTrack = cellList.elements(); eTrack.hasMoreElements();)
-// {
-// int [] cellIndices = (int [])eTrack.nextElement();
-// if(cellIndices[2] == cellLayer)
-// {
-// dPhiTkHit = cellIndices[0]-cellPhi;
-// if(Math.abs(cellIndices[0] - cellPhi) >=(nPhi-2*phiNNCut)) // C.M.
-// {
-
-// int sgnDPhi= ((cellIndices[0] - cellPhi)>0)?1:-1;
-// dPhiTkHit=(Math.abs(cellIndices[0]-cellPhi)>nPhi-phiNNCut)?
-// sgnDPhi*(Math.abs(cellIndices[0]-cellPhi)-nPhi):(cellIndices[0] - cellPhi);
-// }// C.M.
-
-
-// }
-// if ( ((Math.abs(cellIndices[0] - cellPhi) <= phiNNCut)
-// || (Math.abs(cellIndices[0] - cellPhi) >= (nPhi-phiNNCut)) ) // C.M.- Added (nPhi-phiNNCut) in conditions
-// && Math.abs(cellIndices[1] - cellTheta) <= thetaNNCut &&
-// cellIndices[2] == cellLayer)
-// {
-// hitList.add(hit);
-// }
-// }
-// }
-// }
+
protected Detector det;
protected Segmentation segm;
protected IDEncoder encoder;
protected CalHitMapMgr dataMgr;
- protected String EcalName = "EcalBarrHits";
- protected String HcalName = "HcalBarrHits";
+ protected String ecalHitmapName = "EcalBarrHits";
+ protected String hcalHitmapName = "HcalBarrHits";
+ protected String ecalSubdetName = "EMBarrel";
+ protected String hcalSubdetName = "HADBarrel";
+ protected String muSubdetName = "MuonBarrel";
protected int nc=0;
protected double BField ;// C.M.-Jan03
protected DeDx dedx;
protected double []rpVect = new double[8] ; // C.M.-4-Vector
- protected String aname; // C.M.-Jan03
+ protected String subdetName; // C.M.-Jan03
protected double rmin;
protected double rmax;
protected double zmax;