SiPhase2Digitizer/plugins/SSDigitizerAlgorithm.cc

0001 #include <iostream>
0002 #include <cmath>
0003
0004 #include "SimTracker/SiPhase2Digitizer/plugins/SSDigitizerAlgorithm.h"
0005
0006 #include "FWCore/Framework/interface/ConsumesCollector.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009
0010 // Geometry
0011 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0012 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0013
0014 #include "CondFormats/DataRecord/interface/SiPhase2OuterTrackerCondDataRecords.h"
0015
0016 using namespace edm;
0017
0018 namespace {
0019   double nFactorial(int n);
0020   double aScalingConstant(int N, int i);
0021 }  // namespace
0022 namespace {
0023   double nFactorial(int n) { return std::tgamma(n + 1); }
0024   double aScalingConstant(int N, int i) {
0025     return std::pow(-1, (double)i) * nFactorial(N) * nFactorial(N + 2) /
0026            (nFactorial(N - i) * nFactorial(N + 2 - i) * nFactorial(i));
0027   }
0028 }  // namespace
0029
0030 void SSDigitizerAlgorithm::init(const edm::EventSetup& es) {
0031   if (use_LorentzAngle_DB_)  // Get Lorentz angle from DB record
0032     siPhase2OTLorentzAngle_ = &es.getData(siPhase2OTLorentzAngleToken_);
0033
0034   if (use_deadmodule_DB_)  // Get Bad Channel (SiStripBadStrip) from DB
0035     badChannelPayload_ = &es.getData(badChannelToken_);
0036
0037   geom_ = &es.getData(geomToken_);
0038 }
0039
0040 SSDigitizerAlgorithm::SSDigitizerAlgorithm(const edm::ParameterSet& conf, edm::ConsumesCollector iC)
0041     : Phase2TrackerDigitizerAlgorithm(conf.getParameter<ParameterSet>("AlgorithmCommon"),
0042                                       conf.getParameter<ParameterSet>("SSDigitizerAlgorithm"),
0043                                       iC),
0044       hitDetectionMode_(conf.getParameter<ParameterSet>("SSDigitizerAlgorithm").getParameter<int>("HitDetectionMode")),
0045       pulseShapeParameters_(conf.getParameter<ParameterSet>("SSDigitizerAlgorithm")
0046                                 .getParameter<std::vector<double> >("PulseShapeParameters")),
0047       deadTime_(conf.getParameter<ParameterSet>("SSDigitizerAlgorithm").getParameter<double>("CBCDeadTime")),
0048       geomToken_(iC.esConsumes()) {
0049   if (use_LorentzAngle_DB_)
0050     siPhase2OTLorentzAngleToken_ = iC.esConsumes();
0051
0052   if (use_deadmodule_DB_) {
0053     std::string badChannelLabel_ = conf.getParameter<ParameterSet>("SSDigitizerAlgorithm")
0054                                        .getUntrackedParameter<std::string>("BadChannelLabel", "");
0055     badChannelToken_ = iC.esConsumes(edm::ESInputTag{"", badChannelLabel_});
0056   }
0057
0058   pixelFlag_ = false;
0059   LogDebug("SSDigitizerAlgorithm ") << "SSDigitizerAlgorithm constructed "
0060                                     << "Configuration parameters:"
0061                                     << "Threshold/Gain = "
0062                                     << "threshold in electron Endcap = " << theThresholdInE_Endcap_
0063                                     << "threshold in electron Barrel = " << theThresholdInE_Barrel_ << " "
0064                                     << theElectronPerADC_ << " " << theAdcFullScale_ << " The delta cut-off is set to "
0065                                     << tMax_ << " pix-inefficiency " << addPixelInefficiency_;
0066   storeSignalShape();
0067 }
0068 SSDigitizerAlgorithm::~SSDigitizerAlgorithm() { LogDebug("SSDigitizerAlgorithm") << "SSDigitizerAlgorithm deleted"; }
0069 //
0070 // -- Select the Hit for Digitization
0071 //
0072 bool SSDigitizerAlgorithm::select_hit(const PSimHit& hit, double tCorr, double& sigScale) const {
0073   bool result = false;
0074   if (hitDetectionMode_ == SSDigitizerAlgorithm::SampledMode)
0075     result = select_hit_sampledMode(hit, tCorr, sigScale);
0076   else if (hitDetectionMode_ == SSDigitizerAlgorithm::LatchedMode)
0077     result = select_hit_latchedMode(hit, tCorr, sigScale);
0078   else {
0079     double toa = hit.tof() - tCorr;
0080     result = (toa > theTofLowerCut_ && toa < theTofUpperCut_);
0081   }
0082   return result;
0083 }
0084 //
0085 // -- Select Hits in Sampled Mode
0086 //
0087 bool SSDigitizerAlgorithm::select_hit_sampledMode(const PSimHit& hit, double tCorr, double& sigScale) const {
0088   double toa = hit.tof() - tCorr;
0089   double sampling_time = bx_time;
0090
0091   DetId det_id = DetId(hit.detUnitId());
0092   float theThresholdInE =
0093       (det_id.subdetId() == StripSubdetector::TOB) ? theThresholdInE_Barrel_ : theThresholdInE_Endcap_;
0094
0095   sigScale = getSignalScale(sampling_time - toa);
0096   return (sigScale * hit.energyLoss() / GeVperElectron_ > theThresholdInE);
0097 }
0098 //
0099 // -- Select Hits in Hit Detection Mode
0100 //
0101 bool SSDigitizerAlgorithm::select_hit_latchedMode(const PSimHit& hit, double tCorr, double& sigScale) const {
0102   float toa = hit.tof() - tCorr;
0103   toa -= hit.eventId().bunchCrossing() * bx_time;
0104
0105   float sampling_time = (-1) * (hit.eventId().bunchCrossing() + 1) * bx_time;
0106
0107   DetId det_id = DetId(hit.detUnitId());
0108   float theThresholdInE =
0109       (det_id.subdetId() == StripSubdetector::TOB) ? theThresholdInE_Barrel_ : theThresholdInE_Endcap_;
0110
0111   bool lastPulse = true;
0112   bool aboveThr = false;
0113   for (float i = deadTime_; i <= bx_time; i++) {
0114     sigScale = getSignalScale(sampling_time - toa + i);
0115
0116     aboveThr = (sigScale * hit.energyLoss() / GeVperElectron_ > theThresholdInE);
0117     if (!lastPulse && aboveThr)
0118       return true;
0119
0120     lastPulse = aboveThr;
0121   }
0122   return false;
0123 }
0124 double SSDigitizerAlgorithm::cbc3PulsePolarExpansion(double x) const {
0125   constexpr size_t max_par = 6;
0126   if (pulseShapeParameters_.size() < max_par)
0127     return -1;
0128   double xOffset = pulseShapeParameters_[0];
0129   double tau = pulseShapeParameters_[1];
0130   double r = pulseShapeParameters_[2];
0131   double theta = pulseShapeParameters_[3];
0132   int nTerms = static_cast<int>(pulseShapeParameters_[4]);
0133
0134   double fN = 0;
0135   double xx = x - xOffset;
0136   if (xx < 0)
0137     return 0;
0138
0139   for (int i = 0; i < nTerms; i++) {
0140     double angularTerm = 0;
0141     double temporalTerm = 0;
0142     double rTerm = std::pow(r, i) / (std::pow(tau, 2. * i) * ::nFactorial(i + 2));
0143     for (int j = 0; j <= i; j++) {
0144       angularTerm += std::pow(std::cos(theta), (double)(i - j)) * std::pow(std::sin(theta), (double)j);
0145       temporalTerm += ::aScalingConstant(i, j) * std::pow(xx, (double)(i - j)) * std::pow(tau, (double)j);
0146     }
0147     double fi = rTerm * angularTerm * temporalTerm;
0148
0149     fN += fi;
0150   }
0151   return fN;
0152 }
0153 double SSDigitizerAlgorithm::signalShape(double x) const {
0154   double xOffset = pulseShapeParameters_[0];
0155   double tau = pulseShapeParameters_[1];
0156   double maxCharge = pulseShapeParameters_[5];
0157
0158   double xx = x - xOffset;
0159   return maxCharge * (std::exp(-xx / tau) * std::pow(xx / tau, 2.) * cbc3PulsePolarExpansion(x));
0160 }
0161 void SSDigitizerAlgorithm::storeSignalShape() {
0162   for (size_t i = 0; i < interpolationPoints; i++) {
0163     float val = i / interpolationStep;
0164
0165     pulseShapeVec_.push_back(signalShape(val));
0166   }
0167 }
0168 double SSDigitizerAlgorithm::getSignalScale(double xval) const {
0169   double res = 0.0;
0170   int len = pulseShapeVec_.size();
0171
0172   if (xval < 0.0 || xval * interpolationStep >= len)
0173     return res;
0174
0175   unsigned int lower = std::floor(xval) * interpolationStep;
0176   unsigned int upper = std::ceil(xval) * interpolationStep;
0177   for (size_t i = lower + 1; i < upper * interpolationStep; i++) {
0178     float val = i * 0.1;
0179     if (val > xval) {
0180       res = pulseShapeVec_[i - 1];
0181       break;
0182     }
0183   }
0184   return res;
0185 }
0186 //
0187 // -- Compare Signal with Threshold
0188 //
0189 bool SSDigitizerAlgorithm::isAboveThreshold(const digitizerUtility::SimHitInfo* const hisInfo,
0190                                             float charge,
0191                                             float thr) const {
0192   return (charge >= thr);
0193 }
0194 //
0195 // -- Read Bad Channels from the Condidion DB and kill channels/module accordingly
0196 //
0197 void SSDigitizerAlgorithm::module_killing_DB(const Phase2TrackerGeomDetUnit* pixdet) {
0198   uint32_t detId = pixdet->geographicalId().rawId();
0199
0200   signal_map_type& theSignal = _signal[detId];
0201   signal_map_type signalNew;
0202
0203   SiStripBadStrip::Range range = badChannelPayload_->getRange(detId);
0204   for (std::vector<unsigned int>::const_iterator badChannel = range.first; badChannel != range.second; ++badChannel) {
0205     const auto& firstStrip = badChannelPayload_->decodePhase2(*badChannel).firstStrip;
0206     const auto& channelRange = badChannelPayload_->decodePhase2(*badChannel).range;
0207
0208     for (int index = 0; index < channelRange; index++) {
0209       for (auto& s : theSignal) {
0210         auto& channel = s.first;
0211         if (channel == firstStrip + index)
0212           s.second.set(0.);
0213       }
0214     }
0215   }
0216 }