lcsim/src/org/lcsim/contrib/SteveMagill
diff -N DTPhoAcceptDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DTPhoAcceptDriver.java 20 Aug 2008 20:25:14 -0000 1.1
@@ -0,0 +1,181 @@
+package org.lcsim.contrib.SteveMagill;
+
+import java.util.*;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.geometry.util.CalorimeterIDDecoder;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import org.lcsim.event.Track;
+import org.lcsim.util.hitmap.HitMap;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.spacegeom.*;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.MCReconstructedParticle;
+import org.lcsim.event.ParticleID;
+import org.lcsim.recon.cluster.cheat.CheatCluster;
+import org.lcsim.event.ParticleID;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.base.BaseParticleID;
+import org.lcsim.event.base.BaseReconstructedParticle;
+import hep.physics.particle.properties.ParticlePropertyManager;
+import hep.physics.particle.properties.ParticleType;
+
+/**
+ * Checks candidate clusters for closeness to tracks and adjusts
+ */
+
+public class DTPhoAcceptDriver extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ private int _mincells;
+ private double _dTrcl;
+ private double _minE;
+ private String[] _phnames;
+ private String _outname;
+ private boolean DTD = false;
+
+ public DTPhoAcceptDriver(int mincells, double dTrcl, double minE)
+ {
+ _mincells = mincells;
+ _dTrcl = dTrcl;
+ _minE = minE;
+ }
+
+ public void process(EventHeader event)
+ {
+ // combine all input clusters into a single shower cluster list
+ List<BasicCluster> combclusters = new ArrayList<BasicCluster>();
+// System.out.println(" Num of Clus Lists " + _clusternames.length);
+ for (int i=0; i<_phnames.length; i++)
+ {
+ try
+ {
+ List<BasicCluster> clus = event.get(BasicCluster.class,_phnames[i]);
+// System.out.println(" Num of Clus in list " + clus.size());
+ for (BasicCluster cl : clus)
+ {
+ BasicCluster cclus = new BasicCluster();
+ if (cl.getSize()>0)
+ {
+ // if (_phnames[i]=="DTPhotonClusters") aida.cloud1D("Energy of DT Clusters").fill(cl.getEnergy());
+ // if (_phnames[i]=="LEPhotonClusters") aida.cloud1D("Energy of LE Clusters").fill(cl.getEnergy());
+ if (DTD) aida.cloud1D("Num Hits in Comb Clus").fill(cl.getSize());
+ if (DTD) aida.cloud2D("Num hits Comb Clus vs E").fill(cl.getEnergy(),cl.getSize());
+ cclus.addCluster(cl);
+ combclusters.add(cclus);
+ }
+ }
+ }
+ catch(java.lang.IllegalArgumentException ex)
+ {
+ System.out.println("requested object not found in event " + _phnames[i]);
+ }
+ }
+ event.put("CombinedPhoCandClusters",combclusters);
+
+ // get layer 0 and shower max maps, also mip cluster map
+// Map<Track, BasicCluster> trmipmap = (Map<Track, BasicCluster>) event.get("TrackMipMap");
+ Map<Track, BasicCluster> trkclmap = (Map<Track, BasicCluster>) event.get("TrackClusMap");
+ Map<Track, SpacePoint> tre0map = (Map<Track, SpacePoint>) event.get("TrackXE0Map");
+ Map<Track, SpacePoint> tresmmap = (Map<Track, SpacePoint>) event.get("TrackXEShMaxMap");
+ // create new final cluster list
+ List<BasicCluster> FPhotonclusters = new ArrayList<BasicCluster>(); // final photon clusters
+
+ // check found photon clusters
+ for (BasicCluster combclus : combclusters)
+ {
+ if (combclus.getEnergy()>_minE && combclus.getSize()>_mincells)
+ {
+ int ncltr = 0;
+ double ecp[] = combclus.getPosition();
+ double ecpx = ecp[0];
+ double ecpy = ecp[1];
+ double ecpz = ecp[2];
+ double ecR = Math.sqrt(ecpx*ecpx+ecpy*ecpy);
+ double clth = Math.atan(ecR/ecpz);
+ if (clth<0) clth+=Math.PI;
+ double clph = Math.atan2(ecpy,ecpx);
+ if (clph<0) clph+=2*Math.PI;
+
+ // now check tracks
+ int tmtch = 0;
+ // first, check with first layer of ECAL
+ for (Track itr : tre0map.keySet())
+ {
+ SpacePoint tre0sp = tre0map.get(itr);
+ double tre0th = Math.atan(tre0sp.rxy()/tre0sp.z());
+ if (tre0th<0) tre0th+=Math.PI;
+ double tre0ph = Math.atan2(tre0sp.y(),tre0sp.x());
+ if (tre0ph<0) tre0ph+=2*Math.PI;
+ double dist0th = Math.abs(clth-tre0th);
+ double dist0ph = Math.abs(clph-tre0ph);
+ if (dist0ph>Math.PI) dist0ph = 2*Math.PI-dist0ph;
+ double dist0clmip = Math.sqrt(dist0th*dist0th+dist0ph*dist0ph);
+ if (dist0clmip<0.2)
+ {
+ if (DTD) aida.cloud1D("PhoCand Clus Tre0 Dist").fill(dist0clmip);
+ if (DTD) aida.cloud2D("Theta vs Phi of Tre0").fill(tre0ph,tre0th);
+ }
+ // now check dist at shower max also
+ SpacePoint trsmsp = tresmmap.get(itr);
+ double tresmth = Math.atan(trsmsp.rxy()/trsmsp.z());
+ if (tresmth<0) tresmth+=Math.PI;
+ double tresmph = Math.atan2(trsmsp.y(),trsmsp.x());
+ if (tresmph<0) tresmph+=2*Math.PI;
+ double distmth = Math.abs(clth-tresmth);
+ double distmph = Math.abs(clph-tresmph);
+ if (distmph>Math.PI) distmph = 2*Math.PI-distmph;
+ double distmclmip = Math.sqrt(distmth*distmth+distmph*distmph);
+ if (distmclmip<0.2)
+ {
+ if (DTD) aida.cloud1D("PhoCand Clus TrSM Dist").fill(distmclmip);
+ if (DTD) aida.cloud2D("Theta vs Phi of TrSM").fill(tresmph,tresmth);
+ }
+ double distclmip = distmclmip;
+ if (dist0clmip<distmclmip) distclmip=dist0clmip;
+ // compare to track and check eoverp
+ if (distclmip<_dTrcl)
+ {
+ // get E/p for this track to see if cluster should be merged or not
+ BasicCluster clus = trkclmap.get(itr);
+ double clE = clus.getEnergy();
+ double[] trP = itr.getMomentum();
+ double pmag = Math.sqrt(trP[0]*trP[0]+trP[1]*trP[1]+trP[2]*trP[2]);
+ double eop = clE/pmag;
+ double eopplus = (clE+combclus.getEnergy())/pmag;
+ if (eopplus<2.0)
+ {
+ tmtch++;
+ // add this cluster to track cluster
+ clus.addCluster(combclus);
+ }
+ }
+ }
+ // if no matches, keep as photon
+ if (tmtch==0)
+ {
+ FPhotonclusters.add(combclus);
+ }
+ }
+ }
+ event.put(_outname,FPhotonclusters);
+ }
+
+ public void setInputClusNames(String[] phnames)
+ {
+ _phnames = phnames;
+ }
+
+ public void setOutputClusNames(String outname)
+ {
+ _outname = outname;
+ }
+
+}
