BTL simulation: OOT effects + premix fix

DataFormats/FTLDigi/interface/PMTDSimAccumulator.h

@@ -29,12 +29,12 @@ class PMTDSimAccumulator {
   };
   class Data {
   public:
-    constexpr static unsigned energyOffset = 15;
-    constexpr static unsigned energyMask = 0x1;
-    constexpr static unsigned sampleOffset = 11;
+    constexpr static unsigned energyOffset = 14;
+    constexpr static unsigned energyMask = 0x3;
+    constexpr static unsigned sampleOffset = 10;
     constexpr static unsigned sampleMask = 0xf;
     constexpr static unsigned dataOffset = 0;
-    constexpr static unsigned dataMask = 0x7ff;
+    constexpr static unsigned dataMask = 0x3ff;
 
     Data() : data_(0) {}
     Data(unsigned short ei, unsigned short si, unsigned short d)

SimFastTiming/FastTimingCommon/interface/BTLElectronicsSim.h

@@ -60,7 +60,7 @@ class BTLElectronicsSim {
   const float DarkCountRate_;
   const float SigmaElectronicNoise_;
   const float SigmaClock_;
-
+  const bool smearTimeForOOTtails_;
   const float Npe_to_pC_;
   const float Npe_to_V_;
 
@@ -80,6 +80,7 @@ class BTLElectronicsSim {
   const float sinPhi_;
 
   const float ScintillatorDecayTime2_;
+  const float ScintillatorDecayTimeInv_;
   const float SPTR2_;
   const float DCRxRiseTime_;
   const float SigmaElectronicNoise2_;

SimFastTiming/FastTimingCommon/interface/MTDDigitizer.h

@@ -59,26 +59,26 @@ namespace mtd_digitizer {
                                     const float minCharge,
                                     const float maxCharge) {
     constexpr auto nEnergies = std::tuple_size<decltype(MTDCellInfo().hit_info)>::value;
-    static_assert(nEnergies <= PMTDSimAccumulator::Data::energyMask + 1,
-                  "PMTDSimAccumulator bit pattern needs to updated");
-    static_assert(nSamples <= PMTDSimAccumulator::Data::sampleMask + 1,
-                  "PMTDSimAccumulator bit pattern needs to updated");
+    static_assert(nEnergies == PMTDSimAccumulator::Data::energyMask + 1,
+                  "PMTDSimAccumulator bit pattern needs to be updated");
+    static_assert(nSamples == PMTDSimAccumulator::Data::sampleMask,
+                  "PMTDSimAccumulator bit pattern needs to be updated");
 
     const float minPackChargeLog = minCharge > 0.f ? std::log(minCharge) : -2;
     const float maxPackChargeLog = std::log(maxCharge);
     constexpr uint16_t base = 1 << PMTDSimAccumulator::Data::sampleOffset;
 
     simResult.reserve(simData.size());
-    // mimicing the digitization
+    // mimicking the digitization
     for (const auto& elem : simData) {
       // store only non-zero
       for (size_t iEn = 0; iEn < nEnergies; ++iEn) {
         const auto& samples = elem.second.hit_info[iEn];
         for (size_t iSample = 0; iSample < nSamples; ++iSample) {
           if (samples[iSample] > minCharge) {
             unsigned short packed;
-            if (iEn == 1) {
-              // assuming linear range for tof of 0..26
+            if (iEn == 1 || iEn == 3) {
+              // assuming linear range for tof of 0..25
               packed = samples[iSample] / PREMIX_MAX_TOF * base;
             } else {
               packed = logintpack::pack16log(samples[iSample], minPackChargeLog, maxPackChargeLog, base);
@@ -107,14 +107,18 @@ namespace mtd_digitizer {
       size_t iEn = detIdIndexHitInfo.energyIndex();
       size_t iSample = detIdIndexHitInfo.sampleIndex();
 
+      if (iEn > PMTDSimAccumulator::Data::energyMask + 1 || iSample > PMTDSimAccumulator::Data::sampleMask)
+        throw cms::Exception("MTDDigitixer::loadSimHitAccumulator")
+            << "Index out of range: iEn = " << iEn << " iSample = " << iSample << std::endl;
+
       float value;
-      if (iEn == 1) {
+      if (iEn == 1 || iEn == 3) {
         value = static_cast<float>(detIdIndexHitInfo.data()) / base * PREMIX_MAX_TOF;
       } else {
         value = logintpack::unpack16log(detIdIndexHitInfo.data(), minPackChargeLog, maxPackChargeLog, base);
       }
 
-      if (iEn == 0) {
+      if (iEn == 0 || iEn == 2) {
         hit_info[iEn][iSample] += value;
       } else if (hit_info[iEn][iSample] == 0) {
         // For iEn==1 the digitizers just set the TOF of the first SimHit

SimFastTiming/FastTimingCommon/interface/MTDDigitizerTypes.h

@@ -15,12 +15,15 @@ namespace mtd_digitizer {
   typedef std::array<MTDSimData_t, nSamples> MTDSimHitData;
 
   struct MTDCellInfo {
-    //1st array=energy, 2nd array=time-of-flight
-    std::array<MTDSimHitData, 2> hit_info;
+    // for the BTL tile geometry and ETL:
+    //     1st array=energy, 2nd array=time-of-flight
+    // for the BTL bar geometry:
+    //     3rd array=energy (right side), 4th array=time-of-flight (right side)
+    std::array<MTDSimHitData, 4> hit_info;
   };
 
   // Maximum value of time-of-flight for premixing packing
-  constexpr float PREMIX_MAX_TOF = 26.0f;
+  constexpr float PREMIX_MAX_TOF = 25.0f;
 
   struct MTDCellId {
     MTDCellId() : detid_(0), row_(0), column_(0) {}
@@ -37,6 +40,9 @@ namespace mtd_digitizer {
   // intermediate det id for ETL
   typedef std::unordered_map<MTDCellId, MTDCellInfo> MTDSimHitDataAccumulator;
 
+  constexpr int kNumberOfBX = 15;
+  constexpr int kInTimeBX = 9;
+
 }  // namespace mtd_digitizer
 
 namespace std {

SimFastTiming/FastTimingCommon/python/mtdDigitizer_cfi.py

@@ -33,6 +33,7 @@
         SigmaElectronicNoise       = cms.double(1.),    # [p.e.]
         SigmaClock                 = cms.double(0.015), # [ns]
         CorrelationCoefficient     = cms.double(1.),
+        SmearTimeForOOTtails       = cms.bool(True),
         Npe_to_pC                  = cms.double(0.016), # [pC] 
         Npe_to_V                   = cms.double(0.0064),# [V] 
 

SimFastTiming/FastTimingCommon/src/BTLBarDeviceSim.cc

@@ -84,36 +84,43 @@ void BTLBarDeviceSim::getHitsResponse(const std::vector<std::tuple<int, uint32_t
     // --- Get the simHit energy and convert it from MeV to photo-electrons
     float Npe = 1000. * hit.energyLoss() * LightYield_ * LightCollEff_ * PDE_;
 
-    // --- Get the simHit time of arrival
-    float toa = std::get<2>(hitRef);
-
-    if (toa > bxTime_ || toa < 0)  //just consider BX==0
-      continue;
-
-    // --- Accumulate the energy of simHits in the same crystal for the BX==0
-    // this is to simulate the charge integration in a 25 ns window
-    (simHitIt->second).hit_info[0][0] += Npe;
-    (simHitIt->second).hit_info[0][1] += Npe;
-
+    // --- Calculate the light propagation time to the crystal bases (labeled L and R)
     double distR = 0.5 * topo.pitch().second - 0.1 * hit.localPosition().y();
     double distL = 0.5 * topo.pitch().second + 0.1 * hit.localPosition().y();
 
-    // This is for the layout with bars along phi
+    // This is for the layouts with bars along phi
     if (topo_->getMTDTopologyMode() == (int)BTLDetId::CrysLayout::bar ||
         topo_->getMTDTopologyMode() == (int)BTLDetId::CrysLayout::barphiflat) {
       distR = 0.5 * topo.pitch().first - 0.1 * hit.localPosition().x();
       distL = 0.5 * topo.pitch().first + 0.1 * hit.localPosition().x();
     }
 
-    double tR = toa + LightCollSlopeR_ * distR;
-    double tL = toa + LightCollSlopeL_ * distL;
+    double tR = std::get<2>(hitRef) + LightCollSlopeR_ * distR;
+    double tL = std::get<2>(hitRef) + LightCollSlopeL_ * distL;
+
+    // --- Accumulate in 15 buckets of 25ns (9 pre-samples, 1 in-time, 5 post-samples)
+    const int iBXR = std::floor(tR / bxTime_) + mtd_digitizer::kInTimeBX;
+    const int iBXL = std::floor(tL / bxTime_) + mtd_digitizer::kInTimeBX;
 
-    // --- Store the time of the first SimHit
-    if ((simHitIt->second).hit_info[1][0] == 0 || tR < (simHitIt->second).hit_info[1][0])
-      (simHitIt->second).hit_info[1][0] = tR;
+    // --- Right side
+    if (iBXR > 0 && iBXR < mtd_digitizer::kNumberOfBX) {
+      // Accumulate the energy of simHits in the same crystal
+      (simHitIt->second).hit_info[0][iBXR] += Npe;
 
-    if ((simHitIt->second).hit_info[1][1] == 0 || tL < (simHitIt->second).hit_info[1][1])
-      (simHitIt->second).hit_info[1][1] = tL;
+      // Store the time of the first SimHit in the i-th BX
+      if ((simHitIt->second).hit_info[1][iBXR] == 0 || tR < (simHitIt->second).hit_info[1][iBXR])
+        (simHitIt->second).hit_info[1][iBXR] = tR - (iBXR - mtd_digitizer::kInTimeBX) * bxTime_;
+    }
+
+    // --- Left side
+    if (iBXL > 0 && iBXL < mtd_digitizer::kNumberOfBX) {
+      // Accumulate the energy of simHits in the same crystal
+      (simHitIt->second).hit_info[2][iBXL] += Npe;
+
+      // Store the time of the first SimHit in the i-th BX
+      if ((simHitIt->second).hit_info[3][iBXL] == 0 || tL < (simHitIt->second).hit_info[3][iBXL])
+        (simHitIt->second).hit_info[3][iBXL] = tL - (iBXL - mtd_digitizer::kInTimeBX) * bxTime_;
+    }
 
   }  // hitRef loop
 }

SimFastTiming/FastTimingCommon/src/BTLElectronicsSim.cc

@@ -23,6 +23,7 @@ BTLElectronicsSim::BTLElectronicsSim(const edm::ParameterSet& pset)
       DarkCountRate_(pset.getParameter<double>("DarkCountRate")),
       SigmaElectronicNoise_(pset.getParameter<double>("SigmaElectronicNoise")),
       SigmaClock_(pset.getParameter<double>("SigmaClock")),
+      smearTimeForOOTtails_(pset.getParameter<bool>("SmearTimeForOOTtails")),
       Npe_to_pC_(pset.getParameter<double>("Npe_to_pC")),
       Npe_to_V_(pset.getParameter<double>("Npe_to_V")),
       adcNbits_(pset.getParameter<uint32_t>("adcNbits")),
@@ -37,6 +38,7 @@ BTLElectronicsSim::BTLElectronicsSim(const edm::ParameterSet& pset)
       cosPhi_(0.5 * (sqrt(1. + CorrCoeff_) + sqrt(1. - CorrCoeff_))),
       sinPhi_(0.5 * CorrCoeff_ / cosPhi_),
       ScintillatorDecayTime2_(ScintillatorDecayTime_ * ScintillatorDecayTime_),
+      ScintillatorDecayTimeInv_(1. / ScintillatorDecayTime_),
       SPTR2_(SinglePhotonTimeResolution_ * SinglePhotonTimeResolution_),
       DCRxRiseTime_(DarkCountRate_ * ScintillatorRiseTime_),
       SigmaElectronicNoise2_(SigmaElectronicNoise_ * SigmaElectronicNoise_),
@@ -48,17 +50,20 @@ void BTLElectronicsSim::run(const mtd::MTDSimHitDataAccumulator& input,
   MTDSimHitData chargeColl, toa1, toa2;
 
   for (MTDSimHitDataAccumulator::const_iterator it = input.begin(); it != input.end(); it++) {
+    // --- Digitize only the in-time bucket:
+    const unsigned int iBX = mtd_digitizer::kInTimeBX;
+
     chargeColl.fill(0.f);
     toa1.fill(0.f);
     toa2.fill(0.f);
-    for (size_t i = 0; i < it->second.hit_info[0].size(); i++) {
+    for (size_t iside = 0; iside < 2; iside++) {
       // --- Fluctuate the total number of photo-electrons
-      float Npe = CLHEP::RandPoissonQ::shoot(hre, (it->second).hit_info[0][i]);
+      float Npe = CLHEP::RandPoissonQ::shoot(hre, (it->second).hit_info[2 * iside][iBX]);
       if (Npe < EnergyThreshold_)
         continue;
 
       // --- Get the time of arrival and add a channel time offset
-      float finalToA1 = (it->second).hit_info[1][i] + ChannelTimeOffset_;
+      float finalToA1 = (it->second).hit_info[1 + 2 * iside][iBX] + ChannelTimeOffset_;
 
       if (smearChannelTimeOffset_ > 0.) {
         float timeSmearing = CLHEP::RandGaussQ::shoot(hre, 0., smearChannelTimeOffset_);
@@ -77,6 +82,26 @@ void BTLElectronicsSim::run(const mtd::MTDSimHitDataAccumulator& input,
       float finalToA2 = finalToA1 + times[1];
       finalToA1 += times[0];
 
+      // --- Estimate the time uncertainty due to photons from earlier OOT hits in the current BTL cell
+      if (smearTimeForOOTtails_) {
+        float rate_oot = 0.;
+        // Loop on earlier OOT hits
+        for (int ibx = 0; ibx < mtd_digitizer::kInTimeBX; ++ibx) {
+          if ((it->second).hit_info[2 * iside][ibx] > 0.) {
+            float hit_time = (it->second).hit_info[1 + 2 * iside][ibx] + bxTime_ * (ibx - mtd_digitizer::kInTimeBX);
+            float npe_oot = CLHEP::RandPoissonQ::shoot(hre, (it->second).hit_info[2 * iside][ibx]);
+            rate_oot += npe_oot * exp(hit_time * ScintillatorDecayTimeInv_) * ScintillatorDecayTimeInv_;
+          }
+        }  // ibx loop
+
+        if (rate_oot > 0.) {
+          float sigma_oot = sqrt(rate_oot * ScintillatorRiseTime_) * ScintillatorDecayTime_ / Npe;
+          float smearing_oot = CLHEP::RandGaussQ::shoot(hre, 0., sigma_oot);
+          finalToA1 += smearing_oot;
+          finalToA2 += smearing_oot;
+        }
+      }  // if smearTimeForOOTtails_
+
       // --- Uncertainty due to the fluctuations of the n-th photon arrival time:
       if (testBeamMIPTimeRes_ > 0.) {
         // In this case the time resolution is parametrized from the testbeam.
@@ -115,17 +140,19 @@ void BTLElectronicsSim::run(const mtd::MTDSimHitDataAccumulator& input,
       finalToA1 += cosPhi_ * smearing_thr1_uncorr + sinPhi_ * smearing_thr2_uncorr;
       finalToA2 += sinPhi_ * smearing_thr1_uncorr + cosPhi_ * smearing_thr2_uncorr;
 
-      chargeColl[i] = Npe * Npe_to_pC_;  // the p.e. number is here converted to pC
+      chargeColl[iside] = Npe * Npe_to_pC_;  // the p.e. number is here converted to pC
+
+      toa1[iside] = finalToA1;
+      toa2[iside] = finalToA2;
 
-      toa1[i] = finalToA1;
-      toa2[i] = finalToA2;
-    }
+    }  // iside loop
 
     //run the shaper to create a new data frame
     BTLDataFrame rawDataFrame(it->first.detid_);
     runTrivialShaper(rawDataFrame, chargeColl, toa1, toa2, it->first.row_, it->first.column_);
     updateOutput(output, rawDataFrame);
-  }
+
+  }  // MTDSimHitDataAccumulator loop
 }
 
 void BTLElectronicsSim::runTrivialShaper(BTLDataFrame& dataFrame,

SimFastTiming/FastTimingCommon/src/BTLTileDeviceSim.cc

@@ -84,14 +84,16 @@ void BTLTileDeviceSim::getHitsResponse(const std::vector<std::tuple<int, uint32_
     if (smearLightCollTime_ > 0.)
       toa += CLHEP::RandGaussQ::shoot(hre, 0., smearLightCollTime_);
 
-    if (toa > bxTime_ || toa < 0)  //just consider BX==0
+    // Accumulate in 15 buckets of 25ns (9 pre-samples, 1 in-time, 5 post-samples)
+    const int iBX = std::floor(toa / bxTime_) + mtd_digitizer::kInTimeBX;
+    if (iBX < 0 || iBX >= mtd_digitizer::kNumberOfBX)
       continue;
 
-    (simHitIt->second).hit_info[0][0] += Npe;
+    (simHitIt->second).hit_info[0][iBX] += Npe;
 
-    // --- Store the time of the first SimHit
-    if ((simHitIt->second).hit_info[1][0] == 0 || toa < (simHitIt->second).hit_info[1][0])
-      (simHitIt->second).hit_info[1][0] = toa;
+    // --- Store the time of the first SimHit in the i-th BX
+    if ((simHitIt->second).hit_info[1][iBX] == 0 || toa < (simHitIt->second).hit_info[1][iBX])
+      (simHitIt->second).hit_info[1][iBX] = toa - (iBX - mtd_digitizer::kInTimeBX) * bxTime_;
 
   }  // hitRef loop
 }