Fixing iss887

detector-data/detectors/HPS_Run2021Pass2FEE_1pt92/SamplingFractions/Ecal.properties

@@ -0,0 +1 @@
+samplingFraction: 1.0

detector-data/detectors/HPS_Run2021Pass2FEE_1pt92/SamplingFractions/Hodoscope.properties

@@ -0,0 +1 @@
+samplingFraction: 1.0

detector-data/detectors/HPS_Run2021Pass2FEE_1pt92/detector.properties

@@ -0,0 +1 @@
+name: HPS_Run2021Pass2FEE_1pt92

ecal-recon/src/main/java/org/hps/recon/ecal/IterateGainFactorDriver.java

@@ -52,16 +52,43 @@ public void setOutputCollectionName(final String outputCollectionName) {
     }
 
     private final String ecalReadoutName = "EcalHits";
-
-    /**
+    private int calibYear=0;
+    private ArrayList<Integer> badChannels = null;
+    /*
      * Basic no argument constructor.
+       A.C.: I add a year-specific set of crystals that are broken, with the corresponding energy to be set to zero for both data and MC.
+       For the moment, this is hard-coded
      */
     public IterateGainFactorDriver() {
+        badChannels=new ArrayList<Integer>();
     }
 
     public void setGainFile(String filename) {
         this.gainFileName = filename;
     }
+
+    public void setCalibYear(int year) {
+        this.calibYear=year;
+        System.out.println("IterateGainFactorDriver: setting calib year to "+this.calibYear);  
+
+    }
+
+
+    public void setBadChannels() {
+        if (this.calibYear==2021) {
+            System.out.println("Adding the 2021 bad channels");
+            badChannels.add(6);
+            badChannels.add(15);
+            badChannels.add(26);     
+            badChannels.add(153);
+            badChannels.add(192);
+            badChannels.add(198);
+            badChannels.add(275);
+            badChannels.add(334);
+            badChannels.add(419);       
+        }
+        System.out.println("IterateGainFactorDriver: bad channels are "+badChannels);    
+    }
 
     /**
      * Read in a text file that has multiplicative factors on the original gain
@@ -105,13 +132,15 @@ private void readGainFile() {
         } catch (IOException e) {
             e.printStackTrace();
         }
+        System.out.println("ECAL Gain Factors were read");
     }
 
     @Override
     public void detectorChanged(Detector detector) {
         // ECAL combined conditions object.
         ecalConditions = DatabaseConditionsManager.getInstance().getEcalConditions();
         readGainFile();
+        setBadChannels();
     }
 
     /**
@@ -127,9 +156,14 @@ public List<CalorimeterHit> iterateHits(final List<CalorimeterHit> hits) {
         for (final CalorimeterHit hit : hits) {
             double time = hit.getTime();
             long cellID = hit.getCellID();
-            double energy = hit.getCorrectedEnergy() * gainFileGains.get(findChannelId(cellID));
-            CalorimeterHit newHit = CalorimeterHitUtilities.create(energy, time, cellID);
-            newHits.add(newHit);
+            /*If the channels is not flagged as "bad", re-compute the energy and create a new CalorimterHit
+             * Otherwise ignore the hit
+             * */    
+            if (this.badChannels.contains(findChannelId(cellID))==false) {
+                double energy = hit.getCorrectedEnergy() * gainFileGains.get(findChannelId(cellID));
+                CalorimeterHit newHit = CalorimeterHitUtilities.create(energy, time, cellID);
+                newHits.add(newHit);
+            }
         }
         return newHits;
     }

ecal-recon/src/main/java/org/hps/recon/ecal/cluster/ClusterCorrectionUtilities.java

