Commit in lcsim-contrib/src/main/java/org/lcsim/contrib/SteveMagill on MAIN | |||
TrMipDriver.java | +655 | added 1.1 |
diff -N TrMipDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ TrMipDriver.java 14 Feb 2012 19:21:18 -0000 1.1 @@ -0,0 +1,655 @@
+package org.lcsim.contrib.SteveMagill; + +// This driver extrapolates tracks into CAL, associating mips and determining interaction layer +// output : Mip Clusters and IL per track + +import java.util.*; +import org.lcsim.event.CalorimeterHit; +import org.lcsim.event.Track; +import org.lcsim.event.EventHeader; +import org.lcsim.util.Driver; +import org.lcsim.util.swim.*; +import org.lcsim.recon.cluster.util.*; +import org.lcsim.util.aida.*; +import org.lcsim.geometry.*; +import hep.physics.vec.Hep3Vector; +import hep.physics.vec.BasicHep3Vector; +import org.lcsim.geometry.subdetector.CylindricalCalorimeter; +import hep.aida.*; +import org.lcsim.spacegeom.*; +import org.lcsim.util.hitmap.HitMap; + +public class TrMipDriver extends Driver +{ + ITree tree; + IHistogramFactory histogramFactory; + IAnalysisFactory analysisFactory = IAnalysisFactory.create(); + private AIDA aida = AIDA.defaultInstance(); + private boolean _initialized = false; +// The number of layers in the EM Calorimeter Barrel, Endcap + private int _numbemlayers; + private int _numecemlayers; +// The number of layers in the Hadronic Calorimeter Barrel, Endcap + private int _numbhadlayers; + private int _numechadlayers; +// The EM calorimeter hits, by layer in barrel, endcap + private List[] _emBLayerHits; + private List[] _emECLayerHits; +// The HAD calorimeter hits, by layer in barrel, endcap + private List[] _hadBLayerHits; + private List[] _hadECLayerHits; +// The radii of the barrel calorimeter layers + private double[] _emBRadii; + private double[] _hadBRadii; + private double[] BRadii = new double[100]; + private double[] ECZs = new double[100]; +// modified Barrel helix swimmer for org.lcsim test +// private Helix[] _emBSwimmers; //EM +// private Helix[] _hadBSwimmers; //HAD +// Z extent of the central barrel calorimeters. + private double _embZ; //EM Barrel Z + private double[] _emecZ; //EM Endcap Z + private double _hadbZ; //HAD barrel z + private double[] _hadecZ; // HAD endcap z + private double _dminE; + private double _dminH; + private double _dminTrC; + private double _mincd; + private String _nameExt; + private String[] _hitcollnames; + private String _tmname; + private String _tilname; + private String _intname; + private boolean trmipD = false; +// private boolean trmipRes = false; +// private SpacePoint trackspacept; + private double zField; +// private String[] _hitmapnames; + + public TrMipDriver(double dminE, double dminH, double dminTrC, double mincd) + { + _dminE = dminE; + _dminH = dminH; + _dminTrC = dminTrC; + _mincd = mincd; + } + + protected void process(EventHeader event) + { + super.process(event); // executes all added drivers + +// Initialize things here + if(!_initialized) + { + // setup specific detector stuff here + Detector det = event.getDetector(); +// double[] zero = {0, 0, 0}; +// _fieldStrength = det.getFieldMap().getField(zero); + // try this + Hep3Vector ip = new BasicHep3Vector(); + zField = det.getFieldMap().getField(ip).z(); + if (trmipD) System.out.println("B Field " + zField + " Tesla"); + + // Setup layer counts and distances for barrel, endcap + CylindricalCalorimeter emb = ((CylindricalCalorimeter) det.getSubdetectors().get("EMBarrel")); + _embZ = emb.getZMin(); + _numbemlayers = emb.getLayering().getLayerCount(); + if (trmipD) System.out.println("EM Barrel Layers " +_numbemlayers); + for (int i=0; i<_numbemlayers; ++i) + { + BRadii[i]=emb.getLayering().getDistanceToLayerSensorMid(i); + if (trmipD) System.out.println("EM Barrel Layer Number " +i+ " EM Barrel Radius " + BRadii[i]); + } + + CylindricalCalorimeter emec = ((CylindricalCalorimeter) det.getSubdetectors().get("EMEndcap")); + _numecemlayers = emec.getLayering().getLayerCount(); + if (trmipD) System.out.println("EM EC Layers " +_numecemlayers); + for (int i=0; i<_numecemlayers; ++i) + { + ECZs[i] = emec.getLayering().getDistanceToLayerSensorMid(i); + if (trmipD) System.out.println("EM Endcap Layer Number " +i+ " EM Endcap Z " +ECZs[i]); + } + + CylindricalCalorimeter hadb = ((CylindricalCalorimeter) det.getSubdetectors().get("HADBarrel")); + _hadbZ = hadb.getZMin(); + _numbhadlayers = hadb.getLayering().getLayerCount(); + if (trmipD) System.out.println("HAD Barrel Layers " +_numbhadlayers); + for (int i=0; i<_numbhadlayers; ++i) + { + BRadii[i+_numbemlayers]=hadb.getLayering().getDistanceToLayerSensorMid(i); + if (trmipD) System.out.println("HAD Barrel Layer Number " +i+ " HAD Barrel Radius " + BRadii[i+_numbemlayers]); + } + + CylindricalCalorimeter hadec = ((CylindricalCalorimeter) det.getSubdetectors().get("HADEndcap")); + _numechadlayers = hadec.getLayering().getLayerCount(); + if (trmipD) System.out.println("HAD Endcap Layers " +_numechadlayers); + for (int i=0; i<_numechadlayers; ++i) + { + ECZs[i+_numecemlayers] = hadec.getLayering().getDistanceToLayerSensorMid(i); + if (trmipD) System.out.println("HAD Endcap Layer Number " +i+ " HAD Endcap Z " +ECZs[i]); + } + + _initialized = true; + } // end of initialization section + + // Need cells in CAL, but need to have them arranged in layers with density for each cell + // only done once, and only internal to this piece of code + // variables for barrel + int nhitsBTotCAL=0; int nhitsEMB=0; int nhitsHADB=0; + double[][] celbtheta = new double[10000][100]; + double[][] celbphi = new double[10000][100]; + double[][] celbden = new double[10000][100]; + long[][] celbID = new long[10000][100]; + CalorimeterHit[][] celbhit = new CalorimeterHit[10000][100]; + int[] nhitbLayer = new int[200]; + // variables for endcap + int nhitsECTotCAL=0; int nhitsEMEC=0; int nhitsHADEC=0; + double[][] celectheta = new double[10000][100]; + double[][] celecphi = new double[10000][100]; + double[][] celecden = new double[10000][100]; + long[][] celecID = new long[10000][100]; + CalorimeterHit[][] celechit = new CalorimeterHit[10000][100]; + int[] nhitecLayer = new int[200]; + + for (int i=0; i<_hitcollnames.length; i++) + { + List<CalorimeterHit> chits = event.get(CalorimeterHit.class,_hitcollnames[i]); + IDDecoder caldecoder = event.getMetaData(chits).getIDDecoder(); + if (trmipD) System.out.println("Collection is " +_hitcollnames[i]); + if (_hitcollnames[i].equals("Ceren_EcalBarrDigiHits")) + { + for (CalorimeterHit chit : chits) + { + long cID = chit.getCellID(); + caldecoder.setID(cID); + int nlay = caldecoder.getLayer(); + double[] htpos = chit.getPosition(); + double cellphi = Math.atan2(htpos[1],htpos[0]); + if (cellphi<0) cellphi+=2*Math.PI; + double htr = Math.sqrt(htpos[0]*htpos[0]+htpos[1]*htpos[1]); + double celltheta = Math.atan(htr/htpos[2]); + if (celltheta<0) celltheta+=Math.PI; + celbtheta[nhitbLayer[nlay]][nlay] = celltheta; + celbphi[nhitbLayer[nlay]][nlay] = cellphi; + celbID[nhitbLayer[nlay]][nlay] = cID; + celbhit[nhitbLayer[nlay]][nlay] = chit; + nhitbLayer[nlay]++; + nhitsEMB++; + nhitsBTotCAL++; + } + } else if (_hitcollnames[i].equals("Ceren_HcalBarrDigiHits")) + { + for (CalorimeterHit chit : chits) + { + long cID = chit.getCellID(); + caldecoder.setID(cID); + int nlay = _numbemlayers + caldecoder.getLayer(); + double[] htpos = chit.getPosition(); + double cellphi = Math.atan2(htpos[1],htpos[0]); + if (cellphi<0) cellphi+=2*Math.PI; + double htr = Math.sqrt(htpos[0]*htpos[0]+htpos[1]*htpos[1]); + double celltheta = Math.atan(htr/htpos[2]); + if (celltheta<0) celltheta+=Math.PI; + celbtheta[nhitbLayer[nlay]][nlay] = celltheta; + celbphi[nhitbLayer[nlay]][nlay] = cellphi; + celbID[nhitbLayer[nlay]][nlay] = cID; + celbhit[nhitbLayer[nlay]][nlay] = chit; + nhitbLayer[nlay]++; + nhitsHADB++; + nhitsBTotCAL++; + } + } else if (_hitcollnames[i].equals("Ceren_EcalEndcapDigiHits")) + { + for (CalorimeterHit chit : chits) + { + long cID = chit.getCellID(); + caldecoder.setID(cID); + int nlay = caldecoder.getLayer(); + double[] htpos = chit.getPosition(); + double cellphi = Math.atan2(htpos[1],htpos[0]); + if (cellphi<0) cellphi+=2*Math.PI; + double htr = Math.sqrt(htpos[0]*htpos[0]+htpos[1]*htpos[1]); + double celltheta = Math.atan(htr/htpos[2]); + if (celltheta<0) celltheta+=Math.PI; + celectheta[nhitecLayer[nlay]][nlay] = celltheta; + celecphi[nhitecLayer[nlay]][nlay] = cellphi; + celecID[nhitecLayer[nlay]][nlay] = cID; + celechit[nhitecLayer[nlay]][nlay] = chit; + nhitecLayer[nlay]++; + nhitsEMEC++; + nhitsECTotCAL++; + } + } else if (_hitcollnames[i].equals("Ceren_HcalEndcapDigiHits")) + { + for (CalorimeterHit chit : chits) + { + long cID = chit.getCellID(); + caldecoder.setID(cID); + int nlay = _numecemlayers + caldecoder.getLayer(); + double[] htpos = chit.getPosition(); + double cellphi = Math.atan2(htpos[1],htpos[0]); + if (cellphi<0) cellphi+=2*Math.PI; + double htr = Math.sqrt(htpos[0]*htpos[0]+htpos[1]*htpos[1]); + double celltheta = Math.atan(htr/htpos[2]); + if (celltheta<0) celltheta+=Math.PI; + celectheta[nhitecLayer[nlay]][nlay] = celltheta; + celecphi[nhitecLayer[nlay]][nlay] = cellphi; + celecID[nhitecLayer[nlay]][nlay] = cID; + celechit[nhitecLayer[nlay]][nlay] = chit; + nhitecLayer[nlay]++; + nhitsHADEC++; + nhitsECTotCAL++; + } + } else + { + System.out.println("From TrackMipClusterDriver - IMPROPER COLLECTIONS"); + } + } + if (trmipD) aida.cloud1D("Total number of hits in B CAL").fill(nhitsBTotCAL); + if (trmipD) aida.cloud1D("Total number of hits in EC CAL").fill(nhitsECTotCAL); + // aida.cloud2D("Number of hits B vs EC CAL").fill(nhitsBTotCAL,nhitsECTotCAL); + + // Now calculate densities for hits in barrel +// HitMap lodenshitmap = new HitMap(); +// HitMap hidenshitmap = new HitMap(); + for (int i=0; i<_numbemlayers+_numbhadlayers; ++i) + { + for (int j=0; j<nhitbLayer[i]; ++j) + { + celbden[j][i] = 1; + for (int k=0; k<nhitbLayer[i]; ++k) + { + if (k == j) continue; + double dmin = _dminE; + double deltheta = Math.abs(celbtheta[k][i]-celbtheta[j][i]); + if (trmipD) aida.cloud1D("Delta Theta B Cell").fill(deltheta); + double delphi = Math.abs(celbphi[k][i]-celbphi[j][i]); + if (delphi > Math.PI) delphi = 2*Math.PI-delphi; + if (trmipD) aida.cloud1D("Delta Phi B Cell").fill(delphi); + double distance = Math.sqrt(deltheta*deltheta+delphi*delphi); + if (i >= _numbemlayers) dmin = _dminH; + if (distance < dmin) + { +// if (i < _numbemlayers) celbden[j][i] += 0.0026/distance; // default is 0.0026 +// if (i > _numbemlayers-1) celbden[j][i] += 0.0054/distance; // default was 0.0052, try 0.0062 + celbden[j][i]++; //just add 1 to density per hit in window + if (trmipD) + { + if (i == 0) aida.cloud1D("Distance between EM B CAL hits layer 0").fill(distance); + if (i < _numbemlayers) aida.cloud1D("Distance between EM B CAL hits").fill(distance); + if (i == _numbemlayers) aida.cloud1D("Distance between HAD B CAL hits layer 0").fill(distance); + if (i > _numbemlayers-1) aida.cloud1D("Distance between HAD B CAL hits").fill(distance); + } + } + } + if (trmipD) aida.cloud1D("Cell Density B").fill(celbden[j][i]); + if (trmipD) aida.cloud2D("Cell Density B vs Layer").fill(i,celbden[j][i]); + // aida.cloud2D("Cell Density vs Vis E").fill(celbhit[j][i].getRawEnergy(),celbden[j][i]); +// if (celbden[j][i] > 0 && celbden[j][i] < _mincd) +// { +// lodenshitmap.put(celbID[j][i],celbhit[j][i]); +// if (trmipD) aida.cloud1D("Visible E LowD hits").fill(celbhit[j][i].getRawEnergy()); +// } +// if (celbden[j][i] > _mincd) +// { +// hidenshitmap.put(celbID[j][i],celbhit[j][i]); +// if (trmipD) aida.cloud1D("Visible E HighD hits").fill(celbhit[j][i].getRawEnergy()); +// } + } + } + + // Now calculate densities for hits in endcap + for (int i=0; i<_numecemlayers+_numechadlayers; ++i) + { + for (int j=0; j<nhitecLayer[i]; ++j) + { + celecden[j][i] = 1; + for (int k=0; k<nhitecLayer[i]; ++k) + { + if (k == j) continue; + double dmin = _dminE; + double deltheta = Math.abs(celectheta[k][i]-celectheta[j][i]); + if (trmipD) aida.cloud1D("Delta Theta EC Cell").fill(deltheta); + double delphi = Math.abs(celecphi[k][i]-celecphi[j][i]); + if (delphi > Math.PI) delphi = 2*Math.PI-delphi; + if (trmipD) aida.cloud1D("Delta Phi EC Cell").fill(delphi); + double distance = Math.sqrt(deltheta*deltheta+delphi*delphi); + if (i >= _numecemlayers) dmin = _dminH; + if (distance < dmin) + { +// if (i < _numecemlayers) celecden[j][i] += 0.0028/distance; // default was 0.0017, try 0.0026 +// if (i > _numecemlayers-1) celecden[j][i] += 0.0085/distance; // default was 0.0042, try 0.0062 + celecden[j][i]++; //just add 1 to density per hit in window + if (trmipD) + { + if (i == 0) aida.cloud1D("Distance between EM EC CAL hits layer 0").fill(distance); + if (i < _numecemlayers) aida.cloud1D("Distance between EM EC CAL hits").fill(distance); + if (i > _numecemlayers-1) aida.cloud1D("Distance between HAD EC CAL hits").fill(distance); + if (i == _numecemlayers) aida.cloud1D("Distance between HAD EC CAL hits layer 0").fill(distance); + } + } + } + if (trmipD) aida.cloud1D("Cell Density EC").fill(celecden[j][i]); + if (trmipD) aida.cloud2D("Cell Density EC vs Layer").fill(i,celecden[j][i]); +// if (celecden[j][i] > 0 && celecden[j][i] < _mincd) +// { +// lodenshitmap.put(celecID[j][i],celechit[j][i]); +// if (trmipD) aida.cloud1D("Visible E LowD hits").fill(celechit[j][i].getRawEnergy()); +// } +// if (celecden[j][i] > _mincd) +// { +// hidenshitmap.put(celecID[j][i],celechit[j][i]); +// if (trmipD) aida.cloud1D("Visible E HighD hits").fill(celechit[j][i].getRawEnergy()); +// } + } + } +// List<CalorimeterHit> lodensHits = new Vector<CalorimeterHit>(); +// lodensHits.addAll(lodenshitmap.values()); +// event.put("LowDensityHits",lodensHits); +// List<CalorimeterHit> hidensHits = new Vector<CalorimeterHit>(); +// hidensHits.addAll(hidenshitmap.values()); +// event.put("HighDensityHits",hidensHits); + + // Loop over tracks in event, extrapolate to layers and match hits if mips + List<Track> evtracks = event.get(Track.class,_intname); +// aida.cloud1D("Number of Tracks").fill(evtracks.size()); + List<BasicCluster> mipclusters = new ArrayList<BasicCluster>(); + // also, need a map linking tracks to mip clusters for future use and map linking IL to track + Map<Track, BasicCluster> trkmipmap = new HashMap<Track, BasicCluster>(); + Map<Track, Integer> TrILmap = new HashMap<Track, Integer>(); + Map<BasicCluster, Integer> clusILmap = new HashMap<BasicCluster, Integer>(); + Map<Track, SpacePoint> trkposmap = new HashMap<Track, SpacePoint>(); + Map<Track, SpacePoint> trkmepmap = new HashMap<Track, SpacePoint>(); + // Make hitmaps that contain mip hits + HitMap embhitmap = (HitMap) (event.get("ChEBarrhitmap")); + HitMap hadbhitmap = (HitMap) (event.get("ChHBarrhitmap")); + HitMap emechitmap = (HitMap) (event.get("ChEEndcaphitmap")); + HitMap hadechitmap = (HitMap) (event.get("ChHEndcaphitmap")); + HitMap sembhitmap = (HitMap) (event.get("ScEBarrhitmap")); + HitMap shadbhitmap = (HitMap) (event.get("ScHBarrhitmap")); + HitMap semechitmap = (HitMap) (event.get("ScEEndcaphitmap")); + HitMap shadechitmap = (HitMap) (event.get("ScHEndcaphitmap")); + double TotTrP = 0; + int jt = 0; + int[] IntLay = new int[200]; +// System.out.println(" Starting Track Loop in TMC"); + SpacePoint trackspacept = new SpacePoint(); + SpacePoint lastmipSP = new SpacePoint(); + for (Track itrack : evtracks) + { + BasicCluster mipclus = new BasicCluster(); + double mipE = 0.; + double trpt = Math.sqrt(itrack.getPX()*itrack.getPX()+itrack.getPY()*itrack.getPY()); + double TrP = Math.sqrt(itrack.getPX()*itrack.getPX()+itrack.getPY()*itrack.getPY()+itrack.getPZ()*itrack.getPZ()); + double[] trp = itrack.getMomentum(); + if (TrP<0.75) + { + if (trp[2]>0) trp[2] = trp[2]-.002; // dEdx subtraction approx on z comp + if (trp[2]<0) trp[2] = trp[2]+.002; + } + Hep3Vector trp3 = new BasicHep3Vector(trp); + double[] trrp = itrack.getReferencePoint(); + double[] trpar = itrack.getTrackParameters(); + double[] tror = new double[3]; + tror[0] = -trpar[0]*Math.sin(trpar[1])+trrp[0]; + tror[1] = trpar[0]*Math.cos(trpar[1])+trrp[1]; + tror[2] = trpar[3]+trrp[2]; +// System.out.println(" Track Origin " + tror[0] + " " + tror[1] + " " + tror[2]); + Hep3Vector tror3 = new BasicHep3Vector(tror); + SpacePoint trorsp = new SpacePoint(tror3); // old swimmer was deprecated, uses SP now + int trq = itrack.getCharge(); +// System.out.println("Track Charge " + trq + " Track P " + trp3); + if (trq == 0) continue; // shouldn't need this, but it doesn't hurt + TotTrP += Math.sqrt(itrack.getPX()*itrack.getPX()+itrack.getPY()*itrack.getPY()+itrack.getPZ()*itrack.getPZ()); + if (trmipD) aida.cloud1D("pT of all Tracks").fill(trpt); + HelixSwimmer tswim = new HelixSwimmer(zField); + // swim track to 1st layer of CAL, arguments are momentum, origin, and charge +// tswim.setTrack(trp3, tror3, trq); // deprecated + tswim.setTrack(trp3, trorsp, trq); + double torad = tswim.getDistanceToRadius(BRadii[0]); // this is midpoint of first active layer of EM Barrel Cal + double toz = ECZs[0]; + if (trp[2]>0) toz = tswim.getDistanceToZ(ECZs[0]); // this is midpoint of first active layer of EM Endcap Cal + if (trp[2]<0) toz = tswim.getDistanceToZ(-ECZs[0]); // other endcap +// System.out.println(" Dist to BRad " + torad + " Dist to ECZ " + toz); + // if distance to radius is less than distance to z, use barrel hit collection, else use endcap + if (torad<Math.abs(toz)) // in Barrel + { + if (trmipD) aida.cloud1D("Track pT B Tracks").fill(trpt); + boolean nointeraction = true; + IntLay[jt] = _numbemlayers+_numbhadlayers; // set to last layer + // loop over all barrel layers, swimming track to radius and getting spacepoint + for (int j=0; j<_numbemlayers+_numbhadlayers; ++j) + { + if (nointeraction) + { + double prden = 0; + int ntcmatch = 0; + long[] layerid = new long[200]; + CalorimeterHit[] layerHit = new CalorimeterHit[200]; + double drad = tswim.getDistanceToRadius(BRadii[j]); + SpacePoint trackSP = tswim.getPointAtLength(drad); +// System.out.println("Barrel Layer " + j + "Barrel Radius " + BRadii[j] + "SpacePoint " + trackSP); + double trbtheta = Math.atan(trackSP.rxy()/trackSP.z()); + if (trbtheta<0) trbtheta+=Math.PI; + double trbphi = Math.atan2(trackSP.y(),trackSP.x()); + if (trbphi<0) trbphi+=2*Math.PI; + if (trmipD) aida.cloud1D("Track Theta at Brad").fill(trbtheta); + if (trmipD) aida.cloud1D("Track Phi at Brad").fill(trbphi); +// System.out.println("Track Theta " + trbtheta + "Track phi " + trbphi); + if (trmipD) aida.cloud2D("Theta vs Phi for tracks at Brad").fill(trbphi,trbtheta); + // now loop over all hits in layer, adding mips or determining IL + for (int i=0; i<nhitbLayer[j]; ++i) + { + if (trmipD) aida.cloud2D("Theta vs Phi for CAL B hits at Brad").fill(celbphi[i][j],celbtheta[i][j]); + if (trmipD) aida.cloud1D("Del theta track CAL B cell").fill(Math.abs(trbtheta-celbtheta[i][j])); + double dthe = Math.abs(trbtheta-celbtheta[i][j]); + double dphi = Math.abs(trbphi-celbphi[i][j]); + if (dphi > Math.PI) dphi = 2*Math.PI-dphi; + if (trmipD) aida.cloud1D("Del phi track CAL B cell").fill(dphi); + double delpos = Math.sqrt(dphi*dphi+dthe*dthe); + if (trmipD) aida.cloud1D("Del Pos Track CAL B cell").fill(delpos); + if (trmipD && j<_numbemlayers) aida.cloud1D("Del Pos Track CAL B EM cell").fill(delpos); + if (trmipD && j>_numbemlayers-1) aida.cloud1D("Del Pos Track CAL B HAD cell").fill(delpos); + if (trmipD) aida.cloud2D("Del Pos Track CAL B cell vs cell density").fill(celbden[i][j],delpos); + double trcedmin = _dminTrC; + if (j >= _numbemlayers) trcedmin = 1.5*trcedmin; + if (delpos<trcedmin) + { + prden += celbden[i][j]; + layerid[ntcmatch] = celbID[i][j]; + layerHit[ntcmatch] = celbhit[i][j]; + ntcmatch++; + } + } + if (ntcmatch>0) prden = prden/ntcmatch; + if (prden == 0 || prden>_mincd) + { + IntLay[jt] = j; + trackspacept = trackSP; +// System.out.println("Track SpacePoint at IL" + trackspacept); + if (trmipD) + { + if (j==0) aida.cloud1D("Track pt for B layer 0 IL").fill(trpt); + if (prden==0 && j==0) aida.cloud1D("No Hits in IL B layer 0").fill(prden); + if (prden>_mincd && j==0) aida.cloud1D("Many Hits in IL B layer 0").fill(prden); + } + nointeraction = false; + if (trmipD) aida.cloud1D("Cell Density Sum for B Interactions").fill(prden); + } else + { + lastmipSP = trackSP; + if (trmipD) aida.cloud1D("Cell Density Sum for B Noninteractions").fill(prden); + // add hits in this layer to mip cluster, remove from hitmap + for (int l=0; l<ntcmatch; ++l) + { +// System.out.println("Add this hit to B mip cluster " + l + " Layer " +j); +// System.out.println("Hit ID " + layerid[l] + " Hit " + layerHit[l]); + mipclus.addHit(layerHit[l]); + // increment hit energy using dEdx per layer + if (j<_numbemlayers) mipE += 0.027; // for em section + if (j>_numbemlayers-1) mipE += 0.054; // for had section + if (j<_numbemlayers) + { + embhitmap.remove(layerid[l]); + // also remove this hit from scin hitmap + sembhitmap.remove(layerid[l]); + } else + { + hadbhitmap.remove(layerid[l]); + shadbhitmap.remove(layerid[l]); + } + } + } + } + } + if (mipclus.getCalorimeterHits().size()==0) IntLay[jt]=0; + aida.cloud1D("Interaction Layer in B").fill(IntLay[jt]); + } else // Endcap + { + if (trmipD) aida.cloud1D("Track pT EC Tracks").fill(trpt); + boolean nointeraction = true; + IntLay[jt] = _numecemlayers+_numechadlayers; // set to last layer + // loop over all barrel layers, swimming track to radius and getting spacepoint + for (int j=0; j<_numecemlayers+_numechadlayers; ++j) + { + if (nointeraction) + { + double prden = 0; + int ntcmatch = 0; + long[] layerid = new long[200]; + CalorimeterHit[] layerHit = new CalorimeterHit[200]; + double drad = 0.; + if (trp[2]>0) drad = tswim.getDistanceToZ(ECZs[j]); + if (trp[2]<0) drad = tswim.getDistanceToZ(-ECZs[j]); + SpacePoint trackSP = tswim.getPointAtLength(drad); + double trectheta = Math.atan(trackSP.rxy()/trackSP.z()); + if (trectheta<0) trectheta+=Math.PI; + double trecphi = Math.atan2(trackSP.y(),trackSP.x()); + if (trecphi<0) trecphi+=2*Math.PI; + if (trmipD) aida.cloud1D("Track Theta at ECZ").fill(trectheta); + if (trmipD) aida.cloud1D("Track Phi at ECZ").fill(trecphi); + if (trmipD) aida.cloud2D("Theta vs Phi for tracks at ECZ").fill(trecphi,trectheta); + // have track theta, phi in Endcaps - check for mip hits in Endcap Collections + for (int i=0; i<nhitecLayer[j]; ++i) + { + if (trmipD) aida.cloud2D("Theta vs Phi for CAL EC hits at ECZ").fill(celecphi[i][j],celectheta[i][j]); + if (trmipD) aida.cloud1D("Del theta track CAL EC cell").fill(Math.abs(trectheta-celectheta[i][j])); + double dthe = Math.abs(trectheta-celectheta[i][j]); + double dphi = Math.abs(trecphi-celecphi[i][j]); + if (dphi > Math.PI) dphi = 2*Math.PI-dphi; + if (trmipD) aida.cloud1D("Del phi track CAL EC cell").fill(dphi); + double delpos = Math.sqrt(dphi*dphi+dthe*dthe); + if (trmipD) aida.cloud1D("Del Pos Track CAL EC cell").fill(delpos); + if (trmipD && j<_numecemlayers) aida.cloud1D("Del Pos Track CAL EC EM cell").fill(delpos); + if (trmipD && j>_numecemlayers-1) aida.cloud1D("Del Pos Track CAL EC HAD cell").fill(delpos); + double trcedmin = _dminTrC; + if (j >= _numecemlayers) trcedmin = 1.5*trcedmin; + if (delpos<trcedmin) + { + prden += celecden[i][j]; + layerid[ntcmatch] = celecID[i][j]; + layerHit[ntcmatch] = celechit[i][j]; + ntcmatch++; + } + } + if (ntcmatch>0) prden = prden/ntcmatch; + if (prden == 0 || prden>_mincd) + { + IntLay[jt] = j; + trackspacept = trackSP; + if (trmipD) + { + if (j==0) aida.cloud1D("Track pt for EC layer 0 IL").fill(trpt); + if (prden==0 && j==0) aida.cloud1D("No Hits in IL EC layer 0").fill(prden); + if (prden>_mincd && j==0) aida.cloud1D("Many Hits in IL EC layer 0").fill(prden); + } + nointeraction = false; + if (trmipD) aida.cloud1D("Cell Density Sum for EC Interactions").fill(prden); + } else + { + lastmipSP = trackSP; + if (trmipD) aida.cloud1D("Cell Density Sum for EC Noninteractions").fill(prden); + // add hits in this layer to mip cluster +// System.out.println("Number of matches " + ntcmatch); +// System.out.println("Number of hits in this layer " + nhitecLayer[j]); + for (int l=0; l<ntcmatch; ++l) + { +// System.out.println("Add this hit to EC mip cluster " + l + " Layer " +j); +// System.out.println("Hit ID " + layerid[l] + " Hit " + layerHit[l]); + mipclus.addHit(layerHit[l]); + // increment hit energy using dEdx per layer + if (j<_numecemlayers) mipE += 0.027; + if (j>_numecemlayers-1) mipE += 0.054; + if (j<_numecemlayers) + { + emechitmap.remove(layerid[l]); + semechitmap.remove(layerid[l]); + } else + { + hadechitmap.remove(layerid[l]); + shadechitmap.remove(layerid[l]); + } + } + } + } + } + if (mipclus.getCalorimeterHits().size()==0) IntLay[jt]=0; + aida.cloud1D("Interaction Layer in EC").fill(IntLay[jt]); + } + // add mip cluster to cluster list and fill maps + // first, set energy based on dEdx in various layers + mipclus.setEnergy(mipE); + mipclusters.add(mipclus); + trkmipmap.put(itrack, mipclus); + clusILmap.put(mipclus, IntLay[jt]); + TrILmap.put(itrack, IntLay[jt]); + trkposmap.put(itrack, trackspacept); + trkmepmap.put(itrack, lastmipSP); + jt++; // increment track index + } + if (trmipD) aida.cloud1D("Total Track Momentum").fill(TotTrP); + if (trmipD) aida.cloud1D("Number of Perfect Tracks").fill(evtracks.size()); + if (trmipD) aida.cloud1D("Number of mip clusters").fill(mipclusters.size()); + if (trmipD) aida.cloud2D("No of Tracks vs No of MipClus").fill(mipclusters.size(),evtracks.size()); + int strmiph = 0; + for (BasicCluster mipcluster : mipclusters) + { + if (evtracks.size()==1) strmiph = mipcluster.getCalorimeterHits().size(); + if (trmipD) aida.cloud1D("Number of hits in Mip Cluster").fill(mipcluster.getCalorimeterHits().size()); + } + for (int i=0; i<evtracks.size(); ++i) + { + aida.cloud1D("Interaction Layer all Tracks").fill(IntLay[i]); + if (evtracks.size()==1 && trmipD) aida.cloud1D("Interaction Layer 1 Track only").fill(IntLay[i]); + if (evtracks.size()==1 && trmipD) aida.cloud2D("NumHits MipCLus vs IL 1Tr").fill(IntLay[0],strmiph); +// System.out.println("Interaction Layer for Track " + IntLay[i]); + } + if (mipclusters.size() > 0) event.put(_nameExt, mipclusters); + event.put(_tmname,trkmipmap); + event.put("MipClusILMap",clusILmap); + event.put(_tilname,trkposmap); + event.put("TrackMipEPMap",trkmepmap); + event.put("TrackILMap",TrILmap); + } + public void setMipClusterListName(String ext) + { + _nameExt = ext; + } + public void setCollectionNames(String[] names) + { + _hitcollnames = names; + } + public void setTrackMipClusMap(String tmname) + { + _tmname = tmname; + } + public void setTrackILPosMap(String tilname) + { + _tilname = tilname; + } + public void setInputTrackList(String intname) + { + _intname = intname; + } +} +
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1