lcsim/src/org/lcsim/contrib/SteveMagill
diff -u -r1.1 -r1.2
--- NHCandClusterDriver.java 18 May 2007 19:58:44 -0000 1.1
+++ NHCandClusterDriver.java 20 Aug 2008 20:24:20 -0000 1.2
@@ -14,6 +14,8 @@
import org.lcsim.geometry.*;
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.BasicHep3Vector;
+import org.lcsim.spacegeom.*;
+import org.lcsim.event.Track;
import hep.aida.*;
public class NHCandClusterDriver extends Driver
@@ -23,14 +25,15 @@
private String _outclusname;
private int _mincells;
private double _Emin;
- private double _ClCldist;
+ private double _TrCldist;
+ private boolean nhdb = false;
- public NHCandClusterDriver(int mincells, double Emin, double ClCldist)
+ public NHCandClusterDriver(int mincells, double Emin, double TrCldist)
{
// add arguments if needed
_mincells = mincells;
_Emin = Emin;
- _ClCldist = ClCldist;
+ _TrCldist = TrCldist;
}
protected void process(EventHeader event)
@@ -38,8 +41,10 @@
super.process(event); // executes all added drivers
List<BasicCluster> nhclus = new ArrayList<BasicCluster>();
- List<PerfectTrack> trks = event.get(PerfectTrack.class, "PerfectTracks");
- List<BasicCluster> trshcls = event.get(BasicCluster.class, "TrackCALClusters");
+ Map<Track, SpacePoint> tresmmap = (Map<Track, SpacePoint>) event.get("TrackXEShMaxMap");
+ Map<Track, SpacePoint> trh0map = (Map<Track, SpacePoint>) event.get("TrackXH0Map");
+ Map<Track, BasicCluster> trkclmap = (Map<Track, BasicCluster>) event.get("TrackClusMap");
+
// to keep as neutral, minimum number of hits, distance from nearest track cluster, etc.
try
{
@@ -48,50 +53,131 @@
{
if (cl.getCalorimeterHits().size()>_mincells && cl.getEnergy()>_Emin)
{
+ // if cluster has hit in ECAL layer 0 and no hits in hcal, not NH
+ int nzhits = 0;
+ int nhchits = 0;
+ String emb = "EMBarrel";
+ String emec = "EMEndcap";
+ String hadb = "HADBarrel";
+ String hadec = "HADEndcap";
+ List<CalorimeterHit> clhits = cl.getCalorimeterHits();
+ for(CalorimeterHit clhit : clhits)
+ {
+ int calflag = 0;
+ IDDecoder _decoder = clhit.getIDDecoder();
+ String dname = _decoder.getSubdetector().getName();
+// System.out.println(" name of subdetector for this hit " + dname);
+ if (dname.equals(emb)) calflag = 1;
+ if (dname.equals(emec)) calflag = 1;
+ if (dname.equals(hadb)) calflag = 2;
+ if (dname.equals(hadec)) calflag = 2;
+// System.out.println(" calflag is " + calflag);
+ _decoder.setID(clhit.getCellID());
+ int layer =_decoder.getLayer();
+ if (layer==0 && calflag==1) nzhits++;
+ if (calflag==2) nhchits++;
+ }
+// System.out.println(" num zero layer hits " + nzhits + " HAD hits " + nhchits);
+ if (nzhits==1 && nhchits==0) continue; // not a nh - probably a mip stub
+
// get cluster hit positions
- aida.cloud1D("NHCand Clus E").fill(cl.getEnergy());
+ if (nhdb) aida.cloud1D("NHCand Clus E").fill(cl.getEnergy());
double[] clpos = cl.getPosition();
double cph = Math.atan2(clpos[1],clpos[0]);
if (cph<0) cph+=2*Math.PI;
double cshr = Math.sqrt(clpos[0]*clpos[0]+clpos[1]*clpos[1]);
double cth = Math.atan(cshr/clpos[2]);
if (cth<0) cth+=Math.PI;
+ double clrawE = cl.getRawEnergy();
+
+ // now check distance of cluster from track positions
int nclmtch = 0;
- for (BasicCluster trshcl : trshcls)
+ double dist1 = 999.;
+ double dist2 = 999.;
+ for (Track itr : tresmmap.keySet())
{
- // get this cluster position
- double[] tpos = trshcl.getPosition();
- double tph = Math.atan2(tpos[1],tpos[0]);
- if (tph<0) tph+=2*Math.PI;
- double tshr = Math.sqrt(tpos[0]*tpos[0]+tpos[1]*tpos[1]);
- double tth = Math.atan(tshr/tpos[2]);
- if (tth<0) tth+=Math.PI;
- double dccth = Math.abs(tth-cth);
- double dccph = Math.abs(tph-cph);
- if (dccph>Math.PI) dccph = 2*Math.PI-dccph;
- double dcc = Math.sqrt(dccth*dccth+dccph*dccph);
- aida.cloud1D("ClusClus distance in NHCand").fill(dcc);
- if (dcc<_ClCldist) nclmtch++;
- if (nclmtch>0) continue;
+ // use maps at ESM and H0
+ SpacePoint trSPesm = tresmmap.get(itr);
+ double tresmth = Math.atan(trSPesm.rxy()/trSPesm.z());
+ if (tresmth<0) tresmth+=Math.PI;
+ double tresmph = Math.atan2(trSPesm.y(),trSPesm.x());
+ if (tresmph<0) tresmph+=2*Math.PI;
+ double d1ccth = Math.abs(tresmth-cth);
+ double d1ccph = Math.abs(tresmph-cph);
+ if (d1ccph>Math.PI) d1ccph = 2*Math.PI-d1ccph;
+ dist1 = Math.sqrt(d1ccth*d1ccth+d1ccph*d1ccph);
+
+ SpacePoint trSPh0 = trh0map.get(itr);
+ double trh0th = Math.atan(trSPh0.rxy()/trSPh0.z());
+ if (trh0th<0) trh0th+=Math.PI;
+ double trh0ph = Math.atan2(trSPh0.y(),trSPh0.x());
+ if (trh0ph<0) trh0ph+=2*Math.PI;
+ double d2ccth = Math.abs(trh0th-cth);
+ double d2ccph = Math.abs(trh0ph-cph);
+ if (d2ccph>Math.PI) d2ccph = 2*Math.PI-d2ccph;
+ dist2 = Math.sqrt(d2ccth*d2ccth+d2ccph*d2ccph);
+
+ if (dist1<dist2) dist2 = dist1;
+ if (dist2<_TrCldist)
+ {
+ // 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+cl.getEnergy())/pmag;
+ if (eopplus<2.0)
+ {
+ nclmtch++;
+ // add this cluster to track cluster
+ clus.addCluster(cl);
+ }
+ }
}
- if (nclmtch == 0) nhclus.add(cl);
- }
- }
- if (nhclus.size()>0) event.put(_outclusname,nhclus);
+ if (nclmtch>0) continue;
+
+ // check average energy-weighted cluster hit distances
+ List<CalorimeterHit> chits = cl.getCalorimeterHits();
+ double Efdist = 0;
+ for (CalorimeterHit chit : chits)
+ {
+ double[] htpos = chit.getPosition();
+ double htph = Math.atan2(htpos[1],htpos[0]);
+ if (htph<0) htph+=2*Math.PI;
+ double htr = Math.sqrt(htpos[0]*htpos[0]+htpos[1]*htpos[1]);
+ double htth = Math.atan(htr/htpos[2]);
+ if (htth<0) htth+=Math.PI;
+ // get distance of hit from cluster center
+ double disth = Math.abs(htth-cth);
+ double disph = Math.abs(htph-cph);
+ if (disph>Math.PI) disph = 2*Math.PI-disph;
+ aida.cloud2D("deltheta vs delphi dist").fill(disph,disth);
+ Efdist += chit.getRawEnergy()*Math.sqrt(disth*disth+disph*disph);
+ }
+ double AEdis = Efdist/clrawE;
+ aida.cloud1D("EwgtSum Hit Dist from NH cluster center").fill(AEdis);
+// if (AEdis<0.01) continue;
+
+ nhclus.add(cl);
+ } // test on cluster size and energy
+ } // end of cluster loop
+// if (nhclus.size()>0) event.put(_outclusname,nhclus);
+ event.put(_outclusname,nhclus);
double nhclE = 0;
for (BasicCluster nhcl : nhclus)
{
- nhclE += nhcl.getEnergy()*1.3;
+ nhclE += nhcl.getEnergy();
}
- aida.cloud1D("PFA Neutral Hadron ESum").fill(nhclE);
+ if (nhdb) aida.cloud1D("PFA Neutral Hadron ESum").fill(nhclE);
}
catch(java.lang.IllegalArgumentException ex)
{
System.out.println("requested object not found in event " + _clusname);
}
- } // end of event loop
-
+ } // end of event loop
+
public void setInputClusterName(String name)
{
_clusname = name;