@@ -20,14 +20,18 @@ public final class ClusterCorrectionUtilities {
     static final double NOISE_B = 1.3725E-4;
     static final double NOISE_C = 3.01E-4;
 
-    static final int N_ITERATIONS_2019 = 5;
-
     static final Random random = new Random();
 
+    static final int N_ITERATIONS_2019 = 5;
     static final double deltaClusterEnergy_thr2019 = 0.5 / 100; // relative deltaE/E
     static final double deltaClusterX_thr2019 = 0.5; // abs dX in mm
     static final double deltaClusterY_thr2019 = 0.5; // abs dY in mm
 
+    static final int N_ITERATIONS_2021 = 5;
+    static final double deltaClusterEnergy_thr2021 = 0.5 / 100; // relative deltaE/E
+    static final double deltaClusterX_thr2021 = 0.5; // abs dX in mm
+    static final double deltaClusterY_thr2021 = 0.5; // abs dY in mm
+
     // Calculate the noise factor to smear the Ecal energy by
     public static double calcNoise(double energy) {
         return random.nextGaussian() * Math.sqrt(NOISE_A + NOISE_B * energy + NOISE_C * Math.pow(energy, 2));
@@ -98,7 +102,53 @@ public static void applyCorrections(double beamEnergy, HPSEcal3 ecal, List<Clust
                         }
 
                     }
-                } else {
+                } else if (beamEnergy > 3.0) {
+                    for (int it = 0; it < N_ITERATIONS_2021; it++) {
+                        if (it == 0)
+                            addNoise = true;
+                        else
+                            addNoise = false;
+
+                        // get the cluster energy and cluster position from last iteration
+                        double clusterPosition[] = baseCluster.getPosition();
+                        double clusterEnergy = baseCluster.getEnergy();
+
+                        // To correct the energy, need the correct position - from previous iteration -
+                        // and the raw energy
+                        baseCluster.setPosition(clusterPosition);
+                        baseCluster.setEnergy(clusterEnergy_NC);
+                        ClusterEnergyCorrection2021.setCorrectedEnergy(ecal, baseCluster, isMC, addNoise);
+                        double clusterEnergy_C = baseCluster.getEnergy();
+
+                        // To correct the position, need the correct energy - from previous iteration -
+                        // and the raw position
+                        baseCluster.setEnergy(clusterEnergy);
+                        baseCluster.setPosition(clusterPosition_NC);
+                        ClusterPositionCorrection2021.setCorrectedPosition(baseCluster);
+
+                        // now set back the energy after correction (the position is already set)
+                        baseCluster.setEnergy(clusterEnergy_C);
+
+                        // absolute POS, relative in energy
+                        double deltaClusterEnergy = clusterEnergy - clusterEnergy_C;
+                        double deltaClusterX = clusterPosition[0] - baseCluster.getPosition()[0];
+                        double deltaClusterY = clusterPosition[1] - baseCluster.getPosition()[1];
+
+                        // A.C. : add a condition on the relative variation of the energy before-after
+                        // iteration
+                        // add a condition on the absolute variation of the position before-after
+                        // iteration
+                        // if all are satisfied, break the loop
+                        if ((Math.abs(deltaClusterEnergy / clusterEnergy) < deltaClusterEnergy_thr2021)
+                                && (Math.abs(deltaClusterX) < deltaClusterX_thr2021)
+                                && (Math.abs(deltaClusterY) < deltaClusterY_thr2021)) {
+                            break;
+                        }
+
+                    }
+                }
+
+                else {
                     ClusterPositionCorrection.setCorrectedPosition(baseCluster);
                     ClusterEnergyCorrection.setCorrectedEnergy(ecal, baseCluster, isMC);
                 }
