Hcal tot digitization

Hcal/include/Hcal/HcalDigiProducer.h

@@ -88,6 +88,12 @@ class HcalDigiProducer : public framework::Producer {
   /// Helpful when comparing with test beam data
   bool zeroSuppression_{true};
 
+  /// If true, save digitization emulation in a file
+  bool digiEmulationDebug_{false};
+
+  /// Name of txt file to save digitization emulation
+  std::string digiEmulationFile_;
+
   /// Hgcroc Emulator to digitize analog voltage signals
   std::unique_ptr<ldmx::HgcrocEmulator> hgcroc_;
 

Hcal/python/digi.py

@@ -89,7 +89,7 @@ def __init__(self, instance_name = 'hcalDigis') :
 
         # avg parameters
         self.avgReadoutThreshold = 4. #ADCs - noise config only
-        self.avgGain = 1.2 #noise config only
+        self.avgGain = 5100./1023. #noise config only
         self.avgPedestal = 1. #noise config only
         # avg noise set to 0.02PE
         self.avgNoiseRMS = self.hgcroc.calculateVoltageHcal(0.02)/self.avgGain
@@ -98,6 +98,11 @@ def __init__(self, instance_name = 'hcalDigis') :
         # to e.g. simulate a pedestal measurement
         self.zeroSuppression = True
 
+        # if True, save pulse shapes to a file
+        self.digiEmulationDebug = False
+        # name of the file to save pulse shapes to
+        self.digiEmulationFile = ""
+
         # input and output collection name parameters
         self.inputCollName = 'HcalSimHits'
         self.inputPassName = ''
@@ -155,7 +160,7 @@ def __init__(self, instance_name = 'hcalRecon') :
 
         # avg parameters
         self.avgToaThreshold = 1.6 # mV - correction config only
-        self.avgGain = 1.2 # correction config only 
+        self.avgGain = 5100./1023. # correction config only 
         self.avgPedestal = 1. #noise config only   
 
 class HcalSingleEndRecProducer(Producer) :

Hcal/python/hcal_hardcoded_conditions.py

@@ -14,7 +14,7 @@
 HcalTrigPrimConditionsHardcode.validForAllRows([ 1 , # ADC_PEDESTAL -- should match value from HgcrocEmulator
                                                  5 , # ADC_THRESHOLD -- current noise is
                                                  1,  # TOT_PEDESTAL -- currently set to match ADC pedestal
-                                                 10000,  # TOT_THRESHOLD -- rather large value...
+                                                 5100.,  # TOT_THRESHOLD 
                                                  # Rounding because trigger primitives shouldn't represent floating point operations.
                                                  # See https://github.com/LDMX-Software/Hcal/issues/66#issuecomment-1719663799
                                                  round(2.5) ] # TOT_GAIN, ratio of recon TOT gain over recon ADC gain
@@ -24,17 +24,18 @@
 adc_pedestal.validForAllRows([1.]) # should match HgcrocEmulator
 
 adc_gain = SimpleCSVDoubleTableProvider("hcal_adc_gain",["gain"])
-adc_gain.validForAllRows([1.2]) # 4 ADCs per PE - maxADCRange/readoutPadCapacitance/1024
+adc_gain.validForAllRows([5100./1024.]) 
 
-tot_calib = SimpleCSVDoubleTableProvider("hcal_tot_calibration",
+tot_calib = SimpleCSVDoubleTableProvider("hcal_tot_pedestal",
                                          ["m_adc_i","cut_point_tot","high_slope","high_offset",
                                           "low_slope","low_power","low_offset","tot_not",
-                                          "channel","flagged"])
+                                          "channel","flagged","tot_pedestal","tot_gain"])
 tot_calib.validForAllRows([
     1.,1.,1.,1.,
     1.,1.,1.,1.,
+    1.,1.,
     1.,1.
-]) # dummy value since TOT is not implemented
+]) # needs to be tweaked 
 
 toa_calib =  SimpleCSVDoubleTableProvider("hcal_toa_calibration",
                                           ["bx_shift","mean_shift"])
