5 modified files
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -u -r1.2 -r1.3
--- MuonsBoth.java 23 Jun 2009 01:17:46 -0000 1.2
+++ MuonsBoth.java 13 Aug 2009 20:40:12 -0000 1.3
@@ -9,6 +9,7 @@
import hep.aida.IProfile1D;
import java.io.IOException;
import java.util.ArrayList;
+import java.util.HashSet;
import java.util.List;
import java.util.Set;
import java.util.logging.Level;
@@ -25,6 +26,7 @@
import org.lcsim.event.MCParticle;
import org.lcsim.event.RelationalTable;
import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
import org.lcsim.recon.tracking.seedtracker.*;
import org.lcsim.util.aida.AIDA;
@@ -38,13 +40,15 @@
private AIDA aida = AIDA.defaultInstance();
private IProfile1D efficiencyVSEta,avgEfficiencyVSEta;
private IProfile1D effMuons;
+ public ITupleWrapper itw;
private java.util.Random gtor;
public String outputPlots="myplots.aida";
public MuonsBoth(){
IHistogramFactory hf = aida.histogramFactory();
- efficiencyVSEta = hf.createProfile1D("efficiency vs eta",81,-2.025,2.025);
- avgEfficiencyVSEta= hf.createProfile1D("Avg efficiency vs eta",50,0.,2.5);
+ efficiencyVSEta = hf.createProfile1D("efficiency vs eta",103,-2.825,2.825);
+ avgEfficiencyVSEta= hf.createProfile1D("Avg efficiency vs eta",60,0.,2.8);
gtor = new java.util.Random(10101);
+ itw = new ITupleWrapper();
}
@Override
@@ -97,15 +101,23 @@
MCParticleExtended mcp = new MCParticleExtended(mcpx);
Set<Track> setOfTrack= trktomc.allTo(mcpx);
double ptotal = mcp.getPTotal();
+ aida.cloud1D("all particles").fill(mcp.getEta());
+ itw.tuple("all muons eta").fill(mcp.getEta());
double ptcut=50.;
if( ptotal < ptcut
||(mcp.getPDGID()!=13 && mcp.getPDGID()!=-13)
||mcp.getGeneratorStatus()!= MCParticleExtended.FINAL_STATE
- ||Math.abs(mcp.getEta())>2.5)
- {continue;}
+ //||Math.abs(mcp.getEta())>2.8
+ )
+ {
+ itw.tuple("not final state or not muons").fill(mcp.getEta());
+ continue;
+ }
else{
+ itw.tuple("eta for findable muons").fill(mcp.getEta());
if(!findable.isFindable(mcpx, slist)){
System.out.println("muon non trouvable...");
+ aida.cloud1D("eta for not findable muons").fill(mcp.getEta());
List<Ignore> ignores = new ArrayList<Ignore>();
ignores.add(Ignore.NoZ0Cut);
if (findable.isFindable(mcpx, slist,ignores)){
@@ -141,10 +153,22 @@
aida.cloud1D("track with more than 1 mcp").fill(setOfTrack.size()); continue;
}
if(setOfTrack.size()==0) {
- weight=0.;
+ System.out.println("======================================");
+ weight=0.;
int nhits = findable.LayersHit(mcpx);
aida.cloud1D("mcp/eta for non-finded muons").fill(mcp.getEta());
- aida.cloud1D("mcp/theta for non-finded muons").fill(mcp.getTheta()); aida.cloud1D("nhits for non-finded muons").fill(nhits);
+ itw.tuple("mcp/eta for non-finded muons").fill(mcp.getEta());
+ aida.cloud1D("mcp/theta for non-finded muons").fill(mcp.getTheta());
+ aida.cloud1D("nhits for non-finded muons").fill(nhits);
+ Set<TrackerHit> MCPHits = new HashSet(hittomc.allTo(mcpx));
+ //System.out.println("hit set size "+MCPHits.size());
+ for(TrackerHit hit : MCPHits)
+ {
+ System.out.println("hit :"+hit);
+ //System.out.println("rawhit : "+hit.getRawHits());
+
+ }
+ System.out.println("======================================");
}
else if (setOfTrack.size()==1){
weight=1;
@@ -196,12 +220,13 @@
}
efficiencyVSEta.fill(mcp.getEta(),weight);
aida.cloud1D("mcp/eta for detectable muons").fill(mcp.getEta());
+ itw.tuple("mcp/eta for detectable muons").fill(mcp.getEta());
aida.cloud1D("mcp/theta for detectable muons").fill(mcp.getTheta());
}
}//ENDFOR
- super.printStatistics(System.out);
+ //super.printStatistics(System.out);
aida.cloud1D("mcp/number of detectable muons").fill(numberOfMuons);
aida.cloud1D("mcp/number of detected muons").fill(numberOfMuonsWithTrack);
}
@@ -219,5 +244,11 @@
this.outputPlots=output;
}
+
+ @Override
+ public void suspend(){
+ itw.dump();
+ }
+
}//end class
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -u -r1.4 -r1.5
--- ITupleWrapper.java 20 Jul 2009 21:59:08 -0000 1.4
+++ ITupleWrapper.java 13 Aug 2009 20:40:12 -0000 1.5
@@ -8,16 +8,14 @@
import hep.aida.ITree;
import hep.aida.ITuple;
import hep.aida.ITupleFactory;
-import java.lang.String;
import java.util.HashMap;
import java.util.HashSet;
import java.util.List;
import java.util.Set;
import java.util.Vector;
-//import org.lcsim.util.loop.LCIODriver;
/**
*
- * @author matthiasbussonnier
+ * @author matthias bussonnier [log in to unmask]
*/
public class ITupleWrapper {
private static HashMap entry = new HashMap<String,Entry>();
@@ -68,38 +66,38 @@
for(Number i : setOfSize){
int ii = i.intValue();
-
- String columnString = "int firstcol";
+ String columnString = "";
for(Object key: entry.keySet() ){
Entry ent = (Entry) entry.get(key);
if(ent.size()!=ii){
-
continue;
}
- columnString = columnString+","+ent.type()+key;
+ columnString = columnString+ent.type()+key+",";
}
-
-
+ columnString = columnString.substring(0,columnString.length()-1);
+ System.out.println(columnString);
ITuple tuple = tf.create("tuple"+i, "MyNtule"+i, columnString);
int j;
for(j=0;j<ii;j++ ){
int k = 0;
- tuple.fill(k,0);
+
for(Object key: entry.keySet() ){
Entry ent = (Entry) entry.get(key);
if(ent.size()!=ii)continue;
- k++;
+
if(ent.type().equals("int "))
- {tuple.fill(k,ent.ValueAtIndex(j).intValue());}
+ {
+ if(ii == 0)
+ System.out. println("filling "+k+" with int");
+ tuple.fill(k,ent.ValueAtIndex(j).intValue());
+ }
else
- {tuple.fill(k,ent.ValueAtIndex(j).doubleValue());}
+ {tuple.fill(k,ent.ValueAtIndex(j).doubleValue());}
+ k++;
}
tuple.addRow();
}
-
-
-
}
@@ -118,8 +116,6 @@
private int taille=0;
Entry(){
};
-
-
/**
* lengt of this entry of the Tuple, all entry withe the same size are store in the same
* ITuple when tuple.dump() is called
@@ -172,11 +168,5 @@
System .out.println(" array:"+list.size()+" element 0 : "+list.get(0));
}
}
-
-
-
-
-
-
-
+
}
\ No newline at end of file
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -u -r1.1 -r1.2
--- FindMuonsBoth.java 3 Jun 2009 00:20:09 -0000 1.1
+++ FindMuonsBoth.java 13 Aug 2009 20:40:12 -0000 1.2
@@ -15,12 +15,16 @@
public class FindMuonsBoth extends Driver {
public String outputFile="foobar.slcio";
public String plotsFile="myplots.aida";
+ public static boolean al= false;
private Muons muons;
private MuonsBoth muonsBoth;
public FindMuonsBoth() {
- add(new TrackReconstructionDriver());
- muonsBoth=new MuonsBoth();
- add(muonsBoth);
+ if(!al){
+ add(new TrackReconstructionDriver());
+ muonsBoth=new MuonsBoth();
+ add(muonsBoth);
+ al =true;
+ }
}
public void setOutputFile(String outputFile) {
System.out.println("Will output events to " + outputFile);
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn/HelixTest
diff -u -r1.4 -r1.5
--- HelixTest.java 9 Jul 2009 00:43:49 -0000 1.4
+++ HelixTest.java 13 Aug 2009 20:40:13 -0000 1.5
@@ -5,11 +5,7 @@
package org.lcsim.contrib.Mbussonn.HelixTest;
-import hep.graphics.heprep.HepRepFactory;
-import hep.graphics.heprep.HepRepInstanceTree;
-import hep.graphics.heprep.HepRepType;
import hep.physics.matrix.SymmetricMatrix;
-import org.lcsim.contrib.onoprien.util.swim.Helix;
import hep.physics.vec.Hep3Vector;
import java.util.ArrayList;
import java.util.HashSet;
@@ -20,8 +16,6 @@
import org.freehep.application.studio.Studio;
import org.freehep.util.FreeHEPLookup;
import org.lcsim.fit.helicaltrack.HelixUtils;
-import org.lcsim.contrib.onoprien.util.swim.ZDisk;
-import org.lcsim.contrib.onoprien.util.vector.ConstHep3Vector;
import org.lcsim.detector.*;
import org.lcsim.event.Track;
import org.lcsim.geometry.Detector;
@@ -71,6 +65,8 @@
@Override
protected void process(EventHeader event)
{
+ // just to fix some issues with my computer wen driver are launched twice
+ // some variable assignement to "null" ans clearing empty list are for the same purpose
if(donothing){
System.out.println("doing nothing in here");
return;}
@@ -131,16 +127,13 @@
System.out.println("phi0 = "+hf.phi0());
System.out.println("Z0 = "+hf.z0());
System.out.println("======================");
-// System.out.println("path at 20 = "+HelixUtils.PathToZPlane(hf, 20));
- //System.out.println("path at -20 = "+HelixUtils.PathToZPlane(hf, -20));
- //System.out.println("x at z:20 = "+HelixUtils.isInterceptingZDisk(hf,0.0, 100.0, 20.0));
- //System.out.println("slope = "+hf.slope());
- //System.out.println("===================");
+
+
ISolid trkvol = det.getTrackingVolume().getLogicalVolume().getSolid();
////////////////////////////////////////////////////////////////
//not in use yet... smin is the minimum path length to either://
- // -go outside the dŽtector in the radial direction //
+ // -go outside the detector in the radial direction //
// -go outside the detector in the z direction //
////////////////////////////////////////////////////////////////
//HelixUtils.setMaxPath(2512);
@@ -185,7 +178,7 @@
int i;
//////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
- /// Graph the track path again, to be sure it's the good one ///
+ /// Graph the track path again, to be sure it's the good one (visual debug) ///
// ///
List<Hep3Vector> thehelixpath = new LinkedList<Hep3Vector>(); ///
HelicalTrackFit trackFit = this.TrackToHTF(track); ///
@@ -226,11 +219,17 @@
////////////////////////////////////////////////////////////////////////////////
// create a list of traversed volumes //
////////////////////////////////////////////////////////////////////////////////
- private Set<UniquePV> listOfTraverserVolumes(HelicalTrackFit helix,Set<UniquePV> listeOfVolumes)
+ /**
+ * Find all the volumes traversed by an HelicalTrackFit
+ * @param helix
+ * @param listOfVolumes to be tested
+ * @return an HashSet of the traversed volumes
+ */
+ private Set<UniquePV> listOfTraverserVolumes(HelicalTrackFit helix,Set<UniquePV> listOfVolumes)
{
Set<UniquePV> traversedVolumes = new HashSet<UniquePV>();
- for(UniquePV pv: listeOfVolumes){
+ for(UniquePV pv: listOfVolumes){
try{
if(intersect(helix,pv))
{traversedVolumes.add(pv);}
@@ -462,7 +461,8 @@
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
-
+ ///////////////////// NOT from me, juste a little adapted /////////////////////
+////////////////////////////////////////////////////////////////////////////////
/**
* A "struct" holding geometry information about a single physical volume
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn/JetFinder
diff -u -r1.7 -r1.8
--- JetFinder.java 20 Jul 2009 21:59:08 -0000 1.7
+++ JetFinder.java 13 Aug 2009 20:40:13 -0000 1.8
@@ -17,10 +17,13 @@
import org.lcsim.event.MCParticle;
import org.lcsim.event.ReconstructedParticle;
import org.lcsim.event.RelationalTable;
+import org.lcsim.event.SimTrackerHit;
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;
@@ -43,7 +46,7 @@
private JetFinder()
{
System.out.println("=============================================");
- System.out.println("= creating jetfinderExtended =");
+ System.out.println("= creating jetfinder =");
System.out.println("=============================================");
if(itw == null){
itw = new ITupleWrapper();
@@ -66,16 +69,17 @@
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);
+ RelationalTable helicalTrackHitToMc = 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");
List<Track> tracklist = event.getTracks();
for (LCRelation relation : mcrelations) {
hittomc.add(relation.getFrom(), relation.getTo());
}
+ itw.tuple("hittomc size").fill(hittomc.size());
//******ajoute une ralation entre reconstructed particle and mcp ********//
RelationalTable rc2mc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
List<LCRelation> rcmcfe = event.get(LCRelation.class, "rc2mc");//collection of relation is onyread once here
-
for (LCRelation relation : rcmcfe) { // the collection of rcp is not load anywhere
if(relation.getFrom()==null)
{continue;}
@@ -83,7 +87,7 @@
{continue;}
rc2mc.add(relation.getFrom(), relation.getTo());
}
-
+ itw.tuple("tracklistSize").fill(tracklist.size());
for (Track track : tracklist) {//<editor-fold desc="//construc a relation between track and mcp">
TrackAnalysis tkanal = new TrackAnalysis(track, hittomc);
MCParticle mcp = tkanal.getMCParticle();
@@ -102,6 +106,22 @@
trktomc.add(track, mcp);
}//</editor-fold>
+ ////*********** map the sim tracker hits with the particles********////
+ RelationalTable simhittomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+
+ // Get the collections of SimTrackerHits
+ List<List<SimTrackerHit>> simcols = event.get(SimTrackerHit.class);
+
+ // Loop over the SimTrackerHits and fill in the relational table
+ for (List<SimTrackerHit> simlist : simcols) {
+ for (SimTrackerHit simhit : simlist) {
+ simhittomc.add(simhit, simhit.getMCParticle());
+ }
+ }
+ ////**************************************************************////
+
+
+
int njets=0;
List<MCParticle> mcplist = new Vector<MCParticle>();
FindableTrack findable = new FindableTrack(event);
@@ -135,7 +155,8 @@
int ntracks=0;
int nfindableMCP=0;
MCParticleExtended mcpx = new MCParticleExtended();
-
+ String sfile = StrategyXMLUtils.getDefaultStrategiesPrefix() + "sATLASBarrel-SM08.xml";
+ List<SeedStrategy> slist = StrategyXMLUtils.getStrategyListFromResource(sfile);
for(MCParticle mcp : mcplist){
@@ -145,9 +166,11 @@
double ptotal = mcpx.getPTotal();
double ptcut = 1.0;
if( ptotal < ptcut
- ||findable.LayersHit(mcp)<7
+ //||findable.LayersHit(mcp)<7
||mcp.getGeneratorStatus()!= MCParticleExtended.FINAL_STATE
- ||Math.abs(mcpx.getEta())>2.5)
+ //||Math.abs(mcpx.getEta())>2.5
+ ||!findable.isFindable(mcp, slist)
+ )
{
//System.out.println("skip mcp "+ptotal+" nhits"+findable.LayersHit(mcp));
continue;}
@@ -166,6 +189,7 @@
double pmcp = -10 ;
double trackp = -10;
double mcpp =-10;
+ double trackChiSq = -10000;
@@ -179,15 +203,53 @@
}
else if(setOfTrack.size()==0) {
wgt = 0.0;
+ Set<TrackerHit> MCPHits = hittomc.allTo(mcp);//helicaltrack hits
+ Set<SimTrackerHit> SimHits = simhittomc.allTo(mcp);
+ //loop through mcparticle helical hits,
+ // find the closest simtrackerhit,
+ // find th edifference
+ double chisq=0;
+ for(TrackerHit helicaltrackhit : MCPHits){
+ double chiSq=0;
+ double helicalHitPosition[] = helicaltrackhit.getPosition();
+ SimTrackerHit closerSimHit = SimHits.iterator().next();
+ double firstSimHitPosition[] = closerSimHit.getPoint();
+ double dx = helicalHitPosition[0]-firstSimHitPosition[0];
+ double dy = helicalHitPosition[1]-firstSimHitPosition[1];
+ double dz = helicalHitPosition[2]-firstSimHitPosition[2];
+ double currentDistanceSquared = dx*dx+dy*dy+dz*dz;
+ double[] sigmaMat = helicaltrackhit.getCovMatrix();
+ double sigma = sigmaMat[0]*sigmaMat[0]+sigmaMat[2]*sigmaMat[2]+sigmaMat[5]*sigmaMat[5];
+ for(SimTrackerHit nextSimhit : SimHits){
+ double simHitPosition[] = nextSimhit.getPoint();
+ dx = helicalHitPosition[0]-simHitPosition[0];
+ dy = helicalHitPosition[1]-simHitPosition[1];
+ dz = helicalHitPosition[2]-simHitPosition[2];
+ double nextDistanceSquared = dx*dx+dy*dy+dz*dz;
+ if(nextDistanceSquared<currentDistanceSquared){
+ currentDistanceSquared = nextDistanceSquared;
+ closerSimHit = nextSimhit;
+
+ }
+
+ }
+ chisq += currentDistanceSquared/sigma;
+
+ }
+ itw.Cloud1D("selfChisq").fill(chisq);
+ ////*******************************//
+ itw.tuple("notfoundHelicalHits").fill(MCPHits.size());
itw.tuple("notFoundsNumberOfHits").fillWithInt(findable.LayersHit(mcp));
itw.tuple("notFoundsEta").fillWithDouble(mcpx.getEta());
itw.tuple("notFounds-Angle").fillWithDouble(angle);
}else if (setOfTrack.size()==1){
Set<TrackerHit> MCPHits = new HashSet(hittomc.allTo(mcp));
+
ntracks++;
wgt = 1.0;
Track track = setOfTrack.iterator().next();
+ trackChiSq = track.getChi2();
d0track = track.getTrackParameter(HelicalTrackFit.dcaIndex);
d0mcp = mcpx.getDCA();
z0track = track.getTrackParameter(HelicalTrackFit.z0Index);
@@ -223,7 +285,8 @@
itw.tuple("delta_Eta").fill(deltaEta);
itw.tuple("delta_R").fill(Math.sqrt(deltaPhi*deltaPhi+deltaEta*deltaEta));
- itw.tuple("track_d0").fillWithDouble(d0track);
+ itw.tuple("track_d0").fill(d0track);
+ itw.tuple("track_chiSq").fill(trackChiSq);
itw.tuple("track_z0").fillWithDouble(z0track);
itw.tuple("mcp_d0").fillWithDouble(d0mcp);
itw.tuple("mcp_z0").fillWithDouble(z0mcp);
CVSspam 0.2.8