16 added + 3 modified, total 19 files
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -N SteveTrackMipClusterDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SteveTrackMipClusterDriver.java 5 Dec 2011 18:19:26 -0000 1.1
@@ -0,0 +1,632 @@
+package org.lcsim.contrib.Cassell.recon;
+
+// 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 SteveTrackMipClusterDriver 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 = false;
+// private boolean trmipRes = false;
+// private SpacePoint trackspacept;
+ private double zField;
+// private String[] _hitmapnames;
+
+ public SteveTrackMipClusterDriver(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 " +_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 " + 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 " +_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]);
+ IDDecoder caldecoder = event.getMetaData(chits).getIDDecoder();
+// System.out.println("Collection is " +_hitcollnames[i]);
+ if (_hitcollnames[i].equals("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("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("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("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, "Tracks");
+// 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("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];
+// 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 + " Origin " + tror3);
+ 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 = 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
+ // 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-10) mipE += 0.0055;
+ if (j>_numbemlayers-11 && j<_numbemlayers) mipE += 0.0111;
+ if (j>_numbemlayers-1) mipE += 0.0228;
+ 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-10) mipE += 0.0055;
+ if (j>_numecemlayers-11 && j<_numecemlayers) mipE += 0.0111;
+ if (j>_numecemlayers-1) mipE += 0.0228;
+ 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;
+ }
+}
+
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -N EflowCalibration.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ EflowCalibration.java 5 Dec 2011 18:19:26 -0000 1.1
@@ -0,0 +1 @@
+package org.lcsim.contrib.Cassell.recon;
import org.lcsim.recon.cheater.*;
import java.util.*;
import org.lcsim.geometry.Detector;
import org.lcsim.geometry.Subdetector;
import org.lcsim.util.Driver;
import org.lcsim.event.*;
import org.lcsim.geometry.IDDecoder;
import org.lcsim.util.aida.AIDA;
import hep.physics.vec.*;
import Jama.Matrix;
import org.lcsim.digisim.*;
import hep.aida.*;
import org.lcsim.recon.util.CalorimeterInformation;
import org.lcsim.geometry.Calorimeter.CalorimeterType;
import org.lcsim.recon.util.CalInfoDriver;
import org.lcsim.recon.pfa.identifier.*;
import org.lcsim.recon.cluster.mipfinder.trackxtrap.*;
import org.lcsim.recon.pfa.photonfinder.*;
import org.lcsim.recon.cluster.muonfinder.MuonFinderWrapper3;
import org.lcsim.util.hitmap.*;
import org.lcsim.recon.cluster.util.RemoveHitsFromClusters;
import org.lcsim.recon.cluster.mipfinder.ShowerPointFinderDriver2;
/**
* Calculate sampling fractions from fixed
E
* physics data minimizing event dE/sqrtE
* Ron[...]
\ No newline at end of file
[Note: Some over-long lines of diff output only partialy shown]
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -N FullEWWPrefitRecon.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ FullEWWPrefitRecon.java 5 Dec 2011 18:19:26 -0000 1.1
@@ -0,0 +1,320 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.Cassell.recon;
+import org.lcsim.contrib.Cassell.recon.analysis.*;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.BaseReconstructedParticle;
+import java.util.*;
+import hep.physics.vec.*;
+import hep.physics.jet.*;
+import org.lcsim.util.aida.AIDA;
+import java.io.*;
+
+/**
+ *
+ * @author cassell
+ */
+public class FullEWWPrefitRecon extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ String RPLname = "PandoraPFOCollection";
+ String pjname = "PerfectJets";
+ EventShape es;
+ String FourJetOutN;
+ String RecFourJetOutN = "PandoraPFOCollectionFullE4Jet";
+ FullE4JetFinder fjf;
+ int mx = 10000000;
+ String dsn = "unknown";
+ String FSRPname = "GenFS";
+ PrintStream out;
+ public FullEWWPrefitRecon()
+ {
+ add(new MakeRPlFromMCPl("GenFinalStateParticles",FSRPname));
+ add(new SetEventType2());
+ es = new EventShape();
+ FourJetOutN = FSRPname+"FullE4Jet";
+ fjf = new FullE4JetFinder(FSRPname);
+ add(fjf);
+ }
+
+ public void setDsn(String dsn)
+ {
+ this.dsn = dsn;
+ File file = new File(dsn);
+ try{
+ out = new PrintStream(file);
+ }
+ catch(IOException e){
+ System.out.println("IOException = "+e);
+ }
+ }
+
+ public void setRPLname(String RPLname)
+ {
+ this.RPLname = RPLname;
+ RecFourJetOutN = RPLname+"FullE4Jet";
+ }
+ protected void process(EventHeader event)
+ {
+ super.process(event);
+ Integer itype = (Integer) event.get(Integer.class,"EventType").get(0);
+ int etype = 0;
+ if( (itype > 99)&&(itype < 200) )etype = 1;
+ if( (itype > 199)&&(itype < 300) )etype = 2;
+ List<ReconstructedParticle> genrpl = event.get(ReconstructedParticle.class,FSRPname);
+ List<ReconstructedParticle> orpl = event.get(ReconstructedParticle.class,RPLname);
+ List<ReconstructedParticle> rp4jl = event.get(ReconstructedParticle.class,RecFourJetOutN);
+ List<ReconstructedParticle> prp4jl = event.get(ReconstructedParticle.class,pjname);
+ List<ReconstructedParticle> prp4jl0 = event.get(ReconstructedParticle.class,FourJetOutN);
+ List<ReconstructedParticle> p4jl;
+ double pevtE = 0.;
+ Map<HepLorentzVector,ReconstructedParticle> pltormap = new HashMap<HepLorentzVector,ReconstructedParticle>();
+ for(ReconstructedParticle rp:genrpl)
+ {
+ pltormap.put(new BasicHepLorentzVector(rp.getEnergy(),rp.getMomentum()),rp);
+ pevtE += rp.getEnergy();
+ }
+ es.setEvent(pltormap.keySet());
+ double pobl = es.oblateness();
+ Hep3Vector pthr = es.thrustAxis();
+ double pctthr = Math.abs(pthr.z())/pthr.magnitude();
+ double pM1 = 0.;
+ double pM2 = 0.;
+ p4jl = prp4jl;
+ if(itype == 0)p4jl = prp4jl0;
+ double[] pvals = new double[16];
+ if(p4jl.size() == 4)
+ {
+ double E = p4jl.get(0).getEnergy() + p4jl.get(1).getEnergy();
+ Hep3Vector P = VecOp.add(p4jl.get(0).getMomentum(),p4jl.get(1).getMomentum());
+ pM1 = Math.sqrt(E*E - P.magnitudeSquared());
+ E = p4jl.get(2).getEnergy() + p4jl.get(3).getEnergy();
+ P = VecOp.add(p4jl.get(2).getMomentum(),p4jl.get(3).getMomentum());
+ pM2 = Math.sqrt(E*E - P.magnitudeSquared());
+ for(int i=0;i<4;i++)
+ {
+ pvals[4*i] = p4jl.get(i).getMomentum().x();
+ pvals[4*i+1] = p4jl.get(i).getMomentum().y();
+ pvals[4*i+2] = p4jl.get(i).getMomentum().z();
+ pvals[4*i+3] = p4jl.get(i).getEnergy();
+ }
+ }
+ else
+ {
+ pM1 = -1.;
+ pM2 = -1.;
+ }
+ vars pvars = new vars();
+ pvars.setEevt(pevtE);
+ pvars.setObl(pobl);
+ pvars.setCtthr(pctthr);
+ pvars.setM1(pM1);
+ pvars.setM2(pM2);
+
+ double evtE = 0.;
+ Map<HepLorentzVector,ReconstructedParticle> ltormap = new HashMap<HepLorentzVector,ReconstructedParticle>();
+ for(ReconstructedParticle rp:orpl)
+ {
+ ltormap.put(new BasicHepLorentzVector(rp.getEnergy(),rp.getMomentum()),rp);
+ evtE += rp.getEnergy();
+ }
+ es.setEvent(ltormap.keySet());
+ double robl = es.oblateness();
+ Hep3Vector rthr = es.thrustAxis();
+ double rctthr = Math.abs(rthr.z())/rthr.magnitude();
+ double rM1 = 0.;
+ double rM2 = 0.;
+ double[] rvals = new double[16];
+ if(rp4jl.size() == 4)
+ {
+ List<ReconstructedParticle> crp4jl = new ArrayList<ReconstructedParticle>();
+ for(int i=0;i<4;i++)
+ {
+ crp4jl.add(correctBeta(rp4jl.get(i)));
+ }
+ double E = crp4jl.get(0).getEnergy() + crp4jl.get(1).getEnergy();
+ Hep3Vector P = VecOp.add(crp4jl.get(0).getMomentum(),crp4jl.get(1).getMomentum());
+ rM1 = Math.sqrt(E*E - P.magnitudeSquared());
+ E = crp4jl.get(2).getEnergy() + crp4jl.get(3).getEnergy();
+ P = VecOp.add(crp4jl.get(2).getMomentum(),crp4jl.get(3).getMomentum());
+ rM2 = Math.sqrt(E*E - P.magnitudeSquared());
+ for(int i=0;i<4;i++)
+ {
+ rvals[4*i] = crp4jl.get(i).getMomentum().x();
+ rvals[4*i+1] = crp4jl.get(i).getMomentum().y();
+ rvals[4*i+2] = crp4jl.get(i).getMomentum().z();
+ rvals[4*i+3] = crp4jl.get(i).getEnergy();
+ }
+ }
+ else
+ {
+ rM1 = -1.;
+ rM2 = -1.;
+ }
+ vars rvars = new vars();
+ rvars.setEevt(evtE);
+ rvars.setObl(robl);
+ rvars.setCtthr(rctthr);
+ rvars.setM1(rM1);
+ rvars.setM2(rM2);
+
+
+
+ // Write out the 4 vectors for fitting
+ out4V(etype,pvals,rvals,pvars,rvars);
+ }
+ private void out4V(int type,double[] pvals,double[] rvals,vars pvars,vars rvars)
+ {
+ for(int i=0;i<15;i++)
+ {
+ out.printf("%15.6g", pvals[i]);
+ }
+ out.printf("%15.6g\n", pvals[15]);
+ out.printf("%15.6g", (double) type);
+ out.printf("%15.6g", pvars.getEevt());
+ out.printf("%15.6g", pvars.getObl());
+ out.printf("%15.6g", pvars.getCtthr());
+ out.printf("%15.6g", pvars.getM1());
+ out.printf("%15.6g\n", pvars.getM2());
+ for(int i=0;i<15;i++)
+ {
+ out.printf("%15.6g", rvals[i]);
+ }
+ out.printf("%15.6g\n", rvals[15]);
+ out.printf("%15.6g", rvars.getEevt());
+ out.printf("%15.6g", rvars.getObl());
+ out.printf("%15.6g", rvars.getCtthr());
+ out.printf("%15.6g", rvars.getM1());
+ out.printf("%15.6g\n", rvars.getM2());
+ return;
+ }
+ protected void suspend()
+ {
+ out.close();
+ }
+ public ReconstructedParticle correctBeta(ReconstructedParticle jet)
+ {
+ double Ec = 0;
+ Hep3Vector Pc = new BasicHep3Vector();
+ double jetE = jet.getEnergy();
+ Hep3Vector jetP = jet.getMomentum();
+ for(ReconstructedParticle p:jet.getParticles())
+ {
+ double ecut = nhEcut(p.getEnergy());
+ if( (p.getCharge() == 0.)&&(p.getType() != 22) )
+ {
+ if(p.getEnergy() > ecut)
+ {
+ Ec += p.getEnergy();
+ Pc = VecOp.add(Pc,p.getMomentum());
+ }
+
+ }
+ else
+ {
+ Ec += p.getEnergy();
+ Pc = VecOp.add(Pc,p.getMomentum());
+ }
+ }
+ if(Ec > 10.)
+ {
+ double beta = Pc.magnitude()/Ec;
+ double pmagc = beta*jetE;
+ Pc = VecOp.mult(pmagc/Pc.magnitude(),Pc);
+ return new BaseReconstructedParticle(jetE,Pc);
+ }
+ return jet;
+ }
+ public double nhEcut(double E)
+ {
+ double[] Evals = {25.,75.,125.,175.,225.,275.,325.,375.,425.,475.,525.};
+ double[] nhEcut = {1.,2.,3.,4.,5.,6.,7.,10.,12.,6.,0.};
+ double cut = 0.;
+ if(E > Evals[Evals.length-1])return cut;
+ cut = Evals[0];
+ if(E < Evals[0])return cut;
+ int lb = 0;
+ while(E > Evals[lb+1])lb++;
+ cut = nhEcut[lb] + (nhEcut[lb+1]-nhEcut[lb])*(E-Evals[lb])/(Evals[lb+1]-Evals[lb]);
+ return cut;
+ }
+}
+class vars
+{
+ double Eevt;
+ double obl;
+ double ctthr;
+ double M1;
+ double M2;
+ double ct12 = -2.;
+ vars()
+ {
+
+ }
+
+ public double getEevt()
+ {
+ return Eevt;
+ }
+
+ public void setEevt(double Eevt)
+ {
+ this.Eevt = Eevt;
+ }
+
+ public double getM1()
+ {
+ return M1;
+ }
+
+ public void setM1(double M1)
+ {
+ this.M1 = M1;
+ }
+
+ public double getM2()
+ {
+ return M2;
+ }
+
+ public void setM2(double M2)
+ {
+ this.M2 = M2;
+ }
+
+ public double getCt12()
+ {
+ return ct12;
+ }
+
+ public void setCt12(double ct12)
+ {
+ this.ct12 = ct12;
+ }
+
+ public double getCtthr()
+ {
+ return ctthr;
+ }
+
+ public void setCtthr(double ctthr)
+ {
+ this.ctthr = ctthr;
+ }
+
+ public double getObl()
+ {
+ return obl;
+ }
+
+ public void setObl(double obl)
+ {
+ this.obl = obl;
+ }
+
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -N RemoveHighCtNh.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ RemoveHighCtNh.java 5 Dec 2011 18:19:26 -0000 1.1
@@ -0,0 +1,53 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.Cassell.recon;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.ReconstructedParticle;
+import java.util.*;
+
+/**
+ *
+ * @author cassell
+ */
+public class RemoveHighCtNh extends Driver
+{
+ String RPLname = "PandoraPFOCollection";
+ String outN = RPLname+"NHR";
+ double ctmax = 0.995;
+ public RemoveHighCtNh(String rpln)
+ {
+ RPLname = rpln;
+ outN = RPLname+"NHR";
+ }
+ public RemoveHighCtNh()
+ {
+ }
+
+ public void setCtmax(double ctmax)
+ {
+ this.ctmax = ctmax;
+ }
+
+ public void setRPLname(String RPLname)
+ {
+ this.RPLname = RPLname;
+ outN = RPLname+"NHR";
+ }
+ protected void process(EventHeader event)
+ {
+ List<ReconstructedParticle> outlist = new ArrayList<ReconstructedParticle>(event.get(ReconstructedParticle.class,RPLname));
+ for(ReconstructedParticle p:event.get(ReconstructedParticle.class,RPLname))
+ {
+ if(p.getType() == 2112)
+ {
+ double ct = Math.abs(p.getMomentum().z())/p.getMomentum().magnitude();
+ if(ct > ctmax)outlist.remove(p);
+ }
+ }
+ event.put(outN,outlist,ReconstructedParticle.class,0);
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -N RemoveHighCosThetaNeutralHadrons.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ RemoveHighCosThetaNeutralHadrons.java 5 Dec 2011 18:19:26 -0000 1.1
@@ -0,0 +1,55 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.Cassell.recon;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.ReconstructedParticle;
+import java.util.*;
+
+/**
+ *
+ * @author cassell
+ */
+public class RemoveHighCosThetaNeutralHadrons extends Driver
+{
+ String RPLname = "PandoraPFOCollection";
+ double ctmax = 0.995;
+ public RemoveHighCosThetaNeutralHadrons()
+ {
+ }
+ public RemoveHighCosThetaNeutralHadrons(String name)
+ {
+ RPLname = name;
+ }
+
+ public void setCtmax(double ctmax)
+ {
+ this.ctmax = ctmax;
+ }
+
+ public void setRPLname(String RPLname)
+ {
+ this.RPLname = RPLname;
+ }
+ protected void process(EventHeader event)
+ {
+ List<ReconstructedParticle> rl = event.get(ReconstructedParticle.class,RPLname);
+ List<ReconstructedParticle> rml = new ArrayList<ReconstructedParticle>();
+ for(ReconstructedParticle p:rl)
+ {
+ if(p.getType() == 2112)
+ {
+ double ct = Math.abs(p.getMomentum().z())/p.getMomentum().magnitude();
+ if(ct > ctmax)rml.add(p);
+ }
+ }
+ for(ReconstructedParticle p:rml)
+ {
+ rl.remove(p);
+ }
+
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -N ReconDriver10mmNonProjGap.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ReconDriver10mmNonProjGap.java 5 Dec 2011 18:19:26 -0000 1.1
@@ -0,0 +1,40 @@
+package org.lcsim.contrib.Cassell.recon;
+
+import org.lcsim.digisim.DigiPackageDriver;
+import org.lcsim.recon.pfa.output.FlushReconstructedParticlesDriver;
+import org.lcsim.recon.pfa.structural.SetUpPFA;
+import org.lcsim.recon.util.CalInfoDriver;
+import org.lcsim.util.Driver;
+
+/**
+ * Top-level driver to run UI PFA reconstruction.
+ *
+ * @author cassell
+ * @version $Id: ReconDriver10mmNonProjGap.java,v 1.1 2011/12/05 18:19:26 cassell Exp $
+ */
+
+public class ReconDriver10mmNonProjGap extends Driver
+{
+ /**
+ * Constructor that sets up daughter drivers.
+ */
+ public ReconDriver10mmNonProjGap()
+ {
+ // Cash general calorimeter information
+ add(new CalInfoDriver());
+
+ // Run digisim.
+ add(new DigiPackageDriver());
+
+ add(new RemoveHcalModuleNonProjBorderHits(10.));
+
+ // Run tracking.
+ add(new org.lcsim.recon.tracking.seedtracker.ReconTracking.SiD02ReconTrackingDriver());
+
+ // Set up and run PFA.
+ add(new SetUpPFA("Tracks"));
+
+ // Output collections.
+ add(new FlushReconstructedParticlesDriver("DTreeReclusteredParticles", "ReconstructedParticles", "Clusters"));
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -N RemoveSingleParticleJets.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ RemoveSingleParticleJets.java 5 Dec 2011 18:19:26 -0000 1.1
@@ -0,0 +1,35 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.Cassell.recon;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.ReconstructedParticle;
+import java.util.*;
+
+/**
+ *
+ * @author cassell
+ */
+public class RemoveSingleParticleJets extends Driver
+{
+ String RPLname;
+ String JLname;
+ public RemoveSingleParticleJets(String jname,String rname)
+ {
+ RPLname = rname;
+ JLname = jname;
+ }
+
+ protected void process(EventHeader event)
+ {
+ List<ReconstructedParticle> rl = event.get(ReconstructedParticle.class,RPLname);
+ List<ReconstructedParticle> jl = event.get(ReconstructedParticle.class,JLname);
+ for(ReconstructedParticle jet:jl)
+ {
+ if(jet.getParticles().size() == 1)rl.remove(jet.getParticles().get(0));
+ }
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -N ReconDriver10mmProjGap.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ReconDriver10mmProjGap.java 5 Dec 2011 18:19:26 -0000 1.1
@@ -0,0 +1,40 @@
+package org.lcsim.contrib.Cassell.recon;
+
+import org.lcsim.digisim.DigiPackageDriver;
+import org.lcsim.recon.pfa.output.FlushReconstructedParticlesDriver;
+import org.lcsim.recon.pfa.structural.SetUpPFA;
+import org.lcsim.recon.util.CalInfoDriver;
+import org.lcsim.util.Driver;
+
+/**
+ * Top-level driver to run UI PFA reconstruction.
+ *
+ * @author cassell
+ * @version $Id: ReconDriver10mmProjGap.java,v 1.1 2011/12/05 18:19:26 cassell Exp $
+ */
+
+public class ReconDriver10mmProjGap extends Driver
+{
+ /**
+ * Constructor that sets up daughter drivers.
+ */
+ public ReconDriver10mmProjGap()
+ {
+ // Cash general calorimeter information
+ add(new CalInfoDriver());
+
+ // Run digisim.
+ add(new DigiPackageDriver());
+
+ add(new RemoveHcalModuleProjBorderHits(10.));
+
+ // Run tracking.
+ add(new org.lcsim.recon.tracking.seedtracker.ReconTracking.SiD02ReconTrackingDriver());
+
+ // Set up and run PFA.
+ add(new SetUpPFA("Tracks"));
+
+ // Output collections.
+ add(new FlushReconstructedParticlesDriver("DTreeReclusteredParticles", "ReconstructedParticles", "Clusters"));
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -N MakeCMReconstructedParticles.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MakeCMReconstructedParticles.java 5 Dec 2011 18:19:26 -0000 1.1
@@ -0,0 +1,62 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.Cassell.recon;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.BaseReconstructedParticle;
+import java.util.*;
+import hep.physics.vec.*;
+
+/**
+ *
+ * @author cassell
+ */
+public class MakeCMReconstructedParticles extends Driver
+{
+ String RPLname = "PandoraPFOCollection";
+ String outN = RPLname+"CM";
+ public MakeCMReconstructedParticles(String rpln)
+ {
+ RPLname = rpln;
+ outN = RPLname+"CM";
+ }
+ public MakeCMReconstructedParticles()
+ {
+ }
+
+ public void setRPLname(String RPLname)
+ {
+ this.RPLname = RPLname;
+ outN = RPLname+"CM";
+ }
+ protected void process(EventHeader event)
+ {
+ List<ReconstructedParticle> outlist = new ArrayList<ReconstructedParticle>();
+ double evtE = 0.;
+ Hep3Vector evtP = new BasicHep3Vector();
+ List<ReconstructedParticle> rpl = event.get(ReconstructedParticle.class,RPLname);
+ if(rpl.size() < 2)
+ {
+ event.put(outN,outlist,ReconstructedParticle.class,0);
+ return;
+ }
+ for(ReconstructedParticle rp:rpl)
+ {
+ evtE += rp.getEnergy();
+ evtP = VecOp.add(evtP,rp.getMomentum());
+ }
+ HepLorentzVector evtLV = new BasicHepLorentzVector(evtE,evtP);
+ for(ReconstructedParticle rp:rpl)
+ {
+ HepLorentzVector bp = VecOp.boost(new BasicHepLorentzVector(rp.getEnergy(),rp.getMomentum()),evtLV);
+ BaseReconstructedParticle crp = new BaseReconstructedParticle(bp);
+ crp.addParticle(rp);
+ outlist.add(crp);
+ }
+ event.put(outN,outlist,ReconstructedParticle.class,0);
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -N TestQNeutralHadronClusterEnergyCalculator.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TestQNeutralHadronClusterEnergyCalculator.java 5 Dec 2011 18:19:26 -0000 1.1
@@ -0,0 +1,170 @@
+/*
+ * QNeutralHadronClusterEnergyCalculator.java
+ *
+ */
+package org.lcsim.contrib.Cassell.recon;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.conditions.*;
+import org.lcsim.recon.util.CalorimeterInformation;
+import org.lcsim.geometry.Calorimeter.CalorimeterType;
+
+/**
+ *
+ * @author cassell
+ */
+public class TestQNeutralHadronClusterEnergyCalculator implements ClusterEnergyCalculator
+{
+ protected ConditionsManager _mgr;
+
+ /** Creates a new instance of GenericClusterEnergyCalculator */
+ protected double alpha;
+ String calibrationFilename = "nhQcal-v2r3p10";
+ String analogCalibrationFilename = "nhQcalAnalog-v2r3p10";
+ protected boolean analogHcal = false;
+ protected double[] sf;
+ double[] Em;
+ double[] Ec;
+ protected int nFrontLayersEcal;
+ protected int eBid;
+ protected int eEid;
+ protected int eLid;
+ protected int hBid;
+ protected int hEid;
+ protected int mBid;
+ protected int mEid;
+ protected double EEcal;
+ protected double EHcal;
+ public TestQNeutralHadronClusterEnergyCalculator()
+ {
+ }
+ public TestQNeutralHadronClusterEnergyCalculator(String calFile)
+ {
+ calibrationFilename = calFile;
+ }
+ public TestQNeutralHadronClusterEnergyCalculator(boolean analog)
+ {
+ analogHcal = analog;
+ if(analog)calibrationFilename = analogCalibrationFilename;
+ }
+ public TestQNeutralHadronClusterEnergyCalculator(String calFile,boolean analog)
+ {
+ analogHcal = analog;
+ calibrationFilename = calFile;
+ }
+ protected void init()
+ {
+ _mgr = ConditionsManager.defaultInstance();
+ ConditionsSet cs = _mgr.getConditions("hadronCalibration/"+calibrationFilename);
+ sf = cs.getDoubleArray("SFs");
+ Em = cs.getDoubleArray("Emeas");
+ Ec = cs.getDoubleArray("Ecorr");
+ alpha = cs.getDouble("alpha");
+ nFrontLayersEcal = cs.getInt("nFrontLayersEcal");
+ CalorimeterInformation ci = CalorimeterInformation.instance();
+ eBid = ci.getSystemID(CalorimeterType.EM_BARREL);
+ eEid = ci.getSystemID(CalorimeterType.EM_ENDCAP);
+ eLid = ci.getSystemID(CalorimeterType.LUMI);
+ hBid = ci.getSystemID(CalorimeterType.HAD_BARREL);
+ hEid = ci.getSystemID(CalorimeterType.HAD_ENDCAP);
+ mBid = ci.getSystemID(CalorimeterType.MUON_BARREL);
+ mEid = ci.getSystemID(CalorimeterType.MUON_ENDCAP);
+ }
+ public double getEnergy(Cluster c)
+ {
+ if (_mgr == null)
+ {
+ init();
+ }
+ double EmeasEst = 0.;
+ EEcal = 0.;
+ EHcal = 0.;
+ for(CalorimeterHit hit:c.getCalorimeterHits())
+ {
+ int index = -1;
+ IDDecoder idc = hit.getIDDecoder();
+ idc.setID(hit.getCellID());
+ int detector_index = idc.getValue("system");
+ if(detector_index == eBid)
+ {
+ int layer = idc.getValue("layer");
+ if(layer < nFrontLayersEcal)
+ {
+ index = 0;
+ EmeasEst += hit.getRawEnergy()/sf[index];
+ EEcal += hit.getRawEnergy()/sf[index];
+ }
+ else
+ {
+ index = 1;
+ EmeasEst += hit.getRawEnergy()/sf[index];
+ EEcal += hit.getRawEnergy()/sf[index];
+ }
+ }
+ else if(detector_index == eEid)
+ {
+ int layer = idc.getValue("layer");
+ if(layer < nFrontLayersEcal)
+ {
+ index = 2;
+ EmeasEst += hit.getRawEnergy()/sf[index];
+ EEcal += hit.getRawEnergy()/sf[index];
+ }
+ else
+ {
+ index = 3;
+ EmeasEst += hit.getRawEnergy()/sf[index];
+ EEcal += hit.getRawEnergy()/sf[index];
+ }
+ }
+ else if(detector_index == hBid)
+ {
+ double Ehit = 1.;
+ if(analogHcal)Ehit = hit.getRawEnergy();
+ double[] pos = hit.getPosition();
+ double R = Math.sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
+ double ctheta = Math.abs(pos[2])/R;
+ double st = Math.sqrt(1.-ctheta*ctheta);
+ EmeasEst += Ehit/(1. + alpha*(1./st - 1.))/sf[4];
+ EHcal += Ehit/(1. + alpha*(1./st - 1.))/sf[4];
+ }
+ else if(detector_index == hEid)
+ {
+ double Ehit = 1.;
+ if(analogHcal)Ehit = hit.getRawEnergy();
+ double[] pos = hit.getPosition();
+ double R = Math.sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
+ double st = Math.abs(pos[2])/R;
+ EmeasEst += Ehit/(1. + alpha*(1./st - 1.))/sf[5];
+ EHcal += Ehit/(1. + alpha*(1./st - 1.))/sf[5];
+ }
+ }
+//
+// Now invert
+//
+ double Emeas = linext(EmeasEst);
+ return Emeas;
+ }
+ public double getEEcal(){return EEcal;}
+ public double getEHcal(){return EHcal;}
+ protected double linext(double Emeas)
+ {
+ double E;
+ if(Emeas <= Em[0])E = Emeas*Ec[0]/Em[0];
+ else if(Emeas >= Em[Em.length-1])E = Emeas*Ec[Em.length-1]/Em[Em.length-1];
+ else
+ {
+ int ib = 0;
+ for(int j=1;j<Em.length;j++)
+ {
+ if(Emeas < Em[j])break;
+ ib++;
+ }
+ E = Ec[ib] + (Emeas - Em[ib])*(Ec[ib+1] - Ec[ib])/(Em[ib+1] - Em[ib]);
+ }
+ return E;
+ }
+
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -N RemoveLowMultiplicityNeutralJets.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ RemoveLowMultiplicityNeutralJets.java 5 Dec 2011 18:19:26 -0000 1.1
@@ -0,0 +1,57 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.Cassell.recon;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.ReconstructedParticle;
+import java.util.*;
+
+/**
+ *
+ * @author cassell
+ */
+public class RemoveLowMultiplicityNeutralJets extends Driver
+{
+ String RPLname;
+ String JLname;
+ int npcut;
+ public RemoveLowMultiplicityNeutralJets(String jname,String rname,int cut)
+ {
+ RPLname = rname;
+ JLname = jname;
+ npcut = cut;
+ }
+
+ protected void process(EventHeader event)
+ {
+ List<ReconstructedParticle> rl = event.get(ReconstructedParticle.class,RPLname);
+ List<ReconstructedParticle> jl = event.get(ReconstructedParticle.class,JLname);
+ List<ReconstructedParticle> removel = new ArrayList<ReconstructedParticle>();
+ for(ReconstructedParticle jet:jl)
+ {
+ if(jet.getParticles().size() <= npcut)
+ {
+ int nch = 0;
+ for(ReconstructedParticle p:jet.getParticles())
+ {
+ if(p.getCharge() != 0.)nch++;
+ }
+ if(nch == 0)
+ {
+ removel.add(jet);
+ for(ReconstructedParticle p:jet.getParticles())
+ {
+ rl.remove(p);
+ }
+ }
+ }
+ }
+ for(ReconstructedParticle jet:removel)
+ {
+ jl.remove(jet);
+ }
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -N TwoJetFinder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TwoJetFinder.java 5 Dec 2011 18:19:26 -0000 1.1
@@ -0,0 +1,113 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.Cassell.recon;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.BaseReconstructedParticle;
+import java.util.*;
+import hep.physics.vec.*;
+import hep.physics.jet.*;
+
+/**
+ *
+ * @author cassell
+ */
+public class TwoJetFinder extends Driver
+{
+ String RPLname = "PandoraPFOCollection";
+ EventShape es;
+ String outN = RPLname+"TwoJet";
+ public TwoJetFinder(String rpln)
+ {
+ RPLname = rpln;
+ es = new EventShape();
+ outN = RPLname+"TwoJet";
+ }
+ public TwoJetFinder()
+ {
+ es = new EventShape();
+ }
+
+ public void setRPLname(String RPLname)
+ {
+ this.RPLname = RPLname;
+ outN = RPLname+"TwoJet";
+ }
+ protected void process(EventHeader event)
+ {
+ List<ReconstructedParticle> outlist = new ArrayList<ReconstructedParticle>();
+ List<ReconstructedParticle>[] jet = new List[2];
+ jet[0] = new ArrayList<ReconstructedParticle>();
+ jet[1] = new ArrayList<ReconstructedParticle>();
+ double[] jE = new double[2];
+ Hep3Vector[] jP = new Hep3Vector[2];
+ jP[0] = new BasicHep3Vector();
+ jP[1] = new BasicHep3Vector();
+ double evtE = 0.;
+ Hep3Vector evtP = new BasicHep3Vector();
+ List<ReconstructedParticle> rpl = event.get(ReconstructedParticle.class,RPLname);
+ if(rpl.size() < 2)
+ {
+ event.put(outN,outlist,ReconstructedParticle.class,0);
+ return;
+ }
+ Map<HepLorentzVector,ReconstructedParticle> ltormap = new HashMap<HepLorentzVector,ReconstructedParticle>();
+ for(ReconstructedParticle rp:rpl)
+ {
+ evtE += rp.getEnergy();
+ evtP = VecOp.add(evtP,rp.getMomentum());
+ }
+ HepLorentzVector evtLV = new BasicHepLorentzVector(evtE,evtP);
+ for(ReconstructedParticle rp:rpl)
+ {
+ HepLorentzVector bp = VecOp.boost(new BasicHepLorentzVector(rp.getEnergy(),rp.getMomentum()),evtLV);
+ ltormap.put(bp,rp);
+ }
+ // Find thrust axis, and divide particles with hemisphere cut
+ es.setEvent(ltormap.keySet());
+ Hep3Vector thrusta = es.thrustAxis();
+ for(HepLorentzVector bp:ltormap.keySet())
+ {
+ int id = 1;
+ if(VecOp.dot(thrusta,bp.v3()) > 0)id = 0;
+ ReconstructedParticle rp = ltormap.get(bp);
+ jet[id].add(rp);
+ jE[id] += rp.getEnergy();
+ jP[id] = VecOp.add(jP[id], rp.getMomentum());
+ }
+ // If Beta > 1. abort.
+ // Write empty list to event and return
+ if( (jP[0].magnitude()/jE[0] > 1.)||(jP[1].magnitude()/jE[1] > 1.) )
+ {
+ event.put(outN,outlist,ReconstructedParticle.class,0);
+ return;
+ }
+ // Put the higher energy jet first
+ if(jE[1] > jE[0])
+ {
+ double tE = jE[0];
+ Hep3Vector tP = jP[0];
+ List<ReconstructedParticle> trp = jet[0];
+ jE[0] = jE[1];
+ jP[0] = jP[1];
+ jet[0] = jet[1];
+ jE[1] = tE;
+ jP[1] = tP;
+ jet[1] = trp;
+ }
+ for(int i=0;i<2;i++)
+ {
+ BaseReconstructedParticle jetrp = new BaseReconstructedParticle(jE[i],jP[i]);
+ for(ReconstructedParticle p:jet[i])
+ {
+ jetrp.addParticle(p);
+ }
+ outlist.add(jetrp);
+ }
+ event.put(outN,outlist,ReconstructedParticle.class,0);
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -N ReconDriver5mmNonProjGap.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ReconDriver5mmNonProjGap.java 5 Dec 2011 18:19:26 -0000 1.1
@@ -0,0 +1,40 @@
+package org.lcsim.contrib.Cassell.recon;
+
+import org.lcsim.digisim.DigiPackageDriver;
+import org.lcsim.recon.pfa.output.FlushReconstructedParticlesDriver;
+import org.lcsim.recon.pfa.structural.SetUpPFA;
+import org.lcsim.recon.util.CalInfoDriver;
+import org.lcsim.util.Driver;
+
+/**
+ * Top-level driver to run UI PFA reconstruction.
+ *
+ * @author cassell
+ * @version $Id: ReconDriver5mmNonProjGap.java,v 1.1 2011/12/05 18:19:26 cassell Exp $
+ */
+
+public class ReconDriver5mmNonProjGap extends Driver
+{
+ /**
+ * Constructor that sets up daughter drivers.
+ */
+ public ReconDriver5mmNonProjGap()
+ {
+ // Cash general calorimeter information
+ add(new CalInfoDriver());
+
+ // Run digisim.
+ add(new DigiPackageDriver());
+
+ add(new RemoveHcalModuleNonProjBorderHits(5.));
+
+ // Run tracking.
+ add(new org.lcsim.recon.tracking.seedtracker.ReconTracking.SiD02ReconTrackingDriver());
+
+ // Set up and run PFA.
+ add(new SetUpPFA("Tracks"));
+
+ // Output collections.
+ add(new FlushReconstructedParticlesDriver("DTreeReclusteredParticles", "ReconstructedParticles", "Clusters"));
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -N EflowCalibrationFromSingleQ.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ EflowCalibrationFromSingleQ.java 5 Dec 2011 18:19:26 -0000 1.1
@@ -0,0 +1 @@
+package org.lcsim.contrib.Cassell.recon;
import org.lcsim.recon.cheater.*;
import java.util.*;
import org.lcsim.geometry.Detector;
import org.lcsim.geometry.Subdetector;
import org.lcsim.util.Driver;
import org.lcsim.event.*;
import org.lcsim.geometry.IDDecoder;
import org.lcsim.util.aida.AIDA;
import hep.physics.vec.*;
import Jama.Matrix;
import org.lcsim.digisim.*;
import hep.aida.*;
import org.lcsim.recon.util.CalorimeterInformation;
import org.lcsim.geometry.Calorimeter.CalorimeterType;
import org.lcsim.recon.util.CalInfoDriver;
import org.lcsim.recon.pfa.identifier.*;
import org.lcsim.recon.cluster.mipfinder.trackxtrap.*;
import org.lcsim.recon.pfa.photonfinder.*;
import org.lcsim.recon.cluster.muonfinder.MuonFinderWrapper3;
import org.lcsim.util.hitmap.*;
import org.lcsim.recon.cluster.util.RemoveHitsFromClusters;
import org.lcsim.recon.cluster.mipfinder.ShowerPointFinderDriver2;
/**
* Calculate sampling fractions from fixed
E
* physics data minimizing event dE/sqrtE
* Ron[...]
\ No newline at end of file
[Note: Some over-long lines of diff output only partialy shown]
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -N EflowCalibrationCT.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ EflowCalibrationCT.java 5 Dec 2011 18:19:26 -0000 1.1
@@ -0,0 +1 @@
+package org.lcsim.contrib.Cassell.recon;
import org.lcsim.recon.cheater.*;
import java.util.*;
import org.lcsim.geometry.Detector;
import org.lcsim.geometry.Subdetector;
import org.lcsim.util.Driver;
import org.lcsim.event.*;
import org.lcsim.geometry.IDDecoder;
import org.lcsim.util.aida.AIDA;
import hep.physics.vec.*;
import Jama.Matrix;
import org.lcsim.digisim.*;
import hep.aida.*;
import org.lcsim.recon.util.CalorimeterInformation;
import org.lcsim.geometry.Calorimeter.CalorimeterType;
import org.lcsim.recon.util.CalInfoDriver;
import org.lcsim.recon.pfa.identifier.*;
import org.lcsim.recon.cluster.mipfinder.trackxtrap.*;
import org.lcsim.recon.pfa.photonfinder.*;
import org.lcsim.recon.cluster.muonfinder.MuonFinderWrapper3;
import org.lcsim.util.hitmap.*;
import org.lcsim.recon.cluster.util.RemoveHitsFromClusters;
import org.lcsim.recon.cluster.mipfinder.ShowerPointFinderDriver2;
/**
* Calculate sampling fractions from fixed
E
* physics data minimizing event dE/sqrtE
* Ron[...]
\ No newline at end of file
[Note: Some over-long lines of diff output only partialy shown]
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -N ReconDriver5mmProjGap.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ReconDriver5mmProjGap.java 5 Dec 2011 18:19:26 -0000 1.1
@@ -0,0 +1,40 @@
+package org.lcsim.contrib.Cassell.recon;
+
+import org.lcsim.digisim.DigiPackageDriver;
+import org.lcsim.recon.pfa.output.FlushReconstructedParticlesDriver;
+import org.lcsim.recon.pfa.structural.SetUpPFA;
+import org.lcsim.recon.util.CalInfoDriver;
+import org.lcsim.util.Driver;
+
+/**
+ * Top-level driver to run UI PFA reconstruction.
+ *
+ * @author cassell
+ * @version $Id: ReconDriver5mmProjGap.java,v 1.1 2011/12/05 18:19:26 cassell Exp $
+ */
+
+public class ReconDriver5mmProjGap extends Driver
+{
+ /**
+ * Constructor that sets up daughter drivers.
+ */
+ public ReconDriver5mmProjGap()
+ {
+ // Cash general calorimeter information
+ add(new CalInfoDriver());
+
+ // Run digisim.
+ add(new DigiPackageDriver());
+
+ add(new RemoveHcalModuleProjBorderHits(5.));
+
+ // Run tracking.
+ add(new org.lcsim.recon.tracking.seedtracker.ReconTracking.SiD02ReconTrackingDriver());
+
+ // Set up and run PFA.
+ add(new SetUpPFA("Tracks"));
+
+ // Output collections.
+ add(new FlushReconstructedParticlesDriver("DTreeReclusteredParticles", "ReconstructedParticles", "Clusters"));
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -u -r1.1 -r1.2
--- RemoveHcalModuleNonProjBorderHits.java 16 Jun 2010 16:58:37 -0000 1.1
+++ RemoveHcalModuleNonProjBorderHits.java 5 Dec 2011 18:19:26 -0000 1.2
@@ -2,6 +2,8 @@
import java.util.*;
import org.lcsim.util.Driver;
import org.lcsim.event.EventHeader;
+import org.lcsim.event.Cluster;
+import org.lcsim.recon.cluster.util.BasicCluster;
import org.lcsim.event.CalorimeterHit;
import org.lcsim.geometry.IDDecoder;
import org.lcsim.recon.util.*;
@@ -50,11 +52,16 @@
if(Math.abs(d1) < prox)rm.add(h);
else if(Math.abs(d2) < prox)rm.add(h);
}
+ BasicCluster c = new BasicCluster();
for(CalorimeterHit h:rm)
{
bl.remove(h);
+ c.addHit(h);
}
-// event.put("RemovedHcalBarrelHits", rm);
+ event.put("RemovedHcalBarrelHits", rm);
+ List<Cluster> cl = new ArrayList<Cluster>();
+ cl.add(c);
+ event.put("RemovedHitsCluster", cl);
// System.out.println("Border removal: "+rm.size()+" hits out of "+nh+" removed from HCAL Barrel in evt "+ievt);
ievt++;
}
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -u -r1.3 -r1.4
--- FullE4JetFinder.java 16 Aug 2011 16:48:25 -0000 1.3
+++ FullE4JetFinder.java 5 Dec 2011 18:19:26 -0000 1.4
@@ -131,7 +131,7 @@
BaseReconstructedParticle[] rjet = new BaseReconstructedParticle[4];
int ij = 0;
double maxE = 0.;
- int maxij = -1;
+ int maxij = 0;
for(int i=0;i<2;i++)
{
for(int j=0;j<2;j++)
@@ -141,14 +141,10 @@
{
rjet[ij].addParticle(rp);
}
- if(jE[i][j] > maxE)
- {
- maxE = jE[i][j];
- maxij = ij;
- }
ij++;
}
}
+ if(djP[1].z()/djP[1].magnitude() > djP[0].z()/djP[0].magnitude())maxij = 2;
int[] order = new int[4];
order[0] = maxij;
order[1] = maxij+1;
lcsim-contrib/src/main/java/org/lcsim/contrib/Cassell/recon
diff -u -r1.1 -r1.2
--- DebugReconDriver.java 26 Oct 2010 20:20:14 -0000 1.1
+++ DebugReconDriver.java 5 Dec 2011 18:19:26 -0000 1.2
@@ -31,7 +31,7 @@
clusterername+">5hits"};
String[] rnames = {PPRPflowRname,ClLevelPFlowRname+"cut2",
ClLevelPFlowRname+"cut5"};
- private boolean useNewInitialMipFinding = false;
+ private boolean useNewInitialMipFinding = true;
SetUpPFA setup;
public void setUseNewInitialMipFinding(boolean x)
{
@@ -41,7 +41,7 @@
public DebugReconDriver()
{
ReconDriver rd = new ReconDriver();
- rd.setUseNewInitialMipFinding(true);
+ rd.setUseNewInitialMipFinding(useNewInitialMipFinding);
add(rd);
//
// Do the cheating reconstruction (includes Digisim)
@@ -78,5 +78,9 @@
//
add(new ClusterLevelCheater("E"+clusterername+"recl","H"+clusterername,REClNames,RHClNames,PPRPflowRname,ClLevelPFlowRname,
CheatReconFSname,CheatReconFSTrackedname));
+//
+// Add cheating at UI cluster level
+//
+ add(new UICLCFromRecon());
}
}
CVSspam 0.2.12
