File indexing completed on 2023-01-16 23:37:27
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
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_)
0027 theSiPixelGainCalibrationService_->setESObjects(es);
0028
0029 if (use_deadmodule_DB_)
0030 siPixelBadModule_ = &es.getData(siPixelBadModuleToken_);
0031
0032 if (use_LorentzAngle_DB_)
0033 siPixelLorentzAngle_ = &es.getData(siPixelLorentzAngleToken_);
0034
0035
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
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
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
0099 const float pitch_first = 0.0025;
0100 const float pitch_second = 0.0100;
0101
0102
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();
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
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
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
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];
0282 if (badmodule.errorType == 0) {
0283 for (auto& s : theSignal)
0284 s.second.set(0.);
0285 } else {
0286
0287
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};
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 }