hps-java/src/main/java/org/lcsim/hps/users/phansson
diff -N TruthMomentumResolutionDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TruthMomentumResolutionDriver.java 21 Sep 2012 18:36:34 -0000 1.1
@@ -0,0 +1,329 @@
+package org.lcsim.hps.users.phansson;
+
+import hep.aida.*;
+import java.io.IOException;
+import java.util.*;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.Track;
+import org.lcsim.event.util.ParticleTypeClassifier;
+import org.lcsim.geometry.Detector;
+import org.lcsim.hps.analysis.ecal.HPSMCParticlePlotsDriver;
+import org.lcsim.hps.monitoring.AIDAFrame;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author mgraham
+ */
+public class TruthMomentumResolutionDriver extends Driver {
+
+
+ private String outputPlotFileName="";
+ private boolean hideFrame = false;
+ int totalTracks = 0;
+ private boolean _debug = false;
+ private AIDA aida = AIDA.defaultInstance();
+ private AIDAFrame pFrame;
+ IAnalysisFactory af = aida.analysisFactory();
+ IPlotter pPlotter;
+ IPlotter pPlotter2;
+ IPlotter pPlotter3;
+ ICloud1D hElectronP;
+ ICloud1D hPositronP;
+ ICloud1D hTrackP;
+ ICloud1D hPosTrackP;
+ ICloud1D hTruthMatchedPosTrackP;
+ ICloud1D hTruthMatchedPosTrackPdiff;
+ ICloud1D hNegTrackP;
+ ICloud1D hTruthMatchedNegTrackP;
+ ICloud1D hTruthMatchedNegTrackPdiff;
+ ICloud1D hTruthMatchedPosTrackPdiffPrev[] = new ICloud1D[5];
+ ICloud1D hTruthMatchedNegTrackPdiffPrev[] = new ICloud1D[5];
+ ICloud1D hNTracks;
+ IHistogram1D trkCountVsEventPlot;
+
+ ICloud1D hNPosTracks;
+ ICloud1D hNNegTracks;
+ ICloud1D hNPositronsForTrack;
+ ICloud1D hNElectronsForTrack;
+ ICloud1D hNPositronsForTrackInv;
+ ICloud1D hNElectronsForTrackInv;
+
+ HashMap<Integer,MCParticle> mc_ele_prev = new HashMap<Integer, MCParticle>();
+ HashMap<Integer,MCParticle> mc_pos_prev = new HashMap<Integer, MCParticle>();
+
+
+
+ public void setDebug(boolean v) {
+ this._debug = v;
+ }
+
+ public void setOutputPlotFileName(String filename) {
+ outputPlotFileName = filename;
+ }
+
+ public void setHideFrame(boolean hide) {
+ hideFrame = hide;
+ }
+
+ public TruthMomentumResolutionDriver() {
+
+
+ }
+
+
+
+ public void detectorChanged(Detector detector) {
+
+
+
+ pFrame = new AIDAFrame();
+ pFrame.setTitle("Truth p Plots");
+ makePlots();
+
+
+
+ pFrame.pack();
+ pFrame.setVisible(!hideFrame);
+
+
+ }
+
+ public double getMomentum(Track track) {
+ double[] p_vec = track.getTrackStates().get(0).getMomentum();
+ return Math.sqrt(p_vec[0]*p_vec[0] + p_vec[1]*p_vec[1] +p_vec[2]*p_vec[2]);
+ }
+
+
+ public void process(EventHeader event) {
+
+
+ // Create a map between tracks and the associated MCParticle
+ List<Track> tracklist = event.get(Track.class, "MatchedTracks");
+ if(_debug) System.out.println("Number of Tracks = " + tracklist.size());
+ List<MCParticle> mcparticles = event.get(MCParticle.class).get(0);
+ if(_debug) System.out.println("Number of MC particles = " + mcparticles.size());
+ List<MCParticle> fsParticles = HPSMCParticlePlotsDriver.makeGenFSParticleList(mcparticles);
+ if(_debug) System.out.println("Number of FS MC particles = " + fsParticles.size());
+
+
+
+ MCParticle electron=null;
+ MCParticle positron=null;
+
+ int nele = 0;
+ int nposi = 0;
+ for(MCParticle fs : fsParticles) {
+ if(ParticleTypeClassifier.isElectron(fs.getPDGID())) {
+ ++nele;
+ if(electron==null) electron = fs;
+ else {
+ if(fs.getEnergy()>electron.getEnergy()) {
+ electron = fs;
+ }
+ }
+ } else if(ParticleTypeClassifier.isPositron(fs.getPDGID())) {
+ ++nposi;
+ if(positron==null) positron = fs;
+ else {
+ if(fs.getEnergy()>positron.getEnergy()) {
+ positron = fs;
+ }
+ }
+ }
+ }
+
+ if(electron!=null) this.hElectronP.fill(electron.getMomentum().magnitude());
+ if(positron!=null) this.hPositronP.fill(positron.getMomentum().magnitude());
+ int[] ntrks = {0,0,0};
+
+ for (Track trk : tracklist) {
+ double p = this.getMomentum(trk);
+ this.hTrackP.fill(p);
+ if(trk.getCharge()<0) {
+ this.hPosTrackP.fill(p);
+ if(positron!=null) {
+ if(positron.getCharge()<0) {
+ System.out.println("Error charge!!!");
+ System.exit(1);
+ }
+ this.hTruthMatchedPosTrackP.fill(p);
+ this.hTruthMatchedPosTrackPdiff.fill(p - positron.getMomentum().magnitude());
+ System.out.println("Filling pos for " + mc_pos_prev.size() + " prev ");
+ for(Map.Entry<Integer,MCParticle> prev : mc_pos_prev.entrySet()) {
+ System.out.println("prev " + prev.getKey());
+ this.hTruthMatchedPosTrackPdiffPrev[prev.getKey()].fill(p - prev.getValue().getMomentum().magnitude());
+ }
+ }
+ ++ntrks[1];
+ }
+ else {
+ this.hNegTrackP.fill(p);
+ if(electron!=null) {
+ this.hTruthMatchedNegTrackP.fill(p);
+ this.hTruthMatchedNegTrackPdiff.fill(p - electron.getMomentum().magnitude());
+ System.out.println("Filling ele for " + mc_ele_prev.size() + " prev ");
+ for(Map.Entry<Integer,MCParticle> prev : mc_ele_prev.entrySet()) {
+ System.out.println("prev " + prev.getKey());
+ this.hTruthMatchedNegTrackPdiffPrev[prev.getKey()].fill(p - prev.getValue().getMomentum().magnitude());
+ }
+ }
+ ++ntrks[2];
+ }
+
+ ++totalTracks;
+ ++ntrks[0];
+ //int q = trk.getCharge();
+ }
+ this.hNTracks.fill(ntrks[0]);
+ this.hNPosTracks.fill(ntrks[1]);
+ this.hNNegTracks.fill(ntrks[2]);
+
+ if(ntrks[1]>0) hNPositronsForTrack.fill(nposi);
+ if(ntrks[2]>0) hNElectronsForTrack.fill(nele);
+ if(ntrks[2]>0) hNPositronsForTrackInv.fill(nposi);
+ if(ntrks[1]>0) hNElectronsForTrackInv.fill(nele);
+
+ for(int i=0;i<tracklist.size();++i) trkCountVsEventPlot.fill(event.getEventNumber());
+
+ //Save to list of previous truth particles
+ if(electron!=null) {
+ mc_ele_prev = this.updatePrevMap(mc_ele_prev);
+ mc_ele_prev.put(0, electron);
+ }
+ if(positron!=null) {
+ mc_pos_prev = this.updatePrevMap(mc_pos_prev);
+ mc_pos_prev.put(0, positron);
+ }
+
+ }
+
+ private HashMap<Integer,MCParticle> updatePrevMap(HashMap<Integer,MCParticle> map) {
+ HashMap<Integer,MCParticle> newmap = new HashMap<Integer,MCParticle>();
+ for (Map.Entry<Integer, MCParticle> entry : map.entrySet()) {
+ if(entry.getKey()<4) {
+ System.out.println("Key e = " + entry.getKey() + ", Value = " + entry.getValue());
+ newmap.put(entry.getKey()+1, entry.getValue());
+ }
+ }
+ return newmap;
+ }
+
+
+
+ private MCParticle getHighestEnergyParticle(int pdgId, List<MCParticle> fsParticles) {
+ MCParticle particle = null;
+ for(MCParticle fs : fsParticles) {
+ int fsPdg = fs.getPDGID();
+ if(fsPdg==pdgId) {
+ if(particle==null) {
+ particle = fs;
+ }
+ else {
+ if(fs.getEnergy()>particle.getEnergy()) {
+ particle = fs;
+ }
+ }
+ }
+ }
+ return particle;
+ }
+
+
+ public void endOfData() {
+
+ System.out.println("Total Number of Tracks Found = "+totalTracks);
+
+ if (outputPlotFileName != "")
+ try {
+ aida.saveAs(outputPlotFileName);
+ } catch (IOException ex) {
+ Logger.getLogger(TrigRateDriver.class.getName()).log(Level.SEVERE, "Couldn't save aida plots to file " + outputPlotFileName, ex);
+ }
+
+ }
+
+ private void makePlots() {
+
+
+ pPlotter = af.createPlotterFactory().create("Truth p Plots");
+ pPlotter.setTitle("Truth p Plots");
+ pFrame.addPlotter(pPlotter);
+ IPlotterStyle style0 = pPlotter.style();
+ //style0.dataStyle().fillStyle().setColor("yellow");
+ //style0.dataStyle().errorBarStyle().setVisible(false);
+ pPlotter.createRegions(2, 6);
+ hTrackP = aida.cloud1D("Track p");
+ hPosTrackP = aida.cloud1D("Track p q>0");
+ hNegTrackP = aida.cloud1D("Track p q<0");
+ hTruthMatchedPosTrackP = aida.cloud1D("Track p q>0 e+ match");
+ hTruthMatchedNegTrackP = aida.cloud1D("Track p q<0 e- match");
+ hTruthMatchedPosTrackPdiff = aida.cloud1D("Track p - p(e+) q>0 ");
+ hTruthMatchedNegTrackPdiff = aida.cloud1D("Track p - p(e-) q<0");
+ hNTracks = aida.cloud1D("Ntrks");
+ trkCountVsEventPlot = aida.histogram1D("Number of Tracks vs Event Nr", 501, -0.5, 500.5);
+ trkCountVsEventPlot.annotation().addItem("xAxisLabel", "Event Number");
+ hNPosTracks = aida.cloud1D("Ntrks q>0");
+ hNNegTracks = aida.cloud1D("Ntrks q<0");
+ hElectronP = aida.cloud1D("Electron Momentum");
+ hPositronP = aida.cloud1D("Positron Momentum");
+ hNPositronsForTrack = aida.cloud1D("N positrons given track w/ q>0");
+ hNElectronsForTrack = aida.cloud1D("N electrons given track w/ q<0");
+ hNPositronsForTrackInv = aida.cloud1D("N positrons given track w/ q<0");
+ hNElectronsForTrackInv = aida.cloud1D("N electrons given track w/ q>0");
+
+
+ pPlotter.region(0).plot(hTrackP);
+ pPlotter.region(0).plot(hPosTrackP);
+ pPlotter.region(0).plot(hNegTrackP);
+ pPlotter.region(6).plot(hNTracks);
+ pPlotter.region(6).plot(hNPosTracks);
+ pPlotter.region(6).plot(hNNegTracks);
+ pPlotter.region(1).plot(hElectronP);
+ pPlotter.region(7).plot(hPositronP);
+ pPlotter.region(2).plot(this.hTruthMatchedPosTrackP);
+ pPlotter.region(8).plot(this.hTruthMatchedNegTrackP);
+ pPlotter.region(3).plot(this.hTruthMatchedPosTrackPdiff);
+ pPlotter.region(9).plot(this.hTruthMatchedNegTrackPdiff);
+ pPlotter.region(4).plot(this.hNPositronsForTrack);
+ pPlotter.region(10).plot(this.hNElectronsForTrack);
+ pPlotter.region(5).plot(this.hNPositronsForTrackInv);
+ pPlotter.region(11).plot(this.hNElectronsForTrackInv);
+
+
+ pPlotter2 = af.createPlotterFactory().create("Truth p Plots");
+ pPlotter2.setTitle("Truth p Plots");
+ pFrame.addPlotter(pPlotter2);
+ style0 = pPlotter2.style();
+ //style0.dataStyle().fillStyle().setColor("yellow");
+ //style0.dataStyle().errorBarStyle().setVisible(false);
+ pPlotter2.createRegions(1, 2);
+ pPlotter2.region(1).plot(trkCountVsEventPlot);
+
+
+ pPlotter3 = af.createPlotterFactory().create("Truth p Plots");
+ pPlotter3.setTitle("Prev BS's");
+ pFrame.addPlotter(pPlotter3);
+ pPlotter3.createRegions(2, 5);
+
+ for(int i=0;i<5;++i) {
+ hTruthMatchedPosTrackPdiffPrev[i] = aida.cloud1D("Track p - p(e+) q>0 BC=-" + i);
+ hTruthMatchedNegTrackPdiffPrev[i] = aida.cloud1D("Track p - p(e-) q<0 BC=-" + i);
+ pPlotter3.region(i).plot(hTruthMatchedPosTrackPdiffPrev[i]);
+ pPlotter3.region(5+i).plot(hTruthMatchedNegTrackPdiffPrev[i]);
+ }
+
+ }
+
+
+
+
+
+
+
+}
hps-java/src/main/java/org/lcsim/hps/analysis/ecal
diff -u -r1.4 -r1.5
--- HPSMCParticlePlotsDriver.java 19 Sep 2012 16:24:17 -0000 1.4
+++ HPSMCParticlePlotsDriver.java 21 Sep 2012 18:36:34 -0000 1.5
@@ -1,8 +1,6 @@
package org.lcsim.hps.analysis.ecal;
-import hep.aida.ICloud1D;
-import hep.aida.ICloud2D;
-import hep.aida.IHistogram1D;
+import hep.aida.*;
import java.util.ArrayList;
import java.util.Comparator;
@@ -11,6 +9,7 @@
import org.lcsim.event.EventHeader;
import org.lcsim.event.MCParticle;
import org.lcsim.event.util.ParticleTypeClassifier;
+import org.lcsim.hps.monitoring.AIDAFrame;
import org.lcsim.util.Driver;
import org.lcsim.util.aida.AIDA;
@@ -18,14 +17,18 @@
* Diagnostic plots for HPS ECal.
*
* @author Jeremy McCormick <[log in to unmask]>
- * @version $Id: HPSMCParticlePlotsDriver.java,v 1.4 2012/09/19 16:24:17 phansson Exp $
+ * @version $Id: HPSMCParticlePlotsDriver.java,v 1.5 2012/09/21 18:36:34 phansson Exp $
*/
public class HPSMCParticlePlotsDriver extends Driver {
AIDA aida = AIDA.defaultInstance();
- // MCParticle plots.
+ private AIDAFrame pFrame;
+ IAnalysisFactory af = aida.analysisFactory();
+ public boolean hideFrame = false;
+ // MCParticle plots.
ICloud1D primaryEPlot;
ICloud1D fsCountPlot;
+ IHistogram1D fsCountVsEventPlot;
ICloud1D fsCountTypePlot;
ICloud1D fsCountEventTypePlot;
ICloud1D fsCountEventTypePlot2;
@@ -67,6 +70,9 @@
public void startOfData() {
fsCountPlot = aida.cloud1D("MCParticle: Number of Final State Particles");
fsCountPlot.annotation().addItem("xAxisLabel", "Number of FS Particles");
+
+ fsCountVsEventPlot = aida.histogram1D("MCParticle: Number of Final State Particles vs Event Nr", 501, -0.5, 500.5);
+ fsCountVsEventPlot.annotation().addItem("xAxisLabel", "Event Number");
fsCountTypePlot = aida.cloud1D("MCParticle: Number of Final State Particles Type");
fsCountTypePlot.annotation().addItem("xAxisLabel", "Number of FS Particles of Type");
@@ -138,6 +144,20 @@
eventEPlot = aida.cloud1D("MCParticle: Total Gen FS Electron E in Event");
eventEPlot.annotation().addItem("xAxisLabel", "E [GeV]");
+
+ pFrame = new AIDAFrame();
+ pFrame.setTitle("Truth MC Particle Plots");
+ IPlotter pPlotter = af.createPlotterFactory().create("MCParticle Plots");
+ pPlotter.setTitle("MCParticle Event Types");
+ pPlotter.createRegions(1,2);
+ pPlotter.region(0).plot(fsCountEventTypePlot);
+ pPlotter.region(1).plot(fsCountVsEventPlot);
+ pFrame.addPlotter(pPlotter);
+ pFrame.pack();
+ pFrame.setVisible(!hideFrame);
+
+
+
}
@Override
@@ -152,6 +172,9 @@
//System.out.println("fsParticles="+fsParticles.size());
fsCountPlot.fill(fsParticles.size());
+ for (int i=0;i<fsParticles.size();++i) fsCountVsEventPlot.fill(event.getEventNumber());
+
+
double fsGammaEmax = this.getHighestPhotonE(fsParticles);
//if (fsGammaEmax>0.5) System.out.println("Emax large " + fsGammaEmax);
int[] nelectrons = {0,0};
@@ -279,7 +302,7 @@
return primary;
}
- private static List<MCParticle> makeGenFSParticleList(List<MCParticle> mcparticles) {
+ public static List<MCParticle> makeGenFSParticleList(List<MCParticle> mcparticles) {
List<MCParticle> fsParticles = new ArrayList<MCParticle>();
for (MCParticle mcparticle : mcparticles) {
if (mcparticle.getGeneratorStatus() == MCParticle.FINAL_STATE) {