Back to home page

Project CMSSW displayed by LXR

 
 

    


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 // 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 
0021 using namespace edm;
0022 using namespace sipixelobjects;
0023 
0024 void PixelDigitizerAlgorithm::init(const edm::EventSetup& es) {
0025   if (use_ineff_from_db_)  // load gain calibration service fromdb...
0026     theSiPixelGainCalibrationService_->setESObjects(es);
0027 
0028   if (use_deadmodule_DB_)
0029     siPixelBadModule_ = &es.getData(siPixelBadModuleToken_);
0030 
0031   if (use_LorentzAngle_DB_)  // Get Lorentz angle from DB record
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 // -- Select the Hit for Digitization
0077 //
0078 bool PixelDigitizerAlgorithm::select_hit(const PSimHit& hit, double tCorr, double& sigScale) const {
0079   // in case of signal-shape emulation do not apply [TofLower,TofUpper] selection
0080   double toa = hit.tof() - tCorr;
0081   return apply_timewalk_ || (toa >= theTofLowerCut_ && toa < theTofUpperCut_);
0082 }
0083 
0084 // ======================================================================
0085 //
0086 //  Add  Cross-talk contribution
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   // cross-talk calculation valid for the case of 25x100 pixels
0096   const float pitch_first = 0.0025;
0097   const float pitch_second = 0.0100;
0098 
0099   // 0.5 um tolerance when comparing the pitch to accommodate the small changes in different TK geometrie (temporary fix)
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();  // signal in electrons
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     // subtract the charge which will be shared
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 // -- Compare Signal with Threshold
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 // -- Read Bad Channels from the Condidion DB and kill channels/module accordingly
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 }