File indexing completed on 2024-04-06 12:30:55
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
0021 using namespace edm;
0022 using namespace sipixelobjects;
0023
0024 void PixelDigitizerAlgorithm::init(const edm::EventSetup& es) {
0025 if (use_ineff_from_db_)
0026 theSiPixelGainCalibrationService_->setESObjects(es);
0027
0028 if (use_deadmodule_DB_)
0029 siPixelBadModule_ = &es.getData(siPixelBadModuleToken_);
0030
0031 if (use_LorentzAngle_DB_)
0032 siPixelLorentzAngle_ = &es.getData(siPixelLorentzAngleToken_);
0033
0034 geom_ = &es.getData(geomToken_);
0035 if (useChargeReweighting_) {
0036 theSiPixelChargeReweightingAlgorithm_->init(es);
0037 }
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 geomToken_(iC.esConsumes()) {
0058 if (use_deadmodule_DB_)
0059 siPixelBadModuleToken_ = iC.esConsumes();
0060 if (use_LorentzAngle_DB_)
0061 siPixelLorentzAngleToken_ = iC.esConsumes();
0062 pixelFlag_ = true;
0063 LogDebug("PixelDigitizerAlgorithm") << "Algorithm constructed "
0064 << "Configuration parameters:"
0065 << "Threshold/Gain = "
0066 << "threshold in electron Endcap = " << theThresholdInE_Endcap_
0067 << "threshold in electron Barrel = " << theThresholdInE_Barrel_ << " "
0068 << theElectronPerADC_ << " " << theAdcFullScale_
0069 << " The delta cut-off is set to " << tMax_ << " pix-inefficiency "
0070 << addPixelInefficiency_;
0071 }
0072
0073 PixelDigitizerAlgorithm::~PixelDigitizerAlgorithm() { LogDebug("PixelDigitizerAlgorithm") << "Algorithm deleted"; }
0074
0075
0076
0077
0078 bool PixelDigitizerAlgorithm::select_hit(const PSimHit& hit, double tCorr, double& sigScale) const {
0079
0080 double toa = hit.tof() - tCorr;
0081 return apply_timewalk_ || (toa >= theTofLowerCut_ && toa < theTofUpperCut_);
0082 }
0083
0084
0085
0086
0087
0088
0089 void PixelDigitizerAlgorithm::add_cross_talk(const Phase2TrackerGeomDetUnit* pixdet) {
0090 if (!pixelFlag_)
0091 return;
0092
0093 const Phase2TrackerTopology* topol = &pixdet->specificTopology();
0094
0095
0096 const float pitch_first = 0.0025;
0097 const float pitch_second = 0.0100;
0098
0099
0100 const double pitch_tolerance(0.0005);
0101
0102 if (std::abs(topol->pitch().first - pitch_first) > pitch_tolerance ||
0103 std::abs(topol->pitch().second - pitch_second) > pitch_tolerance)
0104 return;
0105
0106 uint32_t detID = pixdet->geographicalId().rawId();
0107 signal_map_type& theSignal = _signal[detID];
0108 signal_map_type signalNew;
0109
0110 int numRows = topol->nrows();
0111 int numColumns = topol->ncolumns();
0112
0113 for (auto& s : theSignal) {
0114 float signalInElectrons = s.second.ampl();
0115
0116 auto hitChan = PixelDigi::channelToPixel(s.first);
0117
0118 float signalInElectrons_odd_row_Xtalk_next_row = signalInElectrons * odd_row_interchannelCoupling_next_row_;
0119 float signalInElectrons_even_row_Xtalk_next_row = signalInElectrons * even_row_interchannelCoupling_next_row_;
0120 float signalInElectrons_odd_column_Xtalk_next_column =
0121 signalInElectrons * odd_column_interchannelCoupling_next_column_;
0122 float signalInElectrons_even_column_Xtalk_next_column =
0123 signalInElectrons * even_column_interchannelCoupling_next_column_;
0124
0125
0126 s.second.set(signalInElectrons - signalInElectrons_odd_row_Xtalk_next_row -
0127 signalInElectrons_even_row_Xtalk_next_row - signalInElectrons_odd_column_Xtalk_next_column -
0128 signalInElectrons_even_column_Xtalk_next_column);
0129
0130 if (hitChan.first != 0) {
0131 auto XtalkPrev = std::make_pair(hitChan.first - 1, hitChan.second);
0132 int chanXtalkPrev = pixelFlag_ ? PixelDigi::pixelToChannel(XtalkPrev.first, XtalkPrev.second)
0133 : Phase2TrackerDigi::pixelToChannel(XtalkPrev.first, XtalkPrev.second);
0134 if (hitChan.first % 2 == 1)
0135 signalNew.emplace(chanXtalkPrev,
0136 digitizerUtility::Ph2Amplitude(signalInElectrons_even_row_Xtalk_next_row, nullptr, -1.0));
0137 else
0138 signalNew.emplace(chanXtalkPrev,
0139 digitizerUtility::Ph2Amplitude(signalInElectrons_odd_row_Xtalk_next_row, nullptr, -1.0));
0140 }
0141 if (hitChan.first < numRows - 1) {
0142 auto XtalkNext = std::make_pair(hitChan.first + 1, hitChan.second);
0143 int chanXtalkNext = pixelFlag_ ? PixelDigi::pixelToChannel(XtalkNext.first, XtalkNext.second)
0144 : Phase2TrackerDigi::pixelToChannel(XtalkNext.first, XtalkNext.second);
0145 if (hitChan.first % 2 == 1)
0146 signalNew.emplace(chanXtalkNext,
0147 digitizerUtility::Ph2Amplitude(signalInElectrons_odd_row_Xtalk_next_row, nullptr, -1.0));
0148 else
0149 signalNew.emplace(chanXtalkNext,
0150 digitizerUtility::Ph2Amplitude(signalInElectrons_even_row_Xtalk_next_row, nullptr, -1.0));
0151 }
0152
0153 if (hitChan.second != 0) {
0154 auto XtalkPrev = std::make_pair(hitChan.first, hitChan.second - 1);
0155 int chanXtalkPrev = pixelFlag_ ? PixelDigi::pixelToChannel(XtalkPrev.first, XtalkPrev.second)
0156 : Phase2TrackerDigi::pixelToChannel(XtalkPrev.first, XtalkPrev.second);
0157 if (hitChan.second % 2 == 1)
0158 signalNew.emplace(
0159 chanXtalkPrev,
0160 digitizerUtility::Ph2Amplitude(signalInElectrons_even_column_Xtalk_next_column, nullptr, -1.0));
0161 else
0162 signalNew.emplace(
0163 chanXtalkPrev,
0164 digitizerUtility::Ph2Amplitude(signalInElectrons_odd_column_Xtalk_next_column, nullptr, -1.0));
0165 }
0166 if (hitChan.second < numColumns - 1) {
0167 auto XtalkNext = std::make_pair(hitChan.first, hitChan.second + 1);
0168 int chanXtalkNext = pixelFlag_ ? PixelDigi::pixelToChannel(XtalkNext.first, XtalkNext.second)
0169 : Phase2TrackerDigi::pixelToChannel(XtalkNext.first, XtalkNext.second);
0170 if (hitChan.second % 2 == 1)
0171 signalNew.emplace(
0172 chanXtalkNext,
0173 digitizerUtility::Ph2Amplitude(signalInElectrons_odd_column_Xtalk_next_column, nullptr, -1.0));
0174 else
0175 signalNew.emplace(
0176 chanXtalkNext,
0177 digitizerUtility::Ph2Amplitude(signalInElectrons_even_column_Xtalk_next_column, nullptr, -1.0));
0178 }
0179 }
0180 for (auto const& l : signalNew) {
0181 int chan = l.first;
0182 auto iter = theSignal.find(chan);
0183 if (iter != theSignal.end()) {
0184 iter->second += l.second.ampl();
0185 } else {
0186 theSignal.emplace(chan, digitizerUtility::Ph2Amplitude(l.second.ampl(), nullptr, -1.0));
0187 }
0188 }
0189 }
0190
0191 PixelDigitizerAlgorithm::TimewalkCurve::TimewalkCurve(const edm::ParameterSet& pset)
0192 : x_(pset.getParameter<std::vector<double>>("charge")), y_(pset.getParameter<std::vector<double>>("delay")) {
0193 if (x_.size() != y_.size())
0194 throw cms::Exception("Configuration")
0195 << "Timewalk model error: the number of charge values does not match the number of delay values!";
0196 }
0197
0198 double PixelDigitizerAlgorithm::TimewalkCurve::operator()(double x) const {
0199 auto it = std::lower_bound(x_.begin(), x_.end(), x);
0200 if (it == x_.begin())
0201 return y_.front();
0202 if (it == x_.end())
0203 return y_.back();
0204 int index = std::distance(x_.begin(), it);
0205 double x_high = *it;
0206 double x_low = *(--it);
0207 double p = (x - x_low) / (x_high - x_low);
0208 return p * y_[index] + (1 - p) * y_[index - 1];
0209 }
0210
0211 PixelDigitizerAlgorithm::TimewalkModel::TimewalkModel(const edm::ParameterSet& pset) {
0212 threshold_values = pset.getParameter<std::vector<double>>("ThresholdValues");
0213 const auto& curve_psetvec = pset.getParameter<std::vector<edm::ParameterSet>>("Curves");
0214 if (threshold_values.size() != curve_psetvec.size())
0215 throw cms::Exception("Configuration")
0216 << "Timewalk model error: the number of threshold values does not match the number of curves.";
0217 for (const auto& curve_pset : curve_psetvec)
0218 curves.emplace_back(curve_pset);
0219 }
0220
0221 double PixelDigitizerAlgorithm::TimewalkModel::operator()(double q_in, double q_threshold) const {
0222 auto index = find_closest_index(threshold_values, q_threshold);
0223 return curves[index](q_in);
0224 }
0225
0226 std::size_t PixelDigitizerAlgorithm::TimewalkModel::find_closest_index(const std::vector<double>& vec,
0227 double value) const {
0228 auto it = std::lower_bound(vec.begin(), vec.end(), value);
0229 if (it == vec.begin())
0230 return 0;
0231 else if (it == vec.end())
0232 return vec.size() - 1;
0233 else {
0234 auto it_upper = it;
0235 auto it_lower = --it;
0236 auto closest = (value - *it_lower > *it_upper - value) ? it_upper : it_lower;
0237 return std::distance(vec.begin(), closest);
0238 }
0239 }
0240
0241
0242
0243 bool PixelDigitizerAlgorithm::isAboveThreshold(const digitizerUtility::SimHitInfo* hitInfo,
0244 float charge,
0245 float thr) const {
0246 if (charge < thr)
0247 return false;
0248 if (apply_timewalk_ && hitInfo) {
0249 float corrected_time = hitInfo->time();
0250 double time = corrected_time + timewalk_model_(charge, thr);
0251 return (time >= theTofLowerCut_ && time < theTofUpperCut_);
0252 } else
0253 return true;
0254 }
0255
0256
0257
0258 void PixelDigitizerAlgorithm::module_killing_DB(const Phase2TrackerGeomDetUnit* pixdet) {
0259 throw cms::Exception("PixelDigitizerAlgorithm") << "Trying to kill modules from the pixel digitizer."
0260 << " This method is not yet implemented!";
0261 }