 |
 |
| Commit in hps-java/src/main/java/org/lcsim/hps/analysis/ecal on MAIN | |||
| HPSMCParticlePlotsDriver.java | +137 | -12 | 1.2 -> 1.3 |
Added general truth plots.
diff -u -r1.2 -r1.3 --- HPSMCParticlePlotsDriver.java 9 Jan 2012 22:39:08 -0000 1.2 +++ HPSMCParticlePlotsDriver.java 7 Sep 2012 04:35:27 -0000 1.3 @@ -1,6 +1,7 @@
package org.lcsim.hps.analysis.ecal; import hep.aida.ICloud1D;
+import hep.aida.ICloud2D;
import java.util.ArrayList; import java.util.Comparator;
@@ -16,7 +17,7 @@
* Diagnostic plots for HPS ECal. * * @author Jeremy McCormick <[log in to unmask]>
- * @version $Id: HPSMCParticlePlotsDriver.java,v 1.2 2012/01/09 22:39:08 meeg Exp $
+ * @version $Id: HPSMCParticlePlotsDriver.java,v 1.3 2012/09/07 04:35:27 phansson Exp $
*/
public class HPSMCParticlePlotsDriver extends Driver {
@@ -24,14 +25,27 @@
// MCParticle plots. ICloud1D primaryEPlot; ICloud1D fsCountPlot;
+ ICloud1D fsCountTypePlot; + ICloud1D fsCountEventTypePlot; + ICloud1D fsCountEventTypePlot2; + ICloud1D fsCountTypePlot500;
ICloud1D fsEPlot; ICloud1D fsGammaEPlot; ICloud1D fsElectronEPlot; ICloud1D fsPositronEPlot; ICloud1D fsGammaThetaPlot;
+ ICloud1D fsGammaThetayPlot; + ICloud1D fsGammaThetayTrigPlot; + ICloud2D fsGammaThetayEPlot;
ICloud1D fsElectronThetaPlot;
+ ICloud1D fsElectronThetayPlot; + ICloud1D fsElectronThetayTrigPlot; + ICloud2D fsElectronThetayEPlot;
ICloud1D fsPositronThetaPlot;
- ICloud1D eventEPlot;
+ ICloud1D fsPositronThetayPlot; + ICloud1D fsPositronThetayTrigPlot; + ICloud2D fsPositronThetayEPlot; + ICloud1D eventEPlot;
class MCParticleEComparator implements Comparator<MCParticle> {
@@ -53,6 +67,18 @@
fsCountPlot = aida.cloud1D("MCParticle: Number of Final State Particles");
fsCountPlot.annotation().addItem("xAxisLabel", "Number of FS Particles");
+ fsCountTypePlot = aida.cloud1D("MCParticle: Number of Final State Particles Type");
+ fsCountTypePlot.annotation().addItem("xAxisLabel", "Number of FS Particles of Type");
+
+ fsCountTypePlot500 = aida.cloud1D("MCParticle: Number of Final State Particles Type E>0.5GeV");
+ fsCountTypePlot500.annotation().addItem("xAxisLabel", "Number of FS Particles of Type E>0.5GeV");
+
+ fsCountEventTypePlot = aida.cloud1D("MCParticle: Number of Final State Types");
+ fsCountEventTypePlot.annotation().addItem("xAxisLabel", "Number of FS Types");
+
+ fsCountEventTypePlot2 = aida.cloud1D("MCParticle: Number of Final State Types Gamma E>500");
+ fsCountEventTypePlot2.annotation().addItem("xAxisLabel", "Number of FS Types Gamma E>500");
+
fsEPlot = aida.cloud1D("MCParticle: FS Particle E");
fsEPlot.annotation().addItem("xAxisLabel", "Particle E [GeV]");
@@ -66,14 +92,45 @@
fsPositronEPlot.annotation().addItem("xAxisLabel", "Particle E [GeV]");
fsGammaThetaPlot = aida.cloud1D("MCParticle: FS Gamma Theta");
- fsGammaThetaPlot.annotation().addItem("xAxisLabel", "Particle angle [degrees]");
+ fsGammaThetaPlot.annotation().addItem("xAxisLabel", "Particle angle [rad]");
+ fsGammaThetayPlot = aida.cloud1D("MCParticle: FS Gamma Thetay");
+ fsGammaThetayPlot.annotation().addItem("xAxisLabel", "Particle Thetay angle [rad]");
+
+ fsGammaThetayTrigPlot = aida.cloud1D("MCParticle: FS Gamma Thetay Trig");
+ fsGammaThetayTrigPlot.annotation().addItem("xAxisLabel", "Particle Thetay angle [rad]");
+
+ fsGammaThetayEPlot = aida.cloud2D("MCParticle: FS Gamma Thetay vs E");
+ fsGammaThetayEPlot.annotation().addItem("xAxisLabel", "Particle Thetay angle [rad]");
+ fsGammaThetayEPlot.annotation().addItem("yAxisLabel", "Particle Energy [GeV]");
+
fsElectronThetaPlot = aida.cloud1D("MCParticle: FS Electron Theta");
- fsElectronThetaPlot.annotation().addItem("xAxisLabel", "Particle angle [degrees]");
+ fsElectronThetaPlot.annotation().addItem("xAxisLabel", "Particle angle [rad]");
+
+ fsElectronThetayPlot = aida.cloud1D("MCParticle: FS Electron Thetay");
+ fsElectronThetayPlot.annotation().addItem("xAxisLabel", "Particle Thetay angle [rad]");
+
+ fsElectronThetayTrigPlot = aida.cloud1D("MCParticle: FS Electron Thetay Trig");
+ fsElectronThetayTrigPlot.annotation().addItem("xAxisLabel", "Particle Thetay angle [rad]");
+
+ fsElectronThetayEPlot = aida.cloud2D("MCParticle: FS Electron Thetay vs E");
+ fsElectronThetayEPlot.annotation().addItem("xAxisLabel", "Particle Thetay angle [rad]");
+ fsElectronThetayEPlot.annotation().addItem("yAxisLabel", "Particle Energy [GeV]");
fsPositronThetaPlot = aida.cloud1D("MCParticle: FS Positron Theta");
- fsPositronThetaPlot.annotation().addItem("xAxisLabel", "Particle angle [degrees]");
+ fsPositronThetaPlot.annotation().addItem("xAxisLabel", "Particle angle [rad]");
+
+ fsPositronThetayPlot = aida.cloud1D("MCParticle: FS Positron Thetay");
+ fsPositronThetayPlot.annotation().addItem("xAxisLabel", "Particle Thetay angle [rad]");
+
+ fsPositronThetayTrigPlot = aida.cloud1D("MCParticle: FS Positron Thetay Trig");
+ fsPositronThetayTrigPlot.annotation().addItem("xAxisLabel", "Particle Thetay angle [rad]");
+ fsPositronThetayEPlot = aida.cloud2D("MCParticle: FS Positron Thetay vs E");
+ fsPositronThetayEPlot.annotation().addItem("xAxisLabel", "Particle Thetay angle [rad]");
+ fsPositronThetayEPlot.annotation().addItem("yAxisLabel", "Particle Energy [GeV]");
+
+
primaryEPlot = aida.cloud1D("MCParticle: Highest Primary E in Event");
primaryEPlot.annotation().addItem("xAxisLabel", "E [GeV]");
@@ -93,24 +150,57 @@
//System.out.println("fsParticles="+fsParticles.size());
fsCountPlot.fill(fsParticles.size());
-
+
+ double fsGammaEmax = this.getHighestPhotonE(fsParticles);
+ //if (fsGammaEmax>0.5) System.out.println("Emax large " + fsGammaEmax);
+ int[] nelectrons = {0,0};
+ int[] npositrons = {0,0};
+ int[] ngammas = {0,0};
+ int count = 0;
+ double trigThr = 0.5;
for (MCParticle fs : fsParticles) {
+ //System.out.println("Index " + count);
+
double fsE = fs.getEnergy(); double theta = Math.atan2(Math.sqrt(fs.getPX() * fs.getPX() + fs.getPY() * fs.getPY()), fs.getPZ());
+ double thetay = Math.atan2(fs.getPX(), fs.getPZ());
int fsPdg = fs.getPDGID(); fsEPlot.fill(fsE);
+ fsCountTypePlot.fill(fsPdg); + if(fsE>0.5) fsCountTypePlot500.fill(fsPdg);
if (ParticleTypeClassifier.isElectron(fsPdg)) {
fsElectronEPlot.fill(fsE);
- fsElectronThetaPlot.fill(Math.toDegrees(theta));
+ fsElectronThetaPlot.fill(theta); + fsElectronThetayPlot.fill(Math.abs(thetay)); + if(fsE>trigThr) fsElectronThetayTrigPlot.fill(Math.abs(thetay)); + fsElectronThetayEPlot.fill(Math.abs(thetay),fsE); + nelectrons[0]++; + if(fsGammaEmax>0.5) nelectrons[1]++;
} else if (ParticleTypeClassifier.isPositron(fsPdg)) {
fsPositronEPlot.fill(fsE);
- fsPositronThetaPlot.fill(Math.toDegrees(theta));
+ fsPositronThetaPlot.fill(theta); + fsPositronThetayPlot.fill(Math.abs(thetay)); + if(fsE>trigThr) fsPositronThetayTrigPlot.fill(Math.abs(thetay)); + fsPositronThetayEPlot.fill(Math.abs(thetay),fsE); + npositrons[0]++; + if(fsGammaEmax>0.5) npositrons[1]++;
} else if (ParticleTypeClassifier.isPhoton(fsPdg)) {
fsGammaEPlot.fill(fsE);
- fsGammaThetaPlot.fill(Math.toDegrees(theta));
+ fsGammaThetaPlot.fill(theta);
+ fsGammaThetayPlot.fill(Math.abs(thetay));
+ if(fsE>trigThr) fsGammaThetayTrigPlot.fill(Math.abs(thetay));
+ fsGammaThetayEPlot.fill(Math.abs(thetay),fsE);
+ ngammas[0]++;
+ if(fsGammaEmax>0.5) {
+ ngammas[1]++;
+ //System.out.println("Counting high E gamma at count "+ count);
+ }
} }
+ fsCountEventTypePlot.fill(getEventTypeId(nelectrons[0],npositrons[0],ngammas[0])); + fsCountEventTypePlot2.fill(getEventTypeId(nelectrons[1],npositrons[1],ngammas[1])); +
// Sort MCParticles on energy. //Collections.sort(fsParticles, new MCParticleEComparator());
@@ -126,12 +216,47 @@
primaryEPlot.fill(primaryE); }
- // event energy - double eventE = getPrimaryE(fsParticles);
+ // event electron energy + double eventE = getPrimaryElectronE(fsParticles);
eventEPlot.fill(eventE); }
- private double getPrimaryE(List<MCParticle> particles) {
+
+ public int getEventTypeId(int ne, int np, int ng) {
+ //get a unique ID depending on final state particle count
+ if(ne==0 && np==0 & ng==0) return -1;
+ if(ne==1 && np==0 & ng==0) return 1;
+ if(ne==0 && np==1 & ng==0) return 2;
+ if(ne==0 && np==0 & ng==1) return 3;
+ if(ne==1 && np==1 & ng==0) return 4;
+ if(ne==1 && np==0 & ng==1) return 5;
+ if(ne==0 && np==1 & ng==1) return 6;
+ if(ne==1 && np==1 & ng==1) return 7;
+ if(ne>1) return 8;
+ if(np>1) return 9;
+ if(ng>1) return 10;
+
+ return 0;
+ }
+
+ public double getHighestPhotonE(List<MCParticle> particles) {
+ double Emax = -1;
+ double E=0;
+ int count = 0;
+ for (MCParticle particle : particles) {
+ if (ParticleTypeClassifier.isPhoton(particle.getPDGID())) {
+ E = particle.getEnergy();
+ if(E>Emax) {
+ Emax = E;
+ //System.out.println("Emax from photon with index " + count);
+ }
+ count++;
+ }
+ }
+ return Emax;
+ }
+
+ private double getPrimaryElectronE(List<MCParticle> particles) {
double totalE = 0;
for (MCParticle particle : particles) {
if (Math.abs(particle.getPDGID()) == 11) {
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1