Merge pull request #43762 from parbol/PMR_ETL_Digitization_Model

A new digitization model for the ETL detector

DataFormats/FTLDigi/interface/ETLSample.h

@@ -18,23 +18,26 @@ class ETLSample {
     kColumnMask = 0x1f,
     kRowMask = 0x3f,
     kToAMask = 0x7ff,
-    kDataMask = 0xff
+    kDataMask = 0xff,
+    kToTMask = 0x7ff
   };
   enum ETLSampleShifts {
     kThreshShift = 31,
     kModeShift = 30,
     kColumnShift = 25,
     kRowShift = 19,
     kToAShift = 8,
-    kDataShift = 0
+    kDataShift = 0,
+    kToTShift = 0
   };
 
   /**
      @short CTOR
    */
-  ETLSample() : value_(0) {}
-  ETLSample(uint32_t value) : value_(value) {}
-  ETLSample(const ETLSample& o) : value_(o.value_) {}
+  ETLSample() : value_(0), valueToT_(0) {}
+  ETLSample(uint32_t value) : value_(value), valueToT_(0) {}
+  ETLSample(uint32_t value, uint32_t valueToT) : value_(value), valueToT_(valueToT) {}
+  ETLSample(const ETLSample& o) : value_(o.value_), valueToT_(o.valueToT_) {}
   ETLSample& operator=(const ETLSample&) = default;
 
   /**
@@ -45,16 +48,18 @@ class ETLSample {
   void setColumn(uint8_t col) { setWord(col, kColumnMask, kColumnShift); }
   void setRow(uint8_t row) { setWord(row, kRowMask, kRowShift); }
   void setToA(uint16_t toa) { setWord(toa, kToAMask, kToAShift); }
+  void setToT(uint16_t tot) { setWordToT(tot, kToTMask, kToTShift); }
   void setData(uint16_t data) { setWord(data, kDataMask, kDataShift); }
-  void set(bool thr, bool mode, uint16_t toa, uint16_t data, uint8_t row, uint8_t col) {
+  void set(bool thr, bool mode, uint16_t toa, uint16_t tot, uint16_t data, uint8_t row, uint8_t col) {
     value_ = (((uint32_t)thr & kThreshMask) << kThreshShift | ((uint32_t)mode & kModeMask) << kModeShift |
               ((uint32_t)col & kColumnMask) << kColumnShift | ((uint32_t)row & kRowMask) << kRowShift |
               ((uint32_t)toa & kToAMask) << kToAShift | ((uint32_t)data & kDataMask) << kDataShift);
+    valueToT_ = ((uint32_t)tot & kToTMask) << kToTShift;
   }
   void print(std::ostream& out = std::cout) {
     out << "(row,col) : (" << row() << ',' << column() << ") "
-        << "THR: " << threshold() << " Mode: " << mode() << " ToA: " << toa() << " Data: " << data() << " Raw=0x"
-        << std::hex << raw() << std::dec << std::endl;
+        << "THR: " << threshold() << " Mode: " << mode() << " ToA: " << toa() << " ToT: " << tot()
+        << " Data: " << data() << " Raw=0x" << std::hex << raw() << std::dec << std::endl;
   }
 
   /**
@@ -66,6 +71,7 @@ class ETLSample {
   uint32_t column() const { return ((value_ >> kColumnShift) & kColumnMask); }
   uint32_t row() const { return ((value_ >> kRowShift) & kRowMask); }
   uint32_t toa() const { return ((value_ >> kToAShift) & kToAMask); }
+  uint32_t tot() const { return ((valueToT_ >> kToTShift) & kToTMask); }
   uint32_t data() const { return ((value_ >> kDataShift) & kDataMask); }
   uint32_t operator()() { return value_; }
 
@@ -80,8 +86,16 @@ class ETLSample {
     value_ |= ((word & mask) << pos);
   }
 
+  void setWordToT(uint32_t word, uint32_t mask, uint32_t pos) {
+    //clear required bits
+    valueToT_ &= ~(mask << pos);
+    //now set the new value
+    valueToT_ |= ((word & mask) << pos);
+  }
+
   // a 32-bit word
   uint32_t value_;
+  uint32_t valueToT_;
 };
 
 #endif

DataFormats/FTLDigi/src/classes_def.xml

@@ -7,7 +7,10 @@
    <class name="edm::Wrapper<FTLDigiCollection>" />
    <class name="BTLSample"/>
    <class name="std::vector<BTLSample>"/>
-   <class name="ETLSample"/>
+   <class name="ETLSample" ClassVersion="11">
+     <version ClassVersion="11" checksum="1290724074"/>
+     <version ClassVersion="10" checksum="1790811181"/>
+   </class>  
    <class name="std::vector<ETLSample>"/>
    <class name="BTLDataFrame"/>
    <class name="std::vector<BTLDataFrame>"/>

RecoLocalFastTime/FTLCommonAlgos/plugins/ETLUncalibRecHitAlgo.cc

@@ -13,8 +13,11 @@ class ETLUncalibRecHitAlgo : public ETLUncalibratedRecHitAlgoBase {
         adcLSB_(adcSaturation_ / (1 << adcNBits_)),
         toaLSBToNS_(conf.getParameter<double>("toaLSB_ns")),
         tofDelay_(conf.getParameter<double>("tofDelay")),
-        timeError_(conf.getParameter<std::string>("timeResolutionInNs")) {}
-
+        timeError_(conf.getParameter<std::string>("timeResolutionInNs")),
+        timeCorr_p0_(conf.getParameter<double>("timeCorr_p0")),
+        timeCorr_p1_(conf.getParameter<double>("timeCorr_p1")),
+        timeCorr_p2_(conf.getParameter<double>("timeCorr_p2")),
+        timeCorr_p3_(conf.getParameter<double>("timeCorr_p3")) {}
   /// Destructor
   ~ETLUncalibRecHitAlgo() override {}
 
@@ -32,25 +35,44 @@ class ETLUncalibRecHitAlgo : public ETLUncalibratedRecHitAlgoBase {
   const double toaLSBToNS_;
   const double tofDelay_;
   const reco::FormulaEvaluator timeError_;
+  const double timeCorr_p0_;
+  const double timeCorr_p1_;
+  const double timeCorr_p2_;
+  const double timeCorr_p3_;
 };
 
 FTLUncalibratedRecHit ETLUncalibRecHitAlgo::makeRecHit(const ETLDataFrame& dataFrame) const {
   constexpr int iSample = 2;  //only in-time sample
   const auto& sample = dataFrame.sample(iSample);
-
   const std::array<double, 1> amplitudeV = {{double(sample.data()) * adcLSB_}};
-  // NB: Here amplitudeV is defined as an array in order to be used
-  //     below as an input to FormulaEvaluator::evaluate.
+
   double time = double(sample.toa()) * toaLSBToNS_ - tofDelay_;
+  double time_over_threshold = double(sample.tot()) * toaLSBToNS_;
   unsigned char flag = 0;
 
   LogDebug("ETLUncalibRecHit") << "ADC+: set the charge to: " << amplitudeV[0] << ' ' << sample.data() << ' ' << adcLSB_
                                << ' ' << std::endl;
-  LogDebug("ETLUncalibRecHit") << "ADC+: set the time to: " << time << ' ' << sample.toa() << ' ' << toaLSBToNS_ << ' '
-                               << std::endl;
+
+  if (time_over_threshold == 0) {
+    LogDebug("ETLUncalibRecHit") << "ADC+: set the time to: " << time << ' ' << sample.toa() << ' ' << toaLSBToNS_
+                                 << ' ' << std::endl;
+
+  } else {
+    // Time-walk correction for toa
+    double timeWalkCorr = timeCorr_p0_ + timeCorr_p1_ * time_over_threshold +
+                          timeCorr_p2_ * time_over_threshold * time_over_threshold +
+                          timeCorr_p3_ * time_over_threshold * time_over_threshold * time_over_threshold;
+
+    time -= timeWalkCorr;
+
+    LogDebug("ETLUncalibRecHit") << "ADC+: set the time to: " << time << ' ' << sample.toa() << ' ' << toaLSBToNS_
+                                 << " .Timewalk correction: " << timeWalkCorr << std::endl;
+  }
+
   LogDebug("ETLUncalibRecHit") << "Final uncalibrated amplitude : " << amplitudeV[0] << std::endl;
 
   const std::array<double, 1> emptyV = {{0.}};
+
   double timeError = timeError_.evaluate(amplitudeV, emptyV);
 
   return FTLUncalibratedRecHit(dataFrame.id(),
@@ -63,6 +85,5 @@ FTLUncalibratedRecHit ETLUncalibRecHitAlgo::makeRecHit(const ETLDataFrame& dataF
                                -1.f,
                                flag);
 }
-
 #include "FWCore/Framework/interface/MakerMacros.h"
 DEFINE_EDM_PLUGIN(ETLUncalibratedRecHitAlgoFactory, ETLUncalibRecHitAlgo, "ETLUncalibRecHitAlgo");

RecoLocalFastTime/FTLRecProducers/python/mtdUncalibratedRecHits_cfi.py

@@ -22,7 +22,11 @@
     adcSaturation = mtdDigitizer.endcapDigitizer.ElectronicsSimulation.adcSaturation_MIP,
     toaLSB_ns     = mtdDigitizer.endcapDigitizer.ElectronicsSimulation.toaLSB_ns,
     tofDelay      = mtdDigitizer.endcapDigitizer.DeviceSimulation.tofDelay,
-    timeResolutionInNs = cms.string("0.039") # [ns]
+    timeResolutionInNs = cms.string("0.039"), # [ns]
+    timeCorr_p0 = cms.double(0.667116),
+    timeCorr_p1 = cms.double(-0.500665),
+    timeCorr_p2 = cms.double(0.141136),
+    timeCorr_p3 = cms.double(-0.0169169)
 )
 
 

SimFastTiming/FastTimingCommon/interface/ETLDeviceSim.h

@@ -12,6 +12,8 @@
 #include "Geometry/Records/interface/MTDDigiGeometryRecord.h"
 #include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h"
 
+#include "CommonTools/Utils/interface/FormulaEvaluator.h"
+
 #include <tuple>
 
 namespace CLHEP {
@@ -34,10 +36,19 @@ class ETLDeviceSim {
 private:
   const edm::ESGetToken<MTDGeometry, MTDDigiGeometryRecord> geomToken_;
   const MTDGeometry* geom_;
-
   float MIPPerMeV_;
+  const float integratedLum_;
+  const reco::FormulaEvaluator fluence_;
+  const reco::FormulaEvaluator lgadGain_;
+  const reco::FormulaEvaluator lgadGainDegradation_;
+  const bool applyDegradation_;
   float bxTime_;
   float tofDelay_;
+  const reco::FormulaEvaluator MPVMuon_;
+  const reco::FormulaEvaluator MPVPion_;
+  const reco::FormulaEvaluator MPVKaon_;
+  const reco::FormulaEvaluator MPVElectron_;
+  const reco::FormulaEvaluator MPVProton_;
 };
 
 #endif

SimFastTiming/FastTimingCommon/interface/ETLElectronicsSim.h

@@ -16,6 +16,8 @@
 #include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h"
 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
 
+#include "SimFastTiming/FastTimingCommon/interface/ETLPulseShape.h"
+
 namespace mtd = mtd_digitizer;
 
 namespace CLHEP {
@@ -34,7 +36,8 @@ class ETLElectronicsSim {
 
   void runTrivialShaper(ETLDataFrame& dataFrame,
                         const mtd::MTDSimHitData& chargeColl,
-                        const mtd::MTDSimHitData& toa,
+                        const mtd::MTDSimHitData& toa1,
+                        const mtd::MTDSimHitData& toa2,
                         const uint8_t row,
                         const uint8_t column) const;
 
@@ -49,9 +52,8 @@ class ETLElectronicsSim {
   const bool debug_;
   const float bxTime_;
   const float integratedLum_;
-  const reco::FormulaEvaluator fluence_;
-  const reco::FormulaEvaluator lgadGain_;
-  const reco::FormulaEvaluator timeRes2_;
+
+  const ETLPulseShape etlPulseShape_;
 
   // adc/tdc bitwidths
   const uint32_t adcNbits_, tdcNbits_;
@@ -61,8 +63,14 @@ class ETLElectronicsSim {
   const float adcLSB_MIP_;
   const uint32_t adcBitSaturation_;
   const float adcThreshold_MIP_;
+  const float iThreshold_MIP_;
   const float toaLSB_ns_;
   const uint32_t tdcBitSaturation_;
+  const float referenceChargeColl_;
+  const float noiseLevel_;
+  const float sigmaDistorsion_;
+  const float sigmaTDC_;
+  const reco::FormulaEvaluator formulaLandauNoise_;
 };
 
 #endif

SimFastTiming/FastTimingCommon/interface/ETLPulseShape.h

@@ -0,0 +1,16 @@
+#ifndef __SimFastTiming_FastTimingCommon_ETLPulseShape_h__
+#define __SimFastTiming_FastTimingCommon_ETLPulseShape_h__
+
+#include "SimFastTiming/FastTimingCommon/interface/MTDShapeBase.h"
+
+class ETLPulseShape : public MTDShapeBase {
+public:
+  ETLPulseShape();
+
+  ~ETLPulseShape() override;
+
+protected:
+  void fillShape(MTDShapeBase::DVec& aVec) const override;
+};
+
+#endif

SimFastTiming/FastTimingCommon/interface/MTDShapeBase.h

@@ -21,6 +21,8 @@ class MTDShapeBase : public CaloVShape {
   double timeToRise() const override { return 0.; }
 
   std::array<float, 3> timeAtThr(const float scale, const float threshold1, const float threshold2) const;
+  float maximum() const;
+  float fallTime() const;
 
   static constexpr unsigned int kReadoutTimeInterval = 28;  // in nsec
   static constexpr unsigned int kNBinsPerNSec = 100;        // granularity of internal array
@@ -40,6 +42,7 @@ class MTDShapeBase : public CaloVShape {
   const double qNSecPerBin_;
   unsigned int indexOfMax_;
   double timeOfMax_;
+  float fallTime_;
   DVec shape_;
 };
 

SimFastTiming/FastTimingCommon/python/mtdDigitizer_cfi.py

@@ -66,16 +66,23 @@
     premixStage1MinCharge = cms.double(1e-4),
     premixStage1MaxCharge = cms.double(1e6),
     DeviceSimulation  = cms.PSet(
-        bxTime            = cms.double(25),
-        tofDelay          = cms.double(1),
-        meVPerMIP         = cms.double(0.085), # from HGCal
+        bxTime               = cms.double(25),
+        IntegratedLuminosity = cms.double(1000.0),
+        FluenceVsRadius      = cms.string("1.937*TMath::Power(x,-1.706)"),
+        LGADGainVsFluence    = cms.string("TMath::Min(15.,30.-x)"),
+        LGADGainDegradation  = cms.string("TMath::Max(1.0, TMath::Min(x, x + 0.05/0.01 * (x - 1) + y * (1 - x)/0.01))"),
+        applyDegradation     = cms.bool(False),
+        tofDelay             = cms.double(1),
+        meVPerMIP            = cms.double(0.085), #from HGCAL
+        MPVMuon             = cms.string("1.21561e-05 + 8.89462e-07 / (x * x)"),
+        MPVPion             = cms.string("1.24531e-05 + 7.16578e-07 / (x * x)"),
+        MPVKaon             = cms.string("1.20998e-05 + 2.47192e-06 / (x * x * x)"),
+        MPVElectron         = cms.string("1.30030e-05 + 1.55166e-07 / (x * x)"),
+        MPVProton           = cms.string("1.13666e-05 + 1.20093e-05 / (x * x)")
         ),
     ElectronicsSimulation = cms.PSet(
         bxTime               = cms.double(25),
         IntegratedLuminosity = cms.double(1000.),      # [1/fb]
-        FluenceVsRadius      = cms.string("1.937*TMath::Power(x,-1.706)"),
-        LGADGainVsFluence    = cms.string("TMath::Min(15.,30.-x)"),
-        TimeResolution2      = cms.string("0.0225/x"), # [ns^2]
         # n bits for the ADC 
         adcNbits             = cms.uint32(8),
         # n bits for the TDC
@@ -84,8 +91,14 @@
         adcSaturation_MIP  = cms.double(25),
         # for different thickness
         adcThreshold_MIP   = cms.double(0.025),
+        iThreshold_MIP     = cms.double(0.9525),
         # LSB for time of arrival estimate from TDC in ns
         toaLSB_ns          = cms.double(0.013),
+        referenceChargeColl = cms.double(1.0),
+        noiseLevel          = cms.double(0.3554),
+        sigmaDistorsion     = cms.double(0.0),
+        sigmaTDC            = cms.double(0.010),
+        formulaLandauNoise  = cms.string("TMath::Max(0.020, 0.020 * (0.35 * (x - 1.0) + 1.0))") 
         )
 )