lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -u -r1.2 -r1.3
--- ITupleWrapper.java 9 Jul 2009 21:46:17 -0000 1.2
+++ ITupleWrapper.java 14 Jul 2009 17:38:39 -0000 1.3
@@ -8,13 +8,12 @@
import hep.aida.ITree;
import hep.aida.ITuple;
import hep.aida.ITupleFactory;
-import java.lang.Integer;
import java.lang.String;
import java.util.HashMap;
import java.util.HashSet;
-import java.util.LinkedList;
import java.util.List;
import java.util.Set;
+import java.util.Vector;
//import org.lcsim.util.loop.LCIODriver;
/**
*
@@ -22,12 +21,18 @@
*/
public class ITupleWrapper {
private static HashMap entry = new HashMap<String,Entry>();
-
+ private static HashMap sizeMap = new HashMap<Number,Entry>();
public ITupleWrapper() {
}
-
-
+ /**
+ * Just for going easyly from aida to this wrapper
+ * @param name
+ * @return
+ */
+ public Entry Cloud1D(String name){
+ return this.tuple(name);
+ }
/**
* Create a tuple of name "name" if it don't exist, otherwithe return the tuple wit this name
* @param name
@@ -39,7 +44,6 @@
return (Entry) entry.get(name);
}
else{
- //System.out.println("create column"+name);
Entry newEntry = new Entry();
entry.put(name,newEntry);
return newEntry;
@@ -48,43 +52,42 @@
/**
* Write the instance(s) of the ITuple(s) in jas
* Create as many folder as differnts lengts of columns
- *
+ * typically called in "Supend(){}"
*/
public void dump(){
Set<Number> setOfSize = new HashSet<Number>();
for(Object key: entry.keySet() ){
Entry ent = (Entry) entry.get(key);
- // System.out.println("key "+key);
setOfSize.add(ent.size());
+ sizeMap.put(ent.size(), ent);
}
+
IAnalysisFactory af = IAnalysisFactory.create();
ITree tree = af.createTreeFactory().create();
ITupleFactory tf = af.createTupleFactory(tree);
- //System.out.println("set of sizes"+setOfSize);
+
for(Number i : setOfSize){
int ii = i.intValue();
- //System.out.println("size = "+ii);
+
String columnString = "int firstcol";
for(Object key: entry.keySet() ){
Entry ent = (Entry) entry.get(key);
if(ent.size()!=ii){
- // System.out.println("skipping entry "+key+"for number of value "+ii+" because the length is "+ent.size());
+
continue;
}
columnString = columnString+","+ent.type()+key;
}
- // System.out.println("dumping columns:"+columnString);
ITuple tuple = tf.create("tuple"+i, "MyNtule"+i, columnString);
int j;
for(j=0;j<ii;j++ ){
- tuple.addRow();
int k = 0;
tuple.fill(k,0);
for(Object key: entry.keySet() ){
Entry ent = (Entry) entry.get(key);
- if(ent.size()!=i)continue;
+ if(ent.size()!=ii)continue;
k++;
if(ent.type().equals("int "))
{tuple.fill(k,ent.ValueAtIndex(j).intValue());}
@@ -92,7 +95,7 @@
{tuple.fill(k,ent.ValueAtIndex(j).doubleValue());}
}
-
+ tuple.addRow();
}
@@ -102,12 +105,20 @@
}
-
+//****************************************************************************//
+//****************************************************************************//
+//****************************************************************************//
+//****************************************************************************//
+//****************************************************************************//
+//****************************************************************************//
+//****************************************************************************//
public class Entry{
- private List list = new LinkedList();
+ private List list = new Vector(0,100);
private boolean typeDecided = false;
private String type = "";
- Entry(){};
+ private int taille=0;
+ Entry(){
+ };
/**
@@ -116,7 +127,7 @@
* @return
*/
public int size(){
- return list.size();
+ return taille;
}
/**
* localy used to know if the type of value is int or double
@@ -125,7 +136,12 @@
public String type(){
return type;
}
-
+ public void fill(int value){
+ this.fillWithInt(value);
+ }
+ public void fill(double value){
+ this.fillWithDouble(value);
+ }
public void fillWithInt(int value){
this.fillWithValueAndType(value,"int ");
}
@@ -151,6 +167,7 @@
//throw exeption
System.out.println("wrong type");
}
+ taille++;
}
public void description(){
System .out.println(" array:"+list.size()+" element 0 : "+list.get(0));
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn/JetFinder
diff -u -r1.5 -r1.6
--- JetFinder.java 7 Jul 2009 21:20:15 -0000 1.5
+++ JetFinder.java 14 Jul 2009 17:38:39 -0000 1.6
@@ -3,22 +3,13 @@
* To change this template, choose Tools | Templates
* and open the template in the editor.
*/
-import hep.aida.IAnalysisFactory;
import org.lcsim.contrib.Mbussonn.*;
-import hep.aida.IHistogram1D;
-import hep.aida.IHistogramFactory;
import hep.aida.IProfile1D;
-import hep.aida.ITree;
-import hep.aida.ITuple;
-import hep.aida.ITupleFactory;
import hep.physics.vec.Hep3Vector;
-import java.io.IOException;
import java.util.HashSet;
-import java.util.LinkedList;
import java.util.List;
import java.util.Set;
-import java.util.logging.Level;
-import java.util.logging.Logger;
+import java.util.Vector;
import org.lcsim.contrib.Partridge.TrackingTest.TrackAnalysis;
import org.lcsim.contrib.mgraham.sATLASDigi.FindableTrack;
import org.lcsim.event.EventHeader;
@@ -41,50 +32,24 @@
* @author matthiasbussonnier
*/
public class JetFinder extends Driver{
- private AIDA aida = AIDA.defaultInstance();
- IAnalysisFactory af = IAnalysisFactory.create();
- ITree tree = af.createTreeFactory().create();
- ITupleFactory tf = af.createTupleFactory(tree);
- private static int JetZDirection = 0;
- private static int JetTDirection = 1;
- private static int NumberOfParticlePerJet = 2;
- private static int All_Angle = 3;
- private static int NumberOfTrackWitMoreThanOneMCP = 4;
- private static int NotFound_numberOfHits = 5;
- private static int NotFound_Eta = 6;
- private static int NotFound_Angle = 7;
-
-
- //aida.cloud2D("jet direction r z").fill(jpz,jpt);
- //aida.cloud1D("nParticles").fill(listOfParticles.size());
-
- String columnString = "double jetPz,"+
- "double jetPt," +
- "int numberOfParticlesPerJet," +
- "double allparticles_Angle," +
- "int numberOfTrackWithMoreThanOneMcp," +
- "double notfoundnumberofhits," +
- "double notfoudeta," +
- "double notfoundangle";
- ITuple tuple = tf.create("tuple", "MyNtule", columnString);
+ //private AIDA aida = AIDA.defaultInstance();
+
+ private static ITupleWrapper itw= null;
+
private static JetFinder defaultInstance;
public String outputPlots="myplots.aida";
private IProfile1D gx=null;
- private IProfile1D gxnarrow=null;
- private IHistogram1D gy= null;
+ //private IProfile1D gxnarrow=null;
+ //private IHistogram1D gy= null;
private JetFinder()
{
System.out.println("=============================================");
System.out.println("= creating jetfinderExtended =");
System.out.println("=============================================");
- IHistogramFactory hf = aida.histogramFactory();
- if(gx == null){
- 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);
+ if(itw == null){
+ itw = new ITupleWrapper();
+ }
}
public static JetFinder defaultInstance()
{
@@ -97,10 +62,9 @@
@SuppressWarnings("unchecked")
protected void process(EventHeader event)
{
- System.out.println("=============================================");
- System.out.println("=processing event "+event.getEventNumber());
- System.out.println("=============================================");
- System.out.println("processing "+this.getName());
+ if(event.getEventNumber()%100==0)
+ System.out.println("= "+this.getName()+" processing event number "+event.getEventNumber());
+
List<ReconstructedParticle> jets = event.get(ReconstructedParticle.class,"Jets");
//aida.cloud1D("nJets").fill(jets.size());
RelationalTable hittomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
@@ -121,7 +85,6 @@
{continue;}
rc2mc.add(relation.getFrom(), relation.getTo());
}
- System.out.println("nombre de relation :"+rc2mc.size());
for (Track track : tracklist) {//<editor-fold desc="//construc a relation between track and mcp">
TrackAnalysis tkanal = new TrackAnalysis(track, hittomc);
@@ -142,10 +105,10 @@
}//</editor-fold>
int njets=0;
- List<MCParticle> mcplist = new LinkedList<MCParticle>();
+ List<MCParticle> mcplist = new Vector<MCParticle>();
FindableTrack findable = new FindableTrack(event);
for(ReconstructedParticle jet : jets){
- System.out.println("processing jet...");
+
double jpx = jet.getMomentum().x();
double jpy = jet.getMomentum().y();
double jpz = jet.getMomentum().z();
@@ -155,12 +118,13 @@
njets++;
List<ReconstructedParticle> listOfParticles = jet.getParticles();
if(listOfParticles.size()<2){continue;}
- tuple.addRow();
- aida.cloud2D("jet direction r z").fill(jpz,jpt);
- tuple.fill(JetZDirection,jpz);
- tuple.fill(JetTDirection,jpt);
- aida.cloud1D("nParticles").fill(listOfParticles.size());
- tuple.fill(NumberOfParticlePerJet,listOfParticles.size());
+ // tuple.addRow();
+ //aida.cloud2D("jet direction r z").fill(jpz,jpt);
+ itw.tuple("jetZMomentum").fillWithDouble(jpz);
+ itw.tuple("jetTMomentum").fillWithDouble(jpt);
+
+ //aida.cloud1D("nParticles").fill(listOfParticles.size());
+ itw.tuple("numberOfParticlePerJets").fillWithInt(listOfParticles.size());
for(ReconstructedParticle rcpInJet : listOfParticles){
MCParticle mcpp = (MCParticle) rc2mc.to(rcpInJet);
if(mcpp != null)
@@ -170,15 +134,15 @@
int ntracks=0;
int nfindableMCP=0;
MCParticleExtended mcpx = new MCParticleExtended();
- System.out.println("processing "+mcplist.size()+" mcp");
+
for(MCParticle mcp : mcplist){
- tuple.addRow();
+
mcpx.RecycleWithMCParticleAndEvent(mcp, event);//try not to allocate new object inside a tight loop = new MCParticleExtended(mcp,event);
Set<Track> setOfTrack = trktomc.allTo(mcp);
Hep3Vector mmt = mcp.getMomentum();
double ptotal = mcpx.getPTotal();
- double ptcut = 0.0;
+ double ptcut = 1.0;
if( ptotal < ptcut
||findable.LayersHit(mcp)<7
||mcp.getGeneratorStatus()!= MCParticleExtended.FINAL_STATE
@@ -189,68 +153,95 @@
nfindableMCP++;
double angle = getAngle(mmt,jetMomentum);
double wgt = 0.0;
- aida.cloud1D("all/angle").fill(angle);
- tuple.fill(All_Angle,angle);
+ double d0track = -10;
+ double d0mcp = -10;
+ double z0track = -10;
+ double z0mcp = -10;
+ double px = 0;
+ double py = 0;
+ double pz = 0;
+ double pt =-10;
+ double ps =-10;
+ double pmcp = -10 ;
+ double trackp = -10;
+ double mcpp =-10;
+
+
+
+ if(Math.abs(mcpx.getEta())>2.1)continue;
if(setOfTrack.size()>1){
- aida.cloud1D("track with more than 1 mcp").fill(setOfTrack.size());
- tuple.fill(NumberOfTrackWitMoreThanOneMCP,setOfTrack.size());
+
+ itw.tuple("numberOftracksWithMoreThanOneMcp").fillWithInt(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);
- tuple.fill(NotFound_numberOfHits,findable.LayersHit(mcp));
- tuple.fill(NotFound_Eta,mcpx.getEta());
- tuple.fill(NotFound_Angle,angle);
+ itw.tuple("notFoundsNumberOfHits").fillWithInt(findable.LayersHit(mcp));
+ itw.tuple("notFoundsEta").fillWithDouble(mcpx.getEta());
+ itw.tuple("notFounds-Angle").fillWithDouble(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 pmcp= mmt.magnitude();
+ d0track = track.getTrackParameter(HelicalTrackFit.dcaIndex);
+ d0mcp = mcpx.getDCA();
+ z0track = track.getTrackParameter(HelicalTrackFit.z0Index);
+ z0mcp = mcpx.getZ0();
+ px = track.getPX();
+ py = track.getPY();
+ pz = track.getPZ();
+ pt = Math.sqrt(px*px+py*py);
+ ps = Math.sqrt(px*px+py*py+pz*pz);
+ pmcp= mmt.magnitude();
//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);
+ trackp = (double)intersection.size()/(double)trackHits.size();
+ mcpp = (double)intersection.size()/(double)MCPHits.size();
+
}
- aida.cloud1D("angle").fill(angle);
- //tuple.fill(2, angle);
- gy.fill(angle);
- gx.fill(angle,wgt);
- gxnarrow.fill(angle,wgt);
-
+ itw.tuple("track_d0").fillWithDouble(d0track);
+ itw.tuple("track_z0").fillWithDouble(z0track);
+ itw.tuple("mcp_d0").fillWithDouble(d0mcp);
+ itw.tuple("mcp_z0").fillWithDouble(z0mcp);
+ itw.tuple("track_pt").fillWithDouble(pt);
+ itw.tuple("track_pz").fillWithDouble(pz);
+ itw.tuple("mcp_pt").fillWithDouble(mcpx.getPR());
+ itw.tuple("mcp_pz").fillWithDouble(mcpx.getPZ());
+ itw.tuple("mcp_angle").fillWithDouble(angle);
+ itw.tuple("mcp_eta").fill(mcpx.getEta());
+ itw.tuple("mcp_theta").fill(mcpx.getTheta());
+ itw.tuple("mcp_phi").fill(mcpx.getPhi());
+ itw.tuple("track_purity").fill(trackp);
+ itw.tuple("hasThisMcpBeenFounded").fillWithDouble(wgt);
+ itw.tuple("angleForMcpWith0OrOneTrack").fillWithDouble(angle);
+
}//end loop through MCP
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);
+ itw.tuple("nTrackInJet").fill(mcplist.size());
+ itw.tuple("nMCPInJet").fill(ntracks);
+ itw.tuple("nFIndableMCP").fill(nfindableMCP);
+
+
}
mcplist.clear();
}//end looping through jets
- aida.cloud1D("number of jets per event :").fill(njets);
+ itw.tuple("numberOfJetsPerEvent").fill(njets);
+ //aida.cloud1D("number of jets per event :").fill(njets);
}//end proces(event)
-
+ @Override
+ public void suspend() {
+ System.out.println("=================");
+ System.out.println("= suspend =");
+ System.out.println("=================");
+ itw.dump();
+ }
private double getAngle(Hep3Vector h1, Hep3Vector h2){
double ca = VectorArithmetic.dot(h1, h2)/(h1.magnitude()*h2.magnitude());