diff --git a/DataFormats/FTLDigi/interface/PMTDSimAccumulator.h b/DataFormats/FTLDigi/interface/PMTDSimAccumulator.h index 02b1b7a83623c..6c918b42210ad 100644 --- a/DataFormats/FTLDigi/interface/PMTDSimAccumulator.h +++ b/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) diff --git a/SimFastTiming/FastTimingCommon/interface/BTLElectronicsSim.h b/SimFastTiming/FastTimingCommon/interface/BTLElectronicsSim.h index 0de47b0c428bc..fd648bf4cdc88 100644 --- a/SimFastTiming/FastTimingCommon/interface/BTLElectronicsSim.h +++ b/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_; diff --git a/SimFastTiming/FastTimingCommon/interface/MTDDigitizer.h b/SimFastTiming/FastTimingCommon/interface/MTDDigitizer.h index f5dcd9c55a753..12c20345d2956 100644 --- a/SimFastTiming/FastTimingCommon/interface/MTDDigitizer.h +++ b/SimFastTiming/FastTimingCommon/interface/MTDDigitizer.h @@ -59,17 +59,17 @@ namespace mtd_digitizer { const float minCharge, const float maxCharge) { constexpr auto nEnergies = std::tuple_size::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) { @@ -77,8 +77,8 @@ namespace mtd_digitizer { 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(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 diff --git a/SimFastTiming/FastTimingCommon/interface/MTDDigitizerTypes.h b/SimFastTiming/FastTimingCommon/interface/MTDDigitizerTypes.h index 5dc8b89b31803..e38e83a3ef3a1 100644 --- a/SimFastTiming/FastTimingCommon/interface/MTDDigitizerTypes.h +++ b/SimFastTiming/FastTimingCommon/interface/MTDDigitizerTypes.h @@ -15,12 +15,15 @@ namespace mtd_digitizer { typedef std::array MTDSimHitData; struct MTDCellInfo { - //1st array=energy, 2nd array=time-of-flight - std::array 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 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 MTDSimHitDataAccumulator; + constexpr int kNumberOfBX = 15; + constexpr int kInTimeBX = 9; + } // namespace mtd_digitizer namespace std { diff --git a/SimFastTiming/FastTimingCommon/python/mtdDigitizer_cfi.py b/SimFastTiming/FastTimingCommon/python/mtdDigitizer_cfi.py index 92d51a5215226..fa39fe5e31ffe 100644 --- a/SimFastTiming/FastTimingCommon/python/mtdDigitizer_cfi.py +++ b/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] diff --git a/SimFastTiming/FastTimingCommon/src/BTLBarDeviceSim.cc b/SimFastTiming/FastTimingCommon/src/BTLBarDeviceSim.cc index e922ed354437b..b7513b71692b7 100644 --- a/SimFastTiming/FastTimingCommon/src/BTLBarDeviceSim.cc +++ b/SimFastTiming/FastTimingCommon/src/BTLBarDeviceSim.cc @@ -84,36 +84,43 @@ void BTLBarDeviceSim::getHitsResponse(const std::vector(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 } diff --git a/SimFastTiming/FastTimingCommon/src/BTLElectronicsSim.cc b/SimFastTiming/FastTimingCommon/src/BTLElectronicsSim.cc index 4db5c5464b780..c1cfb727bac37 100644 --- a/SimFastTiming/FastTimingCommon/src/BTLElectronicsSim.cc +++ b/SimFastTiming/FastTimingCommon/src/BTLElectronicsSim.cc @@ -23,6 +23,7 @@ BTLElectronicsSim::BTLElectronicsSim(const edm::ParameterSet& pset) DarkCountRate_(pset.getParameter("DarkCountRate")), SigmaElectronicNoise_(pset.getParameter("SigmaElectronicNoise")), SigmaClock_(pset.getParameter("SigmaClock")), + smearTimeForOOTtails_(pset.getParameter("SmearTimeForOOTtails")), Npe_to_pC_(pset.getParameter("Npe_to_pC")), Npe_to_V_(pset.getParameter("Npe_to_V")), adcNbits_(pset.getParameter("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, diff --git a/SimFastTiming/FastTimingCommon/src/BTLTileDeviceSim.cc b/SimFastTiming/FastTimingCommon/src/BTLTileDeviceSim.cc index b2c52331b8732..74706e599818d 100644 --- a/SimFastTiming/FastTimingCommon/src/BTLTileDeviceSim.cc +++ b/SimFastTiming/FastTimingCommon/src/BTLTileDeviceSim.cc @@ -84,14 +84,16 @@ void BTLTileDeviceSim::getHitsResponse(const std::vector 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 }