java/trunk/users/src/main/java/org/hps/users/luca
--- java/trunk/users/src/main/java/org/hps/users/luca/ClusterCount.java 2014-05-13 20:09:17 UTC (rev 575)
+++ java/trunk/users/src/main/java/org/hps/users/luca/ClusterCount.java 2014-05-13 21:35:21 UTC (rev 576)
@@ -1,527 +0,0 @@
-/*
- * 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