java/trunk/users/src/main/java/org/hps/users/luca
--- java/trunk/users/src/main/java/org/hps/users/luca/ClusterCount.java (rev 0)
+++ java/trunk/users/src/main/java/org/hps/users/luca/ClusterCount.java 2014-05-09 10:14:36 UTC (rev 568)
@@ -0,0 +1,527 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.hps.users.luca;
+
+/**/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.hps.users.luca;
+
+
+
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogram2D;
+import java.io.IOException;
+
+import java.io.PrintWriter;
+import java.util.ArrayList;
+import java.util.EnumSet;
+import java.util.LinkedList;
+import java.util.List;
+import java.util.Queue;
+import org.hps.readout.ecal.ClockSingleton;
+import org.hps.readout.ecal.TriggerDriver;
+
+import org.hps.recon.ecal.ECalUtils;
+import org.hps.recon.ecal.HPSEcalCluster;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.geometry.Detector;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ * Reads clusters and makes trigger decision using opposite quadrant criterion.
+ * Prints triggers to file if file path specified.
+ *
+ * @author Omar Moreno <[log in to unmask]>
+ * @author Sho Uemura <[log in to unmask]>
+ * @version $Id: FADCTriggerDriver.java,v 1.4 2013/09/02 21:56:56 phansson Exp $
+ */
+public class ClusterCount extends TriggerDriver {
+
+ int nTriggers;
+ int totalEvents;
+ protected double beamEnergy = -1; //by default, get beam energy from detector name
+ private int minHitCount = 1;
+ private boolean useDefaultCuts = true;
+ private double clusterEnergyHigh = 2.2 * ECalUtils.GeV;
+ private double clusterEnergyLow = .1 * ECalUtils.GeV;
+ private double energySumThreshold = 2.2 * ECalUtils.GeV;
+ private double energyDifferenceThreshold = 2.2 * ECalUtils.GeV;
+ private double maxCoplanarityAngle = 90; // degrees
+// private double energyDistanceDistance = 250; // mm
+// private double energyDistanceThreshold = 0.8 / 2.2;
+ private double energyDistanceDistance = 200; // mm
+ private double energyDistanceThreshold = 0.5; // unitless fraction
+ // maximum time difference between two clusters, in units of readout cycles (4 ns).
+ private int pairCoincidence = 2;
+ private double originX = 1393.0 * Math.tan(0.03052); //ECal midplane, defined by photon beam position (30.52 mrad) at ECal face (z=1393 mm)
+ int allPairs;
+ int oppositeQuadrantCount;
+ int clusterEnergyCount;
+ int energySumCount;
+ int energyDifferenceCount;
+ int energyDistanceCount;
+ int coplanarityCount;
+ AIDA aida = AIDA.defaultInstance();
+ IHistogram2D clusterHitCount2DAll, clusterEnergy2DAll, clusterSumDiff2DAll, energyDistance2DAll, clusterAngles2DAll, clusterCoplanarity2DAll;
+ IHistogram2D clusterHitCount2D, clusterEnergy2D, clusterSumDiff2D, energyDistance2D, clusterAngles2D, clusterCoplanarity2D;
+ IHistogram1D triggerBits1D, triggerTimes1D;
+ IHistogram2D clusterSeeds, trigClusterSeeds;
+ int truthPeriod = 250;
+ private boolean useQuadrants = false;
+ protected String clusterCollectionName = "EcalClusters";
+ // FIFO queues of lists of clusters in each ECal half.
+ // Each list corresponds to one readout cycle.
+ private Queue<List<HPSEcalCluster>> topClusterQueue = null;
+ private Queue<List<HPSEcalCluster>> botClusterQueue = null;
+ PrintWriter pairWriter;
+
+ private enum Flag {
+
+ CLUSTER_HITCOUNT(4), CLUSTER_ENERGY(3), ENERGY_SUM_DIFF(2), ENERGY_DISTANCE(1), COPLANARITY(0);
+ private final int index;
+
+ Flag(int i) {
+ index = i;
+ }
+
+ static int bitmask(EnumSet<Flag> flags) {
+ int mask = 0;
+ for (Flag flag : flags) {
+ mask |= 1 << flag.index;
+ }
+ return mask;
+ }
+ }
+
+ public void setClusterCollectionName(String clusterCollectionName) {
+ this.clusterCollectionName = clusterCollectionName;
+ }
+
+ public void setCutsFromBeamEnergy(double beamEnergy) {
+ if (beamEnergy == 1.1) {
+ System.out.println(this.getClass().getSimpleName() + ": Setting trigger for 1.1 GeV beam");
+ maxCoplanarityAngle = 90;
+ clusterEnergyHigh = .7 * ECalUtils.GeV;
+ clusterEnergyLow = .1 * ECalUtils.GeV;
+ energySumThreshold = 0.8 * ECalUtils.GeV;
+ energyDifferenceThreshold = beamEnergy;
+ } else if (beamEnergy == 2.2) {
+ System.out.println(this.getClass().getSimpleName() + ": Setting trigger for 2.2 GeV beam");
+ maxCoplanarityAngle = 35;
+ clusterEnergyHigh = 1.5 * ECalUtils.GeV;
+ clusterEnergyLow = .1 * ECalUtils.GeV;
+ energySumThreshold = 1.9 * ECalUtils.GeV;
+ energyDifferenceThreshold = beamEnergy;
+ } else if (beamEnergy == 6.6) {
+ System.out.println(this.getClass().getSimpleName() + ": Setting trigger for 6.6 GeV beam");
+ maxCoplanarityAngle = 60;
+ clusterEnergyHigh = 5.0 * ECalUtils.GeV;
+ clusterEnergyLow = .1 * ECalUtils.GeV;
+ energySumThreshold = 5.5 * ECalUtils.GeV;
+ energyDifferenceThreshold = beamEnergy;
+ }
+ }
+
+ protected double getBeamEnergyFromDetector(Detector detector) {
+ if (detector.getName().contains("1pt1")) {
+ return 1.1;
+ } else if (detector.getName().contains("2pt2")) {
+ return 2.2;
+ } else if (detector.getName().contains("6pt6")) {
+ return 6.6;
+ } else {
+ return -1.0;
+ }
+ }
+
+ public void setTruthPeriod(int truthPeriod) {
+ this.truthPeriod = truthPeriod;
+ }
+
+ public void setPairCoincidence(int pairCoincidence) {
+ this.pairCoincidence = pairCoincidence;
+ }
+
+ /**
+ * Set X coordinate used as the origin for cluster coplanarity and distance
+ * calculations. Defaults to the ECal midplane. Units of mm.
+ *
+ * @param originX
+ */
+ public void setOriginX(double originX) {
+ this.originX = originX;
+ }
+
+ /**
+ * Used for plot ranges and cuts that scale with energy. 1.1, 2.2 and 6.6
+ * are associated with default cuts. Units of GeV.
+ *
+ * @param beamEnergy
+ */
+ public void setBeamEnergy(double beamEnergy) {
+ this.beamEnergy = beamEnergy;
+ }
+
+ /**
+ * Use default cuts based on beam energy.
+ *
+ * @param useDefaultCuts
+ */
+ public void setUseDefaultCuts(boolean useDefaultCuts) {
+ this.useDefaultCuts = useDefaultCuts;
+ }
+
+ /**
+ * Minimum hit count for a cluster.
+ *
+ * @param minHitCount
+ */
+ public void setMinHitCount(int minHitCount) {
+ this.minHitCount = minHitCount;
+ }
+
+ /**
+ * Maximum energy for a cluster. Units of GeV.
+ *
+ * @param clusterEnergyHigh
+ */
+ public void setClusterEnergyHigh(double clusterEnergyHigh) {
+ this.clusterEnergyHigh = clusterEnergyHigh * ECalUtils.GeV;
+ }
+
+ /**
+ * Minimum energy for a cluster. Units of GeV.
+ *
+ * @param clusterEnergyLow
+ */
+ public void setClusterEnergyLow(double clusterEnergyLow) {
+ this.clusterEnergyLow = clusterEnergyLow * ECalUtils.GeV;
+ }
+
+ /**
+ * Maximum energy sum of the two clusters in a pair. Units of GeV.
+ *
+ * @param energySumThreshold
+ */
+ public void setEnergySumThreshold(double energySumThreshold) {
+ this.energySumThreshold = energySumThreshold * ECalUtils.GeV;
+ }
+
+ /**
+ * Maximum energy difference between the two clusters in a pair. Units of
+ * GeV.
+ *
+ * @param energyDifferenceThreshold
+ */
+ public void setEnergyDifferenceThreshold(double energyDifferenceThreshold) {
+ this.energyDifferenceThreshold = energyDifferenceThreshold * ECalUtils.GeV;
+ }
+
+ /**
+ * Maximum deviation from coplanarity for the two clusters in a pair. Units
+ * of degrees.
+ *
+ * @param maxCoplanarityAngle
+ */
+ public void setMaxCoplanarityAngle(double maxCoplanarityAngle) {
+ this.maxCoplanarityAngle = maxCoplanarityAngle;
+ }
+
+ /**
+ * Distance threshold for the energy-distance cut. Units of mm.
+ *
+ * @param energyDistanceDistance
+ */
+ public void setEnergyDistanceDistance(double energyDistanceDistance) {
+ this.energyDistanceDistance = energyDistanceDistance;
+ }
+
+ /**
+ * Energy threshold for the energy-distance cut. Units of the beam energy.
+ *
+ * @param energyDistanceThreshold
+ */
+ public void setEnergyDistanceThreshold(double energyDistanceThreshold) {
+ this.energyDistanceThreshold = energyDistanceThreshold;
+ }
+
+ public void setPairFile(String pairFile) {
+ try {
+ pairWriter = new PrintWriter(pairFile);
+ } catch (IOException e) {
+ }
+ }
+
+ @Override
+ public void detectorChanged(Detector detector) {
+ if (beamEnergy < 0) {
+ beamEnergy = this.getBeamEnergyFromDetector(detector);
+ }
+ if (useDefaultCuts) {
+ setCutsFromBeamEnergy(beamEnergy);
+ }
+
+ clusterHitCount2DAll = aida.histogram2D("All cluster pairs: hit count (less energetic vs. more energetic)", 9, 0.5, 9.5, 9, 0.5, 9.5);
+ clusterSumDiff2DAll = aida.histogram2D("All cluster pairs: energy difference vs. sum", 100, 0.0, 2 * beamEnergy, 100, 0.0, beamEnergy);
+ clusterEnergy2DAll = aida.histogram2D("All cluster pairs: energy (less energetic vs. more energetic)", 100, 0.0, 2 * beamEnergy, 100, 0.0, beamEnergy);
+ energyDistance2DAll = aida.histogram2D("All cluster pairs: distance vs. energy (less energetic cluster)", 100, 0.0, 0.5 * beamEnergy, 25, 0.0, 400.0);
+ clusterCoplanarity2DAll = aida.histogram2D("All cluster pairs: cluster angle uncoplanarity vs. less energetic cluster angle", 100, -180.0, 180.0, 100, -180.0, 180.0);
+ clusterAngles2DAll = aida.histogram2D("All cluster pairs: cluster angle (less energetic vs. more energetic)", 100, -180.0, 180.0, 100, -180.0, 180.0);
+
+ clusterHitCount2D = aida.histogram2D("Passed other cuts: hit count (less energetic vs. more energetic)", 9, 0.5, 9.5, 9, 0.5, 9.5);
+ clusterSumDiff2D = aida.histogram2D("Passed other cuts: energy difference vs. sum", 100, 0.0, 2 * beamEnergy, 100, 0.0, beamEnergy);
+ clusterEnergy2D = aida.histogram2D("Passed other cuts: energy (less energetic vs. more energetic)", 100, 0.0, 2 * beamEnergy, 100, 0.0, beamEnergy);
+ energyDistance2D = aida.histogram2D("Passed other cuts: distance vs. energy (less energetic cluster)", 100, 0.0, 0.5 * beamEnergy, 25, 0.0, 400.0);
+ clusterCoplanarity2D = aida.histogram2D("Passed other cuts: cluster angle uncoplanarity vs. less energetic cluster angle", 100, -180.0, 180.0, 100, -180.0, 180.0);
+ clusterAngles2D = aida.histogram2D("Passed other cuts: cluster angle (less energetic vs. more energetic)", 100, -180.0, 180.0, 100, -180.0, 180.0);
+
+ triggerBits1D = aida.histogram1D(detector.getDetectorName() + " : " + clusterCollectionName + " : trigger bits", 33, -1.5, 31.5);
+ triggerTimes1D = aida.histogram1D(detector.getDetectorName() + " : " + clusterCollectionName + " : trigger times", truthPeriod, -0.5, truthPeriod - 0.5);
+
+ clusterSeeds = aida.histogram2D(detector.getDetectorName() + " : " + clusterCollectionName + " : Cluster seeds", 46, -23, 23, 11, -5.5, 5.5);
+ trigClusterSeeds = aida.histogram2D(detector.getDetectorName() + " : " + clusterCollectionName + " : Cluster seeds, with trigger", 46, -23, 23, 11, -5.5, 5.5);
+ }
+
+ @Override
+ public void startOfData() {
+ //initialize queues and fill with empty lists
+ topClusterQueue = new LinkedList<List<HPSEcalCluster>>();
+ botClusterQueue = new LinkedList<List<HPSEcalCluster>>();
+ for (int i = 0; i < 2 * pairCoincidence + 1; i++) {
+ topClusterQueue.add(new ArrayList<HPSEcalCluster>());
+ }
+ for (int i = 0; i < pairCoincidence + 1; i++) {
+ botClusterQueue.add(new ArrayList<HPSEcalCluster>());
+ }
+ super.startOfData();
+ if (clusterCollectionName == null) {
+ throw new RuntimeException("The parameter clusterCollectionName was not set!");
+ }
+
+ allPairs = 0;
+ oppositeQuadrantCount = 0;
+ clusterEnergyCount = 0;
+ energySumCount = 0;
+ energyDifferenceCount = 0;
+ energyDistanceCount = 0;
+ coplanarityCount = 0;
+ }
+
+ @Override
+ public void process(EventHeader event) {
+ if (event.hasCollection(HPSEcalCluster.class, clusterCollectionName)) {
+ // this needs to run every readout cycle whether or not trigger is live
+ updateClusterQueues(event.get(HPSEcalCluster.class, clusterCollectionName));
+
+
+ }
+ super.process(event);
+ }
+
+
+
+
+ @Override
+ public void endOfData() {
+ if (outputStream != null) {
+ printCounts(outputStream);
+ }
+ printCounts(new PrintWriter(System.out));
+ if (pairWriter != null) {
+ pairWriter.close();
+ }
+ super.endOfData();
+ }
+
+ private void printCounts(PrintWriter writer) {
+ writer.printf("Number of pairs: %d\n", allPairs);
+ writer.printf("Number of cluster pairs after successive trigger conditions:\n");
+ if (useQuadrants) {
+ writer.printf("Opposite quadrants: %d\n", oppositeQuadrantCount);
+ }
+ writer.printf("Cluster energy: %d\n", clusterEnergyCount);
+ writer.printf("Energy sum: %d\n", energySumCount);
+ writer.printf("Energy difference: %d\n", energyDifferenceCount);
+ writer.printf("Energy-distance cut: %d\n", energyDistanceCount);
+ writer.printf("Coplanarity: %d\n", coplanarityCount);
+ writer.printf("Trigger count: %d\n", numTriggers);
+ writer.close();
+ }
+
+ protected void updateClusterQueues(List<HPSEcalCluster> ecalClusters) {
+ ArrayList<HPSEcalCluster> topClusterList = new ArrayList<HPSEcalCluster>();
+ ArrayList<HPSEcalCluster> botClusterList = new ArrayList<HPSEcalCluster>();
+ for (HPSEcalCluster ecalCluster : ecalClusters) {
+// System.out.format("add cluster\t%f\t%d\n", ecalCluster.getSeedHit().getTime(), ecalCluster.getSeedHit().getIdentifierFieldValue("iy"));
+ if (ecalCluster.getSeedHit().getIdentifierFieldValue("iy") > 0) {
+ topClusterList.add(ecalCluster);
+ } else {
+ botClusterList.add(ecalCluster);
+ }
+
+ int ix = ecalCluster.getSeedHit().getIdentifierFieldValue("ix");
+ int iy = ecalCluster.getSeedHit().getIdentifierFieldValue("iy");
+ clusterSeeds.fill(ix - 0.5 * Math.signum(ix), iy);
+ }
+
+ topClusterQueue.add(topClusterList);
+ botClusterQueue.add(botClusterList);
+ topClusterQueue.remove();
+ botClusterQueue.remove();
+ }
+
+ /**
+ * Get a list of all unique cluster pairs in the event
+ *
+ * @param ecalClusters : List of ECal clusters
+ * @return list of cluster pairs
+ */
+ protected List<HPSEcalCluster[]> getClusterPairsTopBot() {
+ // Make a list of cluster pairs
+ List<HPSEcalCluster[]> clusterPairs = new ArrayList<HPSEcalCluster[]>();
+
+ // Loop over all top-bottom pairs of clusters; higher-energy cluster goes first in the pair
+ // To apply pair coincidence time, use only bottom clusters from the
+ // readout cycle pairCoincidence readout cycles ago, and top clusters
+ // from all 2*pairCoincidence+1 previous readout cycles
+ for (HPSEcalCluster botCluster : botClusterQueue.element()) {
+ for (List<HPSEcalCluster> topClusters : topClusterQueue) {
+ for (HPSEcalCluster topCluster : topClusters) {
+// System.out.format("%f\t%f\n", topCluster.getSeedHit().getTime(), botCluster.getSeedHit().getTime());
+ if (topCluster.getEnergy() > botCluster.getEnergy()) {
+ HPSEcalCluster[] clusterPair = {topCluster, botCluster};
+ clusterPairs.add(clusterPair);
+ } else {
+ HPSEcalCluster[] clusterPair = {botCluster, topCluster};
+ clusterPairs.add(clusterPair);
+ }
+ }
+ }
+ }
+ return clusterPairs;
+ }
+
+ /**
+ * Checks if the ECal clusters making up a cluster pair lie in opposite
+ * quadrants
+ *
+ * @param clusterPair : pair of clusters
+ * @return true if opposite quadrants, false otherwise
+ */
+ protected boolean oppositeQuadrantsCut(HPSEcalCluster[] clusterPair) {
+ int quad1 = ECalUtils.getQuadrant(clusterPair[0]);
+ int quad2 = ECalUtils.getQuadrant(clusterPair[1]);
+
+ //if clusters are in the same quadrant, they're not opposite quadrants
+ if (quad1 == quad2) {
+ return false;
+ } //opposite pairs of quadrants are either both even (2 and 4) or both odd (1 and 3)
+ else {
+ return ((quad1 & 1) == (quad2 & 1));
+ }
+ }
+
+ /**
+ * Checks if the ECal clusters making up a cluster pair both have at least
+ * the minimum number of hits.
+ *
+ * @param clusterPair: pair of clusters
+ * @return true if pair passes cut, false if fail
+ */
+ protected boolean clusterHitCount(HPSEcalCluster[] clusterPair) {
+ return (clusterPair[0].getCalorimeterHits().size() >= minHitCount
+ && clusterPair[1].getCalorimeterHits().size() >= minHitCount);
+ }
+
+ /**
+ * Checks if the ECal clusters making up a cluster pair lie above the low
+ * energy threshold and below the high energy threshold
+ *
+ * @param clusterPair : pair of clusters
+ * @return true if a pair is found, false otherwise
+ */
+ protected boolean clusterECut(HPSEcalCluster[] clusterPair) {
+ return (clusterPair[0].getEnergy() < clusterEnergyHigh
+ && clusterPair[1].getEnergy() < clusterEnergyHigh
+ && clusterPair[0].getEnergy() > clusterEnergyLow
+ && clusterPair[1].getEnergy() > clusterEnergyLow);
+ }
+
+ /**
+ * Checks if the sum of the energies of ECal clusters making up a cluster
+ * pair is below an energy sum threshold
+ *
+ * @param clusterPair : pair of clusters
+ * @return true if a pair is found, false otherwise
+ */
+ protected boolean energySum(Cluster[] clusterPair) {
+ double clusterESum = clusterPair[0].getEnergy() + clusterPair[1].getEnergy();
+ return (clusterESum < energySumThreshold);
+ }
+
+ /**
+ * Checks if the energy difference between the ECal clusters making up a
+ * cluster pair is below an energy difference threshold
+ *
+ * @param clusterPair : pair of clusters
+ * @return true if pair is found, false otherwise
+ */
+ protected boolean energyDifference(HPSEcalCluster[] clusterPair) {
+ double clusterEDifference = clusterPair[0].getEnergy() - clusterPair[1].getEnergy();
+
+ return (clusterEDifference < energyDifferenceThreshold);
+ }
+
+ /**
+ * Require that the distance from the beam of the lowest energy cluster in a
+ * cluster pair satisfies the following E_low + d_b*.0032 GeV/mm < .8 GeV
+ *
+ * @param clusterPair : pair of clusters
+ * @return true if pair is found, false otherwise
+ */
+ protected boolean energyDistanceCut(HPSEcalCluster[] clusterPair) {
+ HPSEcalCluster lowEnergyCluster = clusterPair[1];
+
+ // Calculate its position
+ double lowEClusterDistance = getClusterDistance(clusterPair[1]);
+ // event passes cut if above the line with X- and Y-intercepts defined by energyDistanceDistance and beamEnergy*energyDistanceThreshold
+ double clusterDistvsE = lowEnergyCluster.getEnergy() + lowEClusterDistance * beamEnergy * energyDistanceThreshold / energyDistanceDistance;
+
+ return (clusterDistvsE > beamEnergy * energyDistanceThreshold);
+ }
+
+ /**
+ * Checks if a cluster pair is coplanar to the beam within a given angle
+ *
+ * @param clusterPair : pair of clusters
+ * @return true if pair is found, false otherwise
+ */
+ protected boolean coplanarityCut(HPSEcalCluster[] clusterPair) {
+ return (Math.abs(pairUncoplanarity(clusterPair)) < maxCoplanarityAngle);
+ }
+
+ protected double pairUncoplanarity(HPSEcalCluster[] clusterPair) { // Find the angle between clusters in the pair
+ double cluster1Angle = (getClusterAngle(clusterPair[0]) + 180.0) % 180.0;
+ double cluster2Angle = (getClusterAngle(clusterPair[1]) + 180.0) % 180.0;
+
+ return cluster2Angle - cluster1Angle;
+ }
+
+ protected double getClusterAngle(HPSEcalCluster cluster) { //returns angle in range of -180 to 180
+ double position[] = cluster.getSeedHit().getPosition();
+ return Math.toDegrees(Math.atan2(position[1], position[0] - originX));
+ }
+
+ protected double getClusterDistance(HPSEcalCluster cluster) {
+ return Math.hypot(cluster.getSeedHit().getPosition()[0] - originX, cluster.getSeedHit().getPosition()[1]);
+ }
+}
\ No newline at end of file