hps-java/src/main/java/org/lcsim/hps/analysis/ecal
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) {