Merge pull request #64 from LDMX-Software/iss40-toa

First pass at double-ended reconstruction for test beam MC and data

exampleConfigs/README.md

@@ -44,7 +44,7 @@ In this configuration file `HcalSingleEndRecProducer` is run by default. A `Doub
 #### Decoded inputs
 A set of decoded April-2022 TB data is available on the SLAC cluster:
 ```
-/nfs/slac/g/ldmx/data/hcal-tb
+/sdf/group/ldmx/data/hcal-tb
 ```
 The directory holds the test beam data from the HCal prototype subsystem.
 
@@ -74,7 +74,7 @@ Event building was not done online, so it needs to be done on raw TB data files
 
 A copy of raw April-2022 TB data is available on the SLAC cluster:
 ```
-/nfs/slac/g/ldmx/CERN-TB-DATA/ldmx/testbeam/data/pf_external/
+/sdf/group/ldmx/CERN-TB-DATA/ldmx/testbeam/data/pf_external/
 ```
 
 - Reformat: 

exampleConfigs/tb_reco.py

@@ -40,7 +40,14 @@
     #     output_name = 'HcalRawDigis'
     # ),
     hcal_digi.HcalSingleEndRecProducer(
-      coll_name = 'HcalRawDigis',
-      pass_name = ''
+        coll_name = 'HcalRawDigis',
+        rec_coll_name = 'HcalSingleEndRecHits',
+        rec_pass_name = '',
+    ),
+    hcal_digi.HcalDoubleEndRecProducer(
+        pass_name = '', 
+        coll_name = 'HcalSingleEndRecHits', 
+        rec_coll_name = 'HcalDoubleEndRecHits',
+        rec_pass_name = ''
     )
 ]

exampleConfigs/tb_reco_pass1.py

@@ -0,0 +1,49 @@
+import argparse, sys, os, pathlib
+
+"""
+Takes raw data file after reformat (with run number)
+and runs single ended reconstruction
+"""
+
+parser = argparse.ArgumentParser(f'ldmx fire {sys.argv[0]}')
+parser.add_argument('input_file')
+parser.add_argument('--pause',action='store_true')
+grp = parser.add_mutually_exclusive_group()
+parser.add_argument('--output_dir',default='.',type=pathlib.Path)
+#parser.add_argument('--max_events',default=-1,type=int)
+arg = parser.parse_args()
+
+from LDMX.Framework import ldmxcfg
+p = ldmxcfg.Process('')
+#p.maxEvents = arg.max_events
+p.termLogLevel = 0
+p.logFrequency = 1000
+
+print(arg.output_dir)
+
+import LDMX.Hcal.HcalGeometry
+import LDMX.Hcal.hcal_testbeam0422_conditions
+import LDMX.Hcal.digi as hcal_digi
+import LDMX.Hcal.hgcrocFormat as hcal_format
+
+base_name = os.path.basename(arg.input_file).replace('.root','')
+dir_name  = os.path.dirname(arg.output_dir)
+
+p.inputFiles = [arg.input_file]
+p.outputFiles = [f'{dir_name}/reco_{base_name}_pass1.root']
+
+# sequence
+tbl = f'{os.environ["LDMX_BASE"]}/ldmx-sw/Hcal/data/testbeam_connections.csv'
+p.sequence = [
+    hcal_digi.HcalSingleEndRecProducer(
+        coll_name = 'HcalRawDigis',
+        rec_coll_name = 'HcalSingleEndRecHits',
+        rec_pass_name = '',
+    ),
+    hcal_digi.HcalDoubleEndRecProducer(
+        pass_name = '',
+        coll_name = 'HcalSingleEndRecHits',
+        rec_coll_name = 'HcalDoubleEndRecHits',
+        rec_pass_name = ''
+    )
+]

exampleConfigs/tb_sim.py

@@ -21,7 +21,7 @@
 p.outputFiles = [
         arg.particle
         +"Sim_%.2fGeV_"%arg.energy
-        + str(p.maxEvents) 
+        + str(p.maxEvents)
         + "_%s.root"%detector
         ]
 
@@ -35,7 +35,7 @@
 from LDMX.SimCore import generators as gen
 myGun = gen.gun('myGun')
 myGun.particle = arg.particle
-myGun.position = [ 0., 0., 0 ] # mm 
+myGun.position = [ 0., 0., -600. ] # mm
 myGun.direction = [ 0., 0., 1] # forward in z
 myGun.energy = arg.energy # GeV
 mySim.generators = [ myGun ]
@@ -44,7 +44,7 @@
 
 # import chip/geometry (hardcoded) conditions
 import LDMX.Hcal.HcalGeometry
-import LDMX.Hcal.hcal_hardcoded_conditions
+import LDMX.Hcal.hcal_testbeamsim_conditions
 import LDMX.Hcal.digi as hcal_digi
 
 # add them to the sequence
@@ -56,6 +56,12 @@
             pass_name = 'sim', coll_name = 'HcalDigis',
             rec_coll_name = 'HcalSingleEndRecHits',
         ),
+        hcal_digi.HcalDoubleEndRecProducer(
+            pass_name = '',
+            coll_name = 'HcalSingleEndRecHits',
+            rec_coll_name = 'HcalDoubleEndRecHits',
+            rec_pass_name = ''
+        )
     ]
 )
 

exampleConfigs/tb_sim_pass1.py

@@ -0,0 +1,68 @@
+import argparse, sys, os, pathlib
+from LDMX.Framework import ldmxcfg
+
+"""
+Simulation of particles through TB prototype
+"""
+
+parser = argparse.ArgumentParser(f'ldmx fire {sys.argv[0]}')
+parser.add_argument('--nevents',default=100,type=int)
+parser.add_argument('--particle',default='neutron') # other options, mu-,e-,pi-,proton
+parser.add_argument('--energy',default=2.0,type=float)
+parser.add_argument('--runnumber',default=1,type=int)
+parser.add_argument('--output_dir',default='.',type=pathlib.Path)
+arg = parser.parse_args()
+
+p = ldmxcfg.Process('sim')
+p.maxEvents = arg.nevents
+p.termLogLevel = 0
+p.logFrequency = 10
+p.run = arg.runnumber
+
+detector = 'ldmx-hcal-prototype-v2.0' # TODO: CHANGE TO FEFIX version
+
+base_name = os.path.basename(arg.particle+"Sim_%.2fGeV_"%arg.energy+str(p.maxEvents)+"_%s"%detector+"_pass1_%i"%arg.runnumber)
+dir_name  = os.path.dirname(arg.output_dir)
+
+p.outputFiles = [f'{dir_name}/{base_name}.root']
+
+from LDMX.SimCore import simulator
+import LDMX.Ecal.EcalGeometry # geometry required by sim
+
+mySim = simulator.simulator('mySim')
+mySim.setDetector(detector)
+
+# Get a pre-written generator
+from LDMX.SimCore import generators as gen
+myGun = gen.gun('myGun')
+myGun.particle = arg.particle
+myGun.position = [ 0., 0., -600. ] # mm
+myGun.direction = [ 0., 0., 1] # forward in z
+myGun.energy = arg.energy # GeV
+mySim.generators = [ myGun ]
+p.sequence.append( mySim )
+# mySim.verbosity = 1
+
+# import chip/geometry (hardcoded) conditions
+import LDMX.Hcal.HcalGeometry
+import LDMX.Hcal.hcal_testbeamsim_conditions
+import LDMX.Hcal.digi as hcal_digi
+
+# add them to the sequence
+p.sequence.extend(
+    [
+        hcal_digi.HcalDigiProducer(),
+        hcal_digi.HcalRecProducer(),
+        hcal_digi.HcalSingleEndRecProducer(
+            pass_name = 'sim', coll_name = 'HcalDigis',
+            rec_coll_name = 'HcalSingleEndRecHits',
+        ),
+        hcal_digi.HcalDoubleEndRecProducer(
+            pass_name = '',
+            coll_name = 'HcalSingleEndRecHits',
+            rec_coll_name = 'HcalDoubleEndRecHits',
+            rec_pass_name = ''
+        )
+    ]
+)
+\

include/Hcal/Event/HcalHit.h

@@ -87,6 +87,30 @@ class HcalHit : public ldmx::CalorimeterHit {
    */
   int getIsADC() const { return isADC_; }
 
+  /**
+   * Get the toa of the positive end.
+   * @return toaPos
+   */
+  int getToaPos() const { return toaPos_; }
+
+  /**
+   * Get the toa of the negative end.
+   * @return toaNeg
+   */
+  int getToaNeg() const { return toaNeg_; }
+
+  /**
+   * Get the amplitude of the positive end.
+   * @return amplitudePos
+   */
+  int getAmplitudePos() const { return amplitudePos_; }
+
+  /**
+   * Get the amplitude of the negative end.
+   * @return amplitudeNeg
+   */
+  int getAmplitudeNeg() const { return amplitudeNeg_; }
+
   /**
    * Set the number of photoelectrons estimated for this hit.
    * @param pe Number of photoelectrons, including noise which affects the
@@ -131,6 +155,45 @@ class HcalHit : public ldmx::CalorimeterHit {
    */
   void setIsADC(int isADC) { isADC_ = isADC; }
 
+  /**
+   * Set time difference (uncorrected)
+   * @param time
+   */
+  void setTimeDiff(double timeDiff) { timeDiff_ = timeDiff; }
+
+  /**
+   * Set toa of the positive end
+   * @param time
+   */
+  void setToaPos(double toaPos) { toaPos_ = toaPos; }
+
+  /**
+   * Set toa of the negative end
+   * @param time
+   */
+  void setToaNeg(double toaNeg) { toaNeg_ = toaNeg; }
+
+  /**
+   * Set amplitude of the positive end
+   * @param amplitude
+   */
+  void setAmplitudePos(double amplitudePos) { amplitudePos_ = amplitudePos; }
+
+  /**
+   * Set amplitude of the negative end
+   * @param amplitude
+   */
+  void setAmplitudeNeg(double amplitudeNeg) { amplitudeNeg_ = amplitudeNeg; }
+
+  /**
+   * Set original position
+   */
+  void setPositionUnchanged(double position, int isX) { position_ = position; isX_ = isX; }
+
+  double getPosition() const { return position_;}
+  int getIsX() const { return isX_;}
+  double getTimeDiff() const { return timeDiff_;}
+
  private:
   /** The number of PE estimated for this hit. */
   float pe_{0};
@@ -148,6 +211,15 @@ class HcalHit : public ldmx::CalorimeterHit {
   /// isADC
   int isADC_;
 
+  double timeDiff_;
+  double position_;
+  double isX_;
+
+  double toaPos_;
+  double toaNeg_;
+  double amplitudePos_;
+  double amplitudeNeg_;
+
   /**
    * The ROOT class definition.
    */

include/Hcal/HcalReconConditions.h

@@ -45,11 +45,14 @@ class HcalReconConditions : public framework::ConditionsObject {
    *
    * @param[in] adc_ped double table of ADC pedestals
    * @param[in] adc_gain double table of ADC gains
-   * @param[in] tot_ped double table of TOT calibrations
+   * @param[in] tot_calib double table of TOT calibrations
+   * @param[in] toa_calib double table of TOA calibrations 
    */
   HcalReconConditions(const conditions::DoubleTableCondition& adc_ped,
                       const conditions::DoubleTableCondition& adc_gain,
-                      const conditions::DoubleTableCondition& tot_calib);
+                      const conditions::DoubleTableCondition& tot_calib,
+		      const conditions::DoubleTableCondition& toa_calib
+		      );
 
   /**
    * get the ADC pedestal
@@ -112,13 +115,26 @@ class HcalReconConditions : public framework::ConditionsObject {
     return tot_calibs_.get(id.raw(), idx);
   }
 
+  /**
+   * get a TOA calibration value
+   *
+   * @param[in] id HCal Digi ID for specific chip
+   * @param[in] index of column in condition file
+   * @returns the TOA calibration for that i
+   */
+  double toaCalib(const ldmx::HcalDigiID& id, int idx = 0) const {
+    return toa_calibs_.get(id.raw(), idx);
+  }
+
  private:
   /// reference to the table of conditions storing the adc pedestals
   const conditions::DoubleTableCondition& adc_pedestals_;
   /// reference to the table of conditions storing the adc gains
   const conditions::DoubleTableCondition& adc_gains_;
   /// reference to the table of conditions storing the tot calibrations
   const conditions::DoubleTableCondition& tot_calibs_;
+  /// reference to the table of conditions storing the toa calibrations
+  const conditions::DoubleTableCondition& toa_calibs_;
 };  // HcalReconConditions
 
 }  // namespace hcal

python/conditions.py

@@ -21,19 +21,21 @@ class HcalReconConditionsProvider(ldmxcfg.ConditionsObjectProvider) :
         provider for the HCal ADC gains
     tot_calib : framework::ConditionsObjectProvider
         provider for the HCal TOT calibrations
+    toa_calib : framework::ConditionsObjectProvider
+        provider for the HCal TOA calibrations
 
     Examples
     --------
     The hcal_hardcoded_conditions.py file provides a working example where each condition
     wrapped here are constant for all runs and all channels.
     """
 
-    def __init__(self,adc_ped,adc_gain,tot_calib) :
+    def __init__(self,adc_ped,adc_gain,tot_calib,toa_calib) :
         super().__init__("HcalReconConditions","hcal::HcalReconConditionsProvider","Hcal")
 
         # our COP only needs the object names but providing the full parent COPs
         #   ensures that they exist
         self.adc_ped = adc_ped.objectName
         self.adc_gain = adc_gain.objectName
         self.tot_calib = tot_calib.objectName
-
+        self.toa_calib = toa_calib.objectName

python/hcal_hardcoded_conditions.py

@@ -34,9 +34,13 @@
     1.,1.
 ]) # dummy value since TOT is not implemented
 
+toa_calib =  SimpleCSVDoubleTableProvider("hcal_toa_calibration",
+                                          ["bx_shift","mean_shift"])
+toa_calib.validForAllRows([0., 0.]) # dummy values
+
 # wrap our tables in the parent object that is used by the processors
 from .conditions import HcalReconConditionsProvider
-HcalReconConditionsProvider(adc_pedestal, adc_gain, tot_calib)
+HcalReconConditionsProvider(adc_pedestal, adc_gain, tot_calib, toa_calib)
 
 HcalHgcrocConditionsHardcode=SimpleCSVDoubleTableProvider("HcalHgcrocConditions", [
             "PEDESTAL",

python/hcal_testbeam0422_conditions.py

@@ -30,5 +30,10 @@
          "low_slope","low_power","low_offset","tot_not","channel","flagged"])
 tot_calib.validForever(f'file://{os.environ["LDMX_BASE"]}/conditions-data/Hcal/testbeam04-2022/tot_calibration/calibrated_tot_calib_v0_1_0.csv')
 
+toa_calib = SimpleCSVDoubleTableProvider("hcal_toa_calibration",
+                                         ["bx_shift", "mean_shift"])
+toa_calib.conditions_baseURL = f'file://{os.environ["LDMX_BASE"]}/conditions-data/'
+toa_calib.entriesURL = '${LDMX_CONDITION_BASEURL}/Hcal/testbeam04-2022/toa_calibration/index_v1_0_0.csv'
+
 from .conditions import HcalReconConditionsProvider
-HcalReconConditionsProvider(adc_pedestal, adc_gain, tot_calib)
+HcalReconConditionsProvider(adc_pedestal, adc_gain, tot_calib, toa_calib)