SiPhase2Digitizer/plugins/PixelDigitizerAlgorithm.cc

0001 #include <iostream>
0002 #include <cmath>
0003
0004 #include "SimTracker/SiPhase2Digitizer/plugins/PixelDigitizerAlgorithm.h"
0005 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0006
0007 #include "FWCore/Framework/interface/ConsumesCollector.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "CalibTracker/SiPixelESProducers/interface/SiPixelGainCalibrationOfflineSimService.h"
0011
0012 // Geometry
0013 #include "CondFormats/SiPixelObjects/interface/GlobalPixel.h"
0014 #include "CondFormats/SiPixelObjects/interface/SiPixelLorentzAngle.h"
0015 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0016 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0017 #include "CondFormats/SiPixelObjects/interface/SiPixelQuality.h"
0018 #include "CondFormats/SiPixelObjects/interface/PixelROC.h"
0019 #include "CondFormats/SiPixelObjects/interface/LocalPixel.h"
0020 #include "CondFormats/SiPixelObjects/interface/CablingPathToDetUnit.h"
0021
0022 using namespace edm;
0023 using namespace sipixelobjects;
0024
0025 void PixelDigitizerAlgorithm::init(const edm::EventSetup& es) {
0026   if (use_ineff_from_db_)  // load gain calibration service fromdb...
0027     theSiPixelGainCalibrationService_->setESObjects(es);
0028
0029   if (use_deadmodule_DB_)
0030     siPixelBadModule_ = &es.getData(siPixelBadModuleToken_);
0031
0032   if (use_LorentzAngle_DB_)  // Get Lorentz angle from DB record
0033     siPixelLorentzAngle_ = &es.getData(siPixelLorentzAngleToken_);
0034
0035   // gets the map and geometry from the DB (to kill ROCs)
0036   fedCablingMap_ = &es.getData(fedCablingMapToken_);
0037   geom_ = &es.getData(geomToken_);
0038 }
0039
0040 PixelDigitizerAlgorithm::PixelDigitizerAlgorithm(const edm::ParameterSet& conf, edm::ConsumesCollector iC)
0041     : Phase2TrackerDigitizerAlgorithm(conf.getParameter<ParameterSet>("AlgorithmCommon"),
0042                                       conf.getParameter<ParameterSet>("PixelDigitizerAlgorithm"),
0043                                       iC),
0044       odd_row_interchannelCoupling_next_row_(conf.getParameter<ParameterSet>("PixelDigitizerAlgorithm")
0045                                                  .getParameter<double>("Odd_row_interchannelCoupling_next_row")),
0046       even_row_interchannelCoupling_next_row_(conf.getParameter<ParameterSet>("PixelDigitizerAlgorithm")
0047                                                   .getParameter<double>("Even_row_interchannelCoupling_next_row")),
0048       odd_column_interchannelCoupling_next_column_(
0049           conf.getParameter<ParameterSet>("PixelDigitizerAlgorithm")
0050               .getParameter<double>("Odd_column_interchannelCoupling_next_column")),
0051       even_column_interchannelCoupling_next_column_(
0052           conf.getParameter<ParameterSet>("PixelDigitizerAlgorithm")
0053               .getParameter<double>("Even_column_interchannelCoupling_next_column")),
0054       apply_timewalk_(conf.getParameter<ParameterSet>("PixelDigitizerAlgorithm").getParameter<bool>("ApplyTimewalk")),
0055       timewalk_model_(
0056           conf.getParameter<ParameterSet>("PixelDigitizerAlgorithm").getParameter<edm::ParameterSet>("TimewalkModel")),
0057       fedCablingMapToken_(iC.esConsumes()),
0058       geomToken_(iC.esConsumes()) {
0059   if (use_deadmodule_DB_)
0060     siPixelBadModuleToken_ = iC.esConsumes();
0061   if (use_LorentzAngle_DB_)
0062     siPixelLorentzAngleToken_ = iC.esConsumes();
0063   pixelFlag_ = true;
0064   LogDebug("PixelDigitizerAlgorithm") << "Algorithm constructed "
0065                                       << "Configuration parameters:"
0066                                       << "Threshold/Gain = "
0067                                       << "threshold in electron Endcap = " << theThresholdInE_Endcap_
0068                                       << "threshold in electron Barrel = " << theThresholdInE_Barrel_ << " "
0069                                       << theElectronPerADC_ << " " << theAdcFullScale_
0070                                       << " The delta cut-off is set to " << tMax_ << " pix-inefficiency "
0071                                       << addPixelInefficiency_;
0072 }
0073 PixelDigitizerAlgorithm::~PixelDigitizerAlgorithm() { LogDebug("PixelDigitizerAlgorithm") << "Algorithm deleted"; }
0074 //
0075 // -- Select the Hit for Digitization
0076 //
0077 bool PixelDigitizerAlgorithm::select_hit(const PSimHit& hit, double tCorr, double& sigScale) const {
0078   double time = hit.tof() - tCorr;
0079   return (time >= theTofLowerCut_ && time < theTofUpperCut_);
0080 }
0081
0082 // ======================================================================
0083 //
0084 //  Add  Cross-talk contribution
0085 //
0086 // ======================================================================
0087 void PixelDigitizerAlgorithm::add_cross_talk(const Phase2TrackerGeomDetUnit* pixdet) {
0088   if (!pixelFlag_)
0089     return;
0090
0091   const Phase2TrackerTopology* topol = &pixdet->specificTopology();
0092
0093   // cross-talk calculation valid for the case of 25x100 pixels
0094   const float pitch_first = 0.0025;
0095   const float pitch_second = 0.0100;
0096
0097   // 0.5 um tolerance when comparing the pitch to accommodate the small changes in different TK geometrie (temporary fix)
0098   const double pitch_tolerance(0.0005);
0099
0100   if (std::abs(topol->pitch().first - pitch_first) > pitch_tolerance ||
0101       std::abs(topol->pitch().second - pitch_second) > pitch_tolerance)
0102     return;
0103
0104   uint32_t detID = pixdet->geographicalId().rawId();
0105   signal_map_type& theSignal = _signal[detID];
0106   signal_map_type signalNew;
0107
0108   int numRows = topol->nrows();
0109   int numColumns = topol->ncolumns();
0110
0111   for (auto& s : theSignal) {
0112     float signalInElectrons = s.second.ampl();  // signal in electrons
0113
0114     auto hitChan = PixelDigi::channelToPixel(s.first);
0115
0116     float signalInElectrons_odd_row_Xtalk_next_row = signalInElectrons * odd_row_interchannelCoupling_next_row_;
0117     float signalInElectrons_even_row_Xtalk_next_row = signalInElectrons * even_row_interchannelCoupling_next_row_;
0118     float signalInElectrons_odd_column_Xtalk_next_column =
0119         signalInElectrons * odd_column_interchannelCoupling_next_column_;
0120     float signalInElectrons_even_column_Xtalk_next_column =
0121         signalInElectrons * even_column_interchannelCoupling_next_column_;
0122
0123     // subtract the charge which will be shared
0124     s.second.set(signalInElectrons - signalInElectrons_odd_row_Xtalk_next_row -
0125                  signalInElectrons_even_row_Xtalk_next_row - signalInElectrons_odd_column_Xtalk_next_column -
0126                  signalInElectrons_even_column_Xtalk_next_column);
0127
0128     if (hitChan.first != 0) {
0129       auto XtalkPrev = std::make_pair(hitChan.first - 1, hitChan.second);
0130       int chanXtalkPrev = pixelFlag_ ? PixelDigi::pixelToChannel(XtalkPrev.first, XtalkPrev.second)
0131                                      : Phase2TrackerDigi::pixelToChannel(XtalkPrev.first, XtalkPrev.second);
0132       if (hitChan.first % 2 == 1)
0133         signalNew.emplace(chanXtalkPrev,
0134                           DigitizerUtility::Amplitude(signalInElectrons_even_row_Xtalk_next_row, nullptr, -1.0));
0135       else
0136         signalNew.emplace(chanXtalkPrev,
0137                           DigitizerUtility::Amplitude(signalInElectrons_odd_row_Xtalk_next_row, nullptr, -1.0));
0138     }
0139     if (hitChan.first < numRows - 1) {
0140       auto XtalkNext = std::make_pair(hitChan.first + 1, hitChan.second);
0141       int chanXtalkNext = pixelFlag_ ? PixelDigi::pixelToChannel(XtalkNext.first, XtalkNext.second)
0142                                      : Phase2TrackerDigi::pixelToChannel(XtalkNext.first, XtalkNext.second);
0143       if (hitChan.first % 2 == 1)
0144         signalNew.emplace(chanXtalkNext,
0145                           DigitizerUtility::Amplitude(signalInElectrons_odd_row_Xtalk_next_row, nullptr, -1.0));
0146       else
0147         signalNew.emplace(chanXtalkNext,
0148                           DigitizerUtility::Amplitude(signalInElectrons_even_row_Xtalk_next_row, nullptr, -1.0));
0149     }
0150
0151     if (hitChan.second != 0) {
0152       auto XtalkPrev = std::make_pair(hitChan.first, hitChan.second - 1);
0153       int chanXtalkPrev = pixelFlag_ ? PixelDigi::pixelToChannel(XtalkPrev.first, XtalkPrev.second)
0154                                      : Phase2TrackerDigi::pixelToChannel(XtalkPrev.first, XtalkPrev.second);
0155       if (hitChan.second % 2 == 1)
0156         signalNew.emplace(chanXtalkPrev,
0157                           DigitizerUtility::Amplitude(signalInElectrons_even_column_Xtalk_next_column, nullptr, -1.0));
0158       else
0159         signalNew.emplace(chanXtalkPrev,
0160                           DigitizerUtility::Amplitude(signalInElectrons_odd_column_Xtalk_next_column, nullptr, -1.0));
0161     }
0162     if (hitChan.second < numColumns - 1) {
0163       auto XtalkNext = std::make_pair(hitChan.first, hitChan.second + 1);
0164       int chanXtalkNext = pixelFlag_ ? PixelDigi::pixelToChannel(XtalkNext.first, XtalkNext.second)
0165                                      : Phase2TrackerDigi::pixelToChannel(XtalkNext.first, XtalkNext.second);
0166       if (hitChan.second % 2 == 1)
0167         signalNew.emplace(chanXtalkNext,
0168                           DigitizerUtility::Amplitude(signalInElectrons_odd_column_Xtalk_next_column, nullptr, -1.0));
0169       else
0170         signalNew.emplace(chanXtalkNext,
0171                           DigitizerUtility::Amplitude(signalInElectrons_even_column_Xtalk_next_column, nullptr, -1.0));
0172     }
0173   }
0174   for (auto const& l : signalNew) {
0175     int chan = l.first;
0176     auto iter = theSignal.find(chan);
0177     if (iter != theSignal.end()) {
0178       iter->second += l.second.ampl();
0179     } else {
0180       theSignal.emplace(chan, DigitizerUtility::Amplitude(l.second.ampl(), nullptr, -1.0));
0181     }
0182   }
0183 }
0184
0185 PixelDigitizerAlgorithm::TimewalkCurve::TimewalkCurve(const edm::ParameterSet& pset)
0186     : x_(pset.getParameter<std::vector<double>>("charge")), y_(pset.getParameter<std::vector<double>>("delay")) {
0187   if (x_.size() != y_.size())
0188     throw cms::Exception("Configuration")
0189         << "Timewalk model error: the number of charge values does not match the number of delay values!";
0190 }
0191
0192 double PixelDigitizerAlgorithm::TimewalkCurve::operator()(double x) const {
0193   auto it = std::lower_bound(x_.begin(), x_.end(), x);
0194   if (it == x_.begin())
0195     return y_.front();
0196   if (it == x_.end())
0197     return y_.back();
0198   int index = std::distance(x_.begin(), it);
0199   double x_high = *it;
0200   double x_low = *(--it);
0201   double p = (x - x_low) / (x_high - x_low);
0202   return p * y_[index] + (1 - p) * y_[index - 1];
0203 }
0204
0205 PixelDigitizerAlgorithm::TimewalkModel::TimewalkModel(const edm::ParameterSet& pset) {
0206   threshold_values = pset.getParameter<std::vector<double>>("ThresholdValues");
0207   const auto& curve_psetvec = pset.getParameter<std::vector<edm::ParameterSet>>("Curves");
0208   if (threshold_values.size() != curve_psetvec.size())
0209     throw cms::Exception("Configuration")
0210         << "Timewalk model error: the number of threshold values does not match the number of curves.";
0211   for (const auto& curve_pset : curve_psetvec)
0212     curves.emplace_back(curve_pset);
0213 }
0214
0215 double PixelDigitizerAlgorithm::TimewalkModel::operator()(double q_in, double q_threshold) const {
0216   auto index = find_closest_index(threshold_values, q_threshold);
0217   return curves[index](q_in);
0218 }
0219
0220 std::size_t PixelDigitizerAlgorithm::TimewalkModel::find_closest_index(const std::vector<double>& vec,
0221                                                                        double value) const {
0222   auto it = std::lower_bound(vec.begin(), vec.end(), value);
0223   if (it == vec.begin())
0224     return 0;
0225   else if (it == vec.end())
0226     return vec.size() - 1;
0227   else {
0228     auto it_upper = it;
0229     auto it_lower = --it;
0230     auto closest = (value - *it_lower > *it_upper - value) ? it_upper : it_lower;
0231     return std::distance(vec.begin(), closest);
0232   }
0233 }
0234 //
0235 // -- Compare Signal with Threshold
0236 //
0237 bool PixelDigitizerAlgorithm::isAboveThreshold(const DigitizerUtility::SimHitInfo* hitInfo,
0238                                                float charge,
0239                                                float thr) const {
0240   if (charge < thr)
0241     return false;
0242   if (apply_timewalk_ && hitInfo) {
0243     float corrected_time = hitInfo->time();
0244     double time = corrected_time + timewalk_model_(charge, thr);
0245     return (time >= theTofLowerCut_ && time < theTofUpperCut_);
0246   } else
0247     return true;
0248 }
0249 //
0250 // -- Read Bad Channels from the Condidion DB and kill channels/module accordingly
0251 //
0252 void PixelDigitizerAlgorithm::module_killing_DB(const Phase2TrackerGeomDetUnit* pixdet) {
0253   bool isbad = false;
0254   uint32_t detID = pixdet->geographicalId().rawId();
0255   int ncol = pixdet->specificTopology().ncolumns();
0256   if (ncol < 0)
0257     return;
0258   std::vector<SiPixelQuality::disabledModuleType> disabledModules = siPixelBadModule_->getBadComponentList();
0259
0260   SiPixelQuality::disabledModuleType badmodule;
0261   for (const auto& mod : disabledModules) {
0262     if (detID == mod.DetID) {
0263       isbad = true;
0264       badmodule = mod;
0265       break;
0266     }
0267   }
0268
0269   if (!isbad)
0270     return;
0271
0272   signal_map_type& theSignal = _signal[detID];  // check validity
0273   if (badmodule.errorType == 0) {               // this is a whole dead module.
0274     for (auto& s : theSignal)
0275       s.second.set(0.);  // reset amplitude
0276   } else {               // all other module types: half-modules and single ROCs.
0277     // Get Bad ROC position:
0278     // follow the example of getBadRocPositions in CondFormats/SiPixelObjects/src/SiPixelQuality.cc
0279     std::vector<GlobalPixel> badrocpositions;
0280     for (size_t j = 0; j < static_cast<size_t>(ncol); j++) {
0281       if (siPixelBadModule_->IsRocBad(detID, j)) {
0282         std::vector<CablingPathToDetUnit> path = fedCablingMap_->pathToDetUnit(detID);
0283         for (auto const& p : path) {
0284           const PixelROC* myroc = fedCablingMap_->findItem(p);
0285           if (myroc->idInDetUnit() == j) {
0286             LocalPixel::RocRowCol local = {39, 25};  //corresponding to center of ROC row, col
0287             GlobalPixel global = myroc->toGlobal(LocalPixel(local));
0288             badrocpositions.push_back(global);
0289             break;
0290           }
0291         }
0292       }
0293     }
0294
0295     for (auto& s : theSignal) {
0296       std::pair<int, int> ip;
0297       if (pixelFlag_)
0298         ip = PixelDigi::channelToPixel(s.first);
0299       else
0300         ip = Phase2TrackerDigi::channelToPixel(s.first);
0301
0302       for (auto const& p : badrocpositions) {
0303         for (auto& k : badPixels_) {
0304           if (p.row == k.getParameter<int>("row") && ip.first == k.getParameter<int>("row") &&
0305               std::abs(ip.second - p.col) < k.getParameter<int>("col")) {
0306             s.second.set(0.);
0307           }
0308         }
0309       }
0310     }
0311   }
0312 }