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       debug_(pset.getUntrackedParameter<bool>("debug", false)),
0013       bxTime_(pset.getParameter<double>("bxTime")),
0014       integratedLum_(pset.getParameter<double>("IntegratedLuminosity")),
0015       fluence_(pset.getParameter<std::string>("FluenceVsRadius")),
0016       lgadGain_(pset.getParameter<std::string>("LGADGainVsFluence")),
0017       timeRes2_(pset.getParameter<std::string>("TimeResolution2")),
0018       adcNbits_(pset.getParameter<uint32_t>("adcNbits")),
0019       tdcNbits_(pset.getParameter<uint32_t>("tdcNbits")),
0020       adcSaturation_MIP_(pset.getParameter<double>("adcSaturation_MIP")),
0021       adcLSB_MIP_(adcSaturation_MIP_ / std::pow(2., adcNbits_)),
0022       adcBitSaturation_(std::pow(2, adcNbits_) - 1),
0023       adcThreshold_MIP_(pset.getParameter<double>("adcThreshold_MIP")),
0024       toaLSB_ns_(pset.getParameter<double>("toaLSB_ns")),
0025       tdcBitSaturation_(std::pow(2, tdcNbits_) - 1) {}
0026
0027 void ETLElectronicsSim::getEventSetup(const edm::EventSetup& evs) { geom_ = &evs.getData(geomToken_); }
0028
0029 void ETLElectronicsSim::run(const mtd::MTDSimHitDataAccumulator& input,
0030                             ETLDigiCollection& output,
0031                             CLHEP::HepRandomEngine* hre) const {
0032   MTDSimHitData chargeColl, toa;
0033
0034   std::vector<double> emptyV;
0035   std::vector<double> radius(1);
0036   std::vector<double> fluence(1);
0037   std::vector<double> gain(1);
0038
0039   for (MTDSimHitDataAccumulator::const_iterator it = input.begin(); it != input.end(); it++) {
0040     chargeColl.fill(0.f);
0041     toa.fill(0.f);
0042
0043     ETLDetId detId = it->first.detid_;
0044     DetId geoId = detId.geographicalId();
0045     const MTDGeomDet* thedet = geom_->idToDet(geoId);
0046     if (thedet == nullptr)
0047       throw cms::Exception("EtlElectronicsSim") << "GeographicalID: " << std::hex << geoId.rawId() << " ("
0048                                                 << detId.rawId() << ") is invalid!" << std::dec << std::endl;
0049     const PixelTopology& topo = static_cast<const PixelTopology&>(thedet->topology());
0050
0051     Local3DPoint local_point(topo.localX(it->first.row_), topo.localY(it->first.column_), 0.);
0052     const auto& global_point = thedet->toGlobal(local_point);
0053
0054     for (size_t i = 0; i < it->second.hit_info[0].size(); i++) {
0055       if ((it->second).hit_info[0][i] < adcThreshold_MIP_)
0056         continue;
0057
0058       // time of arrival
0059       float finalToA = (it->second).hit_info[1][i];
0060
0061       // calculate the LGAD gain as a function of the fluence at R = radius
0062       radius[0] = global_point.perp();
0063       fluence[0] = integratedLum_ * fluence_.evaluate(radius, emptyV);
0064       gain[0] = lgadGain_.evaluate(fluence, emptyV);
0065       if (gain[0] <= 0.)
0066         throw cms::Exception("EtlElectronicsSim") << "Null or negative LGAD gain!" << std::endl;
0067
0068       // Gaussian smearing of the time of arrival
0069       double sigmaToA = sqrt(timeRes2_.evaluate(gain, emptyV));
0070
0071       if (sigmaToA > 0.)
0072         finalToA += CLHEP::RandGaussQ::shoot(hre, 0., sigmaToA);
0073
0074       // fill the time and charge arrays
0075       const unsigned int ibucket = std::floor(finalToA / bxTime_);
0076       if ((i + ibucket) >= chargeColl.size())
0077         continue;
0078
0079       chargeColl[i + ibucket] += (it->second).hit_info[0][i];
0080
0081       if (toa[i + ibucket] == 0. || (finalToA - ibucket * bxTime_) < toa[i + ibucket])
0082         toa[i + ibucket] = finalToA - ibucket * bxTime_;
0083     }
0084
0085     // run the shaper to create a new data frame
0086     ETLDataFrame rawDataFrame(it->first.detid_);
0087     runTrivialShaper(rawDataFrame, chargeColl, toa, it->first.row_, it->first.column_);
0088     updateOutput(output, rawDataFrame);
0089   }
0090 }
0091
0092 void ETLElectronicsSim::runTrivialShaper(ETLDataFrame& dataFrame,
0093                                          const mtd::MTDSimHitData& chargeColl,
0094                                          const mtd::MTDSimHitData& toa,
0095                                          const uint8_t row,
0096                                          const uint8_t col) const {
0097   bool debug = debug_;
0098 #ifdef EDM_ML_DEBUG
0099   for (int it = 0; it < (int)(chargeColl.size()); it++)
0100     debug |= (chargeColl[it] > adcThreshold_MIP_);
0101 #endif
0102
0103   if (debug)
0104     edm::LogVerbatim("ETLElectronicsSim") << "[runTrivialShaper]" << std::endl;
0105
0106   //set new ADCs
0107   for (int it = 0; it < (int)(chargeColl.size()); it++) {
0108     //brute force saturation, maybe could to better with an exponential like saturation
0109     const uint32_t adc = std::min((uint32_t)std::floor(chargeColl[it] / adcLSB_MIP_), adcBitSaturation_);
0110     const uint32_t tdc_time = std::min((uint32_t)std::floor(toa[it] / toaLSB_ns_), tdcBitSaturation_);
0111     ETLSample newSample;
0112     newSample.set(chargeColl[it] > adcThreshold_MIP_, false, tdc_time, adc, row, col);
0113     dataFrame.setSample(it, newSample);
0114
0115     if (debug)
0116       edm::LogVerbatim("ETLElectronicsSim") << adc << " (" << chargeColl[it] << "/" << adcLSB_MIP_ << ") ";
0117   }
0118
0119   if (debug) {
0120     std::ostringstream msg;
0121     dataFrame.print(msg);
0122     edm::LogVerbatim("ETLElectronicsSim") << msg.str() << std::endl;
0123   }
0124 }
0125
0126 void ETLElectronicsSim::updateOutput(ETLDigiCollection& coll, const ETLDataFrame& rawDataFrame) const {
0127   int itIdx(9);
0128   if (rawDataFrame.size() <= itIdx + 2)
0129     return;
0130
0131   ETLDataFrame dataFrame(rawDataFrame.id());
0132   dataFrame.resize(dfSIZE);
0133   bool putInEvent(false);
0134   for (int it = 0; it < dfSIZE; ++it) {
0135     dataFrame.setSample(it, rawDataFrame[itIdx - 2 + it]);
0136     if (it == 2)
0137       putInEvent = rawDataFrame[itIdx - 2 + it].threshold();
0138   }
0139
0140   if (putInEvent) {
0141     coll.push_back(dataFrame);
0142   }
0143 }