lcsim/src/org/lcsim/recon/muon
diff -N EndcapCalSegmentFinder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ EndcapCalSegmentFinder.java 9 Oct 2006 22:30:49 -0000 1.1
@@ -0,0 +1,424 @@
+package org.lcsim.recon.muon;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import org.lcsim.event.*;
+import java.io.*;
+import java.util.*;
+import org.lcsim.util.step.TrackStepper; // C.M. Stepper
+import org.lcsim.recon.cluster.util.Samplexmlclass;
+import org.lcsim.util.step.DeDx;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.subdetector.*;
+import org.lcsim.geometry.compact.Subdetector;
+import org.lcsim.geometry.segmentation.*;
+import org.lcsim.geometry.segmentation.SegmentationBase;
+import org.lcsim.geometry.util.IDEncoder;
+import org.lcsim.geometry.util.IDDescriptor;
+import org.lcsim.geometry.layer.LayerSlice;
+import org.lcsim.recon.cluster.util.CalHitMapMgr;
+import org.lcsim.geometry.layer.*;
+import org.lcsim.geometry.Layered;
+import org.lcsim.geometry.LayeredSubdetector;
+import org.lcsim.geometry.field.FieldOverlay;
+/**
+ * @authors Judith Odili and Caroline Milstene- July 2006
+ * Expension of the Muon package to the EndCaps
+ *
+ */
+public class EndcapCalSegmentFinder extends SegmentFinder// extends BarrelCalSegmentFinder
+{
+ public EndcapCalSegmentFinder(Detector aDet, String subdetectorName)
+ {
+
+ //super(aDet, subdetectorName);
+ AIDA aida = AIDA.defaultInstance();
+
+ if(dedx==null) dedx = DeDx.instance();
+ subdetName = subdetectorName;
+
+ // Get the parameters that describe this detector component
+ det = aDet;
+ String a_nam = det.getHeader().getDetectorName();
+ endcapCal = (CylindricalEndcapCalorimeter)det.getSubdetectors().get(subdetName);
+ Layering layers = ((Layered) endcapCal).getLayering();
+ rmin = endcapCal.getInnerRadius();
+ rmax = endcapCal.getOuterRadius();
+ zmax = endcapCal.getZMax();
+ zmin = endcapCal.getZMin();
+ nLayers = layers.getLayerCount();
+ zLayer0Mid = endcapCal.getLayering().getDistanceToLayerSensorMid(0);
+ dataMgr = CalHitMapMgr.getInstance();
+ segm = (SegmentationBase)endcapCal.getReadout().getSegmentation();
+
+ if(subdetName.equals(ecalendcapSubdetName)) {
+ phiNNCut = 1; // 2;
+ thetaNNCut = 1; // 2;
+ } else
+ if(subdetName.equals(hcalendcapSubdetName)) {
+ phiNNCut = 1; // 3;
+ thetaNNCut = 1; // 2;
+ } else
+ if (subdetName.equals(muendcapSubdetName))
+ {
+ phiNNCut = 2; // 3;
+ thetaNNCut = 2; // 2;
+ }
+
+ Subdetector subDet = aDet.getSubdetector(subdetName);
+ int nLayers = ((Layered) subDet).getLayering().getLayerCount();
+ SegmentationBase base = ((SegmentationBase) subDet.getIDDecoder());
+ // Set system and barrel/endcap values
+ if(debug)System.out.println("ECSF constructor: subdetName="+subdetName);
+
+
+ }
+
+ //----------------------------------!!!
+ public void findSegment(EventHeader event, TrackStepper inStpr)
+ //----------------------------------!!!
+ {
+ stpr = inStpr;
+ int nz = 0; // layer counter
+ int nLayerStepped=0; // curlback layers included
+ double dz = (zmax - zmin)/nLayers; // layer thickness
+ numLayersHit = 0;
+
+
+ double z;
+ double tmpxyz ; double tmpp;
+ String materialName;
+ z=zmin;
+ // UPDATE Phase-Space
+ rpVect = stpr.getNewRp();
+ if(debug)
+ System.out.println("findSegment: subd="+subdetName);
+ if(stpr.getStopTkELow())
+ return;
+ bField=det.getFieldMap().getField(rpVect);
+ if(subdetName.equals(ecalendcapSubdetName)) {
+ // Then do the Tracker up to r=rmin
+ double dTracker=(z-rpVect[2]);
+ materialName = "Air";
+ stpr.tkSteps( z, rmax,"Endcap",bField[2], materialName );
+ rpVect = stpr.getNewRp();
+ if(z<(zmin-dz))
+ return; // does not enter ECAL
+ }
+ if(subdetName.equals(hcalendcapSubdetName)) {
+ if(debug)System.out.println(" !!!=== zmin="+zmin+" z="+z);
+ // Then air to r=rmin
+ rpVect = stpr.getNewRp();
+ double d_EMHD = Math.abs(z-rpVect[2]);
+ if(debug)System.out.println(" distance between calorimeters EM-HD "+d_EMHD);
+ materialName = "Air";
+ stpr.tkSteps( z, rmax,"Endcap",bField[2], materialName );
+ rpVect = stpr.getNewRp();
+ // if(r<(rmin-dr))return; // does not enter HCAL
+ }
+
+ if(subdetName.equals(muendcapSubdetName)) {
+ if(debug)System.out.println(" !!!=== zmin="+zmin+" z="+z);
+ // Then air to r=rmin
+ rpVect = stpr.getNewRp();
+ double d_HDMU = Math.abs(z-rpVect[2]);
+ if(debug)System.out.println(" distance between calorimeters HD-MU "+d_HDMU);
+ materialName = "Air";
+ stpr.tkSteps( z, rmax,"Endcap",bField[2], materialName );
+ rpVect = stpr.getNewRp();
+ // if(r<(rmin-dr))return; // does not enter HCAL
+ }
+
+ if(debug)System.out.println("Stepper pos ENTRY of:("+subdetName+" "+rpVect[0]+"; "+rpVect[1]+"; "+rpVect[2]+"), r="+stpr.partR(rpVect));
+
+ // 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
+ if(debug)System.out.println("subdet="+subdetName+", sensormid="+rLayer0Mid);
+ z = zLayer0Mid;
+ stpr.tkSteps(z, rmax,det, subdetName, nz);
+ rpVect = stpr.getNewRp();
+ if(debug)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
+ nLayerStepped=0;
+ while ((stpr.partR(rpVect)<= rmax )&& (nz < nLayers)) {
+
+ long virtualId = segm.findCellContainingXYZ(rpVect);
+ cellList.add(virtualId);
+
+ Hep3Vector v = new BasicHep3Vector(rpVect[0],rpVect[1],rpVect[2]);
+ hitPos.add(v);
+ if(stpr.getStopTkELow()) break;
+
+ if(stpr.getCurlBack()) {
+ //if(debug)System.out.println("!!!BCF:CURLINGBACK old r="+r+" New r="+(r-dr));
+ z=z-dz ;
+ } else if(nz!=(nLayers-1)) z = z + dz;
+ else {
+ z = zmin + nLayers * dz;
+ if(debug)System.out.println("ECF: AT Exit of "+subdetName+", z="+z);
+ }
+
+ stpr.tkSteps(z, rmax,det, subdetName, nz);
+ rpVect = stpr.getNewRp();
+ // double dThk=(nr!=(nLayers-1))?dr:(dr+dr-(rLayer0Mid-rmin));
+
+ if(!stpr.getCurlBack())nz++;
+ nLayerStepped++;
+ }//End of : Go on to the next layer
+ // Save number of layers track passed through
+ numLayersHit = nz;
+ if(debug)System.out.println("Stepper pos EXIT of:("+subdetName+" "+rpVect[0]+"; "+rpVect[1]+"; "+rpVect[2]+"), r="+stpr.partR(rpVect));
+ rpVect= stpr.getNewRp();
+ if((rpVect[2]>=zmax)||stpr.getStopTkELow()) {// Record Last Layer and Info at Break
+ long virtualId = segm.findCellContainingXYZ(rpVect);
+ cellList.add(virtualId);
+ Hep3Vector v =new BasicHep3Vector(rpVect[0],rpVect[1],rpVect[2]);
+ hitPos.add(v);
+ }
+
+ double rr=stpr.partR(rpVect);
+ // Now match the cells on this list to calorimeter hits
+ if(subdetName.equals(hcalendcapSubdetName)) //C.M.-27Jan03
+ {
+
+// System.out.println("BCSF: Doing matchHits in HCal");
+ matchHitsFast(segm, dataMgr.getCollHitMap(hcalendcapHitmapName) );
+
+ }
+ else if(subdetName.equals(ecalendcapSubdetName)) //
+ {
+// System.out.println("BCSF: Doing matchHits in ECal");
+ matchHitsFast(segm, dataMgr.getCollHitMap(ecalendcapHitmapName) );
+ }
+
+ else if(subdetName.equals(muendcapSubdetName)) //
+ {
+// System.out.println("BCSF: Doing matchHits in ECal");
+ matchHitsFast(segm, dataMgr.getCollHitMap(muendcapHitmapName) );
+ }
+ }
+
+ public void reset() {
+ // Erase lists of hits and cells
+ hitList.removeAllElements();
+ cellList.removeAllElements();
+
+ }
+
+ public void setPhiNNCut(int cut) {
+
+ // Set the number of nearest neighbors cut in phi. Cells
+ // with abs(phi_track - phi_cell) <= cut will be added
+ // to the list
+
+ phiNNCut = cut;
+ }
+
+ public void setThetaNNCut(int cut) {
+
+ // Set the number of nearest neighbors cut in theta. Cells
+ // with abs(theta_track - theta_cell) <= cut will be added
+ // to the list
+
+ thetaNNCut = cut;
+ }
+
+ public double[] getNNCuts() {
+ // Return the number of nearest neighbors cuts in phi and theta
+
+ return(new double[]
+ {phiNNCut, thetaNNCut});
+ }
+
+ public int getNLayersHit() {
+ // Return the number of layers the track passed through
+
+ return(numLayersHit);
+ }
+
+
+ public int xyzToPhiBin(double[] r) {
+ // Determine phi bin corresponding to this Cartesian coordinate
+
+ double phi;
+ int bin;
+ nc++;
+ double r1=r[1]/10; double r0=r[0]/10;
+ phi = Math.atan2(r1, r0);
+ if (phi < 0.) phi = phi + 2*Math.PI;
+ bin = (int)(phi/phiBinSize);
+ return(bin);
+ }
+
+ public int xyzToThetaBin(double[] r) {
+ // Determine theta bin corresponding to this Cartesian coordinate
+ double rho, theta;
+ int bin;
+ double r1=r[1]/10; double r0=r[0]/10; double r2=r[2]/10;
+ rho = Math.sqrt(r0*r0 + r1*r1);
+ theta = Math.atan2(rho, r2);
+ if (theta < 0.) theta = theta + Math.PI;
+ bin = (int)(theta/thetaBinSize);
+ return(bin);
+ }
+
+
+ // The following two methods implement the CalorimeterHits interface
+
+ public Enumeration getHits() {
+ return(hitList.elements());
+ }
+ public Enumeration getHitsPos() {
+ return(hitPos.elements());
+ }
+
+ public int getNHits() {
+ return(hitList.size());
+ }
+ public int getNHitsPos() {
+ return(hitPos.size());
+ }
+
+
+ protected void matchHitsFast(SegmentationBase cell,
+ Map<Long,CalorimeterHit> calHits) {
+ this.matchHitsFast(cell, calHits, true);
+ }
+
+ protected void matchHitsFast(SegmentationBase cell,
+ Map<Long,CalorimeterHit> calHits,
+ boolean debug) {
+ if(debug)System.out.println("Segm:::: " + cell);
+ AIDA aida = AIDA.defaultInstance();
+ int nhitsTotal = 0;
+ // ProjectiveCylinder projCal = (ProjectiveCylinder)segm;
+ if(debug)System.out.println("matchhits: nhits = "+calHits.size());
+ for( Long cellid : cellList ) {
+// if(segm instanceof ProjectiveCylinder) {
+ // check cell itself
+ //System.out.println("IT ENTERED AND CELL ===" + cell);
+ CalorimeterHit hit = calHits.get(cellid);
+ System.out.println("hit ===" + hit);
+ segm.setID(cellid);
+
+ int ilay = segm.getLayer();
+ System.out.println("ILAY ===" + ilay);
+ if(hit!=null) {
+ // found real hit at virtual cell
+ if(debug) System.out.println("matchHits: Adding hit: "+Long.toHexString(cellid));
+ if(debug) {
+ double dall=0.;
+ for(int j=0;j<nLayers;j=j+4) {
+ if(ilay==j) {
+ aida.cloud1D(subdetName+" dphi layer "+ilay).fill(dall);
+ aida.cloud1D(subdetName+" dtheta layer "+ilay).fill(dall);
+ }
+ }
+ }
+ foundHits.add(new BasicHep3Vector(segm.getX(),segm.getY(),segm.getZ()));
+ hitList.add(hit);
+ nhitsTotal++;
+ }
+
+ // now check for the neighbors
+ segm.setID(cellid);
+ int iphi = segm.getValue("phi");
+ System.out.println("IPHI ===" + iphi);
+ int ithe = -99999;
+ try {
+ ithe = segm.getValue("z");
+ System.out.println("ITHE ===" + ithe);
+ } catch(NullPointerException x) {
+ // System.out.println("IN Catch");
+ ithe = segm.getValue("z");
+ }
+ ilay = segm.getLayer();
+ long[] neighs = segm.getNeighbourIDs(0, thetaNNCut, phiNNCut );
+// System.out.println("VirtualID="+MyTools.printID(cellid)
+// +", #neighbors="+neighs.length+" "+thetaNNCut+" "+phiNNCut);
+ for(int i=0; i<neighs.length; ++i) {
+ 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(neighs[i]));
+ segm.setID( neighs[i] );
+ //===========================
+ int jphi = segm.getValue("phi");
+ int jthe = 0;
+ try {
+ jthe = segm.getValue("z");
+ } catch(NullPointerException x) {
+ //System.out.println(".......in catch 2");
+ jthe = segm.getValue("z");
+ }
+ int jlay = segm.getLayer();
+// // if(debug)
+// //{
+ for(int j=0;j<nLayers;j=j+4) {
+ if(ilay==j) {
+ aida.cloud1D(subdetName+" dphi layer "+ilay).fill(jphi-iphi);
+ aida.cloud1D(subdetName+" dtheta layer "+ilay).fill(jthe-ithe);
+ }
+ }
+// // }
+//====================================================
+ foundHits.add(new BasicHep3Vector(segm.getX(),segm.getY(),segm.getZ()));
+ hitList.add(hit);
+ nhitsTotal++;
+ }
+ }
+// }
+ }
+ if(debug) System.out.println("Total # hits: virtual="+cellList.size()
+ +", real="+nhitsTotal);
+ }
+
+
+
+ protected Detector det;
+ protected CylindricalEndcapCalorimeter endcapCal;
+ protected SegmentationBase segm;
+ protected IDEncoder encoder;
+ protected CalHitMapMgr dataMgr;
+ protected String ecalendcapHitmapName = "EcalEndcapHits";
+ protected String hcalendcapHitmapName = "HcalEndcapHits";
+ protected String muendcapHitmapName = "MuonEndcapHits";
+
+ protected String ecalendcapSubdetName = "EMEndcap";
+ protected String hcalendcapSubdetName = "HADEndcap";
+ protected String muendcapSubdetName = "MuonEndcap";
+
+
+
+ protected boolean debug = false;
+ 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 subdetName; // C.M.-Jan03
+ protected double rmin;
+ protected double rLayer0Mid, zLayer0Mid, totalthickness;
+ protected double rmax;
+ protected double zmax, zmin;
+ protected int nLayers;
+ protected int nPhi;
+ protected int nTheta;
+ protected double phiBinSize;
+ protected double thetaBinSize;
+ protected int numLayersHit;
+ protected Vector<Long> cellList = new Vector<Long>(); // list of cell indices
+ protected TrackStepper stpr = new TrackStepper(); //C.M.-Stepper Instance
+ //protected double [] stepConditions = new double[4];
+ protected Map<String, Double> materials;
+ public static Samplexmlclass xmlReader = new Samplexmlclass();
+ private int phiNNCut;
+ private int thetaNNCut;
+ private List<Hep3Vector> foundHits = new ArrayList<Hep3Vector>();
+ protected FieldOverlay Field = new FieldOverlay();
+ protected double[] bField = new double[3];
+}
+
