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