@@ -59,13 +60,13 @@
 
 HcalHgcrocConditionsHardcode.validForAllRows([
     1. , #PEDESTAL 
-    0.02*5/1.2, #NOISE - 0.02 PE with 1 PE ~ 5mV and gain = 1.2
+    0.02*5/(5100./1024.), #NOISE - 0.02 PE with 1 PE ~ 5mV and gain = 5100./1024.
     12.5, #MEAS_TIME - ns - clock_cycle/2 - defines the point in the BX where an in-time (time=0 in times vector) hit would arrive
-    20., #PAD_CAPACITANCE - pF
+    13/5.1, #PAD_CAPACITANCE - pF
     200., #TOT_MAX - ns - maximum time chip would be in TOT mode
-    10240. / 200., #DRAIN_RATE - fC/ns - dummy value for now
-    1.2, #GAIN - large ADC gain for now - conversion from ADC to mV
+    320000./200., #DRAIN_RATE - fC/ns - to drain maximum charge of 320000fC in 200ns
+    5100./1024., #GAIN - take 160 fC - 13 pC to be the linear range of ADC -> V_{adc_max} = 13pC/13/5.1pF = 5100mV
     1. + 4., #READOUT_THRESHOLD - 4 ADC counts above pedestal
-    1.*1.2 + 1*5, #TOA_THRESHOLD - mV - 1 PE above pedestal ( 1 PE  - 5 mV conversion)
-    10000., #TOT_THRESHOLD - mV - very large for now
+    1.*(5100./1024.) + 1*5, #TOA_THRESHOLD - mV - 1 PE above pedestal ( 1 PE  - 5 mV conversion)
+    5100., #TOT_THRESHOLD - mV - TOT begins at V_{adc_max}
     ])

Hcal/src/Hcal/HcalDigiProducer.cxx

@@ -8,7 +8,10 @@
 
 #include "Hcal/HcalDigiProducer.h"
 
+#include <fstream>
+
 #include "Framework/RandomNumberSeedService.h"
