Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:29:46

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