Merge branch 'trunk' into iss1074-geometry

exampleConfigs/README.md

@@ -0,0 +1,92 @@
+# Example Configuration Files
+
+## Test-Beam reconstruction
+
+### Before you start:
+- Verify that you are working with the latest version of the `Hcal/` submodule.
+```
+cd ldmx # the parent directory of ldmx-sw
+cd ldmx-sw/Hcal
+git switch trunk
+# rebuild ldmx-sw
+cd ../build/
+ldmx make install -j 4
+```
+- Verify that you have the [conditions-data](https://github.com/LDMX-Software/conditions-data) repository installed 
+```
+cd ldmx # the parent directory of ldmx-sw
+git clone [email protected]:ldmx-software/conditions-data
+cd ldmx-sw
+git checkout trunk
+git pull
+git submodule update
+# rebuild ldmx-sw
+cd ../build/
+ldmx make install -j 4
+```
+
+### For simulation:
+- Use simulated events with prototype geometry ```detector='ldmx-hcal-prototype-v2.0'```.
+- Run reconstruction on prototype simulation:
+```
+ldmx fire tb_sim.py hcal_XXX_simevents.root
+```
+
+### For data:
+
+- Run reconstruction on prototype simulation:
+```
+ldmx fire tb_reco.py hcal_run_XXX_decoded.root
+```
+
+In this configuration file `HcalSingleEndRecProducer` is run by default. A `DoubleEndRecProducer` will be available too.
+
+#### Decoded inputs
+A set of decoded April-2022 TB data is available on the SLAC cluster:
+```
+/nfs/slac/g/ldmx/data/hcal-tb
+```
+The directory holds the test beam data from the HCal prototype subsystem.
+
+Since event alignment is a tricky business (and sometimes is broken),
+the data has two copies: one with event alignment and one without.
+Please look at the README in that directory for more details.
+
+The following is the table of contents of this folder:
+```
+./
+|-- unaligned/ : unification of the two halves of the HCal was not attempted
+    |-- reformat/ : unpacking of raw data into event objects, no decoding
+    |-- decoded/  : products decoded into HgcrocDigiCollection objects
+    |-- decoded-ntuples/ : HgcrocDigiCollection objects ntuplized for easier analysis
+    |-- fail-decode.list : list of files that failed the decoding step
+    |-- fail-reformat.list : list of files that failed the reformat/unpack step
+|-- aligned/ : two polarfires of the HCal were aligned based on timestamps
+    |-- reformat/ : simple unpacking and alignemnt, no decoding
+    |-- decoded/  : decoded into HgcrocDigiCollection objects
+    |-- decoded-ntuples/ : HgcrocDigiCollection objects ntuplized for easier analysis
+|-- reformat_cfg.py : configuration script for `ldmx reformat ` used for unpacking and alignment
+|-- decode_cfg.py   : configuration script for `ldmx fire ` used for decoding and ntuplization
+````
+
+#### Re-doing data reformatting and decoding
+Event building was not done online, so it needs to be done on raw TB data files (if not using the existing set of decoded test-beam data). 
+
+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/
+```
+
+- Reformat: 
+  - [Build and install reformat](https://github.com/LDMX-Software/ldmx-tb-online/blob/main/reformat/README.md#Building)
+  - An [example config: hcal_cfg.py](https://github.com/LDMX-Software/ldmx-tb-online/blob/main/reformat/TestBeam/hcal_cfg.py) exists. A common run example would be:
+  ```
+  ldmx reformat TestBeam/hcal_cfg.py \
+  --pf0 path/to/fpga_0_run_XXX.raw \
+  --pf1 path/to/fpga_1_run_XXX.raw \
+  --output_filename hcal_run_XXX_reformat.root
+  ```
+- Decode:
+  - One can use the decoding producer in ldmx-sw. 
+  - A common example is available in [decode_cfg.py](https://github.com/LDMX-Software/Hcal/blob/trunk/exampleConfigs/decode_cfg.py)
+  `ldmx-sw fire decode_cfg.py hcal_run_XXX_reformat.root`. By default this decoding does not pedestal subtraction. 

exampleConfigs/decode_cfg.py

@@ -0,0 +1,89 @@
+"""Decoding configuration for raw testbeam data                                                                                                                                                                                                                                                                                                                      
+                                                                                                                                                                                                                                                                                                                                                                     
+Decoding **DOES NOT** attempt                                                                                                                                                                                                                                                                                                       
+to align the two halves of the HCal. We assume a local home path of ldmx-sw installation e.g.                                                                                                                                                                                                                                                                                                                  
+                                                                                                                                                                                                                                                                                                                                                                     
+  LDMX_BASE=/local/cms/user/eichl008/ldmx/                                                                                                                                                                                                                                                                                                                           
+                                                                                                                                                                                                                                                                                                                                                                     
+so that we can construct the output path correctly.                                                                                                                                                                                                                                                                                                                  
+                                                                                                                                                                                                                                                                                                                                                                     
+The run number is deduced from the file name.
+
+@author Tom Eichlersmith, University of Minnesota
+"""
+
+import argparse, sys
+
+parser = argparse.ArgumentParser(f'ldmx fire {sys.argv[0]}',
+    description=__doc__)
+
+parser.add_argument('input_file')
+parser.add_argument('--pause',action='store_true')
+parser.add_argument('--max_events',type=int)
+parser.add_argument('--pedestals',default=None,type=str)
+
+arg = parser.parse_args()
+
+from LDMX.Framework import ldmxcfg
+
+p = ldmxcfg.Process('decode')
+if arg.max_events is not None :
+    p.maxEvents = arg.max_events
+p.termLogLevel = 0
+p.logFrequency = 1000
+
+import LDMX.Hcal.hgcrocFormat as hcal_format
+import LDMX.Hcal.HcalGeometry
+import LDMX.Hcal.hcal_hardcoded_conditions
+from LDMX.DQM import dqm
+import os
+from LDMX.Hcal.DetectorMap import HcalDetectorMap
+detmap = HcalDetectorMap(f'{os.environ["LDMX_BASE"]}/ldmx-sw/Hcal/data/testbeam_connections.csv')
+detmap.want_d2e = True # helps quicken the det -> elec translation                                                                                                                                                                                                                                                                                                   
+
+# extract and deduce parameters from input file name                                                                                                                                                                                                                                                                                                                 
+params = os.path.basename(arg.input_file).replace('.root','').split('_')
+run = params[params.index("run")+1]
+day = params[-2]
+time = params[-1]
+if 'fpga' in params :
+    alignment = 'unaligned'
+    pf = params[params.index("fpga")+1]
+    provided = f'fpga_{pf}'
+    input_names = [ f'Polarfire{pf}Raw' ]
+    out_name = input_names[0]+'Digis'
+elif 'hcal' in params :
+    alignment = 'aligned'
+    provided = 'hcal'
+    input_names = [ f'Polarfire{pf}Raw' for pf in [0,1] ]
+    out_name = 'HcalRawDigis'
+else :
+    raise KeyError(f'Unable to deduce alignment from {fp}, need either "fpga" or "hcal" in base file name.')
+
+dir_name  = f'{os.environ["LDMX_BASE"]}/testbeam/{alignment}/v2-decoded'
+os.makedirs(dir_name, exist_ok=True)
+os.makedirs(dir_name+'-ntuple', exist_ok=True)
+
+file_name = f'decoded_{provided}_run_{run}_{day}_{time}'
+
+p.inputFiles = [arg.input_file]
+p.outputFiles = [f'{dir_name}/{file_name}.root']
+p.histogramFile = f'{dir_name}-ntuple/ntuple_{file_name}.root'
+
+# sequence                                                                                                                                                                                                                                                                                                                                                           
+#   1. decode event packet into digi collection                                                                                                                                                                                                                                                                                                                      
+#   2. ntuplize digi collection                                                                                                                                                                                                                                                                                                                                      
+p.sequence = [
+        hcal_format.HcalRawDecoder(
+            input_names = input_names,
+            output_name = out_name
+            ),
+        dqm.NtuplizeHgcrocDigiCollection(
+            input_name = out_name,
+            pedestal_table = arg.pedestals
+            )
+        ]
+
+if arg.pause :
+    p.pause()
+

exampleConfigs/tb_reco.py

@@ -34,11 +34,11 @@
 # sequence
 tbl = f'{os.environ["LDMX_BASE"]}/ldmx-sw/Hcal/data/testbeam_connections.csv'
 p.sequence = [
-    hcal_format.HcalRawDecoder(
-        input_names = ["Polarfire0Raw","Polarfire1Raw"],
-        connections_table = tbl,
-        output_name = 'HcalRawDigis'
-    ),
+    # hcal_format.HcalRawDecoder(
+    #     input_names = ["Polarfire0Raw","Polarfire1Raw"],
+    #     connections_table = tbl,
+    #     output_name = 'HcalRawDigis'
+    # ),
     hcal_digi.HcalSingleEndRecProducer(
       coll_name = 'HcalRawDigis',
       pass_name = ''

exampleConfigs/tb_sim.py

@@ -0,0 +1,63 @@
+import argparse,sys
+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)
+arg = parser.parse_args()
+
+p = ldmxcfg.Process('sim')
+p.maxEvents = arg.nevents
+p.termLogLevel = 0
+p.logFrequency = 10
+
+detector = 'ldmx-hcal-prototype-v2.0' # TODO: CHANGE TO FEFIX version
+
+p.outputFiles = [
+        arg.particle
+        +"Sim_%.2fGeV_"%arg.energy
+        + str(p.maxEvents) 
+        + "_%s.root"%detector
+        ]
+
+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., 0 ] # 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_hardcoded_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',
+        ),
+    ]
+)
+
+# View configuration before actually running
+p.pause()

python/hcal_hardcoded_conditions.py

@@ -24,8 +24,15 @@
 adc_gain = SimpleCSVDoubleTableProvider("hcal_adc_gain",["gain"])
 adc_gain.validForAllRows([1.2]) # 4 ADCs per PE - maxADCRange/readoutPadCapacitance/1024
 
-tot_calib = SimpleCSVDoubleTableProvider("hcal_tot_calibration",["pedestal"])
-tot_calib.validForAllRows([1.]) # dummy value since TOT is not implemented
+tot_calib = SimpleCSVDoubleTableProvider("hcal_tot_calibration",
+                                         ["m_adc_i","cut_point_tot","high_slope","high_offset",
+                                          "low_slope","low_power","low_offset","tot_not",
+                                          "channel","flagged"])
+tot_calib.validForAllRows([
+    1.,1.,1.,1.,
+    1.,1.,1.,1.,
+    1.,1.
+]) # dummy value since TOT is not implemented
 
 # wrap our tables in the parent object that is used by the processors
 from .conditions import HcalReconConditionsProvider

src/Hcal/HcalDoubleEndRecProducer.cxx

@@ -91,6 +91,7 @@ void HcalDoubleEndRecProducer::produce(framework::Event& event) {
       if(!digi_id.isNegativeEnd() && indices.first==-1) {
 	indices.first = iHit;
       }
+      iHit++;
     }
     indicesByID[id] = indices;
   }
