lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn/JetFinder
diff -u -r1.1 -r1.2
--- JetDriverExtended.java 29 Jun 2009 18:13:08 -0000 1.1
+++ JetDriverExtended.java 30 Jun 2009 18:44:56 -0000 1.2
@@ -26,7 +26,7 @@
* A simple driver which can be used to find jets from ReconstructedParticles.
* The resuslting jets are stored in a new collection of ReconstructedParticles.
* @author tonyj
- * @version $Id: JetDriverExtended.java,v 1.1 2009/06/29 18:13:08 mbussonn Exp $
+ * @version $Id: JetDriverExtended.java,v 1.2 2009/06/30 18:44:56 mbussonn Exp $
*/
public class JetDriverExtended extends Driver
{
@@ -145,7 +145,7 @@
//if (hist)
{
AIDA aida = AIDA.defaultInstance();
- aida.cloud1D("JetDriver/sumE").fill(totalEnergy);
+ aida.cloud1D("JDE/JetDriver/sumE").fill(totalEnergy);
double tje =0;
for(int i=0;i<finder.njets();i++){
HepLorentzVector hlv = finder.jet(i);
@@ -160,31 +160,31 @@
double phi = Math.atan2(px,py);
tje +=hlv.v3().magnitude() ;
if(1==finder.nParticlesPerJet(i)){
- aida.cloud1D("One particle/jet eta ").fill(eta);
- aida.cloud1D("One particle/jet phi ").fill(phi);
- aida.cloud1D("One particle/jet pt").fill(pt);
- aida.cloud1D("One particle/jet |H4|").fill(finder.jet(i).magnitude());
- aida.cloud1D("One particle/jet pr").fill(pr);
- aida.cloud1D("One particle/jet pz").fill(pz);
+ aida.cloud1D("JDE/One particle/jet eta ").fill(eta);
+ aida.cloud1D("JDE/One particle/jet phi ").fill(phi);
+ aida.cloud1D("JDE/One particle/jet pt").fill(pt);
+ aida.cloud1D("JDE/One particle/jet |H4|").fill(finder.jet(i).magnitude());
+ aida.cloud1D("JDE/One particle/jet pr").fill(pr);
+ aida.cloud1D("JDE/One particle/jet pz").fill(pz);
//aida.cloud2D("One particle/jet direction").fill(pz,pr);
}else{
- aida.cloud1D("Jet particle/jet eta ").fill(eta);
- aida.cloud1D("Jet particle/jet phi ").fill(phi);
- aida.cloud1D("Jet particle/jet pt").fill(pt);
- aida.cloud1D("Jet particle/jet |H4|").fill(finder.jet(i).magnitude());
- aida.cloud1D("Jet particle/jet pr").fill(pr);
- aida.cloud1D("Jet particle/jet pz").fill(pz);
+ aida.cloud1D("JDE/Jet particles/jet eta ").fill(eta);
+ aida.cloud1D("JDE/Jet particles/jet phi ").fill(phi);
+ aida.cloud1D("JDE/Jet particles/jet pt").fill(pt);
+ aida.cloud1D("JDE/Jet particles/jet |H4|").fill(finder.jet(i).magnitude());
+ aida.cloud1D("JDE/Jet particles/jet pr").fill(pr);
+ aida.cloud1D("JDE/Jet particles/jet pz").fill(pz);
//aida.cloud2D("Jet particle/jet direction").fill(pz,pr);
- aida.cloud2D("JetDriver/jet direction ").fill(hlv.v3().x(),hlv.v3().y());
- aida.cloud1D("JetDriver/particlesPerJet").fill(finder.nParticlesPerJet(i));
- aida.cloud1D("JetDriver/nJets").fill(finder.njets());
+ aida.cloud2D("JDE/JetDriver/jet direction ").fill(hlv.v3().x(),hlv.v3().y());
+ aida.cloud1D("JDE/JetDriver/particlesPerJet").fill(finder.nParticlesPerJet(i));
+ aida.cloud1D("JDE/JetDriver/nJets").fill(finder.njets());
}
}
- aida.cloud1D("jetDriver/total jet energy").fill(tje);
+ aida.cloud1D("JDE/jetDriver/total jet energy").fill(tje);
}
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn/JetFinder
diff -u -r1.1 -r1.2
--- JetFinder.java 29 Jun 2009 18:13:08 -0000 1.1
+++ JetFinder.java 30 Jun 2009 18:44:56 -0000 1.2
@@ -7,10 +7,9 @@
import hep.aida.IHistogram1D;
import hep.aida.IHistogramFactory;
import hep.aida.IProfile1D;
-import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
import java.io.IOException;
-import java.util.ArrayList;
+import java.util.HashSet;
import java.util.LinkedList;
import java.util.List;
import java.util.Set;
@@ -18,18 +17,15 @@
import java.util.logging.Logger;
import org.lcsim.contrib.Partridge.TrackingTest.TrackAnalysis;
import org.lcsim.contrib.mgraham.sATLASDigi.FindableTrack;
-import org.lcsim.contrib.mgraham.sATLASDigi.FindableTrack.Ignore;
-import org.lcsim.contrib.sATLAS.TrackReconstructionDriver;
import org.lcsim.event.EventHeader;
import org.lcsim.event.LCRelation;
import org.lcsim.event.MCParticle;
import org.lcsim.event.ReconstructedParticle;
import org.lcsim.event.RelationalTable;
import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
import org.lcsim.event.base.BaseRelationalTable;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
-import org.lcsim.recon.tracking.seedtracker.SeedStrategy;
-import org.lcsim.recon.tracking.seedtracker.StrategyXMLUtils;
import org.lcsim.recon.vertexing.zvtop4.VectorArithmetic;
import org.lcsim.util.Driver;
import org.lcsim.util.aida.AIDA;
@@ -54,7 +50,6 @@
gx = hf.createProfile1D("track finding efficiency vs angle ",35,0.0,.7);
gxnarrow = hf.createProfile1D("track finding efficiency vs angle(near zero) ",70,0.0,.2);
gy = hf.createHistogram1D("nummber of mcp in bin","this", 35, 0.0,.7);
-
}
System.out.println("creating gx "+gx);
}
@@ -71,7 +66,7 @@
{
System.out.println("processing "+this.getName());
List<ReconstructedParticle> jets = event.get(ReconstructedParticle.class,"Jets");
- aida.cloud1D("nJets").fill(jets.size());
+ //aida.cloud1D("nJets").fill(jets.size());
RelationalTable hittomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
RelationalTable trktomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
List<LCRelation> mcrelations = event.get(LCRelation.class, "HelicalTrackMCRelations");
@@ -104,7 +99,7 @@
trktomc.add(track, mcp);
}//</editor-fold>
-
+ int njets=0;
List<MCParticle> mcplist = new LinkedList<MCParticle>();
FindableTrack findable = new FindableTrack(event);
for(ReconstructedParticle jet : jets){
@@ -112,9 +107,9 @@
double jpy = jet.getMomentum().y();
double jpz = jet.getMomentum().z();
double jpt = Math.sqrt(jpx*jpx +jpy*jpy);
- if (jpt <150)
+ if (jpt <250)
continue;
-
+ njets++;
List<ReconstructedParticle> listOfParticles = jet.getParticles();
if(listOfParticles.size()<2){continue;}
aida.cloud2D("jet direction r z").fill(jpz,jpt);
@@ -126,6 +121,7 @@
}
Hep3Vector jetMomentum = jet.getMomentum();
int ntracks=0;
+ int nfindableMCP=0;
MCParticleExtended mcpx = new MCParticleExtended();
for(MCParticle mcp : mcplist){
mcpx.RecycleWithMCParticleAndEvent(mcp, event);//try not to allocate new object inside a tight loop = new MCParticleExtended(mcp,event);
@@ -138,48 +134,65 @@
||mcp.getGeneratorStatus()!= MCParticleExtended.FINAL_STATE
||Math.abs(mcpx.getEta())>2.5)
{continue;}
- else{
- double angle = getAngle(mmt,jetMomentum);
- double wgt = 0.0;
- aida.cloud1D("angle").fill(angle);
- if(setOfTrack.size()>1){
- wgt = 0.0;
- aida.cloud1D("track with more than 1 mcp").fill(setOfTrack.size()); continue;
- }
- else if(setOfTrack.size()==0) {
- wgt = 0.0;
- aida.cloud1D("not found/nhits").fill(findable.LayersHit(mcp));
- aida.cloud1D("not found/eta") .fill(mcpx.getEta());
- }else if (setOfTrack.size()==1){
- ntracks++;
- wgt = 1.0; //<editor-fold desc="//some useless graphs">
- Track track = setOfTrack.iterator().next();
- double d0track = track.getTrackParameter(HelicalTrackFit.dcaIndex);
+ nfindableMCP++;
+ double angle = getAngle(mmt,jetMomentum);
+ double wgt = 0.0;
+ aida.cloud1D("all/angle").fill(angle);
+ if(setOfTrack.size()>1){
+ aida.cloud1D("track with more than 1 mcp").fill(setOfTrack.size());
+ continue;
+ }
+ else if(setOfTrack.size()==0) {
+ wgt = 0.0;
+ aida.cloud1D("not found/nhits").fill(findable.LayersHit(mcp));
+ aida.cloud1D("not found/eta") .fill(mcpx.getEta());
+ aida.cloud1D("not found/angle").fill(angle);
+ }else if (setOfTrack.size()==1){
+ ntracks++;
+ wgt = 1.0;
+ Track track = setOfTrack.iterator().next();
+ double d0track = track.getTrackParameter(HelicalTrackFit.dcaIndex);
double d0mcp = mcpx.getDCA();
double z0track = track.getTrackParameter(HelicalTrackFit.z0Index);
double z0mcp = mcpx.getZ0();
-
- double px = track.getPX();
- double py = track.getPY();
- double pz = track.getPZ();
- double ps = Math.sqrt(px*px+py*py+pz*pz);
+ double px = track.getPX();
+ double py = track.getPY();
+ double pz = track.getPZ();
+ double ps = Math.sqrt(px*px+py*py+pz*pz);
double pmcp= mmt.magnitude();
-
- aida.cloud2D("residu/d0-z0").fill(d0track-d0mcp,z0mcp-z0track);
- aida.cloud1D("residu/momentum").fill(pmcp-ps);
- aida.cloud1D("residu/number of hits").fill(findable.LayersHit(mcp)-track.getTrackerHits().size());
- }
- gy.fill(angle);
- gx.fill(angle,wgt);
- gxnarrow.fill(angle,wgt);
- aida.cloud2D("nTrack vs nMCP in jet").fill(listOfParticles.size(), ntracks);
+ //let's get the purity
+ Set<TrackerHit> trackHits = new HashSet(track.getTrackerHits()) ;
+ Set<TrackerHit> MCPHits = new HashSet(hittomc.allTo(mcp));
+ Set intersection = new HashSet(trackHits);
+ intersection.retainAll(MCPHits);
+ double trackp = (double)intersection.size()/(double)trackHits.size();
+ double mcpp = (double)intersection.size()/(double)MCPHits.size();
+
+ aida.cloud2D("residual/d0-z0").fill(d0track-d0mcp,z0mcp-z0track);
+ aida.cloud1D("residual/momentum").fill(pmcp-ps);
+ aida.cloud1D("residual/number of hits").fill(findable.LayersHit(mcp)-track.getTrackerHits().size());
+ aida.cloud1D("residual/chi squared").fill(track.getChi2());
+ aida.cloud1D("track/purity (intersect track)").fill(trackp);
+ aida.cloud1D("track/purity (intersect MCP)") .fill(mcpp);
+
+ }
+ aida.cloud1D("angle").fill(angle);
+ gy.fill(angle);
+ gx.fill(angle,wgt);
+ gxnarrow.fill(angle,wgt);
- }
}//end loop through MCP
- mcplist.clear();
- }
+ if(nfindableMCP!=0){
+ aida.cloud2D("nTrack vs nMCP in jet").fill(listOfParticles.size(), (double)ntracks/(double)listOfParticles.size());
+ aida.cloud2D("% of finded tracks vs nMCPfindable ").fill(nfindableMCP, (double)ntracks/(double)nfindableMCP);
+ }
+ mcplist.clear();
+ }//end looping through jets
+ aida.cloud1D("number of jets per event :").fill(njets);
+ }//end proces(event)
+
+
- }
private double getAngle(Hep3Vector h1, Hep3Vector h2){
double ca = VectorArithmetic.dot(h1, h2)/(h1.magnitude()*h2.magnitude());
return Math.acos(ca);