@@ -119,11 +169,12 @@ public static void applyCorrections(double beamEnergy, HPSEcal3 ecal, Cluster cl
 
             BaseCluster baseCluster = (BaseCluster) cluster;
             boolean addNoise = false;
-            // Apply PID based energy correction.
+            // Apply PID based energy correction. 2019 and 2021 are a special case
             if (beamEnergy > 4.0) {
                 // here we need to play a little bit.
-                // ClusterEnergyCorrection2019 requires the CORRECT position and the RAW energy
-                // ClusterPositionCorrection2019 requires the CORRECT energy and the RAW
+                // ClusterEnergyCorrection2019/2021 requires the CORRECT position and the RAW
+                // energy
+                // ClusterPositionCorrection2019/2021 requires the CORRECT energy and the RAW
                 // position.
                 // Idea: iterate from the non-corrected values
 
@@ -146,13 +197,16 @@ public static void applyCorrections(double beamEnergy, HPSEcal3 ecal, Cluster cl
                     // and the raw energy
                     baseCluster.setPosition(clusterPosition);
                     baseCluster.setEnergy(clusterEnergy_NC);
+
                     ClusterEnergyCorrection2019.setCorrectedEnergy(ecal, baseCluster, isMC, addNoise);
+
                     double clusterEnergy_C = baseCluster.getEnergy();
 
                     // To correct the position, need the correct energy - from previous iteration -
                     // and the raw position
                     baseCluster.setEnergy(clusterEnergy);
                     baseCluster.setPosition(clusterPosition_NC);
+
                     ClusterPositionCorrection2019.setCorrectedPosition(baseCluster);
 
                     // now set back the energy after correction (the position is already set)
@@ -175,8 +229,68 @@ public static void applyCorrections(double beamEnergy, HPSEcal3 ecal, Cluster cl
                     }
 
                 }
+            } else if (beamEnergy > 3.0) {
+                // here we need to play a little bit.
+                // ClusterEnergyCorrection2019/2021 requires the CORRECT position and the RAW
+                // energy
+                // ClusterPositionCorrection2019/2021 requires the CORRECT energy and the RAW
+                // position.
+                // Idea: iterate from the non-corrected values
 
-            } else {
+                baseCluster.setNeedsPropertyCalculation(false); // should have been set already before calling this -
+                                                                // just in case.
+                double clusterPosition_NC[] = baseCluster.getPosition();
+                double clusterEnergy_NC = baseCluster.getEnergy();
+
+                for (int it = 0; it < N_ITERATIONS_2021; it++) {
+                    if (it == 0)
+                        addNoise = true;
+                    else
+                        addNoise = false;
+
+                    // get the cluster energy and cluster position from last iteration
+                    double clusterPosition[] = baseCluster.getPosition();
+                    double clusterEnergy = baseCluster.getEnergy();
+
+                    // To correct the energy, need the correct position - from previous iteration -
+                    // and the raw energy
+                    baseCluster.setPosition(clusterPosition);
+                    baseCluster.setEnergy(clusterEnergy_NC);
+
+                    ClusterEnergyCorrection2021.setCorrectedEnergy(ecal, baseCluster, isMC, addNoise);
+
+                    double clusterEnergy_C = baseCluster.getEnergy();
+
+                    // To correct the position, need the correct energy - from previous iteration -
+                    // and the raw position
+                    baseCluster.setEnergy(clusterEnergy);
+                    baseCluster.setPosition(clusterPosition_NC);
+
+                    ClusterPositionCorrection2021.setCorrectedPosition(baseCluster);
+
+                    // now set back the energy after correction (the position is already set)
+                    baseCluster.setEnergy(clusterEnergy_C);
+
+                    // absolute POS, relative in energy
+                    double deltaClusterEnergy = clusterEnergy - clusterEnergy_C;
+                    double deltaClusterX = clusterPosition[0] - baseCluster.getPosition()[0];
+                    double deltaClusterY = clusterPosition[1] - baseCluster.getPosition()[1];
+                    // A.C. : add a condition on the relative variation of the energy before-after
+                    // iteration
+                    // add a condition on the absolute variation of the position before-after
+                    // iteration
+                    // if all are satisfied, break the loop
+                    if ((Math.abs(deltaClusterEnergy / clusterEnergy) < deltaClusterEnergy_thr2021)
+                            && (Math.abs(deltaClusterX) < deltaClusterX_thr2021)
+                            && (Math.abs(deltaClusterY) < deltaClusterY_thr2021)) {
+                        break;
+                    }
+
+                }
+
+            }
+
+            else {
                 // Apply PID based position correction, which should happen before final energy
                 // correction.
                 ClusterPositionCorrection.setCorrectedPosition(baseCluster);
@@ -204,7 +318,17 @@ public static void applyCorrections(double beamEnergy, HPSEcal3 ecal, Cluster cl
                 // is different from below, since 2019 corrections are based on the CORRECTED
                 // ENERGY
                 ClusterPositionCorrection2019.setCorrectedPosition(baseCluster);
-            } else {
+            } else if (beamEnergy > 3.0) {
+                // Apply energy correction - this depends on the non-corrected energy (from
+                // baseCluster) and ypos, from tracking.
+                ClusterEnergyCorrection2021.setCorrectedEnergy(ecal, baseCluster, ypos, isMC);
+                // Now the energy is correct, can use the 2021 correction. Note that the order
+                // is different from below, since 2021 corrections are based on the CORRECTED
+                // ENERGY
+                ClusterPositionCorrection2021.setCorrectedPosition(baseCluster);
+            }
+
+            else {
                 // Apply PID based position correction, which should happen before final energy
                 // correction.
                 ClusterPositionCorrection.setCorrectedPosition(baseCluster);