Merge pull request #49221 from casarsa/mc-btlnewdigi_squash

Update of the BTL digitization model

Author: cmsbuild

RecoLocalFastTime/FTLCommonAlgos/interface/MTDTimeCalib.h

@@ -20,11 +20,8 @@ class MTDTimeCalib {
 private:
   const MTDGeometry* geom_;
   const MTDTopology* topo_;
-  float btlTimeOffset_;
-  float etlTimeOffset_;
 
   //specific paramters from BTL simulation
-  float btlLightCollTime_;
   float btlLightCollSlope_;
 };
 

RecoLocalFastTime/FTLCommonAlgos/plugins/BTLUncalibRecHitAlgo.cc

@@ -10,15 +10,14 @@ class BTLUncalibRecHitAlgo : public BTLUncalibratedRecHitAlgoBase {
   /// Constructor
   BTLUncalibRecHitAlgo(const edm::ParameterSet& conf, edm::ConsumesCollector& sumes)
       : MTDUncalibratedRecHitAlgoBase<BTLDataFrame>(conf, sumes),
-        adcNBits_(conf.getParameter<uint32_t>("adcNbits")),
-        adcSaturation_(conf.getParameter<double>("adcSaturation")),
-        adcLSB_(adcSaturation_ / (1 << adcNBits_)),
-        toaLSBToNS_(conf.getParameter<double>("toaLSB_ns")),
+        invLightSpeedLYSO_(conf.getParameter<double>("invLightSpeedLYSO")),
+        c_LYSO_(1. / invLightSpeedLYSO_),
+        npeToADC_(conf.getParameter<std::vector<double>>("npeToADC")),
+        npePerMeV_(conf.getParameter<double>("npePerMeV")),
+        invADCPerMeV_(1. / (npeToADC_[1] * npePerMeV_)),
+        tdc_to_ns_(conf.getParameter<double>("tdcLSB_ns")),
         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")),
-        c_LYSO_(conf.getParameter<double>("c_LYSO")) {}
+        timeWalkCorr_(conf.getParameter<std::string>("timeWalkCorrection")) {}
 
   /// Destructor
   ~BTLUncalibRecHitAlgo() override {}
@@ -31,26 +30,20 @@ class BTLUncalibRecHitAlgo : public BTLUncalibratedRecHitAlgoBase {
   FTLUncalibratedRecHit makeRecHit(const BTLDataFrame& dataFrame) const final;
 
 private:
-  const uint32_t adcNBits_;
-  const double adcSaturation_;
-  const double adcLSB_;
-  const double toaLSBToNS_;
-  const reco::FormulaEvaluator timeError_;
-  const double timeCorr_p0_;
-  const double timeCorr_p1_;
-  const double timeCorr_p2_;
+  const double invLightSpeedLYSO_;
   const double c_LYSO_;
+  const std::vector<double> npeToADC_;
+  const double npePerMeV_;
+  const double invADCPerMeV_;
+  const double tdc_to_ns_;
+  const reco::FormulaEvaluator timeError_;
+  const reco::FormulaEvaluator timeWalkCorr_;
 };
 
 FTLUncalibratedRecHit BTLUncalibRecHitAlgo::makeRecHit(const BTLDataFrame& dataFrame) const {
-  // The reconstructed amplitudes and times are saved in a std::pair
-  //    BTL tile geometry (1 SiPM): only the first value of the amplitude
-  //                                and time pairs is used.
-  //    BTL bar geometry (2 SiPMs): both values of the amplitude and
-  //                                time pairs are filled.
-
-  std::pair<float, float> amplitude(0., 0.);
-  std::pair<float, float> time(0., 0.);
+  // The reconstructed amplitudes and times of the right and left hits are saved in a std::pair
+  std::pair<double, double> amplitude(0., 0.);
+  std::pair<double, double> time(0., 0.);
 
   unsigned char flag = 0;
 
@@ -59,30 +52,37 @@ FTLUncalibratedRecHit BTLUncalibRecHitAlgo::makeRecHit(const BTLDataFrame& dataF
 
   double nHits = 0.;
 
-  LogDebug("BTLUncalibRecHit") << "Original input time t1,t2 " << float(sampleRight.toa()) * toaLSBToNS_ << ", "
-                               << float(sampleLeft.toa()) * toaLSBToNS_ << std::endl;
+  LogDebug("BTLUncalibRecHit") << "Original input time t1, t2 " << double(sampleRight.toa()) * tdc_to_ns_ << ", "
+                               << double(sampleLeft.toa()) * tdc_to_ns_ << std::endl;
 
+  // --- Reconstruct amplitude and time of the crystal's right channel
   if (sampleRight.data() > 0) {
-    amplitude.first = float(sampleRight.data()) * adcLSB_;
-    time.first = float(sampleRight.toa()) * toaLSBToNS_;
+    // Correct the time of the right SiPM for the time-walk
+    amplitude.first = double(sampleRight.data());
+    time.first = double(sampleRight.toa()) -
+                 timeWalkCorr_.evaluate(std::array<double, 1>{{amplitude.first}}, std::array<double, 1>{{0.0}});
 
-    nHits += 1.;
+    // Convert ADC counts to MeV and TDC counts to ns
+    amplitude.first = (double(sampleRight.data()) - npeToADC_[0]) * invADCPerMeV_;
+    time.first *= tdc_to_ns_;
 
-    // Correct the time of the left SiPM for the time-walk
-    time.first -= timeCorr_p0_ * pow(amplitude.first, timeCorr_p1_) + timeCorr_p2_;
     flag |= 0x1;
+    nHits += 1.;
   }
 
-  // --- If available, reconstruct the amplitude and time of the second SiPM
+  // --- Reconstruct amplitude and time of the crystal's left channel
   if (sampleLeft.data() > 0) {
-    amplitude.second = float(sampleLeft.data()) * adcLSB_;
-    time.second = float(sampleLeft.toa()) * toaLSBToNS_;
+    // Correct the time of the left SiPM for the time-walk
+    amplitude.second = double(sampleLeft.data());
+    time.second = double(sampleLeft.toa()) -
+                  timeWalkCorr_.evaluate(std::array<double, 1>{{amplitude.second}}, std::array<double, 1>{{0.0}});
 
-    nHits += 1.;
+    // Convert ADC counts to MeV and TDC counts to ns
+    amplitude.second = (double(sampleLeft.data()) - npeToADC_[0]) * invADCPerMeV_;
+    time.second *= tdc_to_ns_;
 
-    // Correct the time of the right SiPM for the time-walk
-    time.second -= timeCorr_p0_ * pow(amplitude.second, timeCorr_p1_) + timeCorr_p2_;
     flag |= (0x1 << 1);
+    nHits += 1.;
   }
 
   // --- Calculate the error on the hit time using the provided parameterization
@@ -95,16 +95,16 @@ FTLUncalibratedRecHit BTLUncalibRecHitAlgo::makeRecHit(const BTLDataFrame& dataF
   // Calculate the position
   // Distance from center of bar to hit
 
-  float position = 0.5f * (c_LYSO_ * (time.second - time.first));
-  float positionError = BTLRecHitsErrorEstimatorIM::positionError();
+  double position = 0.5f * (c_LYSO_ * (time.second - time.first));
+  double positionError = BTLRecHitsErrorEstimatorIM::positionError();
 
   LogDebug("BTLUncalibRecHit") << "DetId: " << dataFrame.id().rawId() << " x position = " << position << " +/- "
                                << positionError;
   LogDebug("BTLUncalibRecHit") << "ADC+: set the charge to: (" << amplitude.first << ", " << amplitude.second << ")  ("
-                               << sampleRight.data() << ", " << sampleLeft.data() << ") " << adcLSB_ << ' '
+                               << sampleRight.data() << ", " << sampleLeft.data() << ") " << invADCPerMeV_ << ' '
                                << std::endl;
   LogDebug("BTLUncalibRecHit") << "TDC+: set the time to: (" << time.first << ", " << time.second << ")  ("
-                               << sampleRight.toa() << ", " << sampleLeft.toa() << ") " << toaLSBToNS_ << ' '
+                               << sampleRight.toa() << ", " << sampleLeft.toa() << ") " << tdc_to_ns_ << ' '
                                << std::endl;
 
   return FTLUncalibratedRecHit(

RecoLocalFastTime/FTLCommonAlgos/plugins/MTDRecHitAlgo.cc

@@ -74,7 +74,7 @@ FTLRecHit MTDRecHitAlgo::makeRecHit(const FTLUncalibratedRecHit& uRecHit, uint32
     }
   }
 
-  // --- Energy calibration: for the time being this is just a conversion pC --> MeV
+  // --- Energy calibration
   energy *= calibration_;
 
   // --- Time calibration: for the time being just removes a time offset in BTL

RecoLocalFastTime/FTLCommonAlgos/src/MTDTimeCalib.cc

@@ -9,12 +9,7 @@
 #include "Geometry/MTDCommonData/interface/MTDTopologyMode.h"
 
 MTDTimeCalib::MTDTimeCalib(edm::ParameterSet const& conf, const MTDGeometry* geom, const MTDTopology* topo)
-    : geom_(geom),
-      topo_(topo),
-      btlTimeOffset_(conf.getParameter<double>("BTLTimeOffset")),
-      etlTimeOffset_(conf.getParameter<double>("ETLTimeOffset")),
-      btlLightCollTime_(conf.getParameter<double>("BTLLightCollTime")),
-      btlLightCollSlope_(conf.getParameter<double>("BTLLightCollSlope")) {}
+    : geom_(geom), topo_(topo), btlLightCollSlope_(conf.getParameter<double>("BTLLightCollSlope")) {}
 
 float MTDTimeCalib::getTimeCalib(const MTDDetId& id) const {
   if (id.subDetector() != MTDDetId::FastTime) {
@@ -25,7 +20,6 @@ float MTDTimeCalib::getTimeCalib(const MTDDetId& id) const {
   float time_calib = 0.;
 
   if (id.mtdSubDetector() == MTDDetId::BTL) {
-    time_calib += btlTimeOffset_;
     BTLDetId hitId(id);
     //for BTL topology gives different layout id
     DetId geoId = hitId.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topo_->getMTDTopologyMode()));
@@ -47,7 +41,7 @@ float MTDTimeCalib::getTimeCalib(const MTDDetId& id) const {
           << "BTL topology mode " << static_cast<int>(btlL) << " unsupported! Aborting";
     }
   } else if (id.mtdSubDetector() == MTDDetId::ETL) {
-    time_calib += etlTimeOffset_;
+    // no correction applied for ETL
   } else {
     throw cms::Exception("MTDTimeCalib") << "MTDDetId: " << std::hex << id.rawId() << " is invalid!" << std::dec
                                          << std::endl;

RecoLocalFastTime/FTLRecProducers/plugins/MTDTimeCalibESProducer.cc

@@ -44,10 +44,7 @@ MTDTimeCalibESProducer::MTDTimeCalibESProducer(const edm::ParameterSet& p) {
 // Configuration descriptions
 void MTDTimeCalibESProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
   edm::ParameterSetDescription desc;
-  desc.add<double>("BTLTimeOffset", 0.)->setComment("Time offset (additive) to all the BTL RecHits [ns]");
-  desc.add<double>("ETLTimeOffset", 0.)->setComment("Time offset (additive) to all the ETL RecHits [ns]");
-  desc.add<double>("BTLLightCollTime", 0.2)->setComment("Light collection time for BTL tile geometry [ns]");
-  desc.add<double>("BTLLightCollSlope", 0.075)
+  desc.add<double>("BTLLightCollSlope", 0.095)
       ->setComment("Light collection slope for bar for BTL bar tile geometry [ns/cm]");
   descriptions.add("MTDTimeCalibESProducer", desc);
 }

RecoLocalFastTime/FTLRecProducers/python/MTDTimeCalibESProducers_cff.py

@@ -2,7 +2,5 @@
 
 from RecoLocalFastTime.FTLRecProducers.MTDTimeCalibESProducer_cfi import *
 
-# the following numbers are obtained on single pions 0.7-10 GeV noPU
-# to have backpropagated time average at 0
-MTDTimeCalibESProducer.BTLTimeOffset = cms.double(0.0115)
-MTDTimeCalibESProducer.ETLTimeOffset = cms.double(0.0066)
+from SimFastTiming.FastTimingCommon.mtdDigitizer_cfi import mtdDigitizer
+MTDTimeCalibESProducer.BTLLightCollSlope = mtdDigitizer.barrelDigitizer.DeviceSimulation.LightCollectionSlope

RecoLocalFastTime/FTLRecProducers/python/mtdRecHits_cfi.py

@@ -2,8 +2,8 @@
 
 _barrelAlgo = cms.PSet(
     algoName = cms.string("MTDRecHitAlgo"),
-    thresholdToKeep = cms.double(1.),          # MeV
-    calibrationConstant = cms.double(0.03125), # MeV/pC
+    thresholdToKeep = cms.double(1.), # [MeV]
+    calibrationConstant = cms.double(1.)
 )
 
 

RecoLocalFastTime/FTLRecProducers/python/mtdUncalibratedRecHits_cfi.py

@@ -2,34 +2,29 @@
 
 from SimFastTiming.FastTimingCommon.mtdDigitizer_cfi import mtdDigitizer
 
-
 _barrelAlgo = cms.PSet(
     algoName = cms.string("BTLUncalibRecHitAlgo"),
-    adcNbits = mtdDigitizer.barrelDigitizer.ElectronicsSimulation.adcNbits,
-    adcSaturation = mtdDigitizer.barrelDigitizer.ElectronicsSimulation.adcSaturation_MIP,
-    toaLSB_ns = mtdDigitizer.barrelDigitizer.ElectronicsSimulation.toaLSB_ns,
-    timeResolutionInNs = cms.string("0.308*pow(x,-0.4175)"), # [ns]
-    timeCorr_p0 = cms.double( 2.21103),
-    timeCorr_p1 = cms.double(-0.933552),
-    timeCorr_p2 = cms.double( 0.),
-    c_LYSO = cms.double(13.846235)     # in unit cm/ns
+    invLightSpeedLYSO = mtdDigitizer.barrelDigitizer.DeviceSimulation.LightCollectionSlope, # [ns/cm]
+    npeToADC = mtdDigitizer.barrelDigitizer.ElectronicsSimulation.PulseQParam, # Npe to ADC counts conversion
+    npePerMeV = mtdDigitizer.barrelDigitizer.DeviceSimulation.LightOutput, # [Npe/MeV]
+    tdcLSB_ns = cms.double(0.020), # [ns]
+    timeResolutionInNs = cms.string("0.0593858*pow(x,-1.02826)+0.0156719"), # [ns]
+    timeWalkCorrection = cms.string("1.9e6/0.020*pow(9.389e5/0.0348*(x+22.5),-0.663)-7.5e-4*x-3.5e-3") # linear ad hoc correction for bias from global delay removal
 )
 
-
 _endcapAlgo = cms.PSet(
     algoName      = cms.string("ETLUncalibRecHitAlgo"),
     adcNbits      = mtdDigitizer.endcapDigitizer.ElectronicsSimulation.adcNbits,
     adcSaturation = mtdDigitizer.endcapDigitizer.ElectronicsSimulation.adcSaturation_MIP,
     toaLSB_ns     = mtdDigitizer.endcapDigitizer.ElectronicsSimulation.toaLSB_ns,
     tofDelay      = mtdDigitizer.endcapDigitizer.DeviceSimulation.tofDelay,
     timeResolutionInNs = cms.string("0.0370"), # [ns]
-    timeCorr_p0 = cms.double(0.974683),
+    timeCorr_p0 = cms.double(0.967683), # 0.974683 - 0.007, ad hoc correction for bias from global delay removal
     timeCorr_p1 = cms.double(-0.237274),
     timeCorr_p2 = cms.double(0.021455),
     timeCorr_p3 = cms.double(-0.000727429)
 )
 
-
 mtdUncalibratedRecHits = cms.EDProducer(
     "MTDUncalibratedRecHitProducer",
     barrel = _barrelAlgo,