lcsim/src/org/lcsim/contrib/SteveMagill
diff -u -r1.1 -r1.2
--- TrackShowerClusterDriver.java 23 Apr 2007 20:08:53 -0000 1.1
+++ TrackShowerClusterDriver.java 18 May 2007 19:58:19 -0000 1.2
@@ -23,6 +23,7 @@
import org.lcsim.geometry.util.CalorimeterIDDecoder;
import org.lcsim.event.ReconstructedParticle;
import org.lcsim.event.base.MCReconstructedParticle;
+import org.lcsim.recon.cluster.cheat.CheatCluster;
public class TrackShowerClusterDriver extends Driver
{
@@ -61,24 +62,24 @@
private double _dminE;
private double _dminH;
private double _dminTrC;
+ private double _mindist;
+ private int _nloops;
+ private double _mineop;
private String[] _clusternames;
// private String[] _hitmapnames;
- public TrackShowerClusterDriver()
+ public TrackShowerClusterDriver(double mindist, int nloops, double mineop)
{
// add arguments if needed
+ _mindist = mindist; // minimum distance between objects to merge (track, clusters) (0.015)
+ _nloops = nloops; // number of iterations to merge clusters (until min E/p is met) (4)
+ _mineop = mineop; // minimum E/p value (0.65)
}
protected void process(EventHeader event)
{
super.process(event); // executes all added drivers
- boolean mcfltr = true;
-// mcfltr = (Boolean)event.get("MCFilt");
-// System.out.println("MCFilt is " +mcfltr);
- if (mcfltr)
- {
-
// Initialize things here
if(!_initialized)
{
@@ -155,19 +156,24 @@
}
}
// System.out.println("Number Shower Clusters before Matching " + showclus.size());
-
+ aida.cloud1D("Num Shower Clus before TrackMatch").fill(showclus.size());
// create a list of track-cal associated clusters for display
List<BasicCluster> trkcalclus = new ArrayList<BasicCluster>();
+ List<BasicCluster> trkshoclus = new ArrayList<BasicCluster>();
// Loop over tracks in event, extrapolate to layers and match hits if mips
// List<PerfectTrack> evtracks = event.get(PerfectTrack.class, "PerfectTracks");
// Get maps linking IL and mip clusters to tracks
Map<PerfectTrack, BasicCluster> trmipmap = (Map<PerfectTrack, BasicCluster>) event.get("TrackMipMap");
Map<BasicCluster, Integer> clusilmap = (Map<BasicCluster, Integer>) event.get("MipClusILMap");
double TotTrP = 0;
+ double[] trmclth = new double[50];
+ double[] trmclph = new double[50];
for (PerfectTrack itr : trmipmap.keySet())
{
double TrClE = 0;
+ int nshclmtch = 0;
BasicCluster trclus = new BasicCluster();
+ BasicCluster trshclus = new BasicCluster();
int niter = 0;
// get track and associated mip cluster
// if (trmipmap.isEmpty()) break;
@@ -175,7 +181,7 @@
BasicCluster trmclus = trmipmap.get(itr);
trclus.addCluster(trmclus);
TrClE += trmclus.getEnergy();
- Integer ClIL = clusilmap.get(trmclus);
+ int ClIL = clusilmap.get(trmclus);
// System.out.println("InteractionLayer for track " + ClIL);
double[] trmpos = trmclus.getPosition();
Hep3Vector mclpos = new BasicHep3Vector(trmpos);
@@ -196,22 +202,24 @@
double trphi = 0;
double ctrth = 0;
double ctrph = 0;
+ double clth = 0;
+ double clph = 0;
if (tobrad<toecz) // in barrel
{
SpacePoint trSP = tswim.getPointAtDistance(tobrad);
- trtheta += trSP.theta();
- trphi += trSP.phi();
- ctrth += Math.atan(trSP.rxy()/trSP.z());
+ trtheta = trSP.theta();
+ trphi = trSP.phi();
+ ctrth = Math.atan(trSP.rxy()/trSP.z());
if (ctrth<0) ctrth+=Math.PI;
- ctrph += Math.atan2(trSP.y(),trSP.x());
+ ctrph = Math.atan2(trSP.y(),trSP.x());
if (ctrph<0) ctrph+=2*Math.PI;
// aida.cloud1D("Theta of Track at IL B").fill(trtheta);
aida.cloud1D("Theta of Track c at IL B").fill(ctrth); // this is correct value
double clr = Math.sqrt(trmpos[0]*trmpos[0]+trmpos[1]*trmpos[1]);
- double clth = Math.atan(clr/trmpos[2]);
+ clth = Math.atan(clr/trmpos[2]);
if (clth<0) clth+=Math.PI;
- double clph = Math.atan2(trmpos[1],trmpos[0]);
+ clph = Math.atan2(trmpos[1],trmpos[0]);
if (clph<0) clph+=2*Math.PI;
aida.cloud1D("Theta of mip cluster c at clus B").fill(Math.atan(clr/trmpos[2]));
aida.cloud1D("Phi of Track at IL B").fill(trphi);
@@ -222,17 +230,17 @@
} else // in endcap
{
SpacePoint trSP = tswim.getPointAtDistance(toecz);
- trtheta += trSP.theta();
- trphi += trSP.phi();
- ctrth += Math.atan(trSP.rxy()/trSP.z());
+ trtheta = trSP.theta();
+ trphi = trSP.phi();
+ ctrth = Math.atan(trSP.rxy()/trSP.z());
if (ctrth<0) ctrth+=Math.PI;
- ctrph += Math.atan2(trSP.y(),trSP.x());
+ ctrph = Math.atan2(trSP.y(),trSP.x());
if (ctrph<0) ctrph+=2*Math.PI;
aida.cloud1D("Theta of Track c at IL EC").fill(ctrth);
double clr = Math.sqrt(trmpos[0]*trmpos[0]+trmpos[1]*trmpos[1]);
- double clth = Math.atan(clr/trmpos[2]);
+ clth = Math.atan(clr/trmpos[2]);
if (clth<0) clth+=Math.PI;
- double clph = Math.atan2(trmpos[1],trmpos[0]);
+ clph = Math.atan2(trmpos[1],trmpos[0]);
if (clph<0) clph+=2*Math.PI;
aida.cloud1D("Theta of mip cluster c at clus EC").fill(Math.atan(clr/trmpos[2]));
aida.cloud1D("Phi of Track at IL EC").fill(trphi);
@@ -241,9 +249,10 @@
aida.cloud2D("Theta vs Phi of Track at IL EC").fill(ctrph,ctrth);
aida.cloud2D("Theta vs Phi of mip cluster at IL EC").fill(clph,clth);
}
- aida.cloud2D("TrackMip Theta vs Phi Extrap").fill(trphi,trtheta);
+ trmclth[0] = ctrth; // 0 element is track
+ trmclph[0] = ctrph;
// preselection of clusters near this track - first make a list of clusters that
- // are candidates to be associated with this track
+ // are candidates to be associated with this track - not used
List<BasicCluster> candshclus = new ArrayList<BasicCluster>();
for (BasicCluster ibcl : showclus)
{
@@ -269,6 +278,12 @@
}
}
aida.cloud1D("Number of Cand Clusters added to track").fill(candshclus.size());
+ double cclE = 0;
+ for (BasicCluster iccl : candshclus)
+ {
+ cclE += iccl.getEnergy();
+ }
+ if (cclE>0) aida.cloud1D("E over p of all Cand Clusters added to Track").fill(cclE/TrP);
event.put("CandidateShowerClusters",candshclus);
// check all clusters for matches using spacepoint of extrapolated track and cluster coords
// for each cluster, loop over all hits in the cluster, checking if any hits are close to track
@@ -297,7 +312,24 @@
double dtrhph = Math.abs(ctrph-htph);
if (dtrhph>Math.PI) dtrhph = 2*Math.PI-dtrhph;
double distrh = Math.sqrt(dtrhth*dtrhth+dtrhph*dtrhph);
- if (distrh<niter*0.0075) nhmtch++;
+ // first dip in dist at 0.075(6*0.0125, second dip at 0.15 (6*.025)
+ double[] distrcl = new double[50];
+ // calculate all distances
+ for (int i=0; i<nshclmtch+1; ++i)
+ {
+ double dth = Math.abs(trmclth[i]-htth);
+ double dph = Math.abs(trmclph[i]-htph);
+ if (dph>Math.PI) dph = 2*Math.PI-dph;
+ distrcl[i] = Math.sqrt(dth*dth+dph*dph);
+ }
+ for (int j=0; j<nshclmtch+1; ++j)
+ {
+ if (distrcl[j]<_mindist)
+ {
+ nhmtch++;
+ continue;
+ }
+ }
}
double[] clpos = ishcl.getPosition();
double cph = Math.atan2(clpos[1],clpos[0]);
@@ -317,113 +349,38 @@
aida.cloud1D("Track Shower dist").fill(dist);
aida.cloud1D("Number of hits matched to Track").fill(nhmtch);
// if (dist < niter*0.015)
- if (nhmtch>5)
+ if (nhmtch>0)
{
trclus.addCluster(ishcl);
+ trshclus.addCluster(ishcl);
aida.cloud2D("Track Shower Match Cluster Theta Phi").fill(cph,cth);
aida.cloud2D("Track Shower Match Track Theta Phi").fill(ctrph,ctrth);
- aida.cloud1D("Number of matched hits for matched CLusters").fill(nhmtch);
+ aida.cloud1D("Number of matched hits for matched Clusters").fill(nhmtch);
// try changing theta, phi to average of trackcluster coord if matched
showclE += ishcl.getEnergy();
- ctrph = (TrP*ctrph+showclE*cph)/(TrP+showclE);
- ctrth = (TrP*ctrth+showclE*cth)/(TrP+showclE);
+ nshclmtch++;
+ trmclth[nshclmtch] = cth;
+ trmclph[nshclmtch] = cph;
+// ctrph = (TrP*ctrph+showclE*cph)/(TrP+showclE);
+// ctrth = (TrP*ctrth+showclE*cth)/(TrP+showclE);
shcl.remove();
// System.out.println("Cluster matched to track");
}
}
TrClE += showclE;
if (TrClE>0) aida.cloud1D("E over P Track Shower matches").fill(TrClE/TrP);
- if (niter == 6) break;
- } while (TrClE/TrP<0.55);
- if (TrClE>0) aida.cloud1D("Number of iterations for Track").fill(niter);
- if (TrClE>0) aida.cloud1D("Final E over p Track Shower matches").fill(TrClE/TrP);
+ if (niter == _nloops) break;
+ } while (TrClE/TrP<_mineop);
+ if (TrClE>0) aida.cloud1D("Number of iterations per Track").fill(niter);
+ if (TrClE>0) aida.cloud1D("Number of Shower Clusters added to Track").fill(nshclmtch);
+ if (TrClE>0) aida.cloud1D("Final E over p Track Shower matches").fill(TrClE*1.3/TrP);
trkcalclus.add(trclus);
+ trkshoclus.add(trshclus);
}
if (trkcalclus.size()>0) event.put("TrackCALClusters",trkcalclus);
- // simple check for neutrals - based on min hits in cluster
- List<BasicCluster> neuclus = new ArrayList<BasicCluster>();
- for (BasicCluster shc : showclus)
- {
- aida.cloud1D("Number of hits in NeuHClusters").fill(shc.getCalorimeterHits().size());
- if (shc.getCalorimeterHits().size()>20) neuclus.add(shc);
- }
- if (neuclus.size()>0) event.put("NeuHClusters",neuclus);
-// System.out.println("Number Shower Clusters after Matching " + showclus.size());
-
- // make some plots
- List<ReconstructedParticle> prepars = event.get(ReconstructedParticle.class,"PerfectRecoParticles");
- double perfpho = 0;
- double perfneu = 0;
- double perfchpE = 0;
- double perfchpP = 0;
- for (ReconstructedParticle prepar : prepars)
- {
- double perfpCh = prepar.getCharge();
-// ParticleID perfpID = perfrp.getParticleIDUsed();
- int perfpID = Math.abs(prepar.getParticleIDUsed().getPDG());
-// System.out.println("Particle ID = " + perfpID);
- double perfpE = prepar.getEnergy();
- Hep3Vector perfpP = prepar.getMomentum();
- double pP = perfpP.magnitude();
- if (perfpCh == 0)
- {
- if (perfpID == 22)
- {
- perfpho += perfpE;
- }else
- {
- perfneu += perfpE;
- }
- } else
- {
- perfchpE += perfpE;
- perfchpP += pP;
- }
- }
- List<MCReconstructedParticle> mcpars = event.get(MCReconstructedParticle.class,"CheatReconstructedParticles");
- double trphoE = 0;
- for (MCReconstructedParticle mcpar : mcpars)
- {
- int mcID = Math.abs(mcpar.getParticleIDUsed().getPDG());
- double mcE = mcpar.getEnergy();
- if (mcID == 22) trphoE += mcE;
- }
- double TotTrClE = 0;
- for (BasicCluster trccl : trkcalclus)
- {
- TotTrClE += trccl.getEnergy();
- }
- aida.cloud1D("PFA Charged Hadron ESum").fill(TotTrClE);
- aida.cloud1D("PFA Total Track PSum").fill(TotTrP);
- double TotRemainE = 0;
- for (BasicCluster nc : neuclus)
- {
- TotRemainE += nc.getEnergy()*1.3;
- }
- aida.cloud1D("Neutral Cluster ESum").fill(TotRemainE);
- List<BasicCluster> phobcl = event.get(BasicCluster.class,"PhotonBClusters");
- List<BasicCluster> phoeccl = event.get(BasicCluster.class,"PhotonECClusters");
- double TotPhoE = 0;
- for (BasicCluster phob : phobcl)
- {
- TotPhoE += phob.getEnergy()*1.015*1.013;
- }
- for (BasicCluster phoec : phoeccl)
- {
- TotPhoE += phoec.getEnergy()*1.015*1.013;
- }
- aida.cloud1D("PFA Photon ESum").fill(TotPhoE);
- aida.cloud1D("PFA Track P + PhotonE").fill(TotTrP+TotPhoE);
- double PFAESum = TotTrP+TotPhoE+TotRemainE;
- aida.cloud1D("PFA TrP + PhoE + RemE").fill(TotTrP+TotPhoE+TotRemainE);
- aida.cloud1D("Difference Perf Pho PFA Pho").fill(perfpho-TotPhoE);
- aida.cloud1D("Difference True Pho Perf Pho").fill((trphoE-perfpho)/Math.sqrt(trphoE));
- aida.cloud2D("TrPho PerfPho vs PerfPho PFAPho").fill(trphoE-perfpho,perfpho-TotPhoE);
- aida.cloud1D("Difference Perf Neu PFA Neu").fill(perfneu-TotRemainE);
- aida.cloud1D("Difference Perf PFA PFA").fill(perfpho+perfneu-TotPhoE-TotRemainE);
- aida.cloud1D("Perf Pho Neu PFA Pho Neu").fill(perfpho+perfneu-TotPhoE-TotRemainE);
-
- } // end of filter loop
+ if (trkshoclus.size()>0) event.put("TrackShowerClusters",trkshoclus);
+ if (showclus.size()>0) event.put("ShowClustersLeft",showclus);
+
}
public void setClusterNames(String[] names)
