lcsim-contrib/src/main/java/org/lcsim/contrib/SteveMagill
diff -N GTrackMipClusterDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ GTrackMipClusterDriver.java 12 Aug 2009 18:30:25 -0000 1.1
@@ -0,0 +1,632 @@
+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.geometry.util.CalorimeterIDDecoder;
+import org.lcsim.util.hitmap.HitMap;
+
+public class GTrackMipClusterDriver 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 boolean trmipD = true;
+// private boolean trmipRes = false;
+// private SpacePoint trackspacept;
+ private double zField;
+// private String[] _hitmapnames;
+
+ public GTrackMipClusterDriver(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();
+ 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();
+ 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 " + BRadii[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 " +ECZs[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 " + BRadii[i+_numbemlayers]);
+ }
+
+ 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 " +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]);
+ CalorimeterIDDecoder caldecoder = (CalorimeterIDDecoder) event.getMetaData(chits).getIDDecoder();
+ 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, "PerfectTracks");
+// 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"));
+ 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 = 2.0*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.044; // for em section
+ if (j>_numbemlayers-1) mipE += 0.044; // for had section
+ if (j<_numbemlayers)
+ {
+ embhitmap.remove(layerid[l]);
+ } else
+ {
+ hadbhitmap.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 = 2.0*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.044;
+ if (j>_numecemlayers-1) mipE += 0.044;
+ if (j<_numecemlayers)
+ {
+ emechitmap.remove(layerid[l]);
+ } else
+ {
+ hadechitmap.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("TrackMipMap",trkmipmap);
+ event.put("MipClusILMap",clusILmap);
+ event.put("TrackILPosMap",trkposmap);
+ event.put("TrackMipEPMap",trkmepmap);
+ event.put("TrackILMap",TrILmap);
+ }
+ public void setClusterNameExtension(String ext)
+ {
+ _nameExt = ext;
+ }
+ public void setCollectionNames(String[] names)
+ {
+ _hitcollnames = names;
+ }
+}
+