+#include "Tools/PulseRecord.h"
 
 namespace hcal {
 
@@ -39,6 +42,10 @@ void HcalDigiProducer::configure(framework::config::Parameters& ps) {
   // and generate pedestal noise digis in every empty channel
   zeroSuppression_ = ps.getParameter<bool>("zeroSuppression");
 
+  // If true, save digis in a txt file
+  digiEmulationDebug_ = ps.getParameter<bool>("digiEmulationDebug");
+  digiEmulationFile_ = ps.getParameter<std::string>("digiEmulationFile");
+
   // collection names
   inputCollName_ = ps.getParameter<std::string>("inputCollName");
   inputPassName_ = ps.getParameter<std::string>("inputPassName");
@@ -134,6 +141,10 @@ void HcalDigiProducer::produce(framework::Event& event) {
     std::vector<std::pair<double, double>> pulses_posend;
     std::vector<std::pair<double, double>> pulses_negend;
 
+    // For plotting purposes
+    std::vector<std::tuple<double, int, int>> DigiToPlot;
+    DigiToPlot.clear();
+
     for (auto psimHit : simBar.second) {
       const ldmx::SimCalorimeterHit& simHit = *psimHit;
 
@@ -223,6 +234,10 @@ void HcalDigiProducer::produce(framework::Event& event) {
       for (int iContrib = 0; iContrib < simHit.getNumberOfContribs();
            iContrib++) {
         double voltage = simHit.getContrib(iContrib).edep * MeV_;
+
+        std::cout << "  Energy deposited: " << simHit.getContrib(iContrib).edep
+                  << std::endl;
+
         double time =
             simHit.getContrib(iContrib).time;  // global time (t=0ns at target)
         time -= position.at(2) /
@@ -255,8 +270,37 @@ void HcalDigiProducer::produce(framework::Event& event) {
 
       bool posEndActivity =
           hgcroc_->digitize(posendID.raw(), pulses_posend, digiToAddPosend);
+      if (digiEmulationDebug_) {
+        // Recording digitized pulses in PulseRecord
+        const auto& pulseRecord = hgcroc_->getPulseRecord();
+        for (const auto& record : pulseRecord) {
+          DigiToPlot.emplace_back(record.getVolts(), record.getADC(),
+                                  record.getTOT());
+        }
+      }
+
       bool negEndActivity =
           hgcroc_->digitize(negendID.raw(), pulses_negend, digiToAddNegend);
+      if (digiEmulationDebug_) {
+        // Recording digitized pulses in PulseRecord
+        const auto& pulseRecord = hgcroc_->getPulseRecord();
+        for (const auto& record : pulseRecord) {
+          DigiToPlot.emplace_back(record.getVolts(), record.getADC(),
+                                  record.getTOT());
+        }
+      }
+
+      if (digiEmulationDebug_) {
+        // Print Voltages, adc/tot values to the same txt file
+        std::ofstream dataFile(digiEmulationFile_, std::ios::app);
+        for (const auto& entry : DigiToPlot) {
+          dataFile << std::get<0>(entry) << " " << std::get<1>(entry) << " "
+                   << std::get<2>(entry) << "\n";
+        }
+        // dataFile << "Section: " << section << " Layer: " << layer << " Strip:
+        // " << strip << "\n";
+        dataFile.close();
+      }
 
       if (posEndActivity && negEndActivity && zeroSuppression_) {
         hcalDigis.addDigi(posendID.raw(), digiToAddPosend);
@@ -280,7 +324,6 @@ void HcalDigiProducer::produce(framework::Event& event) {
           hcalDigis.addDigi(negendID.raw(), digi);
         }
       }
-
     } else {
       bool is_posend = false;
       std::vector<ldmx::HgcrocDigiCollection::Sample> digiToAdd;
@@ -300,6 +343,23 @@ void HcalDigiProducer::produce(framework::Event& event) {
               hgcroc_->noiseDigi(digiID.raw(), 0.0);
           hcalDigis.addDigi(digiID.raw(), digi);
         }
+        // Recording digitized pulses in PulseRecord
+        // Access PulseRecord
+        const auto& pulseRecord = hgcroc_->getPulseRecord();
+
+        for (const auto& record : pulseRecord) {
+          DigiToPlot.emplace_back(record.getVolts(), record.getADC(),
+                                  record.getTOT());
+        }
+
+        // Print Voltages, adc/tot values to the same txt file
+        std::ofstream dataFile("Voltage_ADC_TOT.txt", std::ios::app);
+        for (const auto& entry : DigiToPlot) {
+          dataFile << std::get<0>(entry) << " " << std::get<1>(entry) << " "
+                   << std::get<2>(entry) << "\n";
+        }
+        dataFile.close();
+
       } else {
         ldmx::HcalDigiID digiID(section, layer, strip, 1);
         if (hgcroc_->digitize(digiID.raw(), pulses_negend, digiToAdd)) {
@@ -309,13 +369,32 @@ void HcalDigiProducer::produce(framework::Event& event) {
               hgcroc_->noiseDigi(digiID.raw(), 0.0);
           hcalDigis.addDigi(digiID.raw(), digi);
         }
+
+        // Recording digitized pulses in PulseRecord
+        // Access PulseRecord
+        const auto& pulseRecord = hgcroc_->getPulseRecord();
+
+        for (const auto& record : pulseRecord) {
+          DigiToPlot.emplace_back(record.getVolts(), record.getADC(),
+                                  record.getTOT());
+        }
+
+        // Print Voltages, adc/tot values to the same txt file
+        std::ofstream dataFile("Voltage_ADC_TOT.txt", std::ios::app);
+        for (const auto& entry : DigiToPlot) {
+          dataFile << std::get<0>(entry) << " " << std::get<1>(entry) << " "
+                   << std::get<2>(entry) << "\n";
+        }
+        dataFile.close();
       }
     }
   }
 
   /******************************************************************************************
    * Noise Simulation on Empty Channels
    *****************************************************************************************/
+  // bool noise_ = false;
+  // std::cout << "  noise: " << noise_ << std::endl;
   if (noise_) {
     std::vector<ldmx::HcalDigiID> channelMap;
     int numChannels = 0;

Tools/include/Tools/HgcrocEmulator.h

@@ -7,6 +7,7 @@
 #include "Recon/Event/HgcrocDigiCollection.h"
 #include "SimCore/Event/SimCalorimeterHit.h"
 #include "Tools/NoiseGenerator.h"
+#include "Tools/PulseRecord.h"
 
 //----------//
 //   ROOT   //
@@ -441,6 +442,14 @@ class HgcrocEmulator {
    */
   mutable TF1 pulseFunc_;
 
+  mutable std::vector<ldmx::PulseRecord> pulseRecord_;
+
+  // Getter to access pulseRecord_
+ public:
+  const std::vector<ldmx::PulseRecord>& getPulseRecord() const {
+    return pulseRecord_;
+  }
+
 };  // HgcrocEmulator
 
 }  // namespace ldmx

Tools/include/Tools/PulseRecord.h

@@ -0,0 +1,26 @@
+#ifndef TOOLS_PULSERECORD_H_
+#define TOOLS_PULSERECORD_H_
+
+// new class PulseRecord for plotting purposes
+
+namespace ldmx {
+
+class PulseRecord {
+ public:
+  // Constructor to initialize Voltage, ADC, and TOT
+  PulseRecord(double volts, int adc, int tot)
+      : volts_(volts), adc_(adc), tot_(tot) {}
+
+  // Getters for the recorded data
+  double getVolts() const { return volts_; }
+  int getADC() const { return adc_; }
+  int getTOT() const { return tot_; }
+
+ private:
+  double volts_;
+  int adc_;
+  int tot_;
+};
+}  // namespace ldmx
+
+#endif

Tools/src/Tools/HgcrocEmulator.cxx

@@ -48,7 +48,8 @@ bool HgcrocEmulator::digitize(
     std::vector<ldmx::HgcrocDigiCollection::Sample> &digiToAdd) const {
   // step 0: prepare ourselves for emulation
 
-  digiToAdd.clear();  // make sure it is clean
+  digiToAdd.clear();     // make sure it is clean
+  pulseRecord_.clear();  // clear pulseRecord
 
   // Configure chip settings based off of table (that may have been passed)
   double totMax = getCondition(channelID, "TOT_MAX");
@@ -100,8 +101,20 @@ bool HgcrocEmulator::digitize(
         continue;  // if this hit wasn't in the current BX, continue...
 
       double vpeak = pulse(hit.second);
+      double bxvolts = pulse((iADC - iSOI_) * clockCycle_);
+
+      // std::cout << "DEBUG PRINT - COMPARE IN COMING VOLTAGE AND PULSE FUNC
+      // VOLTAGE" << std::endl; std::cout << "    Incoming peak Voltage: " <<
+      // vpeak << std::endl; std::cout << "    Pulse func Voltage: " << bxvolts
+      // << std::endl; std::cout << " " << std::endl;
+
+      // if (vpeak > totThreshold){
+      // startTOT = true;
+      // if (toverTOT < hit.second)
+      // toverTOT = hit.second;  // use the latest time in the window
+      //}
 
-      if (vpeak > totThreshold) {
+      if (bxvolts > totThreshold) {
         startTOT = true;
         if (toverTOT < hit.second)
           toverTOT = hit.second;  // use the latest time in the window
@@ -133,6 +146,8 @@ bool HgcrocEmulator::digitize(
       double charge_deposited =
           (pulse(toverTOT) - gain * pedestal) * padCapacitance;
 
+      double voltaGe = (pulse(toverTOT) - gain * pedestal);
+
       // Measure Time Over Threshold (TOT) by using the drain rate.
       //  1. Use drain rate to see how long it takes for the charge to drain off
       //  2. Translate this into DIGI samples
@@ -149,7 +164,7 @@ bool HgcrocEmulator::digitize(
       //  to measure a maximum of tot Max [ns]
       int tdc_counts = int(tot * 4096 / totMax) + pedestal;
 
-      // were we already over TOA?  TOT is reported in BX where TOA went over
+      // were we already over TOnA?  TOT is reported in BX where TOA went over
       // threshold...
       int toa{0};
       if (wasTOA) {
@@ -172,6 +187,10 @@ bool HgcrocEmulator::digitize(
           toa                     // TOA is third measurement
       );
 
+      // Record in PulseRecord
+      pulseRecord_.clear();
+      pulseRecord_.emplace_back(voltaGe, 0, tdc_counts);
+
       // TODO: properly handle saturation and recovery, eventually.
       // Now just kill everything...
       while (digiToAdd.size() < nADCs_) {
@@ -180,6 +199,7 @@ bool HgcrocEmulator::digitize(
       }
 
       return true;  // always readout
+
     } else {
       // determine the voltage at the sampling time
       double bxvolts = pulse((iADC - iSOI_) * clockCycle_);
@@ -190,6 +210,12 @@ bool HgcrocEmulator::digitize(
       if (adc < 0) adc = 0;
       if (adc > 1023) adc = 1023;
 
+      // Record in PulseRecord
+      if (adc >= readoutThreshold) {
+        // pulseRecord_.clear();
+        pulseRecord_.emplace_back(bxvolts, adc, 0);
+      }
+
       // check for TOA
       int toa(0);
       if (pulse(startBX) < toaThreshold && overTOA) {
@@ -210,6 +236,7 @@ bool HgcrocEmulator::digitize(
           adc,                    // ADC[t] is the second field
           toa                     // TOA is third measurement
       );
+
     }  // TOT or ADC Mode
   }    // sampling baskets