Author: [log in to unmask]
Date: Tue Oct 20 18:57:00 2015
New Revision: 3864
Log:
Added work-in-progress Moller/trident analysis driver.
Added:
java/trunk/users/src/main/java/org/hps/users/kmccarty/TriggerProcessAnalysisDriver.java (with props)
Modified:
java/trunk/users/src/main/java/org/hps/users/kmccarty/plots/formatter/TriggerPlotsFormat.java
Added: java/trunk/users/src/main/java/org/hps/users/kmccarty/TriggerProcessAnalysisDriver.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/kmccarty/TriggerProcessAnalysisDriver.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/kmccarty/TriggerProcessAnalysisDriver.java Tue Oct 20 18:57:00 2015
@@ -0,0 +1,390 @@
+package org.hps.users.kmccarty;
+
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogram2D;
+import hep.physics.vec.VecOp;
+
+import java.util.ArrayList;
+import java.util.HashSet;
+import java.util.List;
+import java.util.Set;
+
+import org.hps.record.triggerbank.TriggerModule;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+public class TriggerProcessAnalysisDriver extends Driver {
+ private int eventsProcessed = 0;
+ private int møllersProcessed = 0;
+ private int tridentsProcessed = 0;
+ private double timeCoincidence = 2.5;
+ private AIDA aida = AIDA.defaultInstance();
+ private String clusterCollectionName = "EcalClustersCorr";
+ private String particleCollectionName = "FinalStateParticles";
+
+ // Define trident cluster-track matched condition plots.
+ private IHistogram1D trctmInvariantMass = aida.histogram1D("Tridents CTMatched/Invariant Mass", 140, 0.0, 0.070);
+ private IHistogram1D trctmInstancesInEvent = aida.histogram1D("Tridents CTMatched/Instances in Event", 9, 0.5, 9.5);
+ private IHistogram1D trctmEnergySum1D = aida.histogram1D("Tridents CTMatched/Cluster Energy Sum", 150, 0.000, 1.500);
+ private IHistogram1D trctmMomentumSum1D = aida.histogram1D("Tridents CTMatched/Track Momentum Sum", 150, 0.000, 1.500);
+ private IHistogram1D trctmElectronEnergy = aida.histogram1D("Tridents CTMatched/Electron Cluster Energy", 150, 0.000, 1.500);
+ private IHistogram1D trctmElectronMomentum = aida.histogram1D("Tridents CTMatched/Electron Track Momentum", 150, 0.000, 1.500);
+ private IHistogram1D trctmPositronEnergy = aida.histogram1D("Tridents CTMatched/Positron Cluster Energy", 150, 0.000, 1.500);
+ private IHistogram1D trctmPositronMomentum = aida.histogram1D("Tridents CTMatched/Positron Track Momentum", 150, 0.000, 1.500);
+ private IHistogram1D trctmTimeCoincidence = aida.histogram1D("Tridents CTMatched/Time Coincidence", 100, -4, 4);
+ private IHistogram2D trctmClusterPosition = aida.histogram2D("Tridents CTMatched/Cluster Seed Position", 46, -23, 23, 11, -5.5, 5.5);
+ private IHistogram2D trctmEnergySum2D = aida.histogram2D("Tridents CTMatched/Cluster Energy Sum 2D", 300, 0.000, 1.500, 300, 0.000, 1.500);
+ private IHistogram2D trctmMomentumSum2D = aida.histogram2D("Tridents CTMatched/Track Momentum Sum 2D", 300, 0.000, 1.500, 300, 0.000, 1.500);
+
+ // Define the Møller cluster-track matched condition plots.
+ private IHistogram1D møctmInvariantMass = aida.histogram1D("Møller CTMatched/Invariant Mass", 140, 0.0, 0.070);
+ private IHistogram1D møctmInstancesInEvent = aida.histogram1D("Møller CTMatched/Instances in Event", 9, 0.5, 9.5);
+ private IHistogram1D møctmEnergySum1D = aida.histogram1D("Møller CTMatched/Cluster Energy Sum", 150, 0.000, 1.500);
+ private IHistogram1D møctmMomentumSum1D = aida.histogram1D("Møller CTMatched/Track Momentum Sum", 150, 0.000, 1.500);
+ private IHistogram1D møctmElectronEnergy = aida.histogram1D("Møller CTMatched/Electron Cluster Energy", 150, 0.000, 1.500);
+ private IHistogram1D møctmElectronMomentum = aida.histogram1D("Møller CTMatched/Electron Track Momentum", 150, 0.000, 1.500);
+ private IHistogram1D møctmTimeCoincidence = aida.histogram1D("Møller CTMatched/Time Coincidence", 100, -4, 4);
+ private IHistogram2D møctmClusterPosition = aida.histogram2D("Møller CTMatched/Cluster Seed Position", 46, -23, 23, 11, -5.5, 5.5);
+ private IHistogram2D møctmEnergySum2D = aida.histogram2D("Møller CTMatched/Cluster Energy Sum 2D", 300, 0.000, 1.500, 300, 0.000, 1.500);
+ private IHistogram2D møctmMomentumSum2D = aida.histogram2D("Møller CTMatched/Track Momentum Sum 2D", 300, 0.000, 1.500, 300, 0.000, 1.500);
+
+ @Override
+ public void endOfData() {
+ // Calculate the scaling factor for Hertz.
+ double scale = 19000.0 / eventsProcessed;
+
+ System.out.println("Processed " + eventsProcessed + " events.");
+ System.out.println("Processed " + møllersProcessed + " Møller events");
+ System.out.println("\tAcceptance :: " + (100.0 * møllersProcessed / eventsProcessed) + "%");
+ System.out.println("\tRate :: " + (møllersProcessed * scale) + "Hz");
+ System.out.println("Processed " + tridentsProcessed + " trident events");
+ System.out.println("\tAcceptance :: " + (100.0 * tridentsProcessed / eventsProcessed) + "%");
+ System.out.println("\tRate :: " + (tridentsProcessed * scale) + "Hz");
+
+ // Scale the cluster-track matched Møller plots.
+ møctmInvariantMass.scale(scale);
+ møctmInstancesInEvent.scale(scale);
+ møctmEnergySum1D.scale(scale);
+ møctmMomentumSum1D.scale(scale);
+ møctmElectronEnergy.scale(scale);
+ møctmElectronMomentum.scale(scale);
+ møctmTimeCoincidence.scale(scale);
+ møctmClusterPosition.scale(scale);
+ møctmEnergySum2D.scale(scale);
+ møctmMomentumSum2D.scale(scale);
+
+ // Scale the cluster-track matched trident plots.
+ trctmInvariantMass.scale(scale);
+ trctmInstancesInEvent.scale(scale);
+ trctmEnergySum1D.scale(scale);
+ trctmMomentumSum1D.scale(scale);
+ trctmElectronEnergy.scale(scale);
+ trctmElectronMomentum.scale(scale);
+ trctmPositronEnergy.scale(scale);
+ trctmPositronMomentum.scale(scale);
+ trctmTimeCoincidence.scale(scale);
+ trctmClusterPosition.scale(scale);
+ trctmEnergySum2D.scale(scale);
+ trctmMomentumSum2D.scale(scale);
+ }
+
+ @Override
+ public void process(EventHeader event) {
+ // Check whether the SVT was active in this event and, if so,
+ // skip it.
+ final String[] flagNames = { "svt_bias_good", "svt_burstmode_noise_good", "svt_position_good" };
+ boolean svtGood = true;
+ for(int i = 0; i < flagNames.length; i++) {
+ int[] flag = event.getIntegerParameters().get(flagNames[i]);
+ if(flag == null || flag[0] == 0) {
+ svtGood = false;
+ }
+ }
+ if(!svtGood) { return; }
+
+ // Track the number of events with good SVT.
+ eventsProcessed++;
+
+ // Check if the event has a collection of tracks. If it exists,
+ // extract it. Otherwise, skip the event.
+ if(!event.hasCollection(ReconstructedParticle.class, particleCollectionName)) {
+ return;
+ }
+ List<ReconstructedParticle> trackList = event.get(ReconstructedParticle.class, particleCollectionName);
+
+ // Check if the event has a collection of clusters. If it
+ // exists, extract it. Otherwise, skip the event.
+ if(!event.hasCollection(Cluster.class, clusterCollectionName)) {
+ return;
+ }
+ List<Cluster> clusterList = event.get(Cluster.class, clusterCollectionName);
+
+ // Get cluster-track matched top/bottom pairs.
+ List<ReconstructedParticle[]> ctMatchedPairs = getTopBottomTracksCTMatched(trackList);
+
+ // Get the trident and Møller tracks for the matched track
+ // and cluster pair condition sets.
+ List<ReconstructedParticle[]> møllers = getMøllerTracksCTMatched(ctMatchedPairs);
+ List<ReconstructedParticle[]> tridents = getTridentTracksCTMatched(ctMatchedPairs);
+
+ // Track how many events had tridents and Møllers.
+ if(!møllers.isEmpty()) { møllersProcessed++; }
+ if(!tridents.isEmpty()) { tridentsProcessed++; }
+
+ // Produce Møller cluster-track matched plots.
+ møctmInstancesInEvent.fill(møllers.size());
+ for(ReconstructedParticle[] pair : møllers) {
+ // Get the track clusters.
+ Cluster[] trackClusters = { pair[0].getClusters().get(0), pair[1].getClusters().get(0) };
+
+ // Populate the cluster plots.
+ møctmElectronEnergy.fill(trackClusters[0].getEnergy());
+ møctmElectronEnergy.fill(trackClusters[1].getEnergy());
+ møctmEnergySum1D.fill(TriggerModule.getValueEnergySum(trackClusters));
+ møctmEnergySum2D.fill(trackClusters[0].getEnergy(), trackClusters[1].getEnergy());
+ møctmTimeCoincidence.fill(TriggerModule.getClusterTime(trackClusters[0]) - TriggerModule.getClusterTime(trackClusters[1]));
+ møctmClusterPosition.fill(TriggerModule.getClusterXIndex(trackClusters[0]), TriggerModule.getClusterYIndex(trackClusters[0]));
+ møctmClusterPosition.fill(TriggerModule.getClusterXIndex(trackClusters[1]), TriggerModule.getClusterYIndex(trackClusters[1]));
+
+ // Populate the momentum plots.
+ møctmInvariantMass.fill(getInvariantMass(pair));
+ møctmElectronMomentum.fill(pair[0].getMomentum().magnitude());
+ møctmElectronMomentum.fill(pair[1].getMomentum().magnitude());
+ møctmMomentumSum1D.fill(VecOp.add(pair[0].getMomentum(), pair[1].getMomentum()).magnitude());
+ møctmMomentumSum2D.fill(pair[0].getMomentum().magnitude(), pair[1].getMomentum().magnitude());
+ }
+
+ // Produce trident cluster-track matched plots.
+ trctmInstancesInEvent.fill(tridents.size());
+ for(ReconstructedParticle[] pair : tridents) {
+ // Get the electron and positron tracks.
+ ReconstructedParticle electronTrack = pair[pair[0].getCharge() < 0 ? 0 : 1];
+ ReconstructedParticle positronTrack = pair[pair[0].getCharge() > 0 ? 0 : 1];
+
+ // Get the track clusters.
+ Cluster electronCluster = electronTrack.getClusters().get(0);
+ Cluster positronCluster = positronTrack.getClusters().get(0);
+ Cluster[] trackClusters = { pair[0].getClusters().get(0), pair[1].getClusters().get(0) };
+
+ // Populate the cluster plots.
+ trctmElectronEnergy.fill(electronCluster.getEnergy());
+ trctmPositronEnergy.fill(positronCluster.getEnergy());
+ trctmEnergySum1D.fill(TriggerModule.getValueEnergySum(trackClusters));
+ trctmEnergySum2D.fill(pair[0].getEnergy(), pair[1].getEnergy());
+ trctmTimeCoincidence.fill(TriggerModule.getClusterTime(trackClusters[0]) - TriggerModule.getClusterTime(trackClusters[1]));
+ trctmClusterPosition.fill(TriggerModule.getClusterXIndex(trackClusters[0]), TriggerModule.getClusterYIndex(trackClusters[0]));
+ trctmClusterPosition.fill(TriggerModule.getClusterXIndex(trackClusters[1]), TriggerModule.getClusterYIndex(trackClusters[1]));
+
+ // Populate the momentum plots.
+ trctmInvariantMass.fill(getInvariantMass(pair));
+ trctmElectronMomentum.fill(electronTrack.getMomentum().magnitude());
+ trctmPositronMomentum.fill(positronTrack.getMomentum().magnitude());
+ trctmMomentumSum1D.fill(VecOp.add(pair[0].getMomentum(), pair[1].getMomentum()).magnitude());
+ trctmMomentumSum2D.fill(pair[0].getMomentum().magnitude(), pair[1].getMomentum().magnitude());
+ }
+ }
+
+ /**
+ * Produces pairs of tracks. The track pairs are required to be
+ * matched to a cluster and the associated clusters must form a
+ * top/bottom pair. If more than one track points to the same
+ * cluster, only the first track is retained.
+ * @param trackList - A list of all tracks.
+ * @return Returns a list of track pairs meeting the aforementioned
+ * conditions.
+ */
+ private static final List<ReconstructedParticle[]> getTopBottomTracksCTMatched(List<ReconstructedParticle> trackList) {
+ // Track clusters that have already been seen to prevent clusters
+ // that have duplicate tracks from reappearing.
+ Set<Cluster> clusterSet = new HashSet<Cluster>();
+
+ // Separate the tracks into top and bottom tracks based on
+ // the track cluster. Filter out tracks with no clusters.
+ List<ReconstructedParticle> topTracks = new ArrayList<ReconstructedParticle>();
+ List<ReconstructedParticle> botTracks = new ArrayList<ReconstructedParticle>();
+ trackLoop:
+ for(ReconstructedParticle track : trackList) {
+ // Check if the track has a cluster. If not, skip it.
+ if(track.getClusters().isEmpty()) {
+ continue trackLoop;
+ }
+
+ // If the track doesn't have actual tracks, skip it.
+ if(track.getTracks().isEmpty()) {
+ continue trackLoop;
+ }
+
+ // Check if the track cluster has already seen.
+ Cluster trackCluster = track.getClusters().get(0);
+ if(clusterSet.contains(trackCluster)) {
+ continue trackLoop;
+ }
+
+ // If the track has a unique cluster, add it to the proper
+ // list based on the cluster y-index.
+ clusterSet.add(trackCluster);
+ if(TriggerModule.getClusterYIndex(trackCluster) > 0) {
+ topTracks.add(track);
+ } else {
+ botTracks.add(track);
+ }
+ }
+
+ // Form all top/bottom pairs with the unique tracks.
+ List<ReconstructedParticle[]> pairList = new ArrayList<ReconstructedParticle[]>();
+ for(ReconstructedParticle topTrack : topTracks) {
+ for(ReconstructedParticle botTrack : botTracks) {
+ pairList.add(new ReconstructedParticle[] { topTrack, botTrack });
+ }
+ }
+
+ // Return the result.
+ return pairList;
+ }
+
+ /**
+ * Gets a list track pairs that meet the trident condition defined
+ * using tracks with matched calorimeter clusters. A pair meets the
+ * cluster/track matched trident condition is it meets the following:
+ * <ul><li>Both tracks have matched clusters.</li>
+ * <li>Has one positive track.</li>
+ * <li>Has one negative track.</li>
+ * <li>Clusters have a time coincidence of 2.5 ns or less.</li>
+ * <li>The electron momentum is below 900 MeV.</li></ul>
+ * @param pairList - A <code>List</code> collection of parameterized
+ * type <code>ReconstructedParticle[]</code> containing all valid
+ * top/bottom pairs of tracks with matched clusters. These will be
+ * tested to see if they meet the process criteria.
+ * @return Returns a list containing pairs of tracks that meet the
+ * trident condition.
+ */
+ private final List<ReconstructedParticle[]> getTridentTracksCTMatched(List<ReconstructedParticle[]> pairList) {
+ // Store the set of track pairs that meet the trident condition.
+ List<ReconstructedParticle[]> tridentTracks = new ArrayList<ReconstructedParticle[]>();
+
+ // Loop over the filtered pair list and apply the trident
+ // condition test.
+ tridentLoop:
+ for(ReconstructedParticle[] pair : pairList) {
+ // There must be one positive and one negative track.
+ ReconstructedParticle electron = null;
+ ReconstructedParticle positron = null;
+ if(pair[0].getCharge() > 0) { positron = pair[0]; }
+ else if(pair[1].getCharge() > 0) { positron = pair[1]; }
+ if(pair[0].getCharge() < 0) { electron = pair[0]; }
+ else if(pair[1].getCharge() < 0) { electron = pair[1]; }
+ if(electron == null || positron == null) {
+ continue tridentLoop;
+ }
+
+ // Make sure that the clusters are not the same. This should
+ // not actually ever be possible...
+ if(pair[0].getClusters().get(0) == pair[1].getClusters().get(0)) {
+ continue tridentLoop;
+ }
+
+ // The clusters must within a limited time window.
+ Cluster[] trackClusters = { pair[0].getClusters().get(0), pair[1].getClusters().get(0) };
+ if(TriggerModule.getValueTimeCoincidence(trackClusters) > timeCoincidence) {
+ continue tridentLoop;
+ }
+
+ // Require that the electron in the pair have an energy
+ // below 900 MeV to exclude elastic electrons.
+ if(electron.getMomentum().magnitude() >= 0.900) {
+ continue tridentLoop;
+ }
+
+ // If all the above conditions are met, the pair is to be
+ // considered a trident pair. Add it to the list.
+ tridentTracks.add(pair);
+ }
+
+ // Return the list of pairs that passed the condition.
+ return tridentTracks;
+ }
+
+ /**
+ * Gets a list track pairs that meet the Møller condition defined
+ * using tracks with matched calorimeter clusters. A pair meets the
+ * cluster/track matched Møller condition is it meets the following:
+ * <ul><li>Both tracks have matched clusters.</li>
+ * <li>Both tracks are negative.</li>
+ * <li>Clusters have a time coincidence of 2.5 ns or less.</li>
+ * <li>The electron momenta are below 900 MeV.</li>
+ * <li>The momentum sum of the tracks is in the range <code>800 MeV
+ * ⤠p1 + p2 ⤠1500 MeV</li></ul>
+ * @param pairList - A <code>List</code> collection of parameterized
+ * type <code>ReconstructedParticle[]</code> containing all valid
+ * top/bottom pairs of tracks with matched clusters. These will be
+ * tested to see if they meet the process criteria.
+ * @return Returns a list containing pairs of tracks that meet the
+ * Møller condition.
+ */
+ private final List<ReconstructedParticle[]> getMøllerTracksCTMatched(List<ReconstructedParticle[]> pairList) {
+ // Store the set of track pairs that meet the Møller condition.
+ List<ReconstructedParticle[]> møllerTracks = new ArrayList<ReconstructedParticle[]>();
+
+ // Loop over the filtered pair list and apply the Møller
+ // condition test.
+ møllerLoop:
+ for(ReconstructedParticle[] pair : pairList) {
+ // Both tracks must be negatively charged.
+ if(pair[0].getCharge() > 0 || pair[1].getCharge() > 0) {
+ continue møllerLoop;
+ }
+
+ // The clusters must within a limited time window.
+ Cluster[] trackClusters = { pair[0].getClusters().get(0), pair[1].getClusters().get(0) };
+ if(TriggerModule.getValueTimeCoincidence(trackClusters) > timeCoincidence) {
+ continue møllerLoop;
+ }
+
+ // Require that the electrons in the pair have energies
+ // below 900 MeV to exclude elastic electrons.
+ if(pair[0].getMomentum().magnitude() > 0.900 || pair[1].getMomentum().magnitude() > 0.900) {
+ continue møllerLoop;
+ }
+
+ // Require that the energy of the pair be within a range
+ // that is sufficiently "Møller-like."
+ double momentumSum = VecOp.add(pair[0].getMomentum(), pair[1].getMomentum()).magnitude();
+ if(momentumSum < 0.800 || momentumSum > 1.500) {
+ continue møllerLoop;
+ }
+
+ // If all the above conditions are met, the pair is to be
+ // considered a trident pair. Add it to the list.
+ møllerTracks.add(pair);
+ }
+
+ // Return the list of pairs that passed the condition.
+ return møllerTracks;
+ }
+
+ /**
+ * Calculates the approximate invariant mass for a pair of tracks
+ * from their momentum. This assumes that the particles are either
+ * electrons or positrons, and thusly have a sufficiently small
+ * mass term that it can be safely excluded.
+ * @param pair - The track pair for which to calculate the invariant
+ * mass.
+ * @return Returns the approximate invariant mass in units of GeV.
+ */
+ private static final double getInvariantMass(ReconstructedParticle[] pair) {
+ // Get the momentum squared.
+ double p2 = Math.pow(pair[0].getMomentum().magnitude() + pair[1].getMomentum().magnitude(), 2);
+
+ // Get the remaining terms.
+ double xPro = pair[0].getMomentum().x() + pair[1].getMomentum().x();
+ double yPro = pair[0].getMomentum().y() + pair[1].getMomentum().y();
+ double zPro = pair[0].getMomentum().z() + pair[1].getMomentum().z();
+
+ // Calculate the invariant mass.
+ return Math.sqrt(p2 - Math.pow(xPro, 2) - Math.pow(yPro, 2) - Math.pow(zPro, 2));
+ }
+}
+// [log in to unmask]
Modified: java/trunk/users/src/main/java/org/hps/users/kmccarty/plots/formatter/TriggerPlotsFormat.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/kmccarty/plots/formatter/TriggerPlotsFormat.java (original)
+++ java/trunk/users/src/main/java/org/hps/users/kmccarty/plots/formatter/TriggerPlotsFormat.java Tue Oct 20 18:57:00 2015
@@ -57,22 +57,23 @@
*/
public static void main(String[] args) throws IOException {
// Define the root directory for the plots.
- String rootDir = "D:\\cygwin64\\home\\Kyle\\aprime-plots\\";
+ String rootDir = "D:\\cygwin64\\home\\Kyle\\beam-plots\\base\\";
+ //String rootDir = "D:\\cygwin64\\home\\Kyle\\aprime-plots\\base\\readout-plots\\";
// Define the new name of the file containing the trigger plots.
String[] plotFile = {
- //rootDir + "background-ana_triggerPlots.aida"
- rootDir + "15-MeV\\compiled-plots.aida",
- rootDir + "20-MeV\\compiled-plots.aida",
- rootDir + "30-MeV\\compiled-plots.aida",
- rootDir + "40-MeV\\compiled-plots.aida",
- rootDir + "50-MeV\\compiled-plots.aida"
+ rootDir + "compiled-plots.aida"
+ //rootDir + "15-MeV\\compiled-plots.aida",
+ //rootDir + "20-MeV\\compiled-plots.aida",
+ //rootDir + "30-MeV\\compiled-plots.aida",
+ //rootDir + "40-MeV\\compiled-plots.aida",
+ //rootDir + "50-MeV\\compiled-plots.aida"
};
// Define the names of each plot. This will be used for the
// legend in the case of multiple plots.
String[] treeName = {
- //"Background",
+ "Background",
"15 MeV A'",
"20 MeV A'",
"30 MeV A'",
@@ -82,13 +83,13 @@
// Define the color style for the plots.
ColorStyle[] dataColorStyle = {
- //ColorStyle.GREY,
+ ColorStyle.GREY,
ColorStyle.MS_GREEN,
ColorStyle.MS_BLUE,
ColorStyle.MS_ORANGE,
ColorStyle.MS_RED,
+ ColorStyle.TEAL,
ColorStyle.CRIMSON,
- ColorStyle.TEAL,
ColorStyle.FOREST
};
