Hcal tot digitization

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
@@ -151,7 +151,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
@@ -59,13 +59,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./1023.), #NOISE - 0.02 PE with 1 PE ~ 5mV and gain = 5100./1023.
     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./1023., #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./1023.) + 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 {
 
@@ -130,6 +133,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;
 
@@ -213,6 +220,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) /
@@ -246,7 +257,44 @@ void HcalDigiProducer::produce(framework::Event& event) {
           hgcroc_->digitize(negendID.raw(), pulses_negend, digiToAddNegend)) {
         hcalDigis.addDigi(posendID.raw(), digiToAddPosend);
         hcalDigis.addDigi(negendID.raw(), digiToAddNegend);
-      }  // Back Hcal needs to digitize both pulses or none
+
+        // for(const auto& digi: hcalDigis){
+        // std::cout << "    DEBUG PRINT 2: " << digi.at(2) << std::endl;
+        // }
+      }
+
+      // Recording digitized pulses in PulseRecord
+      // Since two pulses are digitized, using existing loop is bit untidy.
+      // Access PulseRecord
+      if (hgcroc_->digitize(posendID.raw(), pulses_posend, digiToAddPosend) &&
+          hgcroc_->digitize(negendID.raw(), pulses_negend, digiToAddNegend)) {
+        const auto& pulseRecord = hgcroc_->getPulseRecord();
+
+        for (const auto& record : pulseRecord) {
+          DigiToPlot.emplace_back(record.getVolts(), record.getADC(),
+                                  record.getTOT());
+        }
+      }
+
+      if (hgcroc_->digitize(posendID.raw(), pulses_negend, digiToAddNegend) &&
+          hgcroc_->digitize(posendID.raw(), pulses_posend, digiToAddPosend)) {
+        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();
+
+      // Back Hcal needs to digitize both pulses or none
     } else {
       bool is_posend = false;
       std::vector<ldmx::HgcrocDigiCollection::Sample> digiToAdd;
@@ -262,18 +310,54 @@ void HcalDigiProducer::produce(framework::Event& event) {
         if (hgcroc_->digitize(digiID.raw(), pulses_posend, digiToAdd)) {
           hcalDigis.addDigi(digiID.raw(), digiToAdd);
         }
+        // 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)) {
           hcalDigis.addDigi(digiID.raw(), digiToAdd);
         }
+
+        // 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_) {
     int numChannels = 0;
     for (int section = 0; section < hcalGeometry.getNumSections(); section++) {
@@ -349,7 +433,8 @@ void HcalDigiProducer::produce(framework::Event& event) {
         }
       }
     }  // loop over noise amplitudes
-  }    // if we should add noise
+
+  }  // if we should add noise
 
   event.add(digiCollName_, hcalDigis);
 

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");
@@ -101,8 +102,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
@@ -134,6 +147,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
@@ -153,7 +168,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) {
@@ -176,14 +191,19 @@ 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_) {
         digiToAdd.emplace_back(true, false,  // flags to mark type of sample
                                0x3FF, 0x3FF, 0);
       }
 
-      return true;  // always readout
+      for (const auto &digi : digiToAdd) return true;  // always readout
+
     } else {
       // determine the voltage at the sampling time
       double bxvolts = pulse((iADC - iSOI_) * clockCycle_);
@@ -194,6 +214,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) {
@@ -214,6 +240,7 @@ bool HgcrocEmulator::digitize(
           adc,                    // ADC[t] is the second field
           toa                     // TOA is third measurement
       );
+
     }  // TOT or ADC Mode
   }    // sampling baskets