Author: [log in to unmask]
Date: Mon Aug 3 14:06:45 2015
New Revision: 3328
Log:
Added driver to analyze Moller Monte Carlo particle data.
Added:
java/trunk/users/src/main/java/org/hps/users/kmccarty/ParticleMCAnalysisDriver.java
Added: java/trunk/users/src/main/java/org/hps/users/kmccarty/ParticleMCAnalysisDriver.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/kmccarty/ParticleMCAnalysisDriver.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/kmccarty/ParticleMCAnalysisDriver.java Mon Aug 3 14:06:45 2015
@@ -0,0 +1,121 @@
+package org.hps.users.kmccarty;
+
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogram2D;
+import hep.physics.vec.Hep3Vector;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+public class ParticleMCAnalysisDriver extends Driver {
+ // Store collection names.
+ private String particleCollectionName = "MCParticle";
+
+ // Declare plots.
+ private AIDA aida = AIDA.defaultInstance();
+ private IHistogram1D chargedTracksPlot = aida.histogram1D("MC Analysis/Event Tracks", 10, -0.5, 9.5);
+ private IHistogram1D allPlot = aida.histogram1D("MC Analysis/Electron Energy Distribution", 110, 0, 2.2);
+ private IHistogram1D electronPlot = aida.histogram1D("MC Analysis/Electron Energy Distribution", 110, 0, 2.2);
+ private IHistogram1D positronPlot = aida.histogram1D("MC Analysis/Positron Energy Distribution", 110, 0, 2.2);
+ private IHistogram1D anglePlot = aida.histogram1D("MC Analysis/Positron\\Electron Pair Angle Distribution", 90, 0, 180);
+ private IHistogram1D momentumSumPlot = aida.histogram1D("MC Analysis/Positron\\Electron Momentum Sum Distribution", 220, 0, 2.2);
+ private IHistogram2D momentumPlot = aida.histogram2D("MC Analysis/Particle Momentum Distribution", 110, 0.0, 1.1, 110, 0.0, 1.1);
+ private IHistogram2D momentumSum2DPlot = aida.histogram2D("MC Analysis/Positron\\Electron 2D Momentum Distribution", 55, 0, 1.1, 55, 0, 1.1);
+
+ @Override
+ public void process(EventHeader event) {
+ // Skip the event if there is no Monte Carlo collection.
+ if(!event.hasCollection(MCParticle.class, particleCollectionName)) {
+ return;
+ }
+
+ // Get the list of Monte Carlo particles.
+ List<MCParticle> particleList = event.get(MCParticle.class, particleCollectionName);
+
+ // Track the positive and negative particles.
+ List<MCParticle> electronList = new ArrayList<MCParticle>();
+ List<MCParticle> positronList = new ArrayList<MCParticle>();
+
+ // Count the number of particles in the event.
+ int chargedParticles = 0;
+
+ // Iterate through the particles.
+ for(MCParticle particle : particleList) {
+ // Look at only t = 0 particles.
+ if(particle.getProductionTime() == 0) {
+ // Plot the x/y momentum of each particle.
+ momentumPlot.fill(particle.getMomentum().x(), particle.getMomentum().y());
+
+ // If the particle is charged, increment the charged
+ // particle count.
+ if(particle.getCharge() > 0) {
+ chargedParticles++;
+ }
+
+ // Populate the general momentum plot.
+ allPlot.fill(particle.getMomentum().magnitude());
+ momentumPlot.fill(particle.getMomentum().x(), particle.getMomentum().y());
+
+ // Store each particle based on its PID and populate
+ // the appropriate plot.
+ if(particle.getPDGID() == 11) {
+ electronList.add(particle);
+ electronPlot.fill(particle.getMomentum().magnitude());
+ } else if(particle.getPDGID() == -11) {
+ positronList.add(particle);
+ positronPlot.fill(particle.getMomentum().magnitude());
+ }
+ }
+ }
+
+ // Populate the charged particles plot.
+ chargedTracksPlot.fill(chargedParticles);
+
+ // Form all electron/positron pairs.
+ List<MCParticle[]> pairList = new ArrayList<MCParticle[]>();
+ for(MCParticle electron : electronList) {
+ for(MCParticle positron : positronList) {
+ pairList.add(new MCParticle[] { electron, positron });
+ }
+ }
+
+ // Plot the positron/electron pair distributions.
+ for(MCParticle[] pair : pairList) {
+ anglePlot.fill(getVectorAngle(pair[0].getMomentum(), pair[1].getMomentum()));
+ momentumSumPlot.fill(getVectorSum(pair[0].getMomentum(), pair[1].getMomentum()));
+ momentumSum2DPlot.fill(pair[0].getMomentum().magnitude(), pair[1].getMomentum().magnitude());
+ }
+ }
+
+ private static final double getVectorSum(Hep3Vector v1, Hep3Vector v2) {
+ // Calculate the sum of the sum of the vector components, squared.
+ double sum = 0;
+ for(int i = 0; i < 3; i++) {
+ sum += ((v1.v()[i] + v2.v()[i]) * (v1.v()[i] + v2.v()[i]));
+ }
+
+ // Return the square root of the sum.
+ return Math.sqrt(sum);
+ }
+
+ private static final double getVectorAngle(Hep3Vector v1, Hep3Vector v2) {
+ // The vector angle is defined as Acos[(v1 · v2) / (âv1â à âv2â)]
+ return Math.acos(getDotProduct(v1, v2) / (v1.magnitude() * v2.magnitude())) / Math.PI * 180.0;
+ }
+
+ private static final double getDotProduct(Hep3Vector v1, Hep3Vector v2) {
+ // Calculate the sum of the vector element products.
+ int product = 0;
+ for(int i = 0; i < 3; i++) {
+ product += (v1.v()[i] * v2.v()[i]);
+ }
+
+ // Return the result.
+ return product;
+ }
+}
|