lcsim/src/org/lcsim/contrib/SteveMagill
diff -N TrackMipClusterDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackMipClusterDriver.java 23 Apr 2007 20:09:53 -0000 1.1
@@ -0,0 +1,581 @@
+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.util.lcio.LCIOConstants;
+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 org.lcsim.recon.ztracking.cheater.*;
+import org.lcsim.mc.fast.tracking.*;
+import hep.aida.*;
+import org.lcsim.spacegeom.*;
+import org.lcsim.geometry.util.CalorimeterIDDecoder;
+import org.lcsim.util.hitmap.HitMap;
+
+public class TrackMipClusterDriver 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
+// The central magnetic field strength
+ private double[] _fieldStrength;
+ private double _dminE;
+ private double _dminH;
+ private double _dminTrC;
+ private double _mincd;
+ private String _nameExt;
+ private String[] _hitcollnames;
+ private boolean trmipD = false;
+ private boolean trmipRes = false;
+// private String[] _hitmapnames;
+
+ public TrackMipClusterDriver(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
+
+ boolean mcfltr = true;
+// mcfltr = (Boolean)event.get("MCFilt");
+// System.out.println("MCFilt is " +mcfltr);
+ if (mcfltr)
+ {
+
+// Initialize things here
+ if(!_initialized)
+ {
+ // setup specific detector stuff here
+ Detector det = event.getDetector();
+ double[] zero = {0, 0, 0};
+ _fieldStrength = det.getFieldMap().getField(zero);
+// System.out.println("B Field " +_fieldStrength[2]+ " Tesla");
+
+ // Setup layer counts and distances for barrel, endcap
+ CylindricalCalorimeter emb = ((CylindricalCalorimeter) det.getSubdetectors().get("EMBarrel"));
+ _embZ = emb.getZMin();
+ _numbemlayers = emb.getLayering().getLayerCount();
+// System.out.println("EM Barrel Layers " +_numbemlayers);
+ for (int i=0; i<_numbemlayers; ++i)
+ {
+ BRadii[i]=emb.getLayering().getDistanceToLayerSensorMid(i);
+// System.out.println("EM Barrel Layer Number " +i+ " EM Barrel Radius " +_emBRadii[i]);
+ }
+
+ CylindricalCalorimeter emec = ((CylindricalCalorimeter) det.getSubdetectors().get("EMEndcap"));
+ _numecemlayers = emec.getLayering().getLayerCount();
+// System.out.println("EM EC Layers " +_numecemlayers);
+ for (int i=0; i<_numecemlayers; ++i)
+ {
+ ECZs[i] = emec.getLayering().getDistanceToLayerSensorMid(i);
+// System.out.println("EM Endcap Layer Number " +i+ " EM Endcap Z " +_emecZ[i]);
+ }
+
+ CylindricalCalorimeter hadb = ((CylindricalCalorimeter) det.getSubdetectors().get("HADBarrel"));
+ _hadbZ = hadb.getZMin();
+ _numbhadlayers = hadb.getLayering().getLayerCount();
+// System.out.println("HAD Barrel Layers " +_numbhadlayers);
+ for (int i=0; i<_numbhadlayers; ++i)
+ {
+ BRadii[i+_numbemlayers]=hadb.getLayering().getDistanceToLayerSensorMid(i);
+// System.out.println("HAD Barrel Layer Number " +i+ " HAD Barrel Radius " +_hadBRadii[i]);
+ }
+
+ CylindricalCalorimeter hadec = ((CylindricalCalorimeter) det.getSubdetectors().get("HADEndcap"));
+ _numechadlayers = hadec.getLayering().getLayerCount();
+// System.out.println("HAD Endcap Layers " +_numechadlayers);
+ for (int i=0; i<_numechadlayers; ++i)
+ {
+ ECZs[i+_numecemlayers] = hadec.getLayering().getDistanceToLayerSensorMid(i);
+// System.out.println("HAD Endcap Layer Number " +i+ " HAD Endcap Z " +_hadecZ[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]);
+ CalorimeterIDDecoder caldecoder = (CalorimeterIDDecoder) event.getMetaData(chits).getIDDecoder();
+// System.out.println("Collection is " +_hitcollnames[i]);
+ if (_hitcollnames[i] == "EcalBarrDigiHits")
+ {
+ for (CalorimeterHit chit : chits)
+ {
+ long cID = chit.getCellID();
+ caldecoder.setID(cID);
+ int nlay = caldecoder.getLayer();
+ celbtheta[nhitbLayer[nlay]][nlay] = caldecoder.getTheta();
+ double cellphi = caldecoder.getPhi();
+ if (cellphi>Math.PI) cellphi=cellphi-2*Math.PI;
+ celbphi[nhitbLayer[nlay]][nlay] = cellphi;
+ celbID[nhitbLayer[nlay]][nlay] = cID;
+ celbhit[nhitbLayer[nlay]][nlay] = chit;
+ nhitbLayer[nlay]++;
+ nhitsEMB++;
+ nhitsBTotCAL++;
+ if (trmipD) aida.cloud1D("Raw Cell E in EM CAL").fill(chit.getRawEnergy());
+ if (trmipD) aida.cloud1D("Corrected Cell E in EM CAL").fill(chit.getCorrectedEnergy());
+ }
+ } else if (_hitcollnames[i] == "HcalBarrDigiHits")
+ {
+ for (CalorimeterHit chit : chits)
+ {
+ long cID = chit.getCellID();
+ caldecoder.setID(cID);
+ int nlay = _numbemlayers + caldecoder.getLayer();
+ celbtheta[nhitbLayer[nlay]][nlay] = caldecoder.getTheta();
+ double cellphi = caldecoder.getPhi();
+ if (cellphi>Math.PI) cellphi=cellphi-2*Math.PI;
+ celbphi[nhitbLayer[nlay]][nlay] = cellphi;
+ celbID[nhitbLayer[nlay]][nlay] = cID;
+ celbhit[nhitbLayer[nlay]][nlay] = chit;
+ nhitbLayer[nlay]++;
+ nhitsHADB++;
+ nhitsBTotCAL++;
+ }
+ } else if (_hitcollnames[i] == "EcalEndcapDigiHits")
+ {
+ for (CalorimeterHit chit : chits)
+ {
+ long cID = chit.getCellID();
+ caldecoder.setID(cID);
+ int nlay = caldecoder.getLayer();
+ celectheta[nhitecLayer[nlay]][nlay] = caldecoder.getTheta();
+ double cellphi = caldecoder.getPhi();
+ if (cellphi>Math.PI) cellphi=cellphi-2*Math.PI;
+ celecphi[nhitecLayer[nlay]][nlay] = cellphi;
+ celecID[nhitecLayer[nlay]][nlay] = cID;
+ celechit[nhitecLayer[nlay]][nlay] = chit;
+ nhitecLayer[nlay]++;
+ nhitsEMEC++;
+ nhitsECTotCAL++;
+ }
+ } else if (_hitcollnames[i] == "HcalEndcapDigiHits")
+ {
+ for (CalorimeterHit chit : chits)
+ {
+ long cID = chit.getCellID();
+ caldecoder.setID(cID);
+ int nlay = _numecemlayers + caldecoder.getLayer();
+ celectheta[nhitecLayer[nlay]][nlay] = caldecoder.getTheta();
+ double cellphi = caldecoder.getPhi();
+ if (cellphi>Math.PI) cellphi=cellphi-2*Math.PI;
+ 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
+// Map<Long, CalorimeterHit> lodenshitmap = new HashMap<Long, CalorimeterHit>();
+// Map<Long, CalorimeterHit> hidenshitmap = new HashMap<Long, CalorimeterHit>();
+ 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 < 30) celbden[j][i] += 0.0026/distance;
+ if (i > 29) celbden[j][i] += 0.0052/distance;
+ if (trmipD)
+ {
+ if (i == 0) aida.cloud1D("Distance between EM B CAL hits layer 0").fill(distance);
+ if (i < 30) aida.cloud1D("Distance between EM B CAL hits").fill(distance);
+ if (i == 30) aida.cloud1D("Distance between HAD B CAL hits layer 0").fill(distance);
+ if (i > 29) aida.cloud1D("Distance between HAD B CAL hits").fill(distance);
+ }
+ }
+ }
+ if (trmipD) aida.cloud1D("Cell Density B").fill(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 = 2*_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 = 2*_dminH;
+ if (distance < dmin)
+ {
+ if (i < 30) celecden[j][i] += 0.0017/distance;
+ if (i > 29) celecden[j][i] += 0.0042/distance;
+ if (trmipD)
+ {
+ if (i == 0) aida.cloud1D("Distance between EM EC CAL hits layer 0").fill(distance);
+ if (i < 30) aida.cloud1D("Distance between EM EC CAL hits").fill(distance);
+ if (i > 29) aida.cloud1D("Distance between HAD EC CAL hits").fill(distance);
+ if (i == 30) aida.cloud1D("Distance between HAD EC CAL hits layer 0").fill(distance);
+ }
+ }
+ }
+ if (trmipD) aida.cloud1D("Cell Density EC").fill(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);
+ if (trmipD) aida.cloud2D("LowD vs HighD hits").fill(hidensHits.size(),lodensHits.size());
+
+ // Loop over tracks in event, extrapolate to layers and match hits if mips
+// List<CheatTrack> evtracks = event.get(CheatTrack.class, "CombinedTracks");
+ List<PerfectTrack> evtracks = event.get(PerfectTrack.class, "PerfectTracks");
+// aida.cloud1D("Number of Combined Tracks").fill(evtracks.size());
+// System.out.println("Number of CheatTracks " +evtracks.size());
+ List<Track> ctracks = new ArrayList<Track>();
+ 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<PerfectTrack, BasicCluster> trkmipmap = new HashMap<PerfectTrack, BasicCluster>();
+ Map<BasicCluster, Integer> clusILmap = new HashMap<BasicCluster, Integer>();
+ // Make hitmaps that contain mip hits
+ HitMap embhitmap = (HitMap) (event.get("EMBarrhitmap"));
+ HitMap hadbhitmap = (HitMap) (event.get("HADBarrhitmap"));
+ HitMap emechitmap = (HitMap) (event.get("EMEndcaphitmap"));
+ HitMap hadechitmap = (HitMap) (event.get("HADEndcaphitmap"));
+ double TotTrP = 0;
+ int jt = 0;
+ int[] IntLay = new int[200];
+ for (PerfectTrack itrack : evtracks)
+ {
+// IntLay[jt] = 0; // set to first layer as default if wanted
+ BasicCluster mipclus = new BasicCluster();
+ double trpt = Math.sqrt(itrack.getPX()*itrack.getPX()+itrack.getPY()*itrack.getPY());
+ double[] trp = itrack.getMomentum();
+ Hep3Vector trp3 = new BasicHep3Vector(trp);
+// double[] tror = itrack.getReferencePoint();
+// Hep3Vector tror3 = new BasicHep3Vector(tror);
+ Hep3Vector tror3 = itrack.getOrigin();
+ int trq = itrack.getCharge();
+// System.out.println("Track Charge " + trq + " Track P " + trp3 + " Origin " + tror3);
+ if (trq == 0) continue; // sometimes these are photons!
+ ctracks.add(itrack);
+ // ctracks.set(jt, itrack);
+ TotTrP += Math.sqrt(itrack.getPX()*itrack.getPX()+itrack.getPY()*itrack.getPY()+itrack.getPZ()*itrack.getPZ());
+ if (trmipD) aida.cloud1D("pT of all Combined Tracks").fill(trpt);
+ HelixSwimmer tswim = new HelixSwimmer(_fieldStrength[2]);
+ // swim track to 1st layer of CAL, arguments are momentum, origin, and charge
+ tswim.setTrack(trp3, tror3, trq);
+// System.out.println("Radius of first EM Barrel layer" + BRadii[0]);
+ double torad = tswim.getDistanceToRadius(BRadii[0]); // this is midpoint of first active layer of EM Barrel Cal
+ double toz = tswim.getDistanceToZ(ECZs[0]); // this is midpoint of first active layer of EM Endcap Cal
+// System.out.println("Track pT " +trpt);
+// System.out.println("Track momentum " +trp3);
+// System.out.println("Track origin " +tror3);
+// System.out.println("Rad Dist " +torad+ "Z Dist " +toz);
+ // if distance to radius is less than distance to z, use barrel hit collection, else use endcap
+ if (torad < toz) // in Barrel
+ {
+ if (trmipD) aida.cloud1D("Track pT Barrel 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.getPointAtDistance(drad);
+// System.out.println("Barrel Layer " + j + "Barrel Radius " + BRadii[j] + "SpacePoint " + trackSP);
+ double trbtheta = trackSP.theta();
+ if (trmipD) aida.cloud1D("Track Theta at EMBrad").fill(trbtheta);
+ double trbphi = trackSP.phi();
+ if (trmipD) aida.cloud1D("Track Phi at EMBrad").fill(trbphi);
+// System.out.println("Track Theta " + trbtheta + "Track phi " + trbphi);
+ if (trmipD) aida.cloud2D("Theta vs Phi for tracks at EMBrad").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 EMBrad").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.cloud2D("Del theta vs Del phi track CAL B cell").fill(dphi,dthe);
+ if (trmipD) aida.cloud1D("Del Pos Track CAL B cell").fill(delpos);
+ // if (trpt>5 && j<30 && delpos<0.0075) aida.cloud1D("Cell Density for close B cells").fill(celbden[i][j]);
+ // if (trpt>5 && j<30 && delpos>0.0075) aida.cloud1D("Cell Density for far B cells").fill(celbden[i][j]);
+ // if (trpt>5 && j>29 && delpos<0.025) aida.cloud1D("Cell Density for close B cells").fill(celbden[i][j]);
+ // if (trpt>5 && j>29 && delpos>0.025) aida.cloud1D("Cell Density for far B cells").fill(celbden[i][j]);
+ double trcedmin = _dminTrC;
+ if (j >= _numbemlayers) trcedmin = 2*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;
+ nointeraction = false;
+ if (trmipD) aida.cloud1D("Cell Density Sum for B Interactions").fill(prden);
+ } else
+ {
+ 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]);
+ if (j<_numbemlayers)
+ {
+ embhitmap.remove(layerid[l]);
+ } else
+ {
+ hadbhitmap.remove(layerid[l]);
+ }
+ }
+ }
+// if (trpt>5) System.out.println("Layer number " + j + " Number of matches " + ntcmatch + " Match density " + prden);
+ }
+ }
+ if (mipclus.getCalorimeterHits().size()==0) IntLay[jt]=0;
+ if (trmipD) aida.cloud1D("Interaction Layer in B").fill(IntLay[jt]);
+ } else // Endcap
+ {
+ if (trmipD) aida.cloud1D("Track pT Endcap 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 = tswim.getDistanceToZ(ECZs[j]);
+ SpacePoint trackSP = tswim.getPointAtDistance(drad);
+// System.out.println("Spacepoint in EMEndcap Layer 0 " +trackSP);
+ double trectheta = trackSP.theta();
+ if (trmipD) aida.cloud1D("Track Theta at ECZ").fill(trectheta);
+ double trecphi = trackSP.phi();
+ 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.cloud2D("Del theta vs Del phi track CAL EC cell").fill(dphi,dthe);
+ if (trmipD) aida.cloud1D("Del Pos Track CAL EC cell").fill(delpos);
+ double trcedmin = 2*_dminTrC;
+ if (j >= _numecemlayers) trcedmin = 2*trcedmin;
+ if (delpos<2*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;
+ nointeraction = false;
+ if (trmipD) aida.cloud1D("Cell Density Sum for EC Interactions").fill(prden);
+ } else
+ {
+ 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]);
+ if (j<_numecemlayers)
+ {
+ emechitmap.remove(layerid[l]);
+ } else
+ {
+ hadechitmap.remove(layerid[l]);
+ }
+ }
+ }
+ }
+ }
+ if (mipclus.getCalorimeterHits().size()==0) IntLay[jt]=0;
+ if (trmipD) aida.cloud1D("Interaction Layer in EC").fill(IntLay[jt]);
+ }
+ // add mip cluster to cluster list and fill maps
+ mipclusters.add(mipclus);
+ trkmipmap.put(itrack, mipclus);
+ clusILmap.put(mipclus, IntLay[jt]);
+
+ jt++; // increment track index
+ }
+ if (trmipD) aida.cloud1D("Total Track Momentum").fill(TotTrP);
+ if (trmipRes) aida.cloud1D("Number of Perfect Tracks").fill(evtracks.size());
+ if (trmipRes) aida.cloud1D("Number of mip clusters").fill(mipclusters.size());
+ if (trmipRes) aida.cloud2D("No of Tracks vs No of MipClus").fill(mipclusters.size(),evtracks.size());
+ for (BasicCluster mipcluster : mipclusters)
+ {
+ if (trmipRes) aida.cloud1D("Number of hits in Mip Cluster").fill(mipcluster.getCalorimeterHits().size());
+ }
+ for (int i=0; i<evtracks.size(); ++i)
+ {
+ if (trmipRes) aida.cloud1D("Interaction Layer all Tracks").fill(IntLay[i]);
+// System.out.println("Interaction Layer for Track " + IntLay[i]);
+ }
+ if (ctracks.size() > 0) event.put("Cheated Tracks", ctracks); // cheated tracks now in event
+ if (mipclusters.size() > 0) event.put(_nameExt, mipclusters);
+ event.put("TrackMipMap",trkmipmap);
+ event.put("MipClusILMap",clusILmap);
+// System.out.println("EMB mip hitmap size " + embhitmap.size());
+// System.out.println("HADB mip hitmap size " + hadbhitmap.size());
+// System.out.println("EMEC mip hitmap size " + emechitmap.size());
+// System.out.println("HADEC mip hitmap size " + hadechitmap.size());
+
+ } // end of filter loop
+ }
+ public void setClusterNameExtension(String ext)
+ {
+ _nameExt = ext;
+ }
+ public void setCollectionNames(String[] names)
+ {
+ _hitcollnames = names;
+ }
+}
+