@@ -103,7 +104,7 @@ void HcalDoubleEndRecProducer::produce(framework::Event& event) {
     auto position = hcalGeometry.getStripCenterPosition(id);
 
     // skip non-double-ended layers
-    if (id.layer() != ldmx::HcalID::HcalSection::BACK)
+    if (id.section() != ldmx::HcalID::HcalSection::BACK)
       continue;
 
     // get two hits to reconstruct

src/Hcal/HcalRecProducer.cxx

@@ -204,8 +204,8 @@ void HcalRecProducer::produce(framework::Event& event) {
         }
 
         // set voltage
-        voltage_posend = amplT_posend * the_conditions.adcGain(id_posend);
-        voltage_negend = amplT_negend * the_conditions.adcGain(id_negend);
+	voltage_posend = amplT_posend * the_conditions.adcGain(id_posend, 0);
+        voltage_negend = amplT_negend * the_conditions.adcGain(id_negend, 0);
       }
 
       // get TOA
@@ -351,6 +351,9 @@ void HcalRecProducer::produce(framework::Event& event) {
     recHit.setXPos(position.X());
     recHit.setYPos(position.Y());
     recHit.setZPos(position.Z());
+    recHit.setSection(id.section());
+    recHit.setStrip(id.strip());
+    recHit.setLayer(id.layer());
     recHit.setPE(PEs);
     recHit.setMinPE(minPEs);
     recHit.setAmplitude((amplT / voltage_per_mip_) * mip_energy_);

src/Hcal/HcalSingleEndRecProducer.cxx

@@ -40,7 +40,7 @@ class HcalSingleEndRecProducer : public framework::Producer {
    * with C++17 structured bindings, this tuple return can be bound to separate
    * variables:
    * ```cpp
-   * auto [ toa, sum_adc, sum_tot ] = extract_measurments(digi,pedestal,isoi);
+   * auto [ toa, sum_adc, sum_tot ] = extract_measurements(digi,pedestal,isoi);
    * ```
    * giving us the dual benefit of separate variable names while only having to
    * loop over the samples within a single digi once
@@ -52,7 +52,7 @@ class HcalSingleEndRecProducer : public framework::Producer {
    * @param[in] pedestal pedestal for this channel
    * @return tuple of (toa [ns since SOI], sum_adc, sum_tot)
    */
-  std::tuple<double, double, int> extract_measurments(
+  std::tuple<double, double, int> extract_measurements(
       const ldmx::HgcrocDigiCollection::HgcrocDigi& digi, double pedestal);
 
  public:
@@ -64,7 +64,7 @@ class HcalSingleEndRecProducer : public framework::Producer {
 
 };  // HcalSingleEndRecProducer
 
-std::tuple<double, double, int> HcalSingleEndRecProducer::extract_measurments(
+std::tuple<double, double, int> HcalSingleEndRecProducer::extract_measurements(
     const ldmx::HgcrocDigiCollection::HgcrocDigi& digi, double pedestal) {
   // sum_adc = total of all but first in-time adc measurements
   double sum_adc{0};
@@ -133,7 +133,7 @@ void HcalSingleEndRecProducer::produce(framework::Event& event) {
     // amplitude/TOT reconstruction
     double num_mips_equivalent{0};
     auto [toa, sum_adc, sum_tot] =
-        extract_measurments(digi, conditions.adcPedestal(digi.id()));
+        extract_measurements(digi, conditions.adcPedestal(digi.id()));
     auto is_adc = conditions.is_adc(digi.id(), sum_tot);
     if (is_adc) {
       double adc_calib = sum_adc / conditions.adcGain(digi.id(), 0);