9 added files
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -N TextFileToStdhep.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TextFileToStdhep.java 3 Jun 2009 00:20:09 -0000 1.1
@@ -0,0 +1,149 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.contrib.Mbussonn;
+
+/**
+ *
+ * @author richp
+ * modified by mbussonn to take arguments from command line
+ */
+import java.io.FileReader;
+import java.io.IOException;
+import java.io.StreamTokenizer;
+import java.util.ArrayList;
+import java.util.List;
+
+import hep.io.stdhep.StdhepBeginRun;
+import hep.io.stdhep.StdhepEndRun;
+import hep.io.stdhep.StdhepEvent;
+import hep.io.stdhep.StdhepWriter;
+
+public class TextFileToStdhep {
+
+ /**
+ * @param args the command line arguments
+ * @throws IOException
+ */
+ public static void main(String[] args) throws IOException {
+ String file = args[0];
+ String infile = file;
+ String outfile = file.substring(0, file.length()-3)+"stdhep";//replace .txt by .stdhep
+ System.out.println("==== processing "+infile+" into "+outfile+" ====");
+ process(infile, outfile);
+ }
+
+ private static void process(String infile, String outfile) throws IOException {
+// try {
+ FileReader fr = new FileReader(infile);
+ StreamTokenizer tok = new StreamTokenizer(fr);
+ tok.resetSyntax();
+ tok.wordChars(33, 255);
+ tok.parseNumbers();
+ tok.whitespaceChars(0, ' ');
+ tok.eolIsSignificant(true);
+ List<Double> beg = getNumbersInLine(tok);
+ if (beg.size() != 7) {
+ throw new RuntimeException("Unexpected entry in begin run record");
+ }
+ int nreq = beg.get(0).intValue();
+ int ngen = beg.get(1).intValue();
+ int nwrit = beg.get(2).intValue();
+ float ecm = beg.get(3).floatValue();
+ float xsec = beg.get(4).floatValue();
+ double rn1 = beg.get(5);
+ double rn2 = beg.get(6);
+ StdhepWriter sw = new StdhepWriter(outfile, "Imported Stdhep Events",
+ "From file " + infile, nwrit);
+ sw.setCompatibilityMode(false);
+ StdhepBeginRun sb = new StdhepBeginRun(nreq, ngen, nwrit, ecm, xsec, rn1, rn2);
+ sw.writeRecord(sb);
+ for (int icross = 0; icross < nwrit; icross++) {
+ List<Double> nums = getNumbersInLine(tok);
+ if (nums.size() != 2) {
+ throw new RuntimeException("Unexpected entry for number of particles");
+ }
+ int nevhep = nums.get(0).intValue();
+ int nhep = nums.get(1).intValue();
+ System.out.println("Number of particles for event " + nevhep + ": " + nhep);
+ int isthep[] = new int[nhep];
+ int idhep[] = new int[nhep];
+ int jmohep[] = new int[2 * nhep];
+ int jdahep[] = new int[2 * nhep];
+ double phep[] = new double[5 * nhep];
+ double vhep[] = new double[4 * nhep];
+ for (int i = 0; i < nhep; i++) {
+ List<Double> vals = getNumbersInLine(tok);
+ if (vals.size() != 16) {
+ throw new RuntimeException("Unexpected entry for a particle");
+ }
+ int ip = vals.get(0).intValue();
+ if (ip != i + 1) {
+ throw new RuntimeException("Particle numbering mismatch");
+ }
+ isthep[i] = vals.get(1).intValue();
+ idhep[i] = vals.get(2).intValue();
+ jmohep[2 * i] = vals.get(3).intValue();
+ jmohep[2 * i + 1] = vals.get(4).intValue();
+ jdahep[2 * i] = vals.get(5).intValue();
+ jdahep[2 * i + 1] = vals.get(6).intValue();
+ for (int j = 0; j < 5; j++) {
+ phep[5 * i + j] = vals.get(j + 7);
+ }
+ for (int j = 0; j < 4; j++) {
+ vhep[4 * i + j] = vals.get(j + 12);
+ }
+
+// if (i == 0 || i == nhep - 1) {
+// System.out.println(ip + " st: " + isthep[i] + " id: " + idhep[i] +
+// " jmo: " + jmohep[2 * i] + " " + jmohep[2 * i + 1] +
+// " jda: " + jdahep[2 * i] + " " + jdahep[2 * i + 1]);
+// System.out.println("p: " + phep[5 * i] + " " + phep[5 * i + 1] + " " +
+// phep[5 * i + 2] + " " + phep[5 * i + 3] + " " + phep[5 * i + 4]);
+// System.out.println("v: " + vhep[4 * i] + " " + vhep[4 * i + 1] + " " +
+// vhep[4 * i + 2] + " " + vhep[4 * i + 3]);
+// }
+ }
+ StdhepEvent ev = new StdhepEvent(nevhep, nhep, isthep, idhep, jmohep, jdahep, phep, vhep);
+ sw.writeRecord(ev);
+ }
+ List<Double> end = getNumbersInLine(tok);
+ if (end.size() != 7) {
+ throw new RuntimeException("Unexpected entry in begin run record");
+ }
+ int nreqe = end.get(0).intValue();
+ int ngene = end.get(1).intValue();
+ int nwrite = end.get(2).intValue();
+ float ecme = end.get(3).floatValue();
+ float xsece = end.get(4).floatValue();
+ double rn1e = end.get(5);
+ double rn2e = end.get(6);
+ StdhepEndRun se = new StdhepEndRun(nreqe, ngene, nwrite, ecme, xsece, rn1e, rn2e);
+ sw.writeRecord(se);
+
+ sw.close();
+ fr.close();
+// } catch (IOException e) {
+// System.out.println("Error opening " + infile);
+// }
+ }
+
+ @SuppressWarnings("static-access")
+ private static List<Double> getNumbersInLine(StreamTokenizer tok) throws IOException {
+ List<Double> nums = new ArrayList<Double>();
+ while (tok.nextToken() != tok.TT_EOF) {
+ if (tok.ttype == tok.TT_EOL) {
+ break;
+ }
+
+ if (tok.ttype != tok.TT_NUMBER) {
+ throw new RuntimeException("Non numeric data encountered");
+ }
+
+ nums.add(tok.nval);
+ }
+
+ return nums;
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -N MuonsAnalysis.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MuonsAnalysis.java 3 Jun 2009 00:20:09 -0000 1.1
@@ -0,0 +1,231 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.Mbussonn;
+
+import hep.aida.IHistogramFactory;
+import hep.aida.IProfile1D;
+import java.io.IOException;
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Set;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+//import org.lcsim.contrib.KFFitter.KFTrack;
+//import org.lcsim.contrib.KFFitter.KFTrackParameters;
+//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.event.TrackerHit;
+import org.lcsim.event.base.BaseRelationalTable;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.RelationalTable;
+//import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.Track;
+import org.lcsim.fit.helicaltrack.*;
+import org.lcsim.recon.tracking.seedtracker.*;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.Driver;
+//*/
+/**
+ *
+ * @author matthiasbussonnier
+ */
+public class MuonsAnalysis extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+ private IProfile1D efficiencyVSEta,avgEfficiencyVSEta;
+ private IProfile1D effMuons;
+ private java.util.Random gtor;
+ public String outputPlots="myplots.aida";
+ public MuonsAnalysis(){
+ 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);
+ //gtor = new java.util.Random(10101);
+ }
+
+ @Override
+ @SuppressWarnings({"unchecked", "unchecked","unchecked"})
+ public void process (EventHeader event){
+ super.processChildren(event);
+ int numberOfMuons =0;
+ int numberOfMuonsWithTrack=0;
+ //int nftrack=0;
+
+ // 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");
+ //trktomc = (RelationalTable) event.get("trackToMCRelationTable");"TrackToMCRelation"
+ List<LCRelation> evtTracktoMc = event.get(LCRelation.class, "TrackToMCRelation");
+// List<Track> tracklist = event.getTracks();
+
+//*/
+ //for (LCRelation relation : mcrelations) {
+ // hittomc.add(relation.getFrom(), relation.getTo());
+ // }
+ for (LCRelation relation : evtTracktoMc) {
+ trktomc.add(relation.getFrom(), relation.getTo());
+ }
+
+/*/
+ for (Track track : tracklist) {
+ TrackAnalysis tkanal = new TrackAnalysis(track, hittomc);
+ MCParticle mcp = tkanal.getMCParticle();
+ tkanal=null;
+ double px = track.getPX();
+ double py = track.getPY();
+ double pz = track.getPZ();
+ double ptotal = Math.sqrt(px*px+py*py+pz*pz);
+ double ptcut = 50.;
+ double dca = track.getTrackParameter(HelicalTrackFit.dcaIndex);
+ double z0 = track.getTrackParameter(HelicalTrackFit.z0Index);
+
+ if(track.getCharge()==0)
+ continue;
+ if( ptotal < ptcut)
+ continue;
+
+ aida.cloud2D("other tracks/orrigin ").fill(dca, z0);
+ aida.cloud1D("other tracks/scalar momentum").fill(ptotal);
+ }
+//*/
+
+ //aida.cloud1D("number of track finded per event").fill(nftrack);
+ String sfile = StrategyXMLUtils.getDefaultStrategiesPrefix() + "sATLASFull-matt.xml";
+ List<SeedStrategy> slist = StrategyXMLUtils.getStrategyListFromResource(sfile);
+ FindableTrack findable = new FindableTrack(event);
+
+ List<MCParticle> mcplist = event.getMCParticles();
+ for(MCParticle mcpx : mcplist){
+ MCParticleExtended mcp = new MCParticleExtended(mcpx);
+ Set<Track> setOfTrack= trktomc.allTo(mcpx);
+ double ptotal = mcp.getPTotal();
+ double ptcut=50.;
+ if( ptotal < ptcut
+ ||(mcp.getPDGID()!=13 && mcp.getPDGID()!=-13)
+ ||mcp.getGeneratorStatus()!= MCParticleExtended.FINAL_STATE
+ ||Math.abs(mcp.getEta())>2.5)
+ {
+ System.err.println("pass through processs");
+ continue;}
+ else{
+ if(!findable.isFindable(mcpx, slist)){
+ 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.;
+ numberOfMuons++;
+ 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 muons").fill(mcp.getEta());
+ aida.cloud1D("mcp/theta for non-finded muons").fill(mcp.getTheta()); aida.cloud1D("nhits for non-finded muons").fill(nhits);
+ }
+ else if (setOfTrack.size()==1){
+ weight=1;
+ numberOfMuonsWithTrack++;
+ Track track = setOfTrack.iterator().next();
+ double d0 = track.getTrackParameter(HelicalTrackFit.dcaIndex);
+ double z0 = track.getTrackParameter(HelicalTrackFit.z0Index);
+ double px = track.getPX();
+ 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);
+ aida.cloud1D("track/d0").fill(d0);
+ 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);
+
+ //************************************************//
+ //Lets try the kalmann filter and see what's wrong//
+ //************************************************//
+// KFTrack kft = new KFTrack();
+// kft.defineSeedTrack(track);
+// List<TrackerHit> stHits= (List<TrackerHit>) track.getTrackerHits();
+// for(TrackerHit hit:stHits){
+// kft.addHit(hit,gtor.nextGaussian(),gtor.nextGaussian());
+// }
+// kft.fit();
+// KFTrackParameters kalmanTrack = kft.getKFTrackParameters();
+// aida.cloud2D("kf/d0 kf VS do Htf").fill(kalmanTrack.getD0(),d0);
+// aida.cloud2D("kf/Z0 kf VS Z0 Htf").fill(kalmanTrack.getZ(),z0);
+// aida.cloud2D("kf/Z0 kf VS Z0 Htf").fill(kalmanTrack.getZ(),z0);
+// aida.cloud2D("kf/blobHF").fill(d0,z0);
+// aida.cloud2D("kf/blobKF").fill(kalmanTrack.getD0(),kalmanTrack.getZ());
+
+ }
+ int imax=50;
+ double etamax=2.5;
+ for(int i=0; i <= imax; i++)
+ {
+ double etac = ((double) i / (double) imax)*etamax+0.001;//+0.001 to avoit integer value wich give holes in graph
+
+ double etamcp = Math.abs(mcp.getEta());
+ if(etamcp < etac){
+ avgEfficiencyVSEta.fill(etac,weight);
+ }
+
+ }
+ efficiencyVSEta.fill(mcp.getEta(),weight);
+ aida.cloud1D("mcp/eta for detectable muons").fill(mcp.getEta());
+ aida.cloud1D("mcp/theta for detectable muons").fill(mcp.getTheta());
+ }
+ }//ENDFOR
+
+ aida.cloud1D("mcp/number of detectable muons").fill(numberOfMuons);
+ aida.cloud1D("mcp/number of detected muons").fill(numberOfMuonsWithTrack);
+ }
+ @Override
+ public void endOfData() {
+ try {
+ aida.saveAs(outputPlots);
+ }
+ catch (IOException ex) {
+ Logger.getLogger(MuonsAnalysis.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ //System.out.println("# of reco tracks = " + _nrecTrk + "; # of MC tracks = " + _nmcTrk + "; Efficiency = " + _nrecTrk / _nmcTrk);
+ }
+ public void setOutputPlots(String output){
+ this.outputPlots=output;
+
+ }
+}//end class
+
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -N FindMuonsAnalysis.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ FindMuonsAnalysis.java 3 Jun 2009 00:20:09 -0000 1.1
@@ -0,0 +1,36 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.contrib.Mbussonn;
+
+import org.lcsim.util.Driver;
+import org.lcsim.util.loop.LCIODriver;
+//import org.lcsim.util.loop.LCIODriver;
+/**
+ *
+ * @author matthiasbussonnier
+ */
+public class FindMuonsAnalysis extends Driver {
+ public String outputFile="foobar.slcio";
+ public String plotsFile="myplots.aida";
+ Muons muonsAnalysis;
+ private MuonsAnalysis muons;
+ public FindMuonsAnalysis() {
+
+ muons=new MuonsAnalysis();
+ add(muons);
+ }
+ public void setOutputFile(String outputFile) {
+ System.out.println("Will output events to " + outputFile);
+ add(new LCIODriver(outputFile));
+ }
+ public void setPlotsFile(String plotsFile) {
+ System.out.println("Will output plots to " + plotsFile);
+ muons.setOutputPlots(plotsFile);
+ }
+
+
+
+
+}
\ No newline at end of file
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -N Muons.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ Muons.java 3 Jun 2009 00:20:09 -0000 1.1
@@ -0,0 +1,120 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.Mbussonn;
+
+import java.io.IOException;
+import java.util.ArrayList;
+import java.util.List;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.contrib.Partridge.TrackingTest.TrackAnalysis;
+import org.lcsim.digisim.MyLCRelation;
+import org.lcsim.event.base.BaseRelationalTable;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.RelationalTable;
+import org.lcsim.event.Track;
+import org.lcsim.recon.tracking.seedtracker.*;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.Driver;
+//*/
+/**
+ *
+ * @author matthiasbussonnier
+ */
+public class Muons extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+ //private IProfile1D efficiencyVSEta,avgEfficiencyVSEta;
+ //private IProfile1D effMuons;
+ private java.util.Random gtor;
+ public String outputPlots="myplots.aida";
+ public Muons(){
+ // 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);
+ gtor = new java.util.Random(10101);
+ }
+
+ @Override
+ @SuppressWarnings({"unchecked", "unchecked"})
+ public void process (EventHeader event){
+ int nftrack=0;
+
+ 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");
+ List<LCRelation> trktomc2 = new ArrayList<LCRelation>();
+ /*event.put("trackToMCRelationTable", trktomc);*/
+
+ List<Track> tracklist = event.getTracks();
+
+ for (LCRelation relation : mcrelations) {
+ hittomc.add(relation.getFrom(), relation.getTo());
+ }
+ for (Track track : tracklist) {
+ TrackAnalysis tkanal = new TrackAnalysis(track, hittomc);
+ MCParticle mcp = tkanal.getMCParticle();
+ double px = track.getPX();
+ double py = track.getPY();
+ double pz = track.getPZ();
+ double ptotal = Math.sqrt(px*px+py*py+pz*pz);
+ double ptcut = 50.;
+ if(track.getCharge()==0)
+ continue;
+ if( ptotal < ptcut)
+ continue;
+ nftrack++;
+ if (mcp != null){
+ trktomc2.add(new MyLCRelation(track,mcp));
+ trktomc.add(track, mcp);
+ }
+ }
+
+
+
+ String str="sATLASFull-JeffMarch26.xml";
+ String sfile = StrategyXMLUtils.getDefaultStrategiesPrefix() + str;
+ System.out.println("using strategy in "+str);
+
+ List<MCParticle> mcplist = event.getMCParticles();
+ for(MCParticle mcpx : mcplist){
+ MCParticleExtended mcp = new MCParticleExtended(mcpx);
+ double ptotal = mcp.getPTotal();
+ double ptcut=50.;
+ if( ptotal < ptcut
+ ||(mcp.getPDGID()!=13 && mcp.getPDGID()!=-13)
+ ||mcp.getGeneratorStatus()!= MCParticleExtended.FINAL_STATE
+ ||Math.abs(mcp.getEta())>2.5)
+ {
+ event.getMCParticles().remove(mcp);
+ continue;}
+ }//ENDFOR
+ event.put("TrackToMCRelation", trktomc2, LCRelation.class, 0);
+ return ;
+ }
+
+ /**********************************************//**********************************************/
+ /**********************************************//**********************************************/
+ /**********************************************//**********************************************/
+ /**********************************************//**********************************************/
+ @Override
+ public void endOfData() {
+ try {
+ aida.saveAs(outputPlots);
+ }
+ catch (IOException ex) {
+ Logger.getLogger(Muons.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ //System.out.println("# of reco tracks = " + _nrecTrk + "; # of MC tracks = " + _nmcTrk + "; Efficiency = " + _nrecTrk / _nmcTrk);
+ }
+ public void setOutputPlots(String output){
+ this.outputPlots=output;
+
+ }
+}//end class
+
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -N MuonsBoth.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MuonsBoth.java 3 Jun 2009 00:20:09 -0000 1.1
@@ -0,0 +1,223 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.Mbussonn;
+
+import hep.aida.IHistogramFactory;
+import hep.aida.IProfile1D;
+import java.io.IOException;
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Set;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.contrib.mgraham.sATLASDigi.FindableTrack;
+import org.lcsim.contrib.mgraham.sATLASDigi.FindableTrack.Ignore;
+//import org.lcsim.contrib.Partridge.TrackingTest.FindableTrack.Ignore;
+import org.lcsim.contrib.Partridge.TrackingTest.TrackAnalysis;
+import org.lcsim.contrib.mgraham.sATLASDigi.FindableTrack;
+import org.lcsim.digisim.MyLCRelation;
+import org.lcsim.event.base.BaseRelationalTable;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.RelationalTable;
+import org.lcsim.event.Track;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.recon.tracking.seedtracker.*;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.Driver;
+//*/
+/**
+ *
+ * @author matthiasbussonnier
+ */
+public class MuonsBoth extends Driver {
+private AIDA aida = AIDA.defaultInstance();
+ private IProfile1D efficiencyVSEta,avgEfficiencyVSEta;
+ private IProfile1D effMuons;
+ 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);
+ gtor = new java.util.Random(10101);
+ }
+
+ @Override
+ @SuppressWarnings({"unchecked", "unchecked"})
+ public void process (EventHeader event){
+ int numberOfMuons =0;
+ int numberOfMuonsWithTrack=0;
+ int nftrack=0;
+
+ 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");
+ List<Track> tracklist = event.getTracks();
+
+ for (LCRelation relation : mcrelations) {
+ hittomc.add(relation.getFrom(), relation.getTo());
+ }
+ for (Track track : tracklist) {
+ TrackAnalysis tkanal = new TrackAnalysis(track, hittomc);
+ MCParticle mcp = tkanal.getMCParticle();
+ tkanal=null;
+ double px = track.getPX();
+ double py = track.getPY();
+ double pz = track.getPZ();
+ double ptotal = Math.sqrt(px*px+py*py+pz*pz);
+ double ptcut = 50.;
+ double dca = track.getTrackParameter(HelicalTrackFit.dcaIndex);
+ double z0 = track.getTrackParameter(HelicalTrackFit.z0Index);
+
+ if(track.getCharge()==0)
+ continue;
+ if( ptotal < ptcut)
+ continue;
+
+ aida.cloud2D("other tracks/orrigin ").fill(dca, z0);
+ aida.cloud1D("other tracks/scalar momentum").fill(ptotal);
+ nftrack++;
+ if (mcp != null)
+ trktomc.add(track, mcp);
+ }
+
+
+ aida.cloud1D("number of track finded per event").fill(nftrack);
+ String sfile = StrategyXMLUtils.getDefaultStrategiesPrefix() + "sATLASFull-JeffMarch26.xml";
+ List<SeedStrategy> slist = StrategyXMLUtils.getStrategyListFromResource(sfile);
+ FindableTrack findable = new FindableTrack(event);
+
+ List<MCParticle> mcplist = event.getMCParticles();
+ for(MCParticle mcpx : mcplist){
+ MCParticleExtended mcp = new MCParticleExtended(mcpx);
+ Set<Track> setOfTrack= trktomc.allTo(mcpx);
+ double ptotal = mcp.getPTotal();
+ double ptcut=50.;
+ if( ptotal < ptcut
+ ||(mcp.getPDGID()!=13 && mcp.getPDGID()!=-13)
+ ||mcp.getGeneratorStatus()!= MCParticleExtended.FINAL_STATE
+ ||Math.abs(mcp.getEta())>2.5)
+ {continue;}
+ else{
+ if(!findable.isFindable(mcpx, slist)){
+ 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.;
+ numberOfMuons++;
+ 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 muons").fill(mcp.getEta());
+ aida.cloud1D("mcp/theta for non-finded muons").fill(mcp.getTheta()); aida.cloud1D("nhits for non-finded muons").fill(nhits);
+ }
+ else if (setOfTrack.size()==1){
+ weight=1;
+ numberOfMuonsWithTrack++;
+ Track track = setOfTrack.iterator().next();
+ double d0 = track.getTrackParameter(HelicalTrackFit.dcaIndex);
+ double z0 = track.getTrackParameter(HelicalTrackFit.z0Index);
+ double px = track.getPX();
+ 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);
+ aida.cloud1D("track/d0").fill(d0);
+ 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);
+
+ //************************************************//
+ //Lets try the kalmann filter and see what's wrong//
+ //************************************************//
+// KFTrack kft = new KFTrack();
+// kft.defineSeedTrack(track);
+// List<TrackerHit> stHits= (List<TrackerHit>) track.getTrackerHits();
+// for(TrackerHit hit:stHits){
+// kft.addHit(hit,gtor.nextGaussian(),gtor.nextGaussian());
+// }
+// kft.fit();
+// KFTrackParameters kalmanTrack = kft.getKFTrackParameters();
+// aida.cloud2D("kf/d0 kf VS do Htf").fill(kalmanTrack.getD0(),d0);
+// aida.cloud2D("kf/Z0 kf VS Z0 Htf").fill(kalmanTrack.getZ(),z0);
+// aida.cloud2D("kf/Z0 kf VS Z0 Htf").fill(kalmanTrack.getZ(),z0);
+// aida.cloud2D("kf/blobHF").fill(d0,z0);
+// aida.cloud2D("kf/blobKF").fill(kalmanTrack.getD0(),kalmanTrack.getZ());
+
+ }
+ int imax=50;
+ double etamax=2.5;
+ for(int i=0; i <= imax; i++)
+ {
+ double etac = ((double) i / (double) imax)*etamax+0.001;//+0.001 to avoit integer value wich give holes in graph
+
+ double etamcp = Math.abs(mcp.getEta());
+ if(etamcp < etac){
+ avgEfficiencyVSEta.fill(etac,weight);
+ }
+
+ }
+ efficiencyVSEta.fill(mcp.getEta(),weight);
+ aida.cloud1D("mcp/eta for detectable muons").fill(mcp.getEta());
+ aida.cloud1D("mcp/theta for detectable muons").fill(mcp.getTheta());
+
+
+ }
+ }//ENDFOR
+
+ aida.cloud1D("mcp/number of detectable muons").fill(numberOfMuons);
+ aida.cloud1D("mcp/number of detected muons").fill(numberOfMuonsWithTrack);
+ }
+ @Override
+ public void endOfData() {
+ try {
+ aida.saveAs(outputPlots);
+ }
+ catch (IOException ex) {
+ Logger.getLogger(Muons.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ //System.out.println("# of reco tracks = " + _nrecTrk + "; # of MC tracks = " + _nmcTrk + "; Efficiency = " + _nrecTrk / _nmcTrk);
+ }
+ public void setOutputPlots(String output){
+ this.outputPlots=output;
+
+ }
+}//end class
+
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -N FindMuonsBoth.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ FindMuonsBoth.java 3 Jun 2009 00:20:09 -0000 1.1
@@ -0,0 +1,37 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.contrib.Mbussonn;
+
+import org.lcsim.contrib.sATLAS.TrackReconstructionDriver;
+import org.lcsim.util.Driver;
+import org.lcsim.util.loop.LCIODriver;
+//import org.lcsim.util.loop.LCIODriver;
+/**
+ *
+ * @author matthiasbussonnier
+ */
+public class FindMuonsBoth extends Driver {
+ public String outputFile="foobar.slcio";
+ public String plotsFile="myplots.aida";
+ private Muons muons;
+ private MuonsBoth muonsBoth;
+ public FindMuonsBoth() {
+ add(new TrackReconstructionDriver());
+ muonsBoth=new MuonsBoth();
+ add(muonsBoth);
+ }
+ public void setOutputFile(String outputFile) {
+ System.out.println("Will output events to " + outputFile);
+ add(new LCIODriver(outputFile));
+ }
+ public void setPlotsFile(String plotsFile) {
+ System.out.println("Will output plots to " + plotsFile);
+ muons.setOutputPlots(plotsFile);
+ }
+
+
+
+
+}
\ No newline at end of file
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -N FindMuons.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ FindMuons.java 3 Jun 2009 00:20:10 -0000 1.1
@@ -0,0 +1,44 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.contrib.Mbussonn;
+
+import org.lcsim.contrib.sATLAS.*;
+import org.lcsim.util.Driver;
+import org.lcsim.util.loop.LCIODriver;
+//import org.lcsim.util.loop.LCIODriver;
+/**
+ *
+ * @author matthiasbussonnier
+ */
+public class FindMuons extends Driver {
+ public String outputFile="foobar.slcio";
+ public String plotsFile="myplots.aida";
+ protected String stgfile="chose";
+ //TrackAnalysisDriver tad;
+
+ Muons muons;
+ public FindMuons() {
+ add(new TrackReconstructionDriver());
+ muons=new Muons();
+ add(muons);
+ }
+
+ public void setOutputFile(String outputFile) {
+ System.out.println("Will output events to " + outputFile);
+ LCIODriver lcio = new LCIODriver(outputFile);
+ //lcio.setWriteOnlyCollection("trackToMCRelationTable");
+ String[] ignoreCollections = {"HelicalTrackHitRelation","HelicalTrackHitRelations","VtxBarrHits","SCTShortEndcapHits","SCTShortBarrHits","MCParticleEndPointEnergy","MCParticleTree","VtxEndcapHits"};
+ lcio.setIgnoreCollections(ignoreCollections);
+ add(lcio);
+ }
+ public void setPlotsFile(String plotsFile) {
+ System.out.println("Will output plots to " + plotsFile);
+ muons.setOutputPlots(plotsFile);
+ }
+
+
+
+
+}
\ No newline at end of file
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -N MCParticleExtended.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MCParticleExtended.java 3 Jun 2009 00:20:10 -0000 1.1
@@ -0,0 +1,130 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.Mbussonn;
+
+import hep.physics.particle.properties.ParticleType;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.HepLorentzVector;
+import java.util.List;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.base.BaseMCParticle;
+
+/**
+ *
+ * @author matthiasbussonnier
+ */
+public class MCParticleExtended implements MCParticle {
+ public MCParticle mcp;
+
+
+ public MCParticleExtended(MCParticle newMCP){
+ mcp= newMCP;
+ }
+
+ /**
+ * give the absolute value of the radial momentum
+ * @return
+ */
+ public double getPR()
+ {
+ double px = mcp.getPX();
+ double py = mcp.getPY();
+ return Math.sqrt(px*px+py*py);
+ }
+ public double getPTotal()
+ {
+ double px = mcp.getPX();
+ double py = mcp.getPY();
+ double pz = mcp.getPZ();
+ return Math.sqrt(px*px+py*py+pz*pz);
+ }
+ public double getEta() {return -Math.log(Math.tan(this.getTheta()/2));}
+ public double getTheta(){return Math.acos(this.getCosTheta());}
+ public double getCosTheta(){
+ return mcp.getPZ()/this.getPTotal();
+ }
+
+ public List<MCParticle> getParents() {
+ return mcp.getParents();
+ }
+
+ public List<MCParticle> getDaughters() {
+ return mcp.getDaughters();
+ }
+
+ public SimulatorStatus getSimulatorStatus() {
+ return mcp.getSimulatorStatus();
+ }
+
+ public Hep3Vector getEndPoint() {
+ return mcp.getEndPoint();
+ }
+
+ public double getEnergy() {
+ return mcp.getEnergy();
+ }
+
+ public double getMass() {
+ return mcp.getMass();
+ }
+
+ public ParticleType getType() {
+ return mcp.getType();
+ }
+
+ public int getPDGID() {
+ return mcp.getPDGID();
+ }
+
+ public Hep3Vector getMomentum() {
+ return mcp.getMomentum();
+ }
+
+ public Hep3Vector getOrigin() {
+ return mcp.getOrigin();
+ }
+
+ public double getOriginX() {
+ return mcp.getOriginX();
+ }
+
+ public double getOriginY() {
+ return mcp.getEnergy();
+ }
+
+ public double getOriginZ() {
+ return mcp.getOriginZ();
+ }
+
+ public double getPX() {
+ return mcp.getPX();
+ }
+
+ public double getPY() {
+ return mcp.getPY();
+ }
+
+ public double getPZ() {
+ return mcp.getPZ();
+ }
+
+ public double getProductionTime() {
+ return mcp.getProductionTime();
+ }
+
+ public int getGeneratorStatus() {
+ return mcp.getGeneratorStatus();
+ }
+
+ public double getCharge() {
+ return mcp.getCharge();
+ }
+
+ public HepLorentzVector asFourVector() {
+ return mcp.asFourVector();
+ }
+
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/Mbussonn
diff -N JetFinder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ JetFinder.java 3 Jun 2009 00:20:10 -0000 1.1
@@ -0,0 +1,44 @@
+package org.lcsim.contrib.Mbussonn;
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+import java.util.List;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.util.JetDriver;
+import org.lcsim.mc.fast.MCFast;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+
+
+/**
+ *
+ * @author matthiasbussonnier
+ */
+public class JetFinder extends Driver{
+ private AIDA aida = AIDA.defaultInstance();
+ public JetFinder()
+ {
+ //super.add(new MCFast());
+ super.add(new JetDriver());
+ }
+
+ @Override
+ protected void process(EventHeader event)
+ {
+ super.process(event);
+
+ List<ReconstructedParticle> jets = event.get(ReconstructedParticle.class,"Jets");
+ aida.cloud1D("nJets").fill(jets.size());
+ for (ReconstructedParticle jet : jets)
+ {
+ List<ReconstructedParticle> particlesInJet = jet.getParticles();
+ aida.cloud1D("nParticles").fill(particlesInJet.size());
+ }
+ }
+
+
+
+}
CVSspam 0.2.8
