| Commit in lcsim/src/org/lcsim/recon/pfa on MAIN | |||
| photonfinder/AssociatePhotonFragmentsSid01.java | +124 | added 1.1 | |
| /HMatrixVars.java | +153 | added 1.1 | |
| /IdentifyPhotonClustersSid01.java | +327 | added 1.1 | |
| /RonPhotonFinderSid01.java | +55 | added 1.1 | |
| /Sid01PhotonFinder.java | +35 | added 1.1 | |
| structural/NonTrivialPFA.java | +3 | -3 | 1.3 -> 1.4 |
| /RunAndWriteOutPFA.java | +1 | -1 | 1.7 -> 1.8 |
| +698 | -4 | ||
JM: move Ron's photon finder out of contrib
diff -N AssociatePhotonFragmentsSid01.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ AssociatePhotonFragmentsSid01.java 12 Dec 2008 21:49:52 -0000 1.1 @@ -0,0 +1,124 @@
+package org.lcsim.recon.pfa.photonfinder;
+
+import java.util.*;
+import org.lcsim.event.Cluster;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import hep.physics.vec.*;
+
+/**
+ * Photon identifier
+ *
+ * Code by Ron Cassell, put into CVS in a temporary place by Mat Charles.
+ */
+public class AssociatePhotonFragmentsSid01
+{
+ double[] bv = {-.98,-.94,-.88,-.68,-.32,.02,.28,.48,.62,.72,.78,.82,.86,.90,.94,.98};
+ double[] dcut1 = {41.25,41.25,38.75,38.75,41.25,43.75,46.25,48.75,
+ 51.25,53.75,56.25,58.75,61.25,63.75,68.75,78.75};
+ double[] dcut2 = {38.75,38.75,36.25,36.25,38.75,41.25,43.75,46.25,
+ 48.75,51.25,53.75,56.25,58.75,61.25,66.25,76.25};
+ double[] dcut3 = {36.25,36.25,33.75,31.25,33.75,36.25,36.25,41.25,
+ 43.75,46.25,48.75,51.25,53.75,56.25,63.75,73.75};
+ double[] dcut4 = {36.25,33.75,31.25,28.75,28.75,31.25,33.75,36.25,
+ 38.75,41.25,43.75,46.25,48.75,51.25,58.75,71.25};
+ double[] dcut5 = {38.75,28.75,28.75,26.25,23.75,28.75,28.75,28.75,
+ 33.75,36.25,38.75,41.25,43.75,48.75,53.75,66.25};
+ double[] dcut6 = {33.75,28.75,26.25,23.75,23.75,23.75,26.25,26.25,
+ 28.75,33.75,36.25,38.75,38.75,46.25,46.25,63.75};
+ double[] dcut10 = {33.75,28.75,26.25,21.25,21.25,18.75,18.75,21.25,
+ 26.25,28.75,26.25,31.25,31.25,38.75,38.75,56.25};
+ double[][] dcut = new double[7][bv.length];
+ double gdcut = 90.;
+ public AssociatePhotonFragmentsSid01()
+ {
+ for(int i=0;i<bv.length;i++)
+ {
+ dcut[0][i] = dcut1[i];
+ dcut[1][i] = dcut2[i];
+ dcut[2][i] = dcut3[i];
+ dcut[3][i] = dcut4[i];
+ dcut[4][i] = dcut5[i];
+ dcut[5][i] = dcut6[i];
+ dcut[6][i] = dcut10[i];
+ }
+ }
+ public List<Cluster> associateFragments(List<Cluster> photons,List<Cluster> fragments)
+ {
+//
+// Map photon core clusters to new clusters,used to add fragments
+//
+ Map<Cluster,BasicCluster> phmap = new HashMap<Cluster,BasicCluster>();
+ for(Cluster c:photons)
+ {
+ BasicCluster cc = new BasicCluster();
+ cc.addCluster(c);
+ phmap.put(c,cc);
+ }
+//
+// Loop over nonphoton clusters, find nearest photon cluster, and see if we should associate it
+//
+ for(Cluster c:fragments)
+ {
+ Cluster nc = null;
+ double mind = 99999.;
+ int nh = c.getSize();
+ Hep3Vector p1 = new BasicHep3Vector(c.getPosition());
+ Hep3Vector p3 = null;
+ int ib = nh-1;
+ if(nh > 10)
+ {
+ ib = 6;
+ }
+ else if(nh > 5)
+ {
+ ib = 5;
+ }
+ for(Cluster cc:photons)
+ {
+ if(cc.getSize() >= nh)
+ {
+ Hep3Vector p2 = new BasicHep3Vector(cc.getPosition());
+ double d = Math.sqrt( (p1.x() - p2.x())*(p1.x() - p2.x()) + (p1.y() - p2.y())*(p1.y() - p2.y()) + (p1.z() - p2.z())*(p1.z() - p2.z()) );
+ if(d < mind)
+ {
+ nc = cc;
+ mind = d;
+ p3 = p2;
+ }
+ }
+ }
+ if(nc != null)
+ {
+ Hep3Vector p4 = new BasicHep3Vector(p1.x()-p3.x(),p1.y()-p3.y(),p1.z()-p3.z());
+ double ct = (p1.x()*p4.x()+p1.y()*p4.y()+p1.z()*p4.z())/(p1.magnitude()*p4.magnitude());
+ if(mind < gdcut)
+ {
+ int ctbin = 0;
+ for(int jb=0;jb<bv.length;jb++)
+ {
+ if(ct < bv[jb])break;
+ ctbin++;
+ }
+ double distcut = 0.;
+ if(ctbin == 0)
+ {
+ distcut = dcut[ib][1] + (dcut[ib][0] - dcut[ib][1])*(ct - bv[1])/(bv[0] - bv[1]);
+ }
+ else if(ctbin == bv.length)
+ {
+ distcut = dcut[ib][bv.length-1] + (dcut[ib][bv.length-2] - dcut[ib][bv.length-1])*(ct - bv[bv.length-1])/(bv[bv.length-2] - bv[bv.length-1]);
+ }
+ else
+ {
+ distcut = dcut[ib][ctbin] + (dcut[ib][ctbin-1] - dcut[ib][ctbin])*(ct - bv[ctbin])/(bv[ctbin-1] - bv[ctbin]);
+ }
+ if(mind < distcut)
+ {
+ phmap.get(nc).addCluster(c);
+ }
+ }
+ }
+ }
+ return new ArrayList(phmap.values());
+ }
+}
diff -N HMatrixVars.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ HMatrixVars.java 12 Dec 2008 21:49:52 -0000 1.1 @@ -0,0 +1,153 @@
+package org.lcsim.recon.pfa.photonfinder;
+
+import org.lcsim.recon.cluster.util.*;
+import java.util.*;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.math.chisq.ChisqProb;
+import org.lcsim.recon.emid.hmatrix.HMatrix;
+import org.lcsim.recon.emid.hmatrix.HMatrixConditionsConverter;
+import org.lcsim.recon.cluster.nn.NearestNeighborClusterer;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+
+/**
+ * HMatrixPhotonIdentifier
+ *
+ * Code by Ron Cassell, put into CVS in a temporary place by Mat Charles.
+ */
+public class HMatrixVars extends Driver
+{
+ private boolean _initialized;
+ private boolean usecore;
+ private Clusterer nn111clusterer;
+ private int _nLayers;
+ // the number of variables in the measurement vector
+ private int _nmeas;
+ // we handle the energy dependence by correlating to the log of the energy
+ private static final double _log10inv = 1./Math.log(10.0);
+ // the index of the cluster energy (log base 10)
+ private int _logEIndex;
+ // the vector of measured values
+ double[] _vals;
+ // the H-matrix itself
+ private HMatrix _hmx = null;
+ private double chisq;
+ private double chisqD;
+ private double chisqDProb;
+ public HMatrixVars()
+ {
+ this(false);
+ }
+ public HMatrixVars(boolean core)
+ {
+ _initialized = false;
+ getConditionsManager().registerConditionsConverter(new HMatrixConditionsConverter());
+ usecore = core;
+ if(usecore)
+ {
+ nn111clusterer = new NearestNeighborClusterer(1,1,1,0,0.);
+ }
+ }
+ public void process(EventHeader event)
+ {
+ String[] det = {"EMBarrel"};
+ if(!_initialized)
+ {
+ CylindricalCalorimeter calsub = (CylindricalCalorimeter)event.getDetector().getSubdetectors().get(det[0]);
+ _nLayers = calsub.getLayering().getLayerCount();
+ // the vector of measurements starts as the longitudinal layers
+ _nmeas = _nLayers;
+ // add the logarithm of the energy
+ _logEIndex = _nmeas;
+ _nmeas+=1;
+ // would add any additional measurements (e.g. width) here
+ _vals = new double[_nmeas];
+
+ _hmx = getConditionsManager().getCachedConditions(HMatrix.class, "LongitudinalHMatrix.hmx").getCachedData();
+
+ _initialized = true;
+ }
+ }
+ public void setCluster(Cluster cl)
+ {
+ if(_hmx == null)
+ {
+ System.out.println("HMatrixPhotonIdentifier:: HMatrix not set!!");
+ return;
+ }
+ Cluster cluster = cl;
+ if(usecore)
+ {
+ List<Cluster> nn111cl = nn111clusterer.createClusters(cluster.getCalorimeterHits());
+ double maxE = 0.;
+ for(Cluster c:nn111cl)
+ {
+ if(c.getEnergy() > maxE)
+ {
+ maxE = c.getEnergy();
+ cluster = c;
+ }
+ }
+ }
+ // get the longitudinal shower layer energies
+ // should be able to fetch this from the cluster...
+ double[] layerE = layerEnergies(cluster);
+ // accumulate the longitudinal energy fractions...
+ // cng 050820 chi-squared values reported as being too high for late showering showers.
+ // modify to shift vector of values to start with showering layer.
+ int firstNonZeroLayer = 0;
+ boolean showered = false;
+ for(int i=0; i<layerE.length; ++i)
+ {
+ _vals[i] = layerE[i];
+ if(_vals[i]!=0.) showered = true;
+ if(_vals[i]==0. && !showered)
+ {
+ firstNonZeroLayer = i+1;
+ }
+ }
+
+ // shift the energies down so first showering layer is layer0
+ if(firstNonZeroLayer!=0)
+ {
+ double[] newvals = new double[_nmeas];
+ System.arraycopy(_vals,firstNonZeroLayer, newvals, 0, layerE.length-firstNonZeroLayer-1 );
+ _vals = newvals;
+ }
+ // Add the correlation to the log of the cluster energy
+ //We want the common logarithm (log 10) of the energy
+ _vals[_logEIndex]=Math.log(cluster.getEnergy())*_log10inv;
+ // now calculate the cluster chi-squared
+ chisq = _hmx.chisquared(_vals);
+ chisqD = _hmx.chisquaredDiagonal(_vals);
+ chisqDProb = ChisqProb.gammq((double) (_nmeas),chisqD);
+ if(chisqDProb<.0000000001) chisqDProb = 0.0000000001;
+ }
+ public double getChisq(){return chisq;}
+ public double getChisqD(){return chisqD;}
+ public double getChisqDProb(){return chisqDProb;}
+ private double[] layerEnergies(Cluster clus)
+ {
+ //FIXME could reuse this array
+ double[] layerEnergies = new double[_nLayers];
+ double clusterEnergy = 0.;
+ List<CalorimeterHit> hits = clus.getCalorimeterHits();
+ for(CalorimeterHit hit : hits)
+ {
+ IDDecoder decoder = hit.getIDDecoder();
+ decoder.setID(hit.getCellID());
+ double e = hit.getCorrectedEnergy();
+ int layer = decoder.getLayer();
+ clusterEnergy+=e;
+ layerEnergies[layer]+=e;
+ }
+ for(int i=0; i<_nLayers; ++i)
+ {
+ layerEnergies[i]/=clusterEnergy;
+ }
+ return layerEnergies;
+ }
+}
diff -N IdentifyPhotonClustersSid01.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ IdentifyPhotonClustersSid01.java 12 Dec 2008 21:49:52 -0000 1.1 @@ -0,0 +1,327 @@
+package org.lcsim.recon.pfa.photonfinder;
+
+import java.util.*;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.SimCalorimeterHit;
+import hep.physics.vec.*;
+import org.lcsim.geometry.IDDecoder;
+
+/**
+ * Photon identifier
+ *
+ * Code by Ron Cassell, put into CVS in a temporary place by Mat Charles.
+ */
+public class IdentifyPhotonClustersSid01
+{
+ int nhitcut = 5;
+ int[] sizes = {7,10,15,22,30,35,50,100,200,400};
+ HMatrixVars hmv;
+ double[][][] cuts;
+ int failcut;
+ public IdentifyPhotonClustersSid01(HMatrixVars HMV)
+ {
+ cuts = new double[sizes.length+1][2][10];
+//
+// Separate cuts for cluster size, and whether layer 0 is hit.
+// Cut 0 - Hmatrix chisq
+// Cut 1 - First layer hit
+// Cut 2 - Shape 0 (transverse spread)
+// Cut 3 - Shape2-Shape1 (Asymmetry)
+// Cut 4 - Shape2-Shape0 (Long-Trans Spread)
+// Cut 5 - DOCA
+// Cut 6 - mean # hits in first 5 layers
+// Cut 7 - lower bound on chisq for small clusters
+// Cut 8 - lower bound on Shape 0 for small clusters
+// Cut 9 - upper bound on Shape 1 for small clusters
+//
+ cuts[0][0][0] = 0;
+ cuts[0][1][0] = 800;
+ cuts[1][0][0] = 300;
+ cuts[1][1][0] = 800;
+ cuts[2][0][0] = 150;
+ cuts[2][1][0] = 250;
+ cuts[3][0][0] = 150;
+ cuts[3][1][0] = 180;
+ cuts[4][0][0] = 100;
+ cuts[4][1][0] = 170;
+ cuts[5][0][0] = 90;
+ cuts[5][1][0] = 150;
+ cuts[6][0][0] = 80;
+ cuts[6][1][0] = 260;
+ cuts[7][0][0] = 90.;
+ cuts[7][1][0] = 800;
+ cuts[8][0][0] = 90;
+ cuts[8][1][0] = 9999;
+ cuts[9][0][0] = 96;
+ cuts[9][1][0] = 9999;
+ cuts[10][0][0] = 190;
+ cuts[10][1][0] = 9999;
+ cuts[0][0][1] = 17.;
+ cuts[0][1][1] = 7.;
+ cuts[1][0][1] = 17.;
+ cuts[1][1][1] = 8.;
+ cuts[2][0][1] = 17.;
+ cuts[2][1][1] = 10.;
+ cuts[3][0][1] = 17.;
+ cuts[3][1][1] = 9.;
+ cuts[4][0][1] = 17.;
+ cuts[4][1][1] = 9.;
+ cuts[5][0][1] = 17.;
+ cuts[5][1][1] = 9.;
+ cuts[6][0][1] = 17.;
+ cuts[6][1][1] = 17.;
+ cuts[7][0][1] = 17.;
+ cuts[7][1][1] = 17.;
+ cuts[8][0][1] = 17.;
+ cuts[8][1][1] = 17.;
+ cuts[9][0][1] = 17.;
+ cuts[9][1][1] = 17.;
+ cuts[10][0][1] = 17.;
+ cuts[10][1][1] = 17.;
+ cuts[0][0][2] = 15.;
+ cuts[0][1][2] = 10.;
+ cuts[1][0][2] = 10.;
+ cuts[1][1][2] = 12.;
+ cuts[2][0][2] = 15.;
+ cuts[2][1][2] = 30.;
+ cuts[3][0][2] = 15.;
+ cuts[3][1][2] = 30.;
+ cuts[4][0][2] = 30.;
+ cuts[4][1][2] = 40.;
+ cuts[5][0][2] = 30.;
+ cuts[5][1][2] = 41.;
+ cuts[6][0][2] = 30.;
+ cuts[6][1][2] = 45.;
+ cuts[7][0][2] = 37.;
+ cuts[7][1][2] = 80.;
+ cuts[8][0][2] = 40.;
+ cuts[8][1][2] = 9999.;
+ cuts[9][0][2] = 40.;
+ cuts[9][1][2] = 9999.;
+ cuts[10][0][2] = 9999.;
+ cuts[10][1][2] = 9999.;
+ cuts[0][0][3] = 8.;
+ cuts[0][1][3] = 10.;
+ cuts[1][0][3] = 4.;
+ cuts[1][1][3] = 13.;
+ cuts[2][0][3] = 8.;
+ cuts[2][1][3] = 30.;
+ cuts[3][0][3] = 10.;
+ cuts[3][1][3] = 30.;
+ cuts[4][0][3] = 20.;
+ cuts[4][1][3] = 31.;
+ cuts[5][0][3] = 20.;
+ cuts[5][1][3] = 46.;
+ cuts[6][0][3] = 25.;
+ cuts[6][1][3] = 80.;
+ cuts[7][0][3] = 20.;
+ cuts[7][1][3] = 80.;
+ cuts[8][0][3] = 20.;
+ cuts[8][1][3] = 9999.;
+ cuts[9][0][3] = 35.;
+ cuts[9][1][3] = 9999.;
+ cuts[10][0][3] = 9999.;
+ cuts[10][1][3] = 9999.;
+ cuts[0][0][4] = 90.;
+ cuts[0][1][4] = 60.;
+ cuts[1][0][4] = 90.;
+ cuts[1][1][4] = 80.;
+ cuts[2][0][4] = 90.;
+ cuts[2][1][4] = 100.;
+ cuts[3][0][4] = 200.;
+ cuts[3][1][4] = 200.;
+ cuts[4][0][4] = 200.;
+ cuts[4][1][4] = 280.;
+ cuts[5][0][4] = 300.;
+ cuts[5][1][4] = 580.;
+ cuts[6][0][4] = 400.;
+ cuts[6][1][4] = 510.;
+ cuts[7][0][4] = 1100.;
+ cuts[7][1][4] = 1100.;
+ cuts[8][0][4] = 1100.;
+ cuts[8][1][4] = 9999.;
+ cuts[9][0][4] = 1100.;
+ cuts[9][1][4] = 9999.;
+ cuts[10][0][4] = 1100.;
+ cuts[10][1][4] = 9999.;
+ cuts[0][0][5] = 2000.;
+ cuts[0][1][5] = 1360.;
+ cuts[1][0][5] = 2000.;
+ cuts[1][1][5] = 1300.;
+ cuts[2][0][5] = 2000.;
+ cuts[2][1][5] = 2000.;
+ cuts[3][0][5] = 2000.;
+ cuts[3][1][5] = 1800.;
+ cuts[4][0][5] = 2000.;
+ cuts[4][1][5] = 1800.;
+ cuts[5][0][5] = 2000.;
+ cuts[5][1][5] = 1500.;
+ cuts[6][0][5] = 2000.;
+ cuts[6][1][5] = 850.;
+ cuts[7][0][5] = 1000.;
+ cuts[7][1][5] = 500.;
+ cuts[8][0][5] = 500.;
+ cuts[8][1][5] = 9999.;
+ cuts[9][0][5] = 410.;
+ cuts[9][1][5] = 9999.;
+ cuts[10][0][5] = 9999.;
+ cuts[10][1][5] = 9999.;
+ cuts[0][0][6] = 1.3;
+ cuts[0][1][6] = 0.;
+ cuts[1][0][6] = 1.5;
+ cuts[1][1][6] = 0.;
+ cuts[2][0][6] = 1.3;
+ cuts[2][1][6] = 0.;
+ cuts[3][0][6] = 1.3;
+ cuts[3][1][6] = 1.1;
+ cuts[4][0][6] = 1.3;
+ cuts[4][1][6] = 0.;
+ cuts[5][0][6] = 1.3;
+ cuts[5][1][6] = 0.;
+ cuts[6][0][6] = 1.3;
+ cuts[6][1][6] = 0.;
+ cuts[7][0][6] = 1.3;
+ cuts[7][1][6] = 0.;
+ cuts[8][0][6] = 1.3;
+ cuts[8][1][6] = 0.;
+ cuts[9][0][6] = 0.;
+ cuts[9][1][6] = 0.;
+ cuts[10][0][6] =1.3;
+ cuts[10][1][6] = 0.;
+ cuts[0][0][7] = 0.;
+ cuts[0][1][7] = 46.;
+ cuts[1][0][7] = 32.;
+ cuts[1][1][7] = 24.;
+ cuts[2][0][7] = 0.;
+ cuts[2][1][7] = 0.;
+ cuts[3][0][7] = 0.;
+ cuts[3][1][7] = 0.;
+ cuts[4][0][7] = 0.;
+ cuts[4][1][7] = 0.;
+ cuts[5][0][7] = 0.;
+ cuts[5][1][7] = 0.;
+ cuts[6][0][7] = 0.;
+ cuts[6][1][7] = 0.;
+ cuts[7][0][7] = 0.;
+ cuts[7][1][7] = 0.;
+ cuts[8][0][7] = 0.;
+ cuts[8][1][7] = 0.;
+ cuts[9][0][7] = 0.;
+ cuts[9][1][7] = 0.;
+ cuts[10][0][7] = 0.;
+ cuts[10][1][7] = 0.;
+ cuts[0][0][8] = 0.;
+ cuts[0][1][8] = 0.;
+ cuts[1][0][8] = 2.;
+ cuts[1][1][8] = 0.;
+ cuts[2][0][8] = 0.;
+ cuts[2][1][8] = 0.;
+ cuts[3][0][8] = 0.;
+ cuts[3][1][8] = 0.;
+ cuts[4][0][8] = 0.;
+ cuts[4][1][8] = 0.;
+ cuts[5][0][8] = 0.;
+ cuts[5][1][8] = 0.;
+ cuts[6][0][8] = 0.;
+ cuts[6][1][8] = 0.;
+ cuts[7][0][8] = 0.;
+ cuts[7][1][8] = 0.;
+ cuts[8][0][8] = 0.;
+ cuts[8][1][8] = 0.;
+ cuts[9][0][8] = 0.;
+ cuts[9][1][8] = 0.;
+ cuts[10][0][8] = 0.;
+ cuts[10][1][8] = 0.;
+ cuts[0][0][9] = 9990.;
+ cuts[0][1][9] = 9990.;
+ cuts[1][0][9] = 70.;
+ cuts[1][1][9] = 9990.;
+ cuts[2][0][9] = 9990.;
+ cuts[2][1][9] = 9990.;
+ cuts[3][0][9] = 9990.;
+ cuts[3][1][9] = 9990.;
+ cuts[4][0][9] = 9990.;
+ cuts[4][1][9] = 9990.;
+ cuts[5][0][9] = 9990.;
+ cuts[5][1][9] = 9990.;
+ cuts[6][0][9] = 9990.;
+ cuts[6][1][9] = 9990.;
+ cuts[7][0][9] = 9990.;
+ cuts[7][1][9] = 9990.;
+ cuts[8][0][9] = 9990.;
+ cuts[8][1][9] = 9990.;
+ cuts[9][0][9] = 9990.;
+ cuts[9][1][9] = 9990.;
+ cuts[10][0][9] = 9990.;
+ cuts[10][1][9] = 9990.;
+
+ hmv = HMV;
+ }
+ public boolean isPhoton(Cluster c)
+ {
+ int nhits = c.getSize();
+ if(nhits <= nhitcut)return false;
+ int sbin = 0;
+ for(int i=0;i<sizes.length;i++)
+ {
+ if(nhits <= sizes[i])break;
+ sbin++;
+ }
+ int fl = 30;
+ int ll = -1;
+ int[] hpl = new int[31];
+ for(CalorimeterHit h:c.getCalorimeterHits())
+ {
+ SimCalorimeterHit sh = (SimCalorimeterHit) h;
+ IDDecoder iddc = h.getIDDecoder();
+ iddc.setID(h.getCellID());
+ int lay = iddc.getValue("layer");
+ if(lay < fl)fl = lay;
+ if(lay > ll)ll = lay;
+ hpl[lay]++;
+ }
+ double[] sp = c.getShape();
+ double del0 = sp[1] - sp[0];
+ double del1 = sp[2] - sp[0];
+ double del2 = sp[2] - sp[1];
+ double[] pos = c.getPosition();
+ double phi = c.getIPhi();
+ double theta = c.getITheta();
+ double cp = Math.cos(phi);
+ double ssp = Math.sin(phi);
+ double ct = Math.cos(theta);
+ double st = Math.sin(theta);
+ double dcasq = pos[0]*pos[0]*(1.-st*st*cp*cp) +
+ pos[1]*pos[1]*(1.-st*st*ssp*ssp) +
+ pos[2]*pos[2]*(1.-ct*ct) -
+ 2.*pos[0]*pos[1]*st*st*ssp*cp -
+ 2.*pos[0]*pos[2]*st*ct*cp -
+ 2.*pos[1]*pos[2]*st*ct*ssp;
+ double dca = -1;
+ if(dcasq > 0.)dca = Math.sqrt(dcasq);
+ int xll = ll;
+ if(fl+4 < ll)xll = fl+4;
+ double t5 = 0.;
+ for(int i=fl;i<=xll;i++){t5+=hpl[i];}
+ double mhpl5 = t5/(1. + xll - fl);
+ hmv.setCluster(c);
+ double chi = hmv.getChisq();
+ int flb = 0;
+ if(fl > 0)flb = 1;
+ failcut = 0;
+ if(chi > cuts[sbin][flb][0])failcut = 1;
+ else if(fl > cuts[sbin][flb][1])failcut = 2;
+ else if(sp[0] > cuts[sbin][flb][2])failcut = 3;
+ else if( sp[2]-sp[1] > cuts[sbin][flb][3])failcut = 4;
+ else if( sp[2] - sp[0] > cuts[sbin][flb][4])failcut = 5;
+ else if(dca > cuts[sbin][flb][5])failcut = 6;
+ else if(mhpl5 < cuts[sbin][flb][6])failcut = 7;
+ else if(chi < cuts[sbin][flb][7])failcut = 8;
+ else if(sp[0] < cuts[sbin][flb][8])failcut = 9;
+ else if(sp[1] > cuts[sbin][flb][9])failcut = 10;
+ if(failcut == 0)return true;
+ return false;
+ }
+ public int getFailcut(){return failcut;}
+}
diff -N RonPhotonFinderSid01.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ RonPhotonFinderSid01.java 12 Dec 2008 21:49:52 -0000 1.1 @@ -0,0 +1,55 @@
+package org.lcsim.recon.pfa.photonfinder;
+
+import java.util.List;
+
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.recon.cluster.nn.NearestNeighborClusterer;
+import org.lcsim.util.Driver;
+import org.lcsim.util.hitmap.HitMap;
+import org.lcsim.util.lcio.LCIOConstants;
+
+/**
+ * Driver to find photons in Sid01 detector
+ *
+ * Code by Ron Cassell, put into CVS in a temporary place by Mat Charles.
+ */
+public class RonPhotonFinderSid01 extends Driver
+{
+ String inName;
+ String outHName;
+ String outCName;
+ NearestNeighborClusterer nn111;
+ HMatrixVars hmv;
+ Sid01PhotonFinder finder;
+ public RonPhotonFinderSid01(String in, String outH, String outC)
+ {
+ inName = in;
+ outHName = outH;
+ outCName = outC;
+ nn111 = new NearestNeighborClusterer(1,1,1,0,0.);
+ hmv = new HMatrixVars();
+ add(hmv);
+ finder = new Sid01PhotonFinder(hmv);
+ }
+ protected void process(EventHeader event)
+ {
+ super.process(event);
+ HitMap inmap = (HitMap) event.get(inName);
+ List<Cluster> photons = finder.findPhotons(nn111.createClusters(inmap));
+ int flag = 1<<LCIOConstants.CLBIT_HITS;
+ event.put(outCName, photons, Cluster.class, flag );
+ HitMap outputHitMap = new HitMap(inmap);
+// Remove photon hits
+ for (Cluster clus : photons)
+ {
+ for (CalorimeterHit hit : clus.getCalorimeterHits())
+ {
+ Long cellID = new Long(hit.getCellID());
+ outputHitMap.remove(cellID);
+ }
+ }
+ event.put(outHName, outputHitMap);
+ }
+}
diff -N Sid01PhotonFinder.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ Sid01PhotonFinder.java 12 Dec 2008 21:49:52 -0000 1.1 @@ -0,0 +1,35 @@
+package org.lcsim.recon.pfa.photonfinder;
+
+import java.util.*;
+import org.lcsim.event.Cluster;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import hep.physics.vec.*;
+
+/**
+ * Photon finder
+ *
+ * Code by Ron Cassell, put into CVS in a temporary place by Mat Charles.
+ */
+public class Sid01PhotonFinder
+{
+ IdentifyPhotonClustersSid01 ider;
+ AssociatePhotonFragmentsSid01 ass;
+ HMatrixVars hmv;
+ public Sid01PhotonFinder(HMatrixVars HMV)
+ {
+ hmv = HMV;
+ ider = new IdentifyPhotonClustersSid01(hmv);
+ ass = new AssociatePhotonFragmentsSid01();
+ }
+ public List<Cluster> findPhotons(List<Cluster> cllist)
+ {
+ List<Cluster> photons = new ArrayList<Cluster>();
+ List<Cluster> fragments = new ArrayList<Cluster>();
+ for(Cluster c:cllist)
+ {
+ if(ider.isPhoton(c))photons.add(c);
+ else fragments.add(c);
+ }
+ return ass.associateFragments(photons,fragments);
+ }
+}
diff -u -r1.3 -r1.4 --- NonTrivialPFA.java 6 Sep 2008 23:47:26 -0000 1.3 +++ NonTrivialPFA.java 12 Dec 2008 21:49:52 -0000 1.4 @@ -29,7 +29,7 @@
* a List<ReconstructedParticle>, written to the event as * PFAReconstructedParticles. *
- * @version $Id: NonTrivialPFA.java,v 1.3 2008/09/06 23:47:26 mcharles Exp $
+ * @version $Id: NonTrivialPFA.java,v 1.4 2008/12/12 21:49:52 jeremy Exp $
* @author Mat Charles <[log in to unmask]> */
@@ -162,7 +162,7 @@
// Find tracks (cheating) String ronTrackList = "FSReconTracks"; String ronMCList = "ReconFSParticles";
- add(new org.lcsim.contrib.Cassell.recon.Cheat.CheatReconDriver());
+ add(new org.lcsim.recon.cheater.CheatReconDriver());
// Choose which track list to use boolean useRonTrackList = true;
@@ -541,7 +541,7 @@
protected void addRonPhotonFinder(String prefix, String inputHitMap, String outputPhotonClusterList, String outputHitMap)
{
- org.lcsim.contrib.uiowa.RonPhotonFinder.RonPhotonFinderSid01 driver = new org.lcsim.contrib.uiowa.RonPhotonFinder.RonPhotonFinderSid01(inputHitMap, outputHitMap, outputPhotonClusterList);
+ org.lcsim.recon.pfa.photonfinder.RonPhotonFinderSid01 driver = new org.lcsim.recon.pfa.photonfinder.RonPhotonFinderSid01(inputHitMap, outputHitMap, outputPhotonClusterList);
add(driver);
}
diff -u -r1.7 -r1.8 --- RunAndWriteOutPFA.java 22 Oct 2008 18:06:30 -0000 1.7 +++ RunAndWriteOutPFA.java 12 Dec 2008 21:49:52 -0000 1.8 @@ -31,7 +31,7 @@
public RunAndWriteOutPFA()
{
// Prepare to run PFA: Tracks (includes DigiSim)
- add(new org.lcsim.contrib.Cassell.recon.Cheat.CheatReconDriver());
+ add(new org.lcsim.recon.cheater.CheatReconDriver());
// Set up and run PFA
add(new SetUpPFA("FSReconTracks"));