FastTimingCommon/src/ETLElectronicsSim.cc

0001 #include "SimFastTiming/FastTimingCommon/interface/ETLElectronicsSim.h"
0002
0003 #include "FWCore/Framework/interface/ConsumesCollector.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "CLHEP/Random/RandGaussQ.h"
0006
0007 using namespace mtd;
0008
0009 ETLElectronicsSim::ETLElectronicsSim(const edm::ParameterSet& pset, edm::ConsumesCollector iC)
0010     : geomToken_(iC.esConsumes()),
0011       geom_(nullptr),
0012       bxTime_(pset.getParameter<double>("bxTime")),
0013       integratedLum_(pset.getParameter<double>("IntegratedLuminosity")),
0014       adcNbits_(pset.getParameter<uint32_t>("adcNbits")),
0015       tdcNbits_(pset.getParameter<uint32_t>("tdcNbits")),
0016       adcSaturation_MIP_(pset.getParameter<double>("adcSaturation_MIP")),
0017       adcLSB_MIP_(adcSaturation_MIP_ / std::pow(2., adcNbits_)),
0018       adcBitSaturation_(std::pow(2, adcNbits_) - 1),
0019       adcThreshold_MIP_(pset.getParameter<double>("adcThreshold_MIP")),
0020       iThreshold_MIP_(pset.getParameter<double>("iThreshold_MIP")),
0021       toaLSB_ns_(pset.getParameter<double>("toaLSB_ns")),
0022       tdcBitSaturation_(std::pow(2, tdcNbits_) - 1),
0023       referenceChargeColl_(pset.getParameter<double>("referenceChargeColl")),
0024       noiseLevel_(pset.getParameter<double>("noiseLevel")),
0025       sigmaDistorsion_(pset.getParameter<double>("sigmaDistorsion")),
0026       sigmaTDC_(pset.getParameter<double>("sigmaTDC")),
0027       formulaLandauNoise_(pset.getParameter<std::string>("formulaLandauNoise")) {}
0028
0029 void ETLElectronicsSim::getEventSetup(const edm::EventSetup& evs) { geom_ = &evs.getData(geomToken_); }
0030
0031 void ETLElectronicsSim::run(const mtd::MTDSimHitDataAccumulator& input,
0032                             ETLDigiCollection& output,
0033                             CLHEP::HepRandomEngine* hre) const {
0034   MTDSimHitData chargeColl, toa1, toa2, tot;
0035
0036   std::vector<double> emptyV;
0037   std::vector<double> radius(1);
0038   std::vector<double> fluence(1);
0039   std::vector<double> chOverMPV(1);
0040
0041   for (MTDSimHitDataAccumulator::const_iterator it = input.begin(); it != input.end(); it++) {
0042     chargeColl.fill(0.f);
0043     toa1.fill(0.f);
0044     toa2.fill(0.f);
0045     tot.fill(0.f);
0046
0047     ETLDetId detId = it->first.detid_;
0048     DetId geoId = detId.geographicalId();
0049     const MTDGeomDet* thedet = geom_->idToDet(geoId);
0050     if (thedet == nullptr)
0051       throw cms::Exception("EtlElectronicsSim") << "GeographicalID: " << std::hex << geoId.rawId() << " ("
0052                                                 << detId.rawId() << ") is invalid!" << std::dec << std::endl;
0053
0054     const PixelTopology& topo = static_cast<const PixelTopology&>(thedet->topology());
0055     Local3DPoint local_point(topo.localX(it->first.row_), topo.localY(it->first.column_), 0.);
0056
0057     //Here we are sampling over all the buckets. However for the ETL there should be only one bucket at i = 0
0058     for (size_t i = 0; i < it->second.hit_info[0].size(); i++) {
0059       if ((it->second).hit_info[0][i] == 0) {
0060         continue;
0061       }
0062       if ((it->second).hit_info[0][i] < adcThreshold_MIP_)
0063         continue;
0064       // time of arrival
0065       double finalToA = (it->second).hit_info[1][i];
0066       double finalToC = (it->second).hit_info[1][i];
0067
0068       // fill the time and charge arrays
0069       const unsigned int ibucket = std::floor(finalToA / bxTime_);
0070       if ((i + ibucket) >= chargeColl.size())
0071         continue;
0072
0073       //Calculate the jitter
0074       double SignalToNoise =
0075           etlPulseShape_.maximum() * ((it->second).hit_info[0][i] / referenceChargeColl_) / noiseLevel_;
0076       double sigmaJitter1 = etlPulseShape_.timeOfMax() / SignalToNoise;
0077       double sigmaJitter2 = (etlPulseShape_.fallTime() - etlPulseShape_.timeOfMax()) / SignalToNoise;
0078       //Calculate the distorsion
0079       double sigmaDistorsion = sigmaDistorsion_;
0080       //Calculate the TDC
0081       double sigmaTDC = sigmaTDC_;
0082       //Calculate the Landau Noise
0083       chOverMPV[0] = (it->second).hit_info[0][i] / (it->second).hit_info[2][i];
0084       double sigmaLN = formulaLandauNoise_.evaluate(chOverMPV, emptyV);
0085       double sigmaToA = sqrt(sigmaJitter1 * sigmaJitter1 + sigmaDistorsion * sigmaDistorsion + sigmaTDC * sigmaTDC +
0086                              sigmaLN * sigmaLN);
0087       double sigmaToC = sqrt(sigmaJitter2 * sigmaJitter2 + sigmaDistorsion * sigmaDistorsion + sigmaTDC * sigmaTDC +
0088                              sigmaLN * sigmaLN);
0089       double smearing1 = 0.0;
0090       double smearing2 = 0.0;
0091
0092       if (sigmaToA > 0. && sigmaToC > 0.) {
0093         smearing1 = CLHEP::RandGaussQ::shoot(hre, 0., sigmaToA);
0094         smearing2 = CLHEP::RandGaussQ::shoot(hre, 0., sigmaToC);
0095       }
0096
0097       finalToA += smearing1;
0098       finalToC += smearing1 + smearing2;
0099
0100       std::array<float, 3> times = etlPulseShape_.timeAtThr(
0101           (it->second).hit_info[0][i] / referenceChargeColl_, iThreshold_MIP_, iThreshold_MIP_);
0102
0103       //The signal is below the threshold
0104       if (times[0] == 0 && times[1] == 0 && times[2] == 0) {
0105         continue;
0106       }
0107       //The signal is considered to be below the threshold
0108       finalToA += times[0];
0109       finalToC += times[2];
0110       if (finalToA >= finalToC)
0111         continue;
0112       chargeColl[i + ibucket] += (it->second).hit_info[0][i];
0113
0114       if (toa1[i + ibucket] == 0. || (finalToA - ibucket * bxTime_) < toa1[i + ibucket]) {
0115         toa1[i + ibucket] = finalToA - ibucket * bxTime_;
0116         toa2[i + ibucket] = finalToC - ibucket * bxTime_;
0117       }
0118
0119       tot[i + ibucket] = finalToC - finalToA;
0120     }
0121     // Run the shaper to create a new data frame
0122     ETLDataFrame rawDataFrame(it->first.detid_);
0123     runTrivialShaper(rawDataFrame, chargeColl, toa1, tot, it->first.row_, it->first.column_);
0124     updateOutput(output, rawDataFrame);
0125   }
0126 }
0127
0128 void ETLElectronicsSim::runTrivialShaper(ETLDataFrame& dataFrame,
0129                                          const mtd::MTDSimHitData& chargeColl,
0130                                          const mtd::MTDSimHitData& toa,
0131                                          const mtd::MTDSimHitData& tot,
0132                                          const uint8_t row,
0133                                          const uint8_t col) const {
0134 #ifdef EDM_ML_DEBUG
0135   bool dumpInfo(false);
0136   for (int it = 0; it < (int)(chargeColl.size()); it++) {
0137     if (chargeColl[it] > adcThreshold_MIP_) {
0138       dumpInfo = true;
0139       break;
0140     }
0141   }
0142   if (dumpInfo) {
0143     LogTrace("ETLElectronicsSim") << "[runTrivialShaper]";
0144   }
0145 #endif
0146
0147   //set new ADCs. Notice that we are only interested in the first element of the array for the ETL
0148   for (int it = 0; it < (int)(chargeColl.size()); it++) {
0149     //brute force saturation, maybe could to better with an exponential like saturation
0150     const uint32_t adc = std::min((uint32_t)std::floor(chargeColl[it] / adcLSB_MIP_), adcBitSaturation_);
0151     const uint32_t tdc_time1 = std::min((uint32_t)std::floor(toa[it] / toaLSB_ns_), tdcBitSaturation_);
0152     const uint32_t tdc_time2 = std::min((uint32_t)std::floor(tot[it] / toaLSB_ns_), tdcBitSaturation_);
0153     //If time over threshold is 0 the event is assumed to not pass the threshold
0154     bool thres = true;
0155     if (tdc_time2 == 0 || chargeColl[it] < adcThreshold_MIP_)
0156       thres = false;
0157
0158     ETLSample newSample;
0159     newSample.set(thres, false, tdc_time1, tdc_time2, adc, row, col);
0160
0161     //ETLSample newSample;
0162     dataFrame.setSample(it, newSample);
0163
0164 #ifdef EDM_ML_DEBUG
0165     if (dumpInfo) {
0166       LogTrace("ETLElectronicsSim") << adc << " (" << chargeColl[it] << "/" << adcLSB_MIP_ << ") ";
0167     }
0168 #endif
0169   }
0170
0171 #ifdef EDM_ML_DEBUG
0172   if (dumpInfo) {
0173     std::ostringstream msg;
0174     dataFrame.print(msg);
0175     LogTrace("ETLElectronicsSim") << msg.str();
0176   }
0177 #endif
0178 }
0179
0180 void ETLElectronicsSim::updateOutput(ETLDigiCollection& coll, const ETLDataFrame& rawDataFrame) const {
0181   int itIdx(9);
0182   if (rawDataFrame.size() <= itIdx + 2)
0183     return;
0184
0185   ETLDataFrame dataFrame(rawDataFrame.id());
0186   dataFrame.resize(dfSIZE);
0187   bool putInEvent(false);
0188   for (int it = 0; it < dfSIZE; ++it) {
0189     dataFrame.setSample(it, rawDataFrame[itIdx - 2 + it]);
0190     if (it == 2)
0191       putInEvent = rawDataFrame[itIdx - 2 + it].threshold();
0192   }
0193
0194   if (putInEvent) {
0195     coll.push_back(dataFrame);
0196   }
0197 }