lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -u -r1.3 -r1.4
--- JetFinder.java 8 Jun 2009 18:44:44 -0000 1.3
+++ JetFinder.java 19 Jun 2009 19:15:23 -0000 1.4
@@ -6,6 +6,7 @@
import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
import java.util.ArrayList;
+import java.util.LinkedList;
import java.util.List;
import java.util.Set;
import org.lcsim.contrib.Partridge.TrackingTest.TrackAnalysis;
@@ -19,11 +20,10 @@
import org.lcsim.event.RelationalTable;
import org.lcsim.event.Track;
import org.lcsim.event.base.BaseRelationalTable;
-import org.lcsim.event.util.JetDriver;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
-import org.lcsim.fit.helicaltrack.HelixParamCalculator;
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;
@@ -92,82 +92,55 @@
List<SeedStrategy> slist = StrategyXMLUtils.getStrategyListFromResource(sfile);
FindableTrack findable = new FindableTrack(event);
- List<MCParticle> mcplist = event.getMCParticles();
- System.out.println("mcp in event : "+mcplist.size());
+ List<MCParticle> mcplist = new LinkedList<MCParticle>();//event.getMCParticles();
+
//mcplist.clear();
for(ReconstructedParticle jet : jets){
+ double jpx = jet.getMomentum().x();
+ double jpy = jet.getMomentum().y();
+ double jpz = jet.getMomentum().z();
+ double jpt = Math.sqrt(jpx*jpx +jpy*jpy);
+ if (jpt <150)continue;
+ aida.cloud1D("jet direction r z").fill(jpz,jpt);
List<ReconstructedParticle> listOfParticles = jet.getParticles();
for(ReconstructedParticle rcpInJet : listOfParticles){
MCParticle mcpp = (MCParticle) rc2mc.to(rcpInJet);
System.out.println("one more mcp");
//MCParticleExtended mcpe = new MCParticleExtended(mcpp,event);
- //mcplist.add(mcpp);
+ if(mcpp != null)
+ mcplist.add(mcpp);
}
-
- }
- System.out.println("mcp in jets : "+mcplist.size());
- for(MCParticle mcpx : mcplist){
- MCParticleExtended mcp = new MCParticleExtended(mcpx,event);
- Set<Track> setOfTrack= trktomc.allTo(mcpx);
- Set<Track> setOfTx = trktomc.allTo(mcp);
- double ptotal = mcp.getPTotal();
- double ptcut=50.;
- double nrcp = rc2mc.allTo(mcpx).size();
- aida.cloud1D("rcp per mcp").fill(nrcp);
- if( ptotal < ptcut
+ Hep3Vector jetMomentum = jet.getMomentum();
+ System.out.println("mcp in jet : "+mcplist.size());
+ int ntracks=0;
+ for(MCParticle mcp : mcplist){
+ MCParticleExtended mcpx = new MCParticleExtended(mcp,event);
+ Set<Track> setOfTrack= trktomc.allTo(mcp);
+ Hep3Vector mmt = mcp.getMomentum();
+ double angle = getAngle(jetMomentum,mmt);
+ aida.cloud1D("angle").fill(angle);
+ double ptotal = mcpx.getPTotal();
+ double ptcut=50.;
+ if( ptotal < ptcut
//||(mcp.getPDGID()!=13 && mcp.getPDGID()!=-13)
||mcp.getGeneratorStatus()!= MCParticleExtended.FINAL_STATE
- ||Math.abs(mcp.getEta())>2.5)
+ ||Math.abs(mcpx.getEta())>2.5)
{continue;}
- else{
- //if(!findable.isFindable(mcpx, slist)){
- if(false){
- System.out.println("muon non trouvable...");
- List<Ignore> ignores = new ArrayList<Ignore>();
- ignores.add(Ignore.NoZ0Cut);
- if (findable.isFindable(mcpx, slist,ignores)){
- System.out.println("not findable because of Z0");
- continue;
- }
- ignores.add(Ignore.NoDCACut);
- if (findable.isFindable(mcpx, slist,ignores)){
- System.out.println("not findable because of DCA");
- continue;
- }
- ignores.add(Ignore.NoMinHitCut);
- if (findable.isFindable(mcpx, slist,ignores)){
- System.out.println("not findable because of min hit");
- continue;
- }
- ignores.add(Ignore.NoSeedCheck);
- if (findable.isFindable(mcpx, slist,ignores)){
- System.out.println("not findable because of seed check");
- continue;
- }
- ignores.add(Ignore.NoConfirmCheck);
- if (findable.isFindable(mcpx, slist,ignores)){
- System.out.println("not findable because of no confirm check");
- continue;
- }
- ignores=null;
- System.out.println("not findable at all... why ?");
- }
- double weight=0.;
-
- if(setOfTrack.size()>1){
- aida.cloud1D("track with more than 1 mcp").fill(setOfTrack.size()); continue;
+ else{
+ if(setOfTrack.size()>1){
+ aida.cloud1D("track with more than 1 mcp").fill(setOfTrack.size()); continue;
}
- if(setOfTrack.size()==0) {
- weight=0.;
- int nhits = findable.LayersHit(mcpx);
- aida.cloud1D("mcp/eta for non-finded particle").fill(mcp.getEta());
- aida.cloud1D("mcp/theta for non-finded particle").fill(mcp.getTheta());
- aida.cloud1D("nhits for non-finded muons").fill(nhits);
- }
- else if (setOfTrack.size()==1){
- System.out.println("tracktomcp : "+setOfTx.size()+" Vs "+setOfTrack.size());
+ else if(setOfTrack.size()==0) {
+ double weight=0.;
+ int nhits = findable.LayersHit(mcpx);
+ aida.cloud1D("mcp/eta for non-finded particle").fill(mcpx.getEta());
+ aida.cloud1D("mcp/theta for non-finded particle").fill(mcpx.getTheta());
+ aida.cloud1D("nhits for non-finded muons").fill(nhits);
+ }else if (setOfTrack.size()==1){
+ ntracks++;
+ //System.out.println("tracktomcp : "+setOfTx.size()+" Vs "+setOfTrack.size());
System.out.println("rcp associed with mcp "+rc2mc.allTo(mcpx).size());
- weight=1;
+ double weight=1;
Track track = setOfTrack.iterator().next();
double d0 = track.getTrackParameter(HelicalTrackFit.dcaIndex);
double z0 = track.getTrackParameter(HelicalTrackFit.z0Index);
@@ -175,7 +148,11 @@
double py = track.getPY();
double pz = track.getPZ();
double pt = Math.sqrt(px*px+py*py);
- double ps = Math.sqrt(px*px+py*py+pz*pz);
+ double ps = Math.sqrt(px*px+py*py+pz*pz);
+ double p[] = track.getMomentum();
+ Hep3Vector h3vp = new BasicHep3Vector(p[0],p[1],p[2]);
+ double trackangle = getAngle(h3vp,jetMomentum);
+ aida.cloud1D("track angle").fill(trackangle);
// double bfield = event.getDetector().getFieldMap().getField(new BasicHep3Vector(0., 0., 0.)).z();
// HelixParamCalculator helix = new HelixParamCalculator(mcp, bfield);
// double d0mc = helix.getDCA();
@@ -184,23 +161,27 @@
aida.cloud1D("track/z0").fill(z0);
aida.cloud2D("track/do VS Z0").fill(z0, d0);
aida.cloud1D("track/tr momentum").fill(pt);
- aida.cloud1D("track/scalar momentum").fill(ps);
- aida.cloud1D("track/residual t momentum").fill(mcp.getPR()-pt);
- aida.cloud1D("track/residual d0").fill(d0-mcp.getDCA());
- aida.cloud1D("track/residual z0").fill(z0-mcp.getZ0());
+ //aida.cloud1D("track/scalar momentum").fill(ps);
+ aida.cloud1D("track/residual t momentum").fill(mcpx.getPR()-pt);
+ aida.cloud1D("track/residual d0").fill(d0-mcpx.getDCA());
+ aida.cloud1D("track/residual z0").fill(z0-mcpx.getZ0());
aida.cloud1D("track/number of Tracker hits").fill(track.getTrackerHits().size());
aida.cloud1D("track/number of mcp hits").fill(findable.LayersHit(mcp));
aida.cloud1D("track/chiSquared").fill(track.getChi2());
+ }//ed loop through mcp
+ aida.cloud2D("ntrack vs nMCP in jet").fill(mcplist.size(),ntracks);
+ aida.cloud2D("nTrack vs nRCP in jet").fill(listOfParticles.size(), ntracks);
}
-
-
-
- }
- }//ENDFOR
+ }
+ mcplist.clear();
+ }
super.printStatistics(System.out);
}
-
+ private double getAngle(Hep3Vector h1, Hep3Vector h2){
+ double ca = VectorArithmetic.dot(h1, h2)/(h1.magnitude()*h2.magnitude());
+ return Math.acos(ca);
+ }
}
