Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-04-04 01:27:04

0001 //class SiPixelDigitizerAlgorithm SimTracker/SiPixelDigitizer/src/SiPixelDigitizerAlgoithm.cc
0002 
0003 // Original Author Danek Kotlinski
0004 // Ported in CMSSW by  Michele Pioppi-INFN perugia
0005 // Added DB capabilities by F.Blekman, Cornell University
0006 //         Created:  Mon Sep 26 11:08:32 CEST 2005
0007 // Add tof, change AddNoise to tracked. 4/06
0008 // Change drift direction. 6/06 d.k.
0009 // Add the statuis (non-rate dependent) inefficiency.
0010 //     -1 - no ineffciency
0011 //      0 - static inefficency only
0012 //    1,2 - low-lumi rate dependent inefficency added
0013 //     10 - high-lumi inefficiency added
0014 // Adopt the correct drift sign convetion from Morris Swartz. d.k. 8/06
0015 // Add more complex misscalinbration, change kev/e to 3.61, diff=3.7,d.k.9/06
0016 // Add the readout channel electronic noise. d.k. 3/07
0017 // Lower the pixel noise from 500 to 175elec.
0018 // Change the input threshold from noise units to electrons.
0019 // Lower the amount of static dead pixels from 0.01 to 0.001.
0020 // Modify to the new random number services. d.k. 5/07
0021 // Protect against sigma=0 (delta tracks on the surface). d.k.5/07
0022 // Change the TOF cut to lower and upper limit. d.k. 7/07
0023 //
0024 // July 2008: Split Lorentz Angle configuration in BPix/FPix (V. Cuplov)
0025 // tanLorentzAngleperTesla_FPix=0.0912 and tanLorentzAngleperTesla_BPix=0.106
0026 // Sept. 2008: Disable Pixel modules which are declared dead in the configuration python file. (V. Cuplov)
0027 // Oct. 2008: Accessing/Reading the Lorentz angle from the DataBase instead of the cfg file. (V. Cuplov)
0028 // Accessing dead modules from the DB. Implementation done and tested on a test.db
0029 // Do not use this option for now. The PixelQuality Objects are not in the official DB yet.
0030 // Feb. 2009: Split Fpix and Bpix threshold and use official numbers (V. Cuplov)
0031 // ThresholdInElectrons_FPix = 2870 and ThresholdInElectrons_BPix = 3700
0032 // update the electron to VCAL conversion using: VCAL_electrons = VCAL * 65.5 - 414
0033 // Feb. 2009: Threshold gaussian smearing (V. Cuplov)
0034 // March, 2009: changed DB access to *SimRcd objects (to de-couple the DB objects from reco chain) (F. Blekman)
0035 // May, 2009: Pixel charge VCAL smearing. (V. Cuplov)
0036 // November, 2009: new parameterization of the pixel response. (V. Cuplov)
0037 // December, 2009: Fix issue with different compilers.
0038 // October, 2010: Improvement: Removing single dead ROC (V. Cuplov)
0039 // November, 2010: Bug fix in removing TBMB/A half-modules (V. Cuplov)
0040 // February, 2011: Time improvement in DriftDirection()  (J. Bashir Butt)
0041 // June, 2011: Bug Fix for pixels on ROC edges in module_killing_DB() (J. Bashir Butt)
0042 // February, 2018: Implement cluster charge reweighting (P. Schuetze, with code from A. Hazi)
0043 #include <iostream>
0044 #include <iomanip>
0045 #include <gsl/gsl_sf_erf.h>
0046 
0047 #include "CLHEP/Random/RandFlat.h"
0048 #include "CLHEP/Random/RandGaussQ.h"
0049 #include "CLHEP/Random/RandGeneral.h"
0050 #include "CalibTracker/SiPixelESProducers/interface/SiPixelGainCalibrationOfflineSimService.h"
0051 #include "CondFormats/SiPixelObjects/interface/CablingPathToDetUnit.h"
0052 #include "CondFormats/SiPixelObjects/interface/GlobalPixel.h"
0053 #include "CondFormats/SiPixelObjects/interface/LocalPixel.h"
0054 #include "CondFormats/SiPixelObjects/interface/PixelFEDCabling.h"
0055 #include "CondFormats/SiPixelObjects/interface/PixelFEDLink.h"
0056 #include "CondFormats/SiPixelObjects/interface/PixelIndices.h"
0057 #include "CondFormats/SiPixelObjects/interface/PixelROC.h"
0058 #include "CondFormats/SiPixelObjects/interface/SiPixel2DTemplateDBObject.h"
0059 #include "CondFormats/SiPixelObjects/interface/SiPixelDynamicInefficiency.h"
0060 #include "CondFormats/SiPixelObjects/interface/SiPixelFEDChannelContainer.h"
0061 #include "CondFormats/SiPixelObjects/interface/SiPixelFedCabling.h"
0062 #include "CondFormats/SiPixelObjects/interface/SiPixelFedCablingMap.h"
0063 #include "CondFormats/SiPixelObjects/interface/SiPixelFedCablingTree.h"
0064 #include "CondFormats/SiPixelObjects/interface/SiPixelFrameReverter.h"
0065 #include "CondFormats/SiPixelObjects/interface/SiPixelLorentzAngle.h"
0066 #include "CondFormats/SiPixelObjects/interface/SiPixelQuality.h"
0067 #include "CondFormats/SiPixelObjects/interface/SiPixelQualityProbabilities.h"
0068 #include "DataFormats/DetId/interface/DetId.h"
0069 #include "DataFormats/FEDRawData/interface/FEDNumbering.h"
0070 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0071 #include "DataFormats/SiPixelDigi/interface/PixelDigi.h"
0072 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0073 #include "FWCore/AbstractServices/interface/RandomNumberGenerator.h"
0074 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0075 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0076 #include "FWCore/ServiceRegistry/interface/Service.h"
0077 #include "FWCore/Utilities/interface/Exception.h"
0078 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0079 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0080 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0081 #include "SimDataFormats/PileupSummaryInfo/interface/PileupMixingContent.h"
0082 #include "SimDataFormats/Track/interface/SimTrack.h"
0083 #include "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.h"
0084 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0085 #include "SimGeneral/NoiseGenerators/interface/GaussianTailNoiseGenerator.h"
0086 #include "SimTracker/Common/interface/SiG4UniversalFluctuation.h"
0087 #include "SimTracker/Common/interface/SiPixelChargeReweightingAlgorithm.h"
0088 #include "SimTracker/SiPixelDigitizer/plugins/SiPixelDigitizerAlgorithm.h"
0089 
0090 using namespace edm;
0091 using namespace sipixelobjects;
0092 
0093 void SiPixelDigitizerAlgorithm::init(const edm::EventSetup& es) {
0094   if (use_ineff_from_db_) {  // load gain calibration service fromdb...
0095     theSiPixelGainCalibrationService_->setESObjects(es);
0096   }
0097   if (use_deadmodule_DB_) {
0098     SiPixelBadModule_ = &es.getData(SiPixelBadModuleToken_);
0099   }
0100   if (use_LorentzAngle_DB_) {
0101     // Get Lorentz angle from DB record
0102     SiPixelLorentzAngle_ = &es.getData(SiPixelLorentzAngleToken_);
0103   }
0104   //gets the map and geometry from the DB (to kill ROCs)
0105   map_ = &es.getData(mapToken_);
0106   geom_ = &es.getData(geomToken_);
0107 
0108   if (KillBadFEDChannels) {
0109     scenarioProbability_ = &es.getData(scenarioProbabilityToken_);
0110     quality_map = &es.getData(PixelFEDChannelCollectionMapToken_);
0111 
0112     SiPixelQualityProbabilities::probabilityMap m_probabilities = scenarioProbability_->getProbability_Map();
0113     std::vector<std::string> allScenarios;
0114 
0115     std::transform(quality_map->begin(),
0116                    quality_map->end(),
0117                    std::back_inserter(allScenarios),
0118                    [](const PixelFEDChannelCollectionMap::value_type& pair) { return pair.first; });
0119 
0120     std::vector<std::string> allScenariosInProb;
0121 
0122     for (auto it = m_probabilities.begin(); it != m_probabilities.end(); ++it) {
0123       //int PUbin = it->first;
0124       for (const auto& entry : it->second) {
0125         auto scenario = entry.first;
0126         auto probability = entry.second;
0127         if (probability != 0) {
0128           if (std::find(allScenariosInProb.begin(), allScenariosInProb.end(), scenario) == allScenariosInProb.end()) {
0129             allScenariosInProb.push_back(scenario);
0130           }
0131         }  // if prob!=0
0132       }  // loop on the scenarios for that PU bin
0133     }  // loop on PU bins
0134 
0135     std::vector<std::string> notFound;
0136     std::copy_if(allScenariosInProb.begin(),
0137                  allScenariosInProb.end(),
0138                  std::back_inserter(notFound),
0139                  [&allScenarios](const std::string& arg) {
0140                    return (std::find(allScenarios.begin(), allScenarios.end(), arg) == allScenarios.end());
0141                  });
0142 
0143     if (!notFound.empty()) {
0144       for (const auto& entry : notFound) {
0145         LogError("SiPixelFEDChannelContainer")
0146             << "The requested scenario: " << entry << " is not found in the map!! \n";
0147       }
0148       throw cms::Exception("SiPixelDigitizerAlgorithm") << "Found: " << notFound.size()
0149                                                         << " missing scenario(s) in SiPixelStatusScenariosRcd while "
0150                                                            "present in SiPixelStatusScenarioProbabilityRcd \n";
0151     }
0152   }
0153   LogInfo("PixelDigitizer ") << " PixelDigitizerAlgorithm init \n";
0154   LogInfo("PixelDigitizer ") << " PixelDigitizerAlgorithm  --> UseReweighting " << UseReweighting << "\n";
0155   LogInfo("PixelDigitizer ") << " PixelDigitizerAlgorithm  -->  store_SimHitEntryExitPoints_ "
0156                              << store_SimHitEntryExitPoints_ << "\n";
0157   LogInfo("PixelDigitizer ") << " PixelDigitizerAlgorithm  -->  makeDigiSimLinks_ " << makeDigiSimLinks_ << "\n";
0158 
0159   TheNewSiPixelChargeReweightingAlgorithmClass->init(es);
0160 
0161   int collectionIndex = 0;  // I don't find what are the different collections here
0162   int tofBin = 0;
0163   for (int i1 = 1; i1 < 3; i1++) {
0164     for (int i2 = 0; i2 < 2; i2++) {
0165       if (i2 == 0) {
0166         tofBin = PixelDigiSimLink::LowTof;
0167       } else {
0168         tofBin = PixelDigiSimLink::HighTof;
0169       }
0170       subDetTofBin theSubDetTofBin = std::make_pair(i1, tofBin);
0171       SimHitCollMap[theSubDetTofBin] = collectionIndex;
0172       collectionIndex++;
0173     }
0174   }
0175 }
0176 
0177 //=========================================================================
0178 
0179 SiPixelDigitizerAlgorithm::SiPixelDigitizerAlgorithm(const edm::ParameterSet& conf, edm::ConsumesCollector iC)
0180     : mapToken_(iC.esConsumes()),
0181       geomToken_(iC.esConsumes()),
0182 
0183       _signal(),
0184       makeDigiSimLinks_(conf.getUntrackedParameter<bool>("makeDigiSimLinks", true)),
0185       store_SimHitEntryExitPoints_(
0186           conf.exists("store_SimHitEntryExitPoints") ? conf.getParameter<bool>("store_SimHitEntryExitPoints") : false),
0187       store_SimHitEntryExitPointsLite_(conf.exists("store_SimHitEntryExitPointsLite")
0188                                            ? conf.getParameter<bool>("store_SimHitEntryExitPointsLite")
0189                                            : false),
0190       use_ineff_from_db_(conf.getParameter<bool>("useDB")),
0191       use_module_killing_(conf.getParameter<bool>("killModules")),       // boolean to kill or not modules
0192       use_deadmodule_DB_(conf.getParameter<bool>("DeadModules_DB")),     // boolean to access dead modules from DB
0193       use_LorentzAngle_DB_(conf.getParameter<bool>("LorentzAngle_DB")),  // boolean to access Lorentz angle from DB
0194 
0195       DeadModules(use_deadmodule_DB_ ? Parameters()
0196                                      : conf.getParameter<Parameters>("DeadModules")),  // get dead module from cfg file
0197 
0198       TheNewSiPixelChargeReweightingAlgorithmClass(),
0199 
0200       // Common pixel parameters
0201       // These are parameters which are not likely to be changed
0202       GeVperElectron(3.61E-09),                             // 1 electron(3.61eV, 1keV(277e, mod 9/06 d.k.
0203       Sigma0(0.00037),                                      // Charge diffusion constant 7->3.7
0204       Dist300(0.0300),                                      //   normalized to 300micron Silicon
0205       alpha2Order(conf.getParameter<bool>("Alpha2Order")),  // switch on/off of E.B effect
0206       ClusterWidth(3.),                                     // Charge integration spread on the collection plane
0207 
0208       // get external parameters:
0209       // To account for upgrade geometries do not assume the number
0210       // of layers or disks.
0211       NumberOfBarrelLayers(conf.exists("NumPixelBarrel") ? conf.getParameter<int>("NumPixelBarrel") : 3),
0212       NumberOfEndcapDisks(conf.exists("NumPixelEndcap") ? conf.getParameter<int>("NumPixelEndcap") : 2),
0213 
0214       // ADC calibration 1adc count(135e.
0215       // Corresponds to 2adc/kev, 270[e/kev]/135[e/adc](2[adc/kev]
0216       // Be carefull, this parameter is also used in SiPixelDet.cc to
0217       // calculate the noise in adc counts from noise in electrons.
0218       // Both defaults should be the same.
0219       theElectronPerADC(conf.getParameter<double>("ElectronPerAdc")),
0220 
0221       // ADC saturation value, 255(8bit adc.
0222       //theAdcFullScale(conf.getUntrackedParameter<int>("AdcFullScale",255)),
0223       theAdcFullScale(conf.getParameter<int>("AdcFullScale")),
0224       theAdcFullScLateCR(conf.getParameter<int>("AdcFullScLateCR")),
0225 
0226       // Noise in electrons:
0227       // Pixel cell noise, relevant for generating noisy pixels
0228       theNoiseInElectrons(conf.getParameter<double>("NoiseInElectrons")),
0229 
0230       // Fill readout noise, including all readout chain, relevant for smearing
0231       //theReadoutNoise(conf.getUntrackedParameter<double>("ReadoutNoiseInElec",500.)),
0232       theReadoutNoise(conf.getParameter<double>("ReadoutNoiseInElec")),
0233 
0234       // Pixel threshold in units of noise:
0235       // thePixelThreshold(conf.getParameter<double>("ThresholdInNoiseUnits")),
0236       // Pixel threshold in electron units.
0237       theThresholdInE_FPix(conf.getParameter<double>("ThresholdInElectrons_FPix")),
0238       theThresholdInE_BPix(conf.getParameter<double>("ThresholdInElectrons_BPix")),
0239       theThresholdInE_BPix_L1(conf.exists("ThresholdInElectrons_BPix_L1")
0240                                   ? conf.getParameter<double>("ThresholdInElectrons_BPix_L1")
0241                                   : theThresholdInE_BPix),
0242       theThresholdInE_BPix_L2(conf.exists("ThresholdInElectrons_BPix_L2")
0243                                   ? conf.getParameter<double>("ThresholdInElectrons_BPix_L2")
0244                                   : theThresholdInE_BPix),
0245 
0246       // Add threshold gaussian smearing:
0247       theThresholdSmearing_FPix(conf.getParameter<double>("ThresholdSmearing_FPix")),
0248       theThresholdSmearing_BPix(conf.getParameter<double>("ThresholdSmearing_BPix")),
0249       theThresholdSmearing_BPix_L1(conf.exists("ThresholdSmearing_BPix_L1")
0250                                        ? conf.getParameter<double>("ThresholdSmearing_BPix_L1")
0251                                        : theThresholdSmearing_BPix),
0252       theThresholdSmearing_BPix_L2(conf.exists("ThresholdSmearing_BPix_L2")
0253                                        ? conf.getParameter<double>("ThresholdSmearing_BPix_L2")
0254                                        : theThresholdSmearing_BPix),
0255 
0256       // electrons to VCAL conversion needed in misscalibrate()
0257       electronsPerVCAL(conf.getParameter<double>("ElectronsPerVcal")),
0258       electronsPerVCAL_Offset(conf.getParameter<double>("ElectronsPerVcal_Offset")),
0259       electronsPerVCAL_L1(conf.exists("ElectronsPerVcal_L1") ? conf.getParameter<double>("ElectronsPerVcal_L1")
0260                                                              : electronsPerVCAL),
0261       electronsPerVCAL_L1_Offset(conf.exists("ElectronsPerVcal_L1_Offset")
0262                                      ? conf.getParameter<double>("ElectronsPerVcal_L1_Offset")
0263                                      : electronsPerVCAL_Offset),
0264 
0265       //theTofCut 12.5, cut in particle TOD +/- 12.5ns
0266       //theTofCut(conf.getUntrackedParameter<double>("TofCut",12.5)),
0267       theTofLowerCut(conf.getParameter<double>("TofLowerCut")),
0268       theTofUpperCut(conf.getParameter<double>("TofUpperCut")),
0269 
0270       // Get the Lorentz angle from the cfg file:
0271       tanLorentzAnglePerTesla_FPix(use_LorentzAngle_DB_ ? 0.0
0272                                                         : conf.getParameter<double>("TanLorentzAnglePerTesla_FPix")),
0273       tanLorentzAnglePerTesla_BPix(use_LorentzAngle_DB_ ? 0.0
0274                                                         : conf.getParameter<double>("TanLorentzAnglePerTesla_BPix")),
0275 
0276       // signal response new parameterization: split Fpix and BPix
0277       FPix_p0(conf.getParameter<double>("FPix_SignalResponse_p0")),
0278       FPix_p1(conf.getParameter<double>("FPix_SignalResponse_p1")),
0279       FPix_p2(conf.getParameter<double>("FPix_SignalResponse_p2")),
0280       FPix_p3(conf.getParameter<double>("FPix_SignalResponse_p3")),
0281 
0282       BPix_p0(conf.getParameter<double>("BPix_SignalResponse_p0")),
0283       BPix_p1(conf.getParameter<double>("BPix_SignalResponse_p1")),
0284       BPix_p2(conf.getParameter<double>("BPix_SignalResponse_p2")),
0285       BPix_p3(conf.getParameter<double>("BPix_SignalResponse_p3")),
0286 
0287       // Add noise
0288       addNoise(conf.getParameter<bool>("AddNoise")),
0289 
0290       // Smear the pixel charge with a gaussian which RMS is a function of the
0291       // pixel charge (Danek's study)
0292       addChargeVCALSmearing(conf.getParameter<bool>("ChargeVCALSmearing")),
0293 
0294       // Add noisy pixels
0295       addNoisyPixels(conf.getParameter<bool>("AddNoisyPixels")),
0296 
0297       // Fluctuate charge in track subsegments
0298       fluctuateCharge(conf.getUntrackedParameter<bool>("FluctuateCharge", true)),
0299 
0300       // Control the pixel inefficiency
0301       AddPixelInefficiency(conf.getParameter<bool>("AddPixelInefficiency")),
0302       KillBadFEDChannels(conf.getParameter<bool>("KillBadFEDChannels")),
0303 
0304       // Add threshold gaussian smearing:
0305       addThresholdSmearing(conf.getParameter<bool>("AddThresholdSmearing")),
0306 
0307       // Get the constants for the miss-calibration studies
0308       doMissCalibrate(conf.getParameter<bool>("MissCalibrate")),       // Enable miss-calibration
0309       doMissCalInLateCR(conf.getParameter<bool>("MissCalInLateCR")),   // Enable miss-calibration
0310       theGainSmearing(conf.getParameter<double>("GainSmearing")),      // sigma of the gain smearing
0311       theOffsetSmearing(conf.getParameter<double>("OffsetSmearing")),  //sigma of the offset smearing
0312 
0313       // Add pixel radiation damage for upgrade studies
0314       AddPixelAging(conf.getParameter<bool>("DoPixelAging")),
0315       UseReweighting(conf.getParameter<bool>("UseReweighting")),
0316 
0317       // delta cutoff in MeV, has to be same as in OSCAR(0.030/cmsim=1.0 MeV
0318       //tMax(0.030), // In MeV.
0319       //tMax(conf.getUntrackedParameter<double>("deltaProductionCut",0.030)),
0320       tMax(conf.getParameter<double>("deltaProductionCut")),
0321 
0322       fluctuate(fluctuateCharge ? new SiG4UniversalFluctuation() : nullptr),
0323       theNoiser(addNoise ? new GaussianTailNoiseGenerator() : nullptr),
0324       calmap((doMissCalibrate || doMissCalInLateCR) ? initCal() : std::map<int, CalParameters, std::less<int> >()),
0325       theSiPixelGainCalibrationService_(use_ineff_from_db_ ? new SiPixelGainCalibrationOfflineSimService(conf, iC)
0326                                                            : nullptr),
0327       pixelEfficiencies_(conf, AddPixelInefficiency, NumberOfBarrelLayers, NumberOfEndcapDisks),
0328       pixelAging_(conf, AddPixelAging, NumberOfBarrelLayers, NumberOfEndcapDisks) {
0329   if (use_deadmodule_DB_) {
0330     //string to specify SiPixelQuality label
0331     SiPixelBadModuleToken_ = iC.esConsumes(edm::ESInputTag("", conf.getParameter<std::string>("SiPixelQualityLabel")));
0332   }
0333   if (use_LorentzAngle_DB_) {
0334     SiPixelLorentzAngleToken_ = iC.esConsumes();
0335   }
0336   if (AddPixelInefficiency && !pixelEfficiencies_.FromConfig) {
0337     // TODO: in practice the bunchspacing is known at MixingModule
0338     // construction time, and thus we could declare the consumption of
0339     // the actual product. In principle, however, MixingModule is
0340     // capable of updating (parts of) its configuration from the
0341     // EventSetup, so if that capability is really needed we'd need to
0342     // invent something new (similar to mayConsume in the ESProducer
0343     // side). So for now, let's consume both payloads.
0344     SiPixelDynamicInefficiencyToken_ = iC.esConsumes();
0345   }
0346   if (KillBadFEDChannels) {
0347     scenarioProbabilityToken_ = iC.esConsumes();
0348     PixelFEDChannelCollectionMapToken_ = iC.esConsumes();
0349   }
0350 
0351   LogInfo("PixelDigitizer ") << "SiPixelDigitizerAlgorithm constructed"
0352                              << "Configuration parameters:"
0353                              << "Threshold/Gain = "
0354                              << "threshold in electron FPix = " << theThresholdInE_FPix
0355                              << "threshold in electron BPix = " << theThresholdInE_BPix
0356                              << "threshold in electron BPix Layer1 = " << theThresholdInE_BPix_L1
0357                              << "threshold in electron BPix Layer2 = " << theThresholdInE_BPix_L2 << " "
0358                              << theElectronPerADC << " " << theAdcFullScale << " The delta cut-off is set to " << tMax
0359                              << " pix-inefficiency " << AddPixelInefficiency;
0360 
0361   LogInfo("PixelDigitizer ") << " SiPixelDigitizerAlgorithm constructed  with  UseReweighting " << UseReweighting
0362                              << " and store_SimHitEntryExitPoints_ " << store_SimHitEntryExitPoints_ << " \n";
0363 
0364   TheNewSiPixelChargeReweightingAlgorithmClass = std::make_unique<SiPixelChargeReweightingAlgorithm>(conf, iC);
0365 }
0366 
0367 std::map<int, SiPixelDigitizerAlgorithm::CalParameters, std::less<int> > SiPixelDigitizerAlgorithm::initCal() const {
0368   std::map<int, SiPixelDigitizerAlgorithm::CalParameters, std::less<int> > calmap;
0369   // Prepare for the analog amplitude miss-calibration
0370   LogDebug("PixelDigitizer ") << " miss-calibrate the pixel amplitude \n";
0371 
0372   const bool ReadCalParameters = false;
0373   if (ReadCalParameters) {  // Read the calibration files from file
0374     // read the calibration constants from a file (testing only)
0375     std::ifstream in_file;  // data file pointer
0376     char filename[80] = "phCalibrationFit_C0.dat";
0377 
0378     in_file.open(filename, std::ios::in);  // in C++
0379     if (in_file.bad()) {
0380       LogInfo("PixelDigitizer ") << " File not found \n ";
0381       return calmap;  // signal error
0382     }
0383     LogInfo("PixelDigitizer ") << " file opened : " << filename << "\n";
0384 
0385     char line[500];
0386     for (int i = 0; i < 3; i++) {
0387       in_file.getline(line, 500, '\n');
0388       LogInfo("PixelDigitizer ") << line << "\n";
0389     }
0390 
0391     LogInfo("PixelDigitizer ") << " test map"
0392                                << "\n";
0393 
0394     float par0, par1, par2, par3;
0395     int colid, rowid;
0396     std::string name;
0397     // Read MC tracks
0398     for (int i = 0; i < (52 * 80); i++) {  // loop over tracks
0399       in_file >> par0 >> par1 >> par2 >> par3 >> name >> colid >> rowid;
0400       if (in_file.bad()) {  // check for errors
0401         LogError("PixelDigitizer") << "Cannot read data file for calmap"
0402                                    << "\n";
0403         return calmap;
0404       }
0405       if (in_file.eof() != 0) {
0406         LogError("PixelDigitizer") << "calmap " << in_file.eof() << " " << in_file.gcount() << " " << in_file.fail()
0407                                    << " " << in_file.good() << " end of file "
0408                                    << "\n";
0409         return calmap;
0410       }
0411 
0412       LogDebug("PixelDigitizer ") << " line " << i << " " << par0 << " " << par1 << " " << par2 << " " << par3 << " "
0413                                   << colid << " " << rowid << "\n";
0414 
0415       CalParameters onePix;
0416       onePix.p0 = par0;
0417       onePix.p1 = par1;
0418       onePix.p2 = par2;
0419       onePix.p3 = par3;
0420 
0421       // Convert ROC pixel index to channel
0422       int chan = PixelIndices::pixelToChannelROC(rowid, colid);
0423       calmap.insert(std::pair<int, CalParameters>(chan, onePix));
0424 
0425       // Testing the index conversion, can be skipped
0426       std::pair<int, int> p = PixelIndices::channelToPixelROC(chan);
0427       if (rowid != p.first)
0428         LogInfo("PixelDigitizer ") << " wrong channel row " << rowid << " " << p.first << "\n";
0429       if (colid != p.second)
0430         LogInfo("PixelDigitizer ") << " wrong channel col " << colid << " " << p.second << "\n";
0431 
0432     }  // pixel loop in a ROC
0433 
0434     LogInfo("PixelDigitizer ") << " map size  " << calmap.size() << " max " << calmap.max_size() << " "
0435                                << calmap.empty() << "\n";
0436 
0437     //     LogInfo("PixelDigitizer ") << " map size  " << calmap.size()  << "\n";
0438     //     map<int,CalParameters,std::less<int> >::iterator ix,it;
0439     //     map<int,CalParameters,std::less<int> >::const_iterator ip;
0440     //     for (ix = calmap.begin(); ix != calmap.end(); ++ix) {
0441     //       int i = (*ix).first;
0442     //       std::pair<int,int> p = channelToPixelROC(i);
0443     //       it  = calmap.find(i);
0444     //       CalParameters y  = (*it).second;
0445     //       CalParameters z = (*ix).second;
0446     //       LogInfo("PixelDigitizer ") << i <<" "<<p.first<<" "<<p.second<<" "<<y.p0<<" "<<z.p0<<" "<<calmap[i].p0<<"\n";
0447 
0448     //       //int dummy=0;
0449     //       //cin>>dummy;
0450     //     }
0451 
0452   }  // end if readparameters
0453   return calmap;
0454 }  // end initCal()
0455 
0456 //=========================================================================
0457 SiPixelDigitizerAlgorithm::~SiPixelDigitizerAlgorithm() {
0458   LogDebug("PixelDigitizer") << "SiPixelDigitizerAlgorithm deleted";
0459 }
0460 
0461 // Read DynIneff Scale factors from Configuration
0462 SiPixelDigitizerAlgorithm::PixelEfficiencies::PixelEfficiencies(const edm::ParameterSet& conf,
0463                                                                 bool AddPixelInefficiency,
0464                                                                 int NumberOfBarrelLayers,
0465                                                                 int NumberOfEndcapDisks) {
0466   // pixel inefficiency
0467   // Don't use Hard coded values, read inefficiencies in from DB/python config or don't use any
0468   int NumberOfTotLayers = NumberOfBarrelLayers + NumberOfEndcapDisks;
0469   FPixIndex = NumberOfBarrelLayers;
0470   if (AddPixelInefficiency) {
0471     FromConfig = conf.exists("thePixelColEfficiency_BPix1") && conf.exists("thePixelColEfficiency_BPix2") &&
0472                  conf.exists("thePixelColEfficiency_BPix3") && conf.exists("thePixelColEfficiency_FPix1") &&
0473                  conf.exists("thePixelColEfficiency_FPix2") && conf.exists("thePixelEfficiency_BPix1") &&
0474                  conf.exists("thePixelEfficiency_BPix2") && conf.exists("thePixelEfficiency_BPix3") &&
0475                  conf.exists("thePixelEfficiency_FPix1") && conf.exists("thePixelEfficiency_FPix2") &&
0476                  conf.exists("thePixelChipEfficiency_BPix1") && conf.exists("thePixelChipEfficiency_BPix2") &&
0477                  conf.exists("thePixelChipEfficiency_BPix3") && conf.exists("thePixelChipEfficiency_FPix1") &&
0478                  conf.exists("thePixelChipEfficiency_FPix2");
0479     if (NumberOfBarrelLayers == 3)
0480       FromConfig = FromConfig && conf.exists("theLadderEfficiency_BPix1") && conf.exists("theLadderEfficiency_BPix2") &&
0481                    conf.exists("theLadderEfficiency_BPix3") && conf.exists("theModuleEfficiency_BPix1") &&
0482                    conf.exists("theModuleEfficiency_BPix2") && conf.exists("theModuleEfficiency_BPix3") &&
0483                    conf.exists("thePUEfficiency_BPix1") && conf.exists("thePUEfficiency_BPix2") &&
0484                    conf.exists("thePUEfficiency_BPix3") && conf.exists("theInnerEfficiency_FPix1") &&
0485                    conf.exists("theInnerEfficiency_FPix2") && conf.exists("theOuterEfficiency_FPix1") &&
0486                    conf.exists("theOuterEfficiency_FPix2") && conf.exists("thePUEfficiency_FPix_Inner") &&
0487                    conf.exists("thePUEfficiency_FPix_Outer") && conf.exists("theInstLumiScaleFactor");
0488     if (NumberOfBarrelLayers >= 4)
0489       FromConfig = FromConfig && conf.exists("thePixelColEfficiency_BPix4") &&
0490                    conf.exists("thePixelEfficiency_BPix4") && conf.exists("thePixelChipEfficiency_BPix4");
0491     if (NumberOfEndcapDisks >= 3)
0492       FromConfig = FromConfig && conf.exists("thePixelColEfficiency_FPix3") &&
0493                    conf.exists("thePixelEfficiency_FPix3") && conf.exists("thePixelChipEfficiency_FPix3");
0494     if (FromConfig) {
0495       LogInfo("PixelDigitizer ") << "The PixelDigitizer inefficiency configuration is read from the config file.\n";
0496       theInstLumiScaleFactor = conf.getParameter<double>("theInstLumiScaleFactor");
0497       int i = 0;
0498       thePixelColEfficiency[i++] = conf.getParameter<double>("thePixelColEfficiency_BPix1");
0499       thePixelColEfficiency[i++] = conf.getParameter<double>("thePixelColEfficiency_BPix2");
0500       thePixelColEfficiency[i++] = conf.getParameter<double>("thePixelColEfficiency_BPix3");
0501       if (NumberOfBarrelLayers >= 4) {
0502         thePixelColEfficiency[i++] = conf.getParameter<double>("thePixelColEfficiency_BPix4");
0503       }
0504       //
0505       i = 0;
0506       thePixelEfficiency[i++] = conf.getParameter<double>("thePixelEfficiency_BPix1");
0507       thePixelEfficiency[i++] = conf.getParameter<double>("thePixelEfficiency_BPix2");
0508       thePixelEfficiency[i++] = conf.getParameter<double>("thePixelEfficiency_BPix3");
0509       if (NumberOfBarrelLayers >= 4) {
0510         thePixelEfficiency[i++] = conf.getParameter<double>("thePixelEfficiency_BPix4");
0511       }
0512       //
0513       i = 0;
0514       thePixelChipEfficiency[i++] = conf.getParameter<double>("thePixelChipEfficiency_BPix1");
0515       thePixelChipEfficiency[i++] = conf.getParameter<double>("thePixelChipEfficiency_BPix2");
0516       thePixelChipEfficiency[i++] = conf.getParameter<double>("thePixelChipEfficiency_BPix3");
0517       if (NumberOfBarrelLayers >= 4) {
0518         thePixelChipEfficiency[i++] = conf.getParameter<double>("thePixelChipEfficiency_BPix4");
0519       }
0520       //
0521       if (NumberOfBarrelLayers == 3) {
0522         i = 0;
0523         theLadderEfficiency_BPix[i++] = conf.getParameter<std::vector<double> >("theLadderEfficiency_BPix1");
0524         theLadderEfficiency_BPix[i++] = conf.getParameter<std::vector<double> >("theLadderEfficiency_BPix2");
0525         theLadderEfficiency_BPix[i++] = conf.getParameter<std::vector<double> >("theLadderEfficiency_BPix3");
0526         if (((theLadderEfficiency_BPix[0].size() != 20) || (theLadderEfficiency_BPix[1].size() != 32) ||
0527              (theLadderEfficiency_BPix[2].size() != 44)) &&
0528             (NumberOfBarrelLayers == 3))
0529           throw cms::Exception("Configuration") << "Wrong ladder number in efficiency config!";
0530         //
0531         i = 0;
0532         theModuleEfficiency_BPix[i++] = conf.getParameter<std::vector<double> >("theModuleEfficiency_BPix1");
0533         theModuleEfficiency_BPix[i++] = conf.getParameter<std::vector<double> >("theModuleEfficiency_BPix2");
0534         theModuleEfficiency_BPix[i++] = conf.getParameter<std::vector<double> >("theModuleEfficiency_BPix3");
0535         if (((theModuleEfficiency_BPix[0].size() != 4) || (theModuleEfficiency_BPix[1].size() != 4) ||
0536              (theModuleEfficiency_BPix[2].size() != 4)) &&
0537             (NumberOfBarrelLayers == 3))
0538           throw cms::Exception("Configuration") << "Wrong module number in efficiency config!";
0539         //
0540         thePUEfficiency.push_back(conf.getParameter<std::vector<double> >("thePUEfficiency_BPix1"));
0541         thePUEfficiency.push_back(conf.getParameter<std::vector<double> >("thePUEfficiency_BPix2"));
0542         thePUEfficiency.push_back(conf.getParameter<std::vector<double> >("thePUEfficiency_BPix3"));
0543         if (((thePUEfficiency[0].empty()) || (thePUEfficiency[1].empty()) || (thePUEfficiency[2].empty())) &&
0544             (NumberOfBarrelLayers == 3))
0545           throw cms::Exception("Configuration")
0546               << "At least one PU efficiency (BPix) number is needed in efficiency config!";
0547       }
0548       // The next is needed for Phase2 Tracker studies
0549       if (NumberOfBarrelLayers >= 5) {
0550         if (NumberOfTotLayers > 20) {
0551           throw cms::Exception("Configuration") << "SiPixelDigitizer was given more layers than it can handle";
0552         }
0553         // For Phase2 tracker layers just set the outermost BPix inefficiency to 99.9% THESE VALUES ARE HARDCODED ALSO ELSEWHERE IN THIS FILE
0554         for (int j = 5; j <= NumberOfBarrelLayers; j++) {
0555           thePixelColEfficiency[j - 1] = 0.999;
0556           thePixelEfficiency[j - 1] = 0.999;
0557           thePixelChipEfficiency[j - 1] = 0.999;
0558         }
0559       }
0560       //
0561       i = FPixIndex;
0562       thePixelColEfficiency[i++] = conf.getParameter<double>("thePixelColEfficiency_FPix1");
0563       thePixelColEfficiency[i++] = conf.getParameter<double>("thePixelColEfficiency_FPix2");
0564       if (NumberOfEndcapDisks >= 3) {
0565         thePixelColEfficiency[i++] = conf.getParameter<double>("thePixelColEfficiency_FPix3");
0566       }
0567       i = FPixIndex;
0568       thePixelEfficiency[i++] = conf.getParameter<double>("thePixelEfficiency_FPix1");
0569       thePixelEfficiency[i++] = conf.getParameter<double>("thePixelEfficiency_FPix2");
0570       if (NumberOfEndcapDisks >= 3) {
0571         thePixelEfficiency[i++] = conf.getParameter<double>("thePixelEfficiency_FPix3");
0572       }
0573       i = FPixIndex;
0574       thePixelChipEfficiency[i++] = conf.getParameter<double>("thePixelChipEfficiency_FPix1");
0575       thePixelChipEfficiency[i++] = conf.getParameter<double>("thePixelChipEfficiency_FPix2");
0576       if (NumberOfEndcapDisks >= 3) {
0577         thePixelChipEfficiency[i++] = conf.getParameter<double>("thePixelChipEfficiency_FPix3");
0578       }
0579       // The next is needed for Phase2 Tracker studies
0580       if (NumberOfEndcapDisks >= 4) {
0581         if (NumberOfTotLayers > 20) {
0582           throw cms::Exception("Configuration") << "SiPixelDigitizer was given more layers than it can handle";
0583         }
0584         // For Phase2 tracker layers just set the extra FPix disk inefficiency to 99.9% THESE VALUES ARE HARDCODED ALSO ELSEWHERE IN THIS FILE
0585         for (int j = 4 + FPixIndex; j <= NumberOfEndcapDisks + NumberOfBarrelLayers; j++) {
0586           thePixelColEfficiency[j - 1] = 0.999;
0587           thePixelEfficiency[j - 1] = 0.999;
0588           thePixelChipEfficiency[j - 1] = 0.999;
0589         }
0590       }
0591       //FPix Dynamic Inefficiency
0592       if (NumberOfBarrelLayers == 3) {
0593         i = FPixIndex;
0594         theInnerEfficiency_FPix[i++] = conf.getParameter<double>("theInnerEfficiency_FPix1");
0595         theInnerEfficiency_FPix[i++] = conf.getParameter<double>("theInnerEfficiency_FPix2");
0596         i = FPixIndex;
0597         theOuterEfficiency_FPix[i++] = conf.getParameter<double>("theOuterEfficiency_FPix1");
0598         theOuterEfficiency_FPix[i++] = conf.getParameter<double>("theOuterEfficiency_FPix2");
0599         thePUEfficiency.push_back(conf.getParameter<std::vector<double> >("thePUEfficiency_FPix_Inner"));
0600         thePUEfficiency.push_back(conf.getParameter<std::vector<double> >("thePUEfficiency_FPix_Outer"));
0601         if (((thePUEfficiency[3].empty()) || (thePUEfficiency[4].empty())) && (NumberOfEndcapDisks == 2))
0602           throw cms::Exception("Configuration")
0603               << "At least one (FPix) PU efficiency number is needed in efficiency config!";
0604         pu_scale.resize(thePUEfficiency.size());
0605       }
0606     } else
0607       LogInfo("PixelDigitizer ") << "The PixelDigitizer inefficiency configuration is read from the database.\n";
0608   }
0609   // the first "NumberOfBarrelLayers" settings [0],[1], ... , [NumberOfBarrelLayers-1] are for the barrel pixels
0610   // the next  "NumberOfEndcapDisks"  settings [NumberOfBarrelLayers],[NumberOfBarrelLayers+1], ... [NumberOfEndcapDisks+NumberOfBarrelLayers-1]
0611 }
0612 
0613 // Read DynIneff Scale factors from DB
0614 void SiPixelDigitizerAlgorithm::init_DynIneffDB(const edm::EventSetup& es) {
0615   LogDebug("PixelDigitizer ") << " In SiPixelDigitizerAlgorithm::init_DynIneffDB " << AddPixelInefficiency << "  "
0616                               << pixelEfficiencies_.FromConfig << "\n";
0617   if (AddPixelInefficiency && !pixelEfficiencies_.FromConfig) {
0618     SiPixelDynamicInefficiency_ = &es.getData(SiPixelDynamicInefficiencyToken_);
0619     pixelEfficiencies_.init_from_db(geom_, SiPixelDynamicInefficiency_);
0620   }
0621 }
0622 
0623 void SiPixelDigitizerAlgorithm::PixelEfficiencies::init_from_db(
0624     const TrackerGeometry* geom, const SiPixelDynamicInefficiency* SiPixelDynamicInefficiency) {
0625   theInstLumiScaleFactor = SiPixelDynamicInefficiency->gettheInstLumiScaleFactor();
0626   const std::map<uint32_t, double>& PixelGeomFactorsDBIn = SiPixelDynamicInefficiency->getPixelGeomFactors();
0627   const std::map<uint32_t, double>& ColGeomFactorsDB = SiPixelDynamicInefficiency->getColGeomFactors();
0628   const std::map<uint32_t, double>& ChipGeomFactorsDB = SiPixelDynamicInefficiency->getChipGeomFactors();
0629   const std::map<uint32_t, std::vector<double> >& PUFactors = SiPixelDynamicInefficiency->getPUFactors();
0630   std::vector<uint32_t> DetIdmasks = SiPixelDynamicInefficiency->getDetIdmasks();
0631 
0632   // Loop on all modules, initialize map for easy access
0633   for (const auto& it_module : geom->detUnits()) {
0634     if (dynamic_cast<PixelGeomDetUnit const*>(it_module) == nullptr)
0635       continue;
0636     const DetId detid = it_module->geographicalId();
0637     uint32_t rawid = detid.rawId();
0638     PixelGeomFactors[rawid] = 1;
0639     ColGeomFactors[rawid] = 1;
0640     ChipGeomFactors[rawid] = 1;
0641     PixelGeomFactorsROCStdPixels[rawid] = std::vector<double>(16, 1);
0642     PixelGeomFactorsROCBigPixels[rawid] = std::vector<double>(16, 1);
0643   }
0644 
0645   // ROC level inefficiency for phase 1 (disentangle scale factors for big and std size pixels)
0646   std::map<uint32_t, double> PixelGeomFactorsDB;
0647 
0648   LogDebug("PixelDigitizer ") << " Check PixelEfficiencies -- PixelGeomFactorsDBIn "
0649                               << "\n";
0650   if (geom->isThere(GeomDetEnumerators::P1PXB) || geom->isThere(GeomDetEnumerators::P1PXEC)) {
0651     for (auto db_factor : PixelGeomFactorsDBIn) {
0652       LogDebug("PixelDigitizer ") << "      db_factor  " << db_factor.first << "  " << db_factor.second << "\n";
0653 
0654       int shift = DetId(db_factor.first).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel) ? BPixRocIdShift
0655                                                                                                        : FPixRocIdShift;
0656       unsigned int rocMask = rocIdMaskBits << shift;
0657       unsigned int rocId = (((db_factor.first) & rocMask) >> shift);
0658       if (rocId != 0) {
0659         rocId--;
0660         unsigned int rawid = db_factor.first & (~rocMask);
0661         const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>(geom->idToDet(rawid));
0662         PixelTopology const* topology = &(theGeomDet->specificTopology());
0663         const int nPixelsInROC = topology->rowsperroc() * topology->colsperroc();
0664         const int nBigPixelsInROC = 2 * topology->rowsperroc() + topology->colsperroc() - 2;
0665         double factor = db_factor.second;
0666         double badFraction = 1 - factor;
0667         double bigPixelFraction = static_cast<double>(nBigPixelsInROC) / nPixelsInROC;
0668         double stdPixelFraction = 1. - bigPixelFraction;
0669 
0670         double badFractionBig = std::min(bigPixelFraction, badFraction);
0671         double badFractionStd = std::max(0., badFraction - badFractionBig);
0672         double badFractionBigReNormalized = badFractionBig / bigPixelFraction;
0673         double badFractionStdReNormalized = badFractionStd / stdPixelFraction;
0674         PixelGeomFactorsROCStdPixels[rawid][rocId] *= (1. - badFractionStdReNormalized);
0675         PixelGeomFactorsROCBigPixels[rawid][rocId] *= (1. - badFractionBigReNormalized);
0676       } else {
0677         PixelGeomFactorsDB[db_factor.first] = db_factor.second;
0678       }
0679     }
0680   }  // is Phase 1 geometry
0681   else {
0682     PixelGeomFactorsDB = PixelGeomFactorsDBIn;
0683   }
0684 
0685   LogDebug("PixelDigitizer ")
0686       << " Check PixelEfficiencies -- Loop on all modules and store module level geometrical scale factors  "
0687       << "\n";
0688   // Loop on all modules, store module level geometrical scale factors
0689   for (const auto& it_module : geom->detUnits()) {
0690     if (dynamic_cast<PixelGeomDetUnit const*>(it_module) == nullptr)
0691       continue;
0692     const DetId detid = it_module->geographicalId();
0693     uint32_t rawid = detid.rawId();
0694     for (auto db_factor : PixelGeomFactorsDB) {
0695       LogDebug("PixelDigitizer ") << "      db_factor PixelGeomFactorsDB  " << db_factor.first << "  "
0696                                   << db_factor.second << "\n";
0697       if (matches(detid, DetId(db_factor.first), DetIdmasks))
0698         PixelGeomFactors[rawid] *= db_factor.second;
0699     }
0700     for (auto db_factor : ColGeomFactorsDB) {
0701       LogDebug("PixelDigitizer ") << "      db_factor ColGeomFactorsDB  " << db_factor.first << "  " << db_factor.second
0702                                   << "\n";
0703       if (matches(detid, DetId(db_factor.first), DetIdmasks))
0704         ColGeomFactors[rawid] *= db_factor.second;
0705     }
0706     for (auto db_factor : ChipGeomFactorsDB) {
0707       LogDebug("PixelDigitizer ") << "      db_factor ChipGeomFactorsDB  " << db_factor.first << "  "
0708                                   << db_factor.second << "\n";
0709       if (matches(detid, DetId(db_factor.first), DetIdmasks))
0710         ChipGeomFactors[rawid] *= db_factor.second;
0711     }
0712   }
0713 
0714   // piluep scale factors are calculated once per event
0715   // therefore vector index is stored in a map for each module that matches to a db_id
0716   size_t i = 0;
0717   LogDebug("PixelDigitizer ") << " Check PixelEfficiencies -- PUFactors "
0718                               << "\n";
0719   for (const auto& factor : PUFactors) {
0720     //
0721     LogDebug("PixelDigitizer ") << "      factor  " << factor.first << "  " << factor.second.size() << "\n";
0722     for (size_t i = 0, n = factor.second.size(); i < n; i++) {
0723       LogDebug("PixelDigitizer ") << "     print factor.second for " << i << "   " << factor.second[i] << "\n";
0724     }
0725     //
0726     const DetId db_id = DetId(factor.first);
0727     for (const auto& it_module : geom->detUnits()) {
0728       if (dynamic_cast<PixelGeomDetUnit const*>(it_module) == nullptr)
0729         continue;
0730       const DetId detid = it_module->geographicalId();
0731       if (!matches(detid, db_id, DetIdmasks))
0732         continue;
0733       if (iPU.count(detid.rawId())) {
0734         throw cms::Exception("Database")
0735             << "Multiple db_ids match to same module in SiPixelDynamicInefficiency DB Object";
0736       } else {
0737         iPU[detid.rawId()] = i;
0738       }
0739     }
0740     thePUEfficiency.push_back(factor.second);
0741     ++i;
0742   }
0743   pu_scale.resize(thePUEfficiency.size());
0744 }
0745 
0746 bool SiPixelDigitizerAlgorithm::PixelEfficiencies::matches(const DetId& detid,
0747                                                            const DetId& db_id,
0748                                                            const std::vector<uint32_t>& DetIdmasks) {
0749   if (detid.subdetId() != db_id.subdetId())
0750     return false;
0751   for (size_t i = 0; i < DetIdmasks.size(); ++i) {
0752     DetId maskid = DetId(DetIdmasks.at(i));
0753     if (maskid.subdetId() != db_id.subdetId())
0754       continue;
0755     if ((detid.rawId() & maskid.rawId()) != (db_id.rawId() & maskid.rawId()) &&
0756         (db_id.rawId() & maskid.rawId()) != DetId(db_id.det(), db_id.subdetId()).rawId())
0757       return false;
0758   }
0759   return true;
0760 }
0761 
0762 SiPixelDigitizerAlgorithm::PixelAging::PixelAging(const edm::ParameterSet& conf,
0763                                                   bool AddAging,
0764                                                   int NumberOfBarrelLayers,
0765                                                   int NumberOfEndcapDisks) {
0766   // pixel aging
0767   // Don't use Hard coded values, read aging in from python or don't use any
0768   if (AddAging) {
0769     int NumberOfTotLayers = NumberOfBarrelLayers + NumberOfEndcapDisks;
0770     FPixIndex = NumberOfBarrelLayers;
0771 
0772     int i = 0;
0773     thePixelPseudoRadDamage[i++] = conf.getParameter<double>("thePixelPseudoRadDamage_BPix1");
0774     thePixelPseudoRadDamage[i++] = conf.getParameter<double>("thePixelPseudoRadDamage_BPix2");
0775     thePixelPseudoRadDamage[i++] = conf.getParameter<double>("thePixelPseudoRadDamage_BPix3");
0776     thePixelPseudoRadDamage[i++] = conf.getParameter<double>("thePixelPseudoRadDamage_BPix4");
0777 
0778     // to be removed when Gaelle will have the phase2 digitizer
0779     if (NumberOfBarrelLayers >= 5) {
0780       if (NumberOfTotLayers > 20) {
0781         throw cms::Exception("Configuration") << "SiPixelDigitizer was given more layers than it can handle";
0782       }
0783       // For Phase2 tracker layers just set the outermost BPix aging 0.
0784       for (int j = 5; j <= NumberOfBarrelLayers; j++) {
0785         thePixelPseudoRadDamage[j - 1] = 0.;
0786       }
0787     }
0788     //
0789     i = FPixIndex;
0790     thePixelPseudoRadDamage[i++] = conf.getParameter<double>("thePixelPseudoRadDamage_FPix1");
0791     thePixelPseudoRadDamage[i++] = conf.getParameter<double>("thePixelPseudoRadDamage_FPix2");
0792     thePixelPseudoRadDamage[i++] = conf.getParameter<double>("thePixelPseudoRadDamage_FPix3");
0793 
0794     //To be removed when Phase2 digitizer will be available
0795     if (NumberOfEndcapDisks >= 4) {
0796       if (NumberOfTotLayers > 20) {
0797         throw cms::Exception("Configuration") << "SiPixelDigitizer was given more layers than it can handle";
0798       }
0799       // For Phase2 tracker layers just set the extra FPix disk aging to 0. BE CAREFUL THESE VALUES ARE HARDCODED ALSO ELSEWHERE IN THIS FILE
0800       for (int j = 4 + FPixIndex; j <= NumberOfEndcapDisks + NumberOfBarrelLayers; j++) {
0801         thePixelPseudoRadDamage[j - 1] = 0.;
0802       }
0803     }
0804   }
0805   // the first "NumberOfBarrelLayers" settings [0],[1], ... , [NumberOfBarrelLayers-1] are for the barrel pixels
0806   // the next  "NumberOfEndcapDisks"  settings [NumberOfBarrelLayers],[NumberOfBarrelLayers+1], ... [NumberOfEndcapDisks+NumberOfBarrelLayers-1]
0807 }
0808 
0809 //=========================================================================
0810 void SiPixelDigitizerAlgorithm::accumulateSimHits(std::vector<PSimHit>::const_iterator inputBegin,
0811                                                   std::vector<PSimHit>::const_iterator inputEnd,
0812                                                   const size_t inputBeginGlobalIndex,
0813                                                   const unsigned int tofBin,
0814                                                   const PixelGeomDetUnit* pixdet,
0815                                                   const GlobalVector& bfield,
0816                                                   const TrackerTopology* tTopo,
0817                                                   CLHEP::HepRandomEngine* engine) {
0818   // produce SignalPoint's for all SimHit's in detector
0819   // Loop over hits
0820 
0821   uint32_t detId = pixdet->geographicalId().rawId();
0822   size_t simHitGlobalIndex = inputBeginGlobalIndex;  // This needs to stored to create the digi-sim link later
0823   for (std::vector<PSimHit>::const_iterator ssbegin = inputBegin; ssbegin != inputEnd; ++ssbegin, ++simHitGlobalIndex) {
0824     // skip hits not in this detector.
0825     if ((*ssbegin).detUnitId() != detId) {
0826       continue;
0827     }
0828 
0829 #ifdef TP_DEBUG
0830     LogDebug("Pixel Digitizer") << (*ssbegin).particleType() << " " << (*ssbegin).pabs() << " "
0831                                 << (*ssbegin).energyLoss() << " " << (*ssbegin).tof() << " " << (*ssbegin).trackId()
0832                                 << " " << (*ssbegin).processType() << " " << (*ssbegin).detUnitId()
0833                                 << (*ssbegin).entryPoint() << " " << (*ssbegin).exitPoint();
0834 #endif
0835 
0836     std::vector<EnergyDepositUnit> ionization_points;
0837     std::vector<SignalPoint> collection_points;
0838 
0839     // fill collection_points for this SimHit, indpendent of topology
0840     // Check the TOF cut
0841     if (((*ssbegin).tof() - pixdet->surface().toGlobal((*ssbegin).localPosition()).mag() / 30.) >= theTofLowerCut &&
0842         ((*ssbegin).tof() - pixdet->surface().toGlobal((*ssbegin).localPosition()).mag() / 30.) <= theTofUpperCut) {
0843       primary_ionization(*ssbegin, ionization_points, engine);  // fills _ionization_points
0844       drift(*ssbegin,
0845             pixdet,
0846             bfield,
0847             tTopo,
0848             ionization_points,
0849             collection_points);  // transforms _ionization_points to collection_points
0850       // compute induced signal on readout elements and add to _signal
0851       induce_signal(inputBegin,
0852                     inputEnd,
0853                     *ssbegin,
0854                     simHitGlobalIndex,
0855                     inputBeginGlobalIndex,
0856                     tofBin,
0857                     pixdet,
0858                     collection_points);  // 1st 3 args needed only for SimHit<-->Digi link
0859     }  //  end if
0860   }  // end for
0861 }
0862 
0863 //============================================================================
0864 void SiPixelDigitizerAlgorithm::calculateInstlumiFactor(PileupMixingContent* puInfo) {
0865   //Instlumi scalefactor calculating for dynamic inefficiency
0866 
0867   if (puInfo) {
0868     const std::vector<int>& bunchCrossing = puInfo->getMix_bunchCrossing();
0869     const std::vector<float>& TrueInteractionList = puInfo->getMix_TrueInteractions();
0870     //const int bunchSpacing = puInfo->getMix_bunchSpacing();
0871 
0872     int pui = 0, p = 0;
0873     std::vector<int>::const_iterator pu;
0874     std::vector<int>::const_iterator pu0 = bunchCrossing.end();
0875 
0876     for (pu = bunchCrossing.begin(); pu != bunchCrossing.end(); ++pu) {
0877       if (*pu == 0) {
0878         pu0 = pu;
0879         p = pui;
0880       }
0881       pui++;
0882     }
0883     if (pu0 != bunchCrossing.end()) {
0884       for (size_t i = 0, n = pixelEfficiencies_.thePUEfficiency.size(); i < n; i++) {
0885         double instlumi = TrueInteractionList.at(p) * pixelEfficiencies_.theInstLumiScaleFactor;
0886         double instlumi_pow = 1.;
0887         pixelEfficiencies_.pu_scale[i] = 0;
0888         for (size_t j = 0; j < pixelEfficiencies_.thePUEfficiency[i].size(); j++) {
0889           pixelEfficiencies_.pu_scale[i] += instlumi_pow * pixelEfficiencies_.thePUEfficiency[i][j];
0890           instlumi_pow *= instlumi;
0891         }
0892       }
0893     }
0894   } else {
0895     for (int i = 0, n = pixelEfficiencies_.thePUEfficiency.size(); i < n; i++) {
0896       pixelEfficiencies_.pu_scale[i] = 1.;
0897     }
0898   }
0899 }
0900 
0901 //============================================================================
0902 void SiPixelDigitizerAlgorithm::calculateInstlumiFactor(const std::vector<PileupSummaryInfo>& ps, int bunchSpacing) {
0903   int p = -1;
0904   for (unsigned int i = 0; i < ps.size(); i++)
0905     if (ps[i].getBunchCrossing() == 0)
0906       p = i;
0907 
0908   if (p >= 0) {
0909     for (size_t i = 0, n = pixelEfficiencies_.thePUEfficiency.size(); i < n; i++) {
0910       double instlumi = ps[p].getTrueNumInteractions() * pixelEfficiencies_.theInstLumiScaleFactor;
0911       double instlumi_pow = 1.;
0912       pixelEfficiencies_.pu_scale[i] = 0;
0913       for (size_t j = 0; j < pixelEfficiencies_.thePUEfficiency[i].size(); j++) {
0914         pixelEfficiencies_.pu_scale[i] += instlumi_pow * pixelEfficiencies_.thePUEfficiency[i][j];
0915         instlumi_pow *= instlumi;
0916       }
0917     }
0918   } else {
0919     for (int i = 0, n = pixelEfficiencies_.thePUEfficiency.size(); i < n; i++) {
0920       pixelEfficiencies_.pu_scale[i] = 1.;
0921     }
0922   }
0923 }
0924 
0925 // ==========  StuckTBMs
0926 
0927 bool SiPixelDigitizerAlgorithm::killBadFEDChannels() const { return KillBadFEDChannels; }
0928 
0929 std::unique_ptr<PixelFEDChannelCollection> SiPixelDigitizerAlgorithm::chooseScenario(
0930     const std::vector<PileupSummaryInfo>& ps, CLHEP::HepRandomEngine* engine) {
0931   std::unique_ptr<PixelFEDChannelCollection> PixelFEDChannelCollection_ = nullptr;
0932   pixelEfficiencies_.PixelFEDChannelCollection_ = nullptr;
0933 
0934   std::vector<int> bunchCrossing;
0935   std::vector<float> TrueInteractionList;
0936 
0937   for (unsigned int i = 0; i < ps.size(); i++) {
0938     bunchCrossing.push_back(ps[i].getBunchCrossing());
0939     TrueInteractionList.push_back(ps[i].getTrueNumInteractions());
0940   }
0941 
0942   int pui = 0, p = 0;
0943   std::vector<int>::const_iterator pu;
0944   std::vector<int>::const_iterator pu0 = bunchCrossing.end();
0945 
0946   for (pu = bunchCrossing.begin(); pu != bunchCrossing.end(); ++pu) {
0947     if (*pu == 0) {
0948       pu0 = pu;
0949       p = pui;
0950     }
0951     pui++;
0952   }
0953 
0954   if (pu0 != bunchCrossing.end()) {
0955     unsigned int PUBin = TrueInteractionList.at(p);  // case delta PU=1, fix me
0956     const auto& theProbabilitiesPerScenario = scenarioProbability_->getProbabilities(PUBin);
0957     std::vector<double> probabilities;
0958     probabilities.reserve(theProbabilitiesPerScenario.size());
0959     for (auto it = theProbabilitiesPerScenario.begin(); it != theProbabilitiesPerScenario.end(); it++) {
0960       probabilities.push_back(it->second);
0961     }
0962 
0963     CLHEP::RandGeneral randGeneral(*engine, &(probabilities.front()), probabilities.size());
0964     double x = randGeneral.shoot();
0965     unsigned int index = x * probabilities.size() - 1;
0966     const std::string& scenario = theProbabilitiesPerScenario.at(index).first;
0967 
0968     PixelFEDChannelCollection_ = std::make_unique<PixelFEDChannelCollection>(quality_map->at(scenario));
0969     pixelEfficiencies_.PixelFEDChannelCollection_ =
0970         std::make_unique<PixelFEDChannelCollection>(quality_map->at(scenario));
0971   }
0972 
0973   return PixelFEDChannelCollection_;
0974 }
0975 
0976 std::unique_ptr<PixelFEDChannelCollection> SiPixelDigitizerAlgorithm::chooseScenario(PileupMixingContent* puInfo,
0977                                                                                      CLHEP::HepRandomEngine* engine) {
0978   //Determine scenario to use for the current event based on pileup information
0979 
0980   std::unique_ptr<PixelFEDChannelCollection> PixelFEDChannelCollection_ = nullptr;
0981   pixelEfficiencies_.PixelFEDChannelCollection_ = nullptr;
0982   if (puInfo) {
0983     const std::vector<int>& bunchCrossing = puInfo->getMix_bunchCrossing();
0984     const std::vector<float>& TrueInteractionList = puInfo->getMix_TrueInteractions();
0985 
0986     int pui = 0, p = 0;
0987     std::vector<int>::const_iterator pu;
0988     std::vector<int>::const_iterator pu0 = bunchCrossing.end();
0989 
0990     for (pu = bunchCrossing.begin(); pu != bunchCrossing.end(); ++pu) {
0991       if (*pu == 0) {
0992         pu0 = pu;
0993         p = pui;
0994       }
0995       pui++;
0996     }
0997 
0998     if (pu0 != bunchCrossing.end()) {
0999       unsigned int PUBin = TrueInteractionList.at(p);  // case delta PU=1, fix me
1000       const auto& theProbabilitiesPerScenario = scenarioProbability_->getProbabilities(PUBin);
1001       std::vector<double> probabilities;
1002       probabilities.reserve(theProbabilitiesPerScenario.size());
1003       for (auto it = theProbabilitiesPerScenario.begin(); it != theProbabilitiesPerScenario.end(); it++) {
1004         probabilities.push_back(it->second);
1005       }
1006 
1007       CLHEP::RandGeneral randGeneral(*engine, &(probabilities.front()), probabilities.size());
1008       double x = randGeneral.shoot();
1009       unsigned int index = x * probabilities.size() - 1;
1010       const std::string& scenario = theProbabilitiesPerScenario.at(index).first;
1011 
1012       PixelFEDChannelCollection_ = std::make_unique<PixelFEDChannelCollection>(quality_map->at(scenario));
1013       pixelEfficiencies_.PixelFEDChannelCollection_ =
1014           std::make_unique<PixelFEDChannelCollection>(quality_map->at(scenario));
1015     }
1016   }
1017   return PixelFEDChannelCollection_;
1018 }
1019 
1020 //============================================================================
1021 void SiPixelDigitizerAlgorithm::setSimAccumulator(const std::map<uint32_t, std::map<int, int> >& signalMap) {
1022   for (const auto& det : signalMap) {
1023     auto& theSignal = _signal[det.first];
1024     for (const auto& chan : det.second) {
1025       theSignal[chan.first].set(chan.second *
1026                                 theElectronPerADC);  // will get divided again by theElectronPerAdc in digitize...
1027     }
1028   }
1029 }
1030 
1031 //============================================================================
1032 void SiPixelDigitizerAlgorithm::digitize(const PixelGeomDetUnit* pixdet,
1033                                          std::vector<PixelDigi>& digis,
1034                                          std::vector<PixelDigiSimLink>& simlinks,
1035                                          std::vector<PixelDigiAddTempInfo>& newClass_Digi_extra,
1036                                          const TrackerTopology* tTopo,
1037                                          CLHEP::HepRandomEngine* engine) {
1038   // Pixel Efficiency moved from the constructor to this method because
1039   // the information of the det are not available in the constructor
1040   // Efficiency parameters. 0 - no inefficiency, 1-low lumi, 10-high lumi
1041 
1042   uint32_t detID = pixdet->geographicalId().rawId();
1043   const signal_map_type& theSignal = _signal[detID];
1044 
1045   // Noise already defined in electrons
1046   // thePixelThresholdInE = thePixelThreshold * theNoiseInElectrons ;
1047   // Find the threshold in noise units, needed for the noiser.
1048 
1049   float thePixelThresholdInE = 0.;
1050 
1051   if (theNoiseInElectrons > 0.) {
1052     if (pixdet->type().isTrackerPixel() && pixdet->type().isBarrel()) {  // Barrel modules
1053       int lay = tTopo->layer(detID);
1054       if (addThresholdSmearing) {
1055         if (pixdet->subDetector() == GeomDetEnumerators::SubDetector::PixelBarrel ||
1056             pixdet->subDetector() == GeomDetEnumerators::SubDetector::P1PXB) {
1057           if (lay == 1) {
1058             thePixelThresholdInE = CLHEP::RandGaussQ::shoot(
1059                 engine, theThresholdInE_BPix_L1, theThresholdSmearing_BPix_L1);  // gaussian smearing
1060           } else if (lay == 2) {
1061             thePixelThresholdInE = CLHEP::RandGaussQ::shoot(
1062                 engine, theThresholdInE_BPix_L2, theThresholdSmearing_BPix_L2);  // gaussian smearing
1063           } else {
1064             thePixelThresholdInE =
1065                 CLHEP::RandGaussQ::shoot(engine, theThresholdInE_BPix, theThresholdSmearing_BPix);  // gaussian smearing
1066           }
1067         }
1068       } else {
1069         if (pixdet->subDetector() == GeomDetEnumerators::SubDetector::PixelBarrel ||
1070             pixdet->subDetector() == GeomDetEnumerators::SubDetector::P1PXB) {
1071           if (lay == 1) {
1072             thePixelThresholdInE = theThresholdInE_BPix_L1;
1073           } else if (lay == 2) {
1074             thePixelThresholdInE = theThresholdInE_BPix_L2;
1075           } else {
1076             thePixelThresholdInE = theThresholdInE_BPix;  // no smearing
1077           }
1078         }
1079       }
1080     } else if (pixdet->type().isTrackerPixel()) {  // Forward disks modules
1081       if (addThresholdSmearing) {
1082         thePixelThresholdInE =
1083             CLHEP::RandGaussQ::shoot(engine, theThresholdInE_FPix, theThresholdSmearing_FPix);  // gaussian smearing
1084       } else {
1085         thePixelThresholdInE = theThresholdInE_FPix;  // no smearing
1086       }
1087     } else {
1088       throw cms::Exception("NotAPixelGeomDetUnit") << "Not a pixel geomdet unit" << detID;
1089     }
1090   }
1091 
1092 #ifdef TP_DEBUG
1093   const PixelTopology* topol = &pixdet->specificTopology();
1094   int numColumns = topol->ncolumns();  // det module number of cols&rows
1095   int numRows = topol->nrows();
1096   // full detector thickness
1097   float moduleThickness = pixdet->specificSurface().bounds().thickness();
1098   LogDebug("PixelDigitizer") << " PixelDigitizer " << numColumns << " " << numRows << " " << moduleThickness;
1099 #endif
1100 
1101   if (addNoise)
1102     add_noise(pixdet, thePixelThresholdInE / theNoiseInElectrons, engine);  // generate noise
1103 
1104   // Do only if needed
1105 
1106   if ((AddPixelInefficiency) && (!theSignal.empty()))
1107     pixel_inefficiency(pixelEfficiencies_, pixdet, tTopo, engine);  // Kill some pixels
1108 
1109   if (use_ineff_from_db_ && (!theSignal.empty()))
1110     pixel_inefficiency_db(detID);
1111 
1112   if (use_module_killing_) {
1113     if (use_deadmodule_DB_) {  // remove dead modules using DB
1114       module_killing_DB(detID);
1115     } else {  // remove dead modules using the list in cfg file
1116       module_killing_conf(detID);
1117     }
1118   }
1119 
1120   make_digis(thePixelThresholdInE, detID, pixdet, digis, simlinks, newClass_Digi_extra, tTopo);
1121 
1122 #ifdef TP_DEBUG
1123   LogDebug("PixelDigitizer") << "[SiPixelDigitizerAlgorithm] converted " << digis.size() << " PixelDigis in DetUnit"
1124                              << detID;
1125 #endif
1126 }
1127 
1128 //***********************************************************************/
1129 // Generate primary ionization along the track segment.
1130 // Divide the track into small sub-segments
1131 void SiPixelDigitizerAlgorithm::primary_ionization(const PSimHit& hit,
1132                                                    std::vector<EnergyDepositUnit>& ionization_points,
1133                                                    CLHEP::HepRandomEngine* engine) const {
1134   // Straight line approximation for trajectory inside active media
1135 
1136   const float SegmentLength = 0.0010;  //10microns in cm
1137   float energy;
1138 
1139   // Get the 3D segment direction vector
1140   LocalVector direction = hit.exitPoint() - hit.entryPoint();
1141 
1142   float eLoss = hit.energyLoss();  // Eloss in GeV
1143   float length = direction.mag();  // Track length in Silicon
1144 
1145   int NumberOfSegments = int(length / SegmentLength);  // Number of segments
1146   if (NumberOfSegments < 1)
1147     NumberOfSegments = 1;
1148 
1149 #ifdef TP_DEBUG
1150   LogDebug("Pixel Digitizer") << " enter primary_ionzation " << NumberOfSegments
1151                               << " shift = " << (hit.exitPoint().x() - hit.entryPoint().x()) << " "
1152                               << (hit.exitPoint().y() - hit.entryPoint().y()) << " "
1153                               << (hit.exitPoint().z() - hit.entryPoint().z()) << " " << hit.particleType() << " "
1154                               << hit.pabs();
1155 #endif
1156 
1157   float* elossVector = new float[NumberOfSegments];  // Eloss vector
1158 
1159   if (fluctuateCharge) {
1160     //MP DA RIMUOVERE ASSOLUTAMENTE
1161     int pid = hit.particleType();
1162     //int pid=211;  // assume it is a pion
1163 
1164     float momentum = hit.pabs();
1165     // Generate fluctuated charge points
1166     fluctuateEloss(pid, momentum, eLoss, length, NumberOfSegments, elossVector, engine);
1167   }
1168 
1169   ionization_points.resize(NumberOfSegments);  // set size
1170 
1171   // loop over segments
1172   for (int i = 0; i != NumberOfSegments; i++) {
1173     // Divide the segment into equal length subsegments
1174     Local3DPoint point = hit.entryPoint() + float((i + 0.5) / NumberOfSegments) * direction;
1175 
1176     if (fluctuateCharge)
1177       energy = elossVector[i] / GeVperElectron;  // Convert charge to elec.
1178     else
1179       energy = hit.energyLoss() / GeVperElectron / float(NumberOfSegments);
1180 
1181     EnergyDepositUnit edu(energy, point);  //define position,energy point
1182     ionization_points[i] = edu;            // save
1183 
1184 #ifdef TP_DEBUG
1185     LogDebug("Pixel Digitizer") << i << " " << ionization_points[i].x() << " " << ionization_points[i].y() << " "
1186                                 << ionization_points[i].z() << " " << ionization_points[i].energy();
1187 #endif
1188 
1189   }  // end for loop
1190 
1191   delete[] elossVector;
1192 }
1193 //******************************************************************************
1194 
1195 // Fluctuate the charge comming from a small (10um) track segment.
1196 // Use the G4 routine. For mip pions for the moment.
1197 void SiPixelDigitizerAlgorithm::fluctuateEloss(int pid,
1198                                                float particleMomentum,
1199                                                float eloss,
1200                                                float length,
1201                                                int NumberOfSegs,
1202                                                float elossVector[],
1203                                                CLHEP::HepRandomEngine* engine) const {
1204   // Get dedx for this track
1205   //float dedx;
1206   //if( length > 0.) dedx = eloss/length;
1207   //else dedx = eloss;
1208 
1209   double particleMass = 139.6;  // Mass in MeV, Assume pion
1210   pid = std::abs(pid);
1211   if (pid != 211) {  // Mass in MeV
1212     if (pid == 11)
1213       particleMass = 0.511;
1214     else if (pid == 13)
1215       particleMass = 105.7;
1216     else if (pid == 321)
1217       particleMass = 493.7;
1218     else if (pid == 2212)
1219       particleMass = 938.3;
1220   }
1221   // What is the track segment length.
1222   float segmentLength = length / NumberOfSegs;
1223 
1224   // Generate charge fluctuations.
1225   float de = 0.;
1226   float sum = 0.;
1227   double segmentEloss = (1000. * eloss) / NumberOfSegs;  //eloss in MeV
1228   for (int i = 0; i < NumberOfSegs; i++) {
1229     //       material,*,   momentum,energy,*, *,  mass
1230     //myglandz_(14.,segmentLength,2.,2.,dedx,de,0.14);
1231     // The G4 routine needs momentum in MeV, mass in Mev, delta-cut in MeV,
1232     // track segment length in mm, segment eloss in MeV
1233     // Returns fluctuated eloss in MeV
1234     double deltaCutoff = tMax;  // the cutoff is sometimes redefined inside, so fix it.
1235     de = fluctuate->SampleFluctuations(double(particleMomentum * 1000.),
1236                                        particleMass,
1237                                        deltaCutoff,
1238                                        double(segmentLength * 10.),
1239                                        segmentEloss,
1240                                        engine) /
1241          1000.;  //convert to GeV
1242     elossVector[i] = de;
1243     sum += de;
1244   }
1245 
1246   if (sum > 0.) {  // If fluctuations give eloss>0.
1247     // Rescale to the same total eloss
1248     float ratio = eloss / sum;
1249 
1250     for (int ii = 0; ii < NumberOfSegs; ii++)
1251       elossVector[ii] = ratio * elossVector[ii];
1252   } else {  // If fluctuations gives 0 eloss
1253     float averageEloss = eloss / NumberOfSegs;
1254     for (int ii = 0; ii < NumberOfSegs; ii++)
1255       elossVector[ii] = averageEloss;
1256   }
1257   return;
1258 }
1259 
1260 //*******************************************************************************
1261 // Drift the charge segments to the sensor surface (collection plane)
1262 // Include the effect of E-field and B-field
1263 void SiPixelDigitizerAlgorithm::drift(const PSimHit& hit,
1264                                       const PixelGeomDetUnit* pixdet,
1265                                       const GlobalVector& bfield,
1266                                       const TrackerTopology* tTopo,
1267                                       const std::vector<EnergyDepositUnit>& ionization_points,
1268                                       std::vector<SignalPoint>& collection_points) const {
1269 #ifdef TP_DEBUG
1270   LogDebug("Pixel Digitizer") << " enter drift ";
1271 #endif
1272 
1273   collection_points.resize(ionization_points.size());  // set size
1274 
1275   LocalVector driftDir = DriftDirection(pixdet, bfield, hit.detUnitId());  // get the charge drift direction
1276   if (driftDir.z() == 0.) {
1277     LogWarning("Magnetic field") << " pxlx: drift in z is zero ";
1278     return;
1279   }
1280 
1281   // tangent of Lorentz angle
1282   //float TanLorenzAngleX = driftDir.x()/driftDir.z();
1283   //float TanLorenzAngleY = 0.; // force to 0, driftDir.y()/driftDir.z();
1284 
1285   float TanLorenzAngleX, TanLorenzAngleY, dir_z, CosLorenzAngleX, CosLorenzAngleY;
1286   if (alpha2Order) {
1287     TanLorenzAngleX = driftDir.x();  // tangen of Lorentz angle
1288     TanLorenzAngleY = driftDir.y();
1289     dir_z = driftDir.z();                                                 // The z drift direction
1290     CosLorenzAngleX = 1. / sqrt(1. + TanLorenzAngleX * TanLorenzAngleX);  //cosine
1291     CosLorenzAngleY = 1. / sqrt(1. + TanLorenzAngleY * TanLorenzAngleY);  //cosine;
1292 
1293   } else {
1294     TanLorenzAngleX = driftDir.x();
1295     TanLorenzAngleY = 0.;                                                 // force to 0, driftDir.y()/driftDir.z();
1296     dir_z = driftDir.z();                                                 // The z drift direction
1297     CosLorenzAngleX = 1. / sqrt(1. + TanLorenzAngleX * TanLorenzAngleX);  //cosine to estimate the path length
1298     CosLorenzAngleY = 1.;
1299   }
1300 
1301   float moduleThickness = pixdet->specificSurface().bounds().thickness();
1302 #ifdef TP_DEBUG
1303   LogDebug("Pixel Digitizer") << " Lorentz Tan " << TanLorenzAngleX << " " << TanLorenzAngleY << " " << CosLorenzAngleX
1304                               << " " << CosLorenzAngleY << " " << moduleThickness * TanLorenzAngleX << " " << driftDir;
1305 #endif
1306 
1307   float Sigma_x = 1.;  // Charge spread
1308   float Sigma_y = 1.;
1309   float DriftDistance;  // Distance between charge generation and collection
1310   float DriftLength;    // Actual Drift Lentgh
1311   float Sigma;
1312 
1313   for (unsigned int i = 0; i != ionization_points.size(); i++) {
1314     float SegX, SegY, SegZ;  // position
1315     SegX = ionization_points[i].x();
1316     SegY = ionization_points[i].y();
1317     SegZ = ionization_points[i].z();
1318 
1319     // Distance from the collection plane
1320     //DriftDistance = (moduleThickness/2. + SegZ); // Drift to -z
1321     // Include explixitely the E drift direction (for CMS dir_z=-1)
1322     DriftDistance = moduleThickness / 2. - (dir_z * SegZ);  // Drift to -z
1323 
1324     if (DriftDistance <= 0.)
1325       LogDebug("PixelDigitizer ") << " <=0 " << DriftDistance << " " << i << " " << SegZ << " " << dir_z << " " << SegX
1326                                   << " " << SegY << " " << (moduleThickness / 2) << " " << ionization_points[i].energy()
1327                                   << " " << hit.particleType() << " " << hit.pabs() << " " << hit.energyLoss() << " "
1328                                   << hit.entryPoint() << " " << hit.exitPoint() << "\n";
1329 
1330     if (DriftDistance < 0.) {
1331       DriftDistance = 0.;
1332     } else if (DriftDistance > moduleThickness)
1333       DriftDistance = moduleThickness;
1334 
1335     // Assume full depletion now, partial depletion will come later.
1336     float XDriftDueToMagField = DriftDistance * TanLorenzAngleX;
1337     float YDriftDueToMagField = DriftDistance * TanLorenzAngleY;
1338 
1339     // Shift cloud center
1340     float CloudCenterX = SegX + XDriftDueToMagField;
1341     float CloudCenterY = SegY + YDriftDueToMagField;
1342 
1343     // Calculate how long is the charge drift path
1344     DriftLength = sqrt(DriftDistance * DriftDistance + XDriftDueToMagField * XDriftDueToMagField +
1345                        YDriftDueToMagField * YDriftDueToMagField);
1346 
1347     // What is the charge diffusion after this path
1348     Sigma = sqrt(DriftLength / Dist300) * Sigma0;
1349 
1350     // Project the diffusion sigma on the collection plane
1351     Sigma_x = Sigma / CosLorenzAngleX;
1352     Sigma_y = Sigma / CosLorenzAngleY;
1353 
1354     // Insert a charge loss due to Rad Damage here
1355     float energyOnCollector = ionization_points[i].energy();  // The energy that reaches the collector
1356 
1357     // add pixel aging
1358     if (AddPixelAging) {
1359       float kValue = pixel_aging(pixelAging_, pixdet, tTopo);
1360       energyOnCollector *= exp(-1 * kValue * DriftDistance / moduleThickness);
1361     }
1362 
1363 #ifdef TP_DEBUG
1364     LogDebug("Pixel Digitizer") << " Dift DistanceZ= " << DriftDistance << " module thickness= " << moduleThickness
1365                                 << " Start Energy= " << ionization_points[i].energy()
1366                                 << " Energy after loss= " << energyOnCollector;
1367 #endif
1368     SignalPoint sp(CloudCenterX, CloudCenterY, Sigma_x, Sigma_y, hit.tof(), energyOnCollector);
1369 
1370     // Load the Charge distribution parameters
1371     collection_points[i] = (sp);
1372 
1373   }  // loop over ionization points, i.
1374 
1375 }  // end drift
1376 
1377 //*************************************************************************
1378 // Induce the signal on the collection plane of the active sensor area.
1379 void SiPixelDigitizerAlgorithm::induce_signal(std::vector<PSimHit>::const_iterator inputBegin,
1380                                               std::vector<PSimHit>::const_iterator inputEnd,
1381                                               const PSimHit& hit,
1382                                               const size_t hitIndex,
1383                                               const size_t FirstHitIndex,
1384                                               const unsigned int tofBin,
1385                                               const PixelGeomDetUnit* pixdet,
1386                                               const std::vector<SignalPoint>& collection_points) {
1387   // X  - Rows, Left-Right, 160, (1.6cm)   for barrel
1388   // Y  - Columns, Down-Up, 416, (6.4cm)
1389 
1390   const PixelTopology* topol = &pixdet->specificTopology();
1391   uint32_t detID = pixdet->geographicalId().rawId();
1392   signal_map_type& theSignal = _signal[detID];
1393 
1394 #ifdef TP_DEBUG
1395   LogDebug("Pixel Digitizer") << " enter induce_signal, " << topol->pitch().first << " " << topol->pitch().second;  //OK
1396 #endif
1397 
1398   // local map to store pixels hit by 1 Hit.
1399   typedef std::map<int, float, std::less<int> > hit_map_type;
1400   hit_map_type hit_signal;
1401 
1402   // map to store pixel integrals in the x and in the y directions
1403   std::map<int, float, std::less<int> > x, y;
1404 
1405   // Assign signals to readout channels and store sorted by channel number
1406 
1407   // Iterate over collection points on the collection plane
1408   for (std::vector<SignalPoint>::const_iterator i = collection_points.begin(); i != collection_points.end(); ++i) {
1409     float CloudCenterX = i->position().x();  // Charge position in x
1410     float CloudCenterY = i->position().y();  //                 in y
1411     float SigmaX = i->sigma_x();             // Charge spread in x
1412     float SigmaY = i->sigma_y();             //               in y
1413     float Charge = i->amplitude();           // Charge amplitude
1414 
1415     if (SigmaX == 0 || SigmaY == 0) {
1416       LogDebug("Pixel Digitizer") << SigmaX << " " << SigmaY << " cloud " << i->position().x() << " "
1417                                   << i->position().y() << " " << i->sigma_x() << " " << i->sigma_y() << " "
1418                                   << i->amplitude() << "\n";
1419     }
1420 
1421 #ifdef TP_DEBUG
1422     LogDebug("Pixel Digitizer") << " cloud " << i->position().x() << " " << i->position().y() << " " << i->sigma_x()
1423                                 << " " << i->sigma_y() << " " << i->amplitude();
1424 #endif
1425 
1426     // Find the maximum cloud spread in 2D plane , assume 3*sigma
1427     float CloudRight = CloudCenterX + ClusterWidth * SigmaX;
1428     float CloudLeft = CloudCenterX - ClusterWidth * SigmaX;
1429     float CloudUp = CloudCenterY + ClusterWidth * SigmaY;
1430     float CloudDown = CloudCenterY - ClusterWidth * SigmaY;
1431 
1432     // Define 2D cloud limit points
1433     LocalPoint PointRightUp = LocalPoint(CloudRight, CloudUp);
1434     LocalPoint PointLeftDown = LocalPoint(CloudLeft, CloudDown);
1435 
1436     // This points can be located outside the sensor area.
1437     // The conversion to measurement point does not check for that
1438     // so the returned pixel index might be wrong (outside range).
1439     // We rely on the limits check below to fix this.
1440     // But remember whatever we do here THE CHARGE OUTSIDE THE ACTIVE
1441     // PIXEL AREA IS LOST, it should not be collected.
1442 
1443     // Convert the 2D points to pixel indices
1444     MeasurementPoint mp = topol->measurementPosition(PointRightUp);  //OK
1445 
1446     int IPixRightUpX = int(floor(mp.x()));
1447     int IPixRightUpY = int(floor(mp.y()));
1448 
1449 #ifdef TP_DEBUG
1450     LogDebug("Pixel Digitizer") << " right-up " << PointRightUp << " " << mp.x() << " " << mp.y() << " " << IPixRightUpX
1451                                 << " " << IPixRightUpY;
1452 #endif
1453 
1454     mp = topol->measurementPosition(PointLeftDown);  //OK
1455 
1456     int IPixLeftDownX = int(floor(mp.x()));
1457     int IPixLeftDownY = int(floor(mp.y()));
1458 
1459 #ifdef TP_DEBUG
1460     LogDebug("Pixel Digitizer") << " left-down " << PointLeftDown << " " << mp.x() << " " << mp.y() << " "
1461                                 << IPixLeftDownX << " " << IPixLeftDownY;
1462 #endif
1463 
1464     // Check detector limits to correct for pixels outside range.
1465     int numColumns = topol->ncolumns();  // det module number of cols&rows
1466     int numRows = topol->nrows();
1467 
1468     IPixRightUpX = numRows > IPixRightUpX ? IPixRightUpX : numRows - 1;
1469     IPixRightUpY = numColumns > IPixRightUpY ? IPixRightUpY : numColumns - 1;
1470     IPixLeftDownX = 0 < IPixLeftDownX ? IPixLeftDownX : 0;
1471     IPixLeftDownY = 0 < IPixLeftDownY ? IPixLeftDownY : 0;
1472 
1473     x.clear();  // clear temporary integration array
1474     y.clear();
1475 
1476     // First integrate charge strips in x
1477     int ix;                                               // TT for compatibility
1478     for (ix = IPixLeftDownX; ix <= IPixRightUpX; ix++) {  // loop over x index
1479       float xUB, xLB, UpperBound, LowerBound;
1480 
1481       // Why is set to 0 if ix=0, does it meen that we accept charge
1482       // outside the sensor? CHeck How it was done in ORCA?
1483       //if(ix == 0) LowerBound = 0.;
1484       if (ix == 0 || SigmaX == 0.)  // skip for surface segemnts
1485         LowerBound = 0.;
1486       else {
1487         mp = MeasurementPoint(float(ix), 0.0);
1488         xLB = topol->localPosition(mp).x();
1489         LowerBound = 1 - calcQ((xLB - CloudCenterX) / SigmaX);
1490       }
1491 
1492       if (ix == numRows - 1 || SigmaX == 0.)
1493         UpperBound = 1.;
1494       else {
1495         mp = MeasurementPoint(float(ix + 1), 0.0);
1496         xUB = topol->localPosition(mp).x();
1497         UpperBound = 1. - calcQ((xUB - CloudCenterX) / SigmaX);
1498       }
1499 
1500       float TotalIntegrationRange = UpperBound - LowerBound;  // get strip
1501       x[ix] = TotalIntegrationRange;                          // save strip integral
1502       if (SigmaX == 0 || SigmaY == 0)
1503         LogDebug("Pixel Digitizer") << TotalIntegrationRange << " " << ix << "\n";
1504     }
1505 
1506     // Now integrate strips in y
1507     int iy;                                               // TT for compatibility
1508     for (iy = IPixLeftDownY; iy <= IPixRightUpY; iy++) {  //loope over y ind
1509       float yUB, yLB, UpperBound, LowerBound;
1510 
1511       if (iy == 0 || SigmaY == 0.)
1512         LowerBound = 0.;
1513       else {
1514         mp = MeasurementPoint(0.0, float(iy));
1515         yLB = topol->localPosition(mp).y();
1516         LowerBound = 1. - calcQ((yLB - CloudCenterY) / SigmaY);
1517       }
1518 
1519       if (iy == numColumns - 1 || SigmaY == 0.)
1520         UpperBound = 1.;
1521       else {
1522         mp = MeasurementPoint(0.0, float(iy + 1));
1523         yUB = topol->localPosition(mp).y();
1524         UpperBound = 1. - calcQ((yUB - CloudCenterY) / SigmaY);
1525       }
1526 
1527       float TotalIntegrationRange = UpperBound - LowerBound;
1528       y[iy] = TotalIntegrationRange;  // save strip integral
1529       if (SigmaX == 0 || SigmaY == 0)
1530         LogDebug("Pixel Digitizer") << TotalIntegrationRange << " " << iy << "\n";
1531     }
1532 
1533     // Get the 2D charge integrals by folding x and y strips
1534     int chan;
1535     for (ix = IPixLeftDownX; ix <= IPixRightUpX; ix++) {    // loop over x index
1536       for (iy = IPixLeftDownY; iy <= IPixRightUpY; iy++) {  //loope over y ind
1537 
1538         float ChargeFraction = Charge * x[ix] * y[iy];
1539 
1540         if (ChargeFraction > 0.) {
1541           chan = PixelDigi::pixelToChannel(ix, iy);  // Get index
1542           // Load the amplitude
1543           hit_signal[chan] += ChargeFraction;
1544         }  // endif
1545 
1546 #ifdef TP_DEBUG
1547         mp = MeasurementPoint(float(ix), float(iy));
1548         LocalPoint lp = topol->localPosition(mp);
1549         chan = topol->channel(lp);
1550         LogDebug("Pixel Digitizer") << " pixel " << ix << " " << iy << " - "
1551                                     << " " << chan << " " << ChargeFraction << " " << mp.x() << " " << mp.y() << " "
1552                                     << lp.x() << " " << lp.y() << " "  // givex edge position
1553                                     << chan;                           // edge belongs to previous ?
1554 #endif
1555 
1556       }  // endfor iy
1557     }  //endfor ix
1558 
1559   }  // loop over charge distributions
1560 
1561   // Fill the global map with all hit pixels from this event
1562 
1563   bool reweighted = false;
1564   bool makeDSLinks = store_SimHitEntryExitPoints_ || makeDigiSimLinks_;
1565 
1566   size_t ReferenceIndex4CR = 0;
1567   if (UseReweighting) {
1568     if (hit.processType() == 0) {
1569       ReferenceIndex4CR = hitIndex;
1570       reweighted = TheNewSiPixelChargeReweightingAlgorithmClass->hitSignalReweight<digitizerUtility::Amplitude>(
1571           hit, hit_signal, hitIndex, ReferenceIndex4CR, tofBin, topol, detID, theSignal, hit.processType(), makeDSLinks);
1572     } else {
1573       std::vector<PSimHit>::const_iterator crSimHit = inputBegin;
1574       ReferenceIndex4CR = FirstHitIndex;
1575       // if the first hit in the same detId is not associated to the same trackId, try to find a better match
1576       if ((*inputBegin).trackId() != hit.trackId()) {
1577         // loop over all the hit from the 1st in the same detId to the hit itself to find the primary particle of the same trackId
1578         uint32_t detId = pixdet->geographicalId().rawId();
1579         size_t localIndex = FirstHitIndex;
1580         for (std::vector<PSimHit>::const_iterator ssbegin = inputBegin; localIndex < hitIndex;
1581              ++ssbegin, ++localIndex) {
1582           if ((*ssbegin).detUnitId() != detId) {
1583             continue;
1584           }
1585           if ((*ssbegin).trackId() == hit.trackId() && (*ssbegin).processType() == 0) {
1586             crSimHit = ssbegin;
1587             ReferenceIndex4CR = localIndex;
1588             break;
1589           }
1590         }
1591       }
1592 
1593       reweighted = TheNewSiPixelChargeReweightingAlgorithmClass->hitSignalReweight<digitizerUtility::Amplitude>(
1594           (*crSimHit),
1595           hit_signal,
1596           hitIndex,
1597           ReferenceIndex4CR,
1598           tofBin,
1599           topol,
1600           detID,
1601           theSignal,
1602           hit.processType(),
1603           makeDSLinks);
1604     }
1605   }
1606   if (!reweighted) {
1607     for (hit_map_type::const_iterator im = hit_signal.begin(); im != hit_signal.end(); ++im) {
1608       int chan = (*im).first;
1609       if (ReferenceIndex4CR == 0) {
1610         // no determination has been done previously because !UseReweighting
1611         // we need to determine it now:
1612         if (hit.processType() == 0)
1613           ReferenceIndex4CR = hitIndex;
1614         else {
1615           ReferenceIndex4CR = FirstHitIndex;
1616           // if the first hit in the same detId is not associated to the same trackId, try to find a better match
1617           if ((*inputBegin).trackId() != hit.trackId()) {
1618             // loop on all the hit from the 1st of the collection to the hit itself to find the Primary particle of the same trackId
1619             uint32_t detId = pixdet->geographicalId().rawId();
1620             size_t localIndex = FirstHitIndex;
1621             for (std::vector<PSimHit>::const_iterator ssbegin = inputBegin; localIndex < hitIndex;
1622                  ++ssbegin, ++localIndex) {
1623               if ((*ssbegin).detUnitId() != detId) {
1624                 continue;
1625               }
1626               if ((*ssbegin).trackId() == hit.trackId() && (*ssbegin).processType() == 0) {
1627                 ReferenceIndex4CR = localIndex;
1628                 break;
1629               }
1630             }
1631           }
1632         }
1633       }
1634       theSignal[chan] += (makeDSLinks ? digitizerUtility::Amplitude(
1635                                             (*im).second, &hit, hitIndex, ReferenceIndex4CR, tofBin, (*im).second)
1636                                       : digitizerUtility::Amplitude((*im).second, (*im).second));
1637 
1638 #ifdef TP_DEBUG
1639       std::pair<int, int> ip = PixelDigi::channelToPixel(chan);
1640       LogDebug("Pixel Digitizer") << " pixel " << ip.first << " " << ip.second << " " << theSignal[chan];
1641 #endif
1642     }
1643   }
1644 
1645 }  // end induce_signal
1646 
1647 /***********************************************************************/
1648 
1649 void SiPixelDigitizerAlgorithm::fillSimHitMaps(std::vector<PSimHit> simHits, const unsigned int tofBin) {
1650   // store here the SimHit map for later
1651   int printnum = 0;
1652   for (std::vector<PSimHit>::const_iterator it = simHits.begin(), itEnd = simHits.end(); it != itEnd;
1653        ++it, ++printnum) {
1654     unsigned int detID = (*it).detUnitId();
1655     unsigned int subdetID = DetId(detID).subdetId();
1656     subDetTofBin theSubDetTofBin = std::make_pair(subdetID, tofBin);
1657     SimHitMap[SimHitCollMap[theSubDetTofBin]].push_back(*it);
1658   }
1659 }
1660 
1661 void SiPixelDigitizerAlgorithm::resetSimHitMaps() { SimHitMap.clear(); }
1662 
1663 /***********************************************************************/
1664 
1665 // Build pixels, check threshold, add misscalibration, ...
1666 void SiPixelDigitizerAlgorithm::make_digis(float thePixelThresholdInE,
1667                                            uint32_t detID,
1668                                            const PixelGeomDetUnit* pixdet,
1669                                            std::vector<PixelDigi>& digis,
1670                                            std::vector<PixelDigiSimLink>& simlinks,
1671                                            std::vector<PixelDigiAddTempInfo>& newClass_Digi_extra,
1672                                            const TrackerTopology* tTopo) const {
1673 #ifdef TP_DEBUG
1674   LogDebug("Pixel Digitizer") << " make digis "
1675                               << " "
1676                               << " pixel threshold FPix" << theThresholdInE_FPix << " "
1677                               << " pixel threshold BPix" << theThresholdInE_BPix << " "
1678                               << " pixel threshold BPix Layer1" << theThresholdInE_BPix_L1 << " "
1679                               << " pixel threshold BPix Layer2" << theThresholdInE_BPix_L2 << " "
1680                               << " List pixels passing threshold ";
1681 #endif
1682 
1683   // Loop over hit pixels
1684 
1685   signalMaps::const_iterator it = _signal.find(detID);
1686   if (it == _signal.end()) {
1687     return;
1688   }
1689 
1690   const signal_map_type& theSignal = (*it).second;
1691 
1692   // unsigned long is enough to store SimTrack id and EncodedEventId
1693   using TrackEventId = std::pair<decltype(SimTrack().trackId()), decltype(EncodedEventId().rawId())>;
1694   std::map<TrackEventId, float> simi;  // re-used
1695 
1696   for (signal_map_const_iterator i = theSignal.begin(); i != theSignal.end(); ++i) {
1697     float signalInElectrons = (*i).second;  // signal in electrons
1698 
1699     // Do the miss calibration for calibration studies only.
1700     //if(doMissCalibrate) signalInElectrons = missCalibrate(signalInElectrons)
1701 
1702     // Do only for pixels above threshold
1703 
1704     if (signalInElectrons >= thePixelThresholdInE &&
1705         signalInElectrons > 0.) {  // check threshold, always reject killed (0-charge) digis
1706 
1707       int chan = (*i).first;  // channel number
1708       std::pair<int, int> ip = PixelDigi::channelToPixel(chan);
1709       int adc = 0;  // ADC count as integer
1710 
1711       // Do the miss calibration for calibration studies only.
1712       if (doMissCalibrate) {
1713         int row = ip.first;                                                           // X in row
1714         int col = ip.second;                                                          // Y is in col
1715         adc = int(missCalibrate(detID, tTopo, pixdet, col, row, signalInElectrons));  //full misscalib.
1716       } else {                                             // Just do a simple electron->adc conversion
1717         adc = int(signalInElectrons / theElectronPerADC);  // calibrate gain
1718       }
1719       adc = std::min(adc, theAdcFullScale);  // Check maximum value
1720 #ifdef TP_DEBUG
1721       LogDebug("Pixel Digitizer") << (*i).first << " " << (*i).second << " " << signalInElectrons << " " << adc
1722                                   << ip.first << " " << ip.second;
1723 #endif
1724 
1725       // Load digis
1726       digis.emplace_back(ip.first, ip.second, adc);
1727 
1728       if (makeDigiSimLinks_ && !(*i).second.hitInfos().empty()) {
1729         //digilink
1730         unsigned int il = 0;
1731         for (const auto& info : (*i).second.hitInfos()) {
1732           // note: according to C++ standard operator[] does
1733           // value-initializiation, which for float means initial value of 0
1734           simi[std::make_pair(info.trackId(), info.eventId().rawId())] += (*i).second.individualampl()[il];
1735           il++;
1736         }
1737 
1738         //sum the contribution of the same trackid
1739         for (const auto& info : (*i).second.hitInfos()) {
1740           // skip if track already processed
1741           auto found = simi.find(std::make_pair(info.trackId(), info.eventId().rawId()));
1742           if (found == simi.end())
1743             continue;
1744 
1745           float sum_samechannel = found->second;
1746           float fraction = sum_samechannel / (*i).second;
1747           if (fraction > 1.f)
1748             fraction = 1.f;
1749 
1750           // Approximation: pick hitIndex and tofBin only from the first SimHit
1751           simlinks.emplace_back((*i).first, info.trackId(), info.hitIndex(), info.tofBin(), info.eventId(), fraction);
1752           simi.erase(found);
1753         }
1754         simi.clear();  // although should be empty already
1755       }
1756 
1757       if (store_SimHitEntryExitPoints_ && !(*i).second.hitInfos().empty()) {
1758         // get info stored, like in simlinks...
1759         for (const auto& info : (*i).second.hitInfos()) {
1760           unsigned int CFPostoBeUsed = info.hitIndex4ChargeRew();
1761           // check if the association (chan, index) is already in the newClass_Digi_extra collection
1762           // if yes, don't push a duplicated entry ; if not, push a new entry
1763           std::vector<PixelDigiAddTempInfo>::iterator loopNewClass;
1764           bool already_present = false;
1765           for (loopNewClass = newClass_Digi_extra.begin(); loopNewClass != newClass_Digi_extra.end(); ++loopNewClass) {
1766             if (chan == (int)loopNewClass->channel() && CFPostoBeUsed == loopNewClass->hitIndex()) {
1767               already_present = true;
1768               loopNewClass->addCharge(info.getAmpl());
1769             }
1770           }
1771           if (!already_present) {
1772             unsigned int tofBin = info.tofBin();
1773             // then inspired by https://github.com/cms-sw/cmssw/blob/master/SimTracker/TrackerHitAssociation/src/TrackerHitAssociator.cc#L566 :
1774             subDetTofBin theSubDetTofBin = std::make_pair(DetId(detID).subdetId(), tofBin);
1775             auto it = SimHitCollMap.find(theSubDetTofBin);
1776             if (it != SimHitCollMap.end()) {
1777               auto it2 = SimHitMap.find((it->second));
1778 
1779               if (it2 != SimHitMap.end()) {
1780                 const PSimHit& theSimHit = (it2->second)[CFPostoBeUsed];
1781                 newClass_Digi_extra.emplace_back(chan,
1782                                                  info.hitIndex4ChargeRew(),
1783                                                  theSimHit.entryPoint(),
1784                                                  theSimHit.exitPoint(),
1785                                                  theSimHit.processType(),
1786                                                  theSimHit.trackId(),
1787                                                  detID,
1788                                                  info.getAmpl());
1789               }
1790             }
1791           }
1792         }  // end for
1793       }  // end if store_SimHitEntryExitPoints_
1794     }
1795   }
1796 }
1797 
1798 /***********************************************************************/
1799 
1800 //  Add electronic noise to pixel charge
1801 void SiPixelDigitizerAlgorithm::add_noise(const PixelGeomDetUnit* pixdet,
1802                                           float thePixelThreshold,
1803                                           CLHEP::HepRandomEngine* engine) {
1804 #ifdef TP_DEBUG
1805   LogDebug("Pixel Digitizer") << " enter add_noise " << theNoiseInElectrons;
1806 #endif
1807 
1808   uint32_t detID = pixdet->geographicalId().rawId();
1809   signal_map_type& theSignal = _signal[detID];
1810 
1811   // First add noise to hit pixels
1812   float theSmearedChargeRMS = 0.0;
1813 
1814   if (!theSignal.empty()) {
1815     for (signal_map_iterator i = theSignal.begin(); i != theSignal.end(); i++) {
1816       if (addChargeVCALSmearing) {
1817         if ((*i).second < 3000) {
1818           theSmearedChargeRMS = 543.6 - (*i).second * 0.093;
1819         } else if ((*i).second < 6000) {
1820           theSmearedChargeRMS = 307.6 - (*i).second * 0.01;
1821         } else {
1822           theSmearedChargeRMS = -432.4 + (*i).second * 0.123;
1823         }
1824 
1825         // Noise from Vcal smearing:
1826         float noise_ChargeVCALSmearing = theSmearedChargeRMS * CLHEP::RandGaussQ::shoot(engine, 0., 1.);
1827         // Noise from full readout:
1828         float noise = CLHEP::RandGaussQ::shoot(engine, 0., theReadoutNoise);
1829 
1830         if (((*i).second + digitizerUtility::Amplitude(noise + noise_ChargeVCALSmearing, -1.)) < 0.) {
1831           (*i).second.set(0);
1832         } else {
1833           (*i).second += digitizerUtility::Amplitude(noise + noise_ChargeVCALSmearing, -1.);
1834         }
1835 
1836       }  // End if addChargeVCalSmearing
1837       else {
1838         // Noise: ONLY full READOUT Noise.
1839         // Use here the FULL readout noise, including TBM,ALT,AOH,OPT-REC.
1840         float noise = CLHEP::RandGaussQ::shoot(engine, 0., theReadoutNoise);
1841 
1842         if (((*i).second + digitizerUtility::Amplitude(noise, -1.)) < 0.) {
1843           (*i).second.set(0);
1844         } else {
1845           (*i).second += digitizerUtility::Amplitude(noise, -1.);
1846         }
1847       }  // end if only Noise from full readout
1848     }
1849   }
1850   if (!addNoisyPixels)  // Option to skip noise in non-hit pixels
1851     return;
1852 
1853   const PixelTopology* topol = &pixdet->specificTopology();
1854   int numColumns = topol->ncolumns();  // det module number of cols&rows
1855   int numRows = topol->nrows();
1856 
1857   // Add noise on non-hit pixels
1858   // Use here the pixel noise
1859   int numberOfPixels = (numRows * numColumns);
1860   std::map<int, float, std::less<int> > otherPixels;
1861   std::map<int, float, std::less<int> >::iterator mapI;
1862 
1863   theNoiser->generate(numberOfPixels,
1864                       thePixelThreshold,    //thr. in un. of nois
1865                       theNoiseInElectrons,  // noise in elec.
1866                       otherPixels,
1867                       engine);
1868 
1869 #ifdef TP_DEBUG
1870   LogDebug("Pixel Digitizer") << " Add noisy pixels " << numRows << " " << numColumns << " " << theNoiseInElectrons
1871                               << " " << theThresholdInE_FPix << theThresholdInE_BPix << " " << numberOfPixels << " "
1872                               << otherPixels.size();
1873 #endif
1874 
1875   // Add noisy pixels
1876   if (otherPixels.empty()) {
1877     return;
1878   }
1879   for (mapI = otherPixels.begin(); mapI != otherPixels.end(); mapI++) {
1880     int iy = ((*mapI).first) / numRows;
1881     int ix = ((*mapI).first) - (iy * numRows);
1882 
1883     // Keep for a while for testing.
1884     if (iy < 0 || iy > (numColumns - 1))
1885       LogWarning("Pixel Geometry") << " error in iy " << iy;
1886     if (ix < 0 || ix > (numRows - 1))
1887       LogWarning("Pixel Geometry") << " error in ix " << ix;
1888 
1889     int chan = PixelDigi::pixelToChannel(ix, iy);
1890 
1891 #ifdef TP_DEBUG
1892     LogDebug("Pixel Digitizer") << " Storing noise = " << (*mapI).first << " " << (*mapI).second << " " << ix << " "
1893                                 << iy << " " << chan;
1894 #endif
1895 
1896     if (theSignal[chan] == 0) {
1897       int inoise = float((*mapI).second);
1898       float noise = (float)inoise;
1899       theSignal[chan] = digitizerUtility::Amplitude(noise, -1.);
1900     }
1901   }
1902 }
1903 
1904 /***********************************************************************/
1905 
1906 // Simulate the readout inefficiencies.
1907 // Delete a selected number of single pixels, dcols and rocs.
1908 void SiPixelDigitizerAlgorithm::pixel_inefficiency(const PixelEfficiencies& eff,
1909                                                    const PixelGeomDetUnit* pixdet,
1910                                                    const TrackerTopology* tTopo,
1911                                                    CLHEP::HepRandomEngine* engine) {
1912   uint32_t detID = pixdet->geographicalId().rawId();
1913   signal_map_type& theSignal = _signal[detID];
1914   const PixelTopology* topol = &pixdet->specificTopology();
1915   int numColumns = topol->ncolumns();  // det module number of cols&rows
1916   int numRows = topol->nrows();
1917   bool isPhase1 = pixdet->subDetector() == GeomDetEnumerators::SubDetector::P1PXB ||
1918                   pixdet->subDetector() == GeomDetEnumerators::SubDetector::P1PXEC;
1919   // Predefined efficiencies
1920   double pixelEfficiency = 1.0;
1921   double columnEfficiency = 1.0;
1922   double chipEfficiency = 1.0;
1923   std::vector<double> pixelEfficiencyROCStdPixels(16, 1);
1924   std::vector<double> pixelEfficiencyROCBigPixels(16, 1);
1925 
1926   auto pIndexConverter = PixelIndices(numColumns, numRows);
1927 
1928   std::vector<int> badRocsFromFEDChannels(16, 0);
1929   if (eff.PixelFEDChannelCollection_ != nullptr) {
1930     PixelFEDChannelCollection::const_iterator it = eff.PixelFEDChannelCollection_->find(detID);
1931 
1932     if (it != eff.PixelFEDChannelCollection_->end()) {
1933       const std::vector<CablingPathToDetUnit>& path = map_->pathToDetUnit(detID);
1934       for (const auto& ch : *it) {
1935         for (unsigned int i_roc = ch.roc_first; i_roc <= ch.roc_last; ++i_roc) {
1936           for (const auto p : path) {
1937             const PixelROC* myroc = map_->findItem(p);
1938             if (myroc->idInDetUnit() == static_cast<unsigned int>(i_roc)) {
1939               LocalPixel::RocRowCol local = {39, 25};  //corresponding to center of ROC row,col
1940               GlobalPixel global = myroc->toGlobal(LocalPixel(local));
1941               int chipIndex(0), colROC(0), rowROC(0);
1942               pIndexConverter.transformToROC(global.col, global.row, chipIndex, colROC, rowROC);
1943               badRocsFromFEDChannels.at(chipIndex) = 1;
1944             }
1945           }
1946         }
1947       }  // loop over channels
1948     }  // detID in PixelFEDChannelCollection_
1949   }  // has PixelFEDChannelCollection_
1950 
1951   if (eff.FromConfig) {
1952     // setup the chip indices conversion
1953     if (pixdet->subDetector() == GeomDetEnumerators::SubDetector::PixelBarrel ||
1954         pixdet->subDetector() == GeomDetEnumerators::SubDetector::P1PXB) {  // barrel layers
1955       int layerIndex = tTopo->layer(detID);
1956       pixelEfficiency = eff.thePixelEfficiency[layerIndex - 1];
1957       columnEfficiency = eff.thePixelColEfficiency[layerIndex - 1];
1958       chipEfficiency = eff.thePixelChipEfficiency[layerIndex - 1];
1959       LogDebug("Pixel Digitizer") << "Using BPix columnEfficiency = " << columnEfficiency
1960                                   << " for layer = " << layerIndex << "\n";
1961       // This should never happen, but only check if it is not an upgrade geometry
1962       if (NumberOfBarrelLayers == 3) {
1963         if (numColumns > 416)
1964           LogWarning("Pixel Geometry") << " wrong columns in barrel " << numColumns;
1965         if (numRows > 160)
1966           LogWarning("Pixel Geometry") << " wrong rows in barrel " << numRows;
1967 
1968         int ladder = tTopo->pxbLadder(detID);
1969         int module = tTopo->pxbModule(detID);
1970         if (module <= 4)
1971           module = 5 - module;
1972         else
1973           module -= 4;
1974 
1975         columnEfficiency *= eff.theLadderEfficiency_BPix[layerIndex - 1][ladder - 1] *
1976                             eff.theModuleEfficiency_BPix[layerIndex - 1][module - 1] * eff.pu_scale[layerIndex - 1];
1977       }
1978     } else if (pixdet->subDetector() == GeomDetEnumerators::SubDetector::PixelEndcap ||
1979                pixdet->subDetector() == GeomDetEnumerators::SubDetector::P1PXEC ||
1980                pixdet->subDetector() == GeomDetEnumerators::SubDetector::P2PXEC) {  // forward disks
1981 
1982       unsigned int diskIndex =
1983           tTopo->layer(detID) + eff.FPixIndex;  // Use diskIndex-1 later to stay consistent with BPix
1984       unsigned int panelIndex = tTopo->pxfPanel(detID);
1985       unsigned int moduleIndex = tTopo->pxfModule(detID);
1986       //if (eff.FPixIndex>diskIndex-1){throw cms::Exception("Configuration") <<"SiPixelDigitizer is using the wrong efficiency value. index = "
1987       //                                                                       <<diskIndex-1<<" , MinIndex = "<<eff.FPixIndex<<" ... "<<tTopo->pxfDisk(detID);}
1988       pixelEfficiency = eff.thePixelEfficiency[diskIndex - 1];
1989       columnEfficiency = eff.thePixelColEfficiency[diskIndex - 1];
1990       chipEfficiency = eff.thePixelChipEfficiency[diskIndex - 1];
1991       LogDebug("Pixel Digitizer") << "Using FPix columnEfficiency = " << columnEfficiency
1992                                   << " for Disk = " << tTopo->pxfDisk(detID) << "\n";
1993       // Sometimes the forward pixels have wrong size,
1994       // this crashes the index conversion, so exit, but only check if it is not an upgrade geometry
1995       if (NumberOfBarrelLayers ==
1996           3) {  // whether it is the present or the phase 1 detector can be checked using GeomDetEnumerators::SubDetector
1997         if (numColumns > 260 || numRows > 160) {
1998           if (numColumns > 260)
1999             LogWarning("Pixel Geometry") << " wrong columns in endcaps " << numColumns;
2000           if (numRows > 160)
2001             LogWarning("Pixel Geometry") << " wrong rows in endcaps " << numRows;
2002           return;
2003         }
2004         if ((panelIndex == 1 && (moduleIndex == 1 || moduleIndex == 2)) ||
2005             (panelIndex == 2 && moduleIndex == 1)) {  //inner modules
2006           columnEfficiency *= eff.theInnerEfficiency_FPix[diskIndex - 1] * eff.pu_scale[3];
2007         } else {  //outer modules
2008           columnEfficiency *= eff.theOuterEfficiency_FPix[diskIndex - 1] * eff.pu_scale[4];
2009         }
2010       }  // current detector, forward
2011     } else if (pixdet->subDetector() == GeomDetEnumerators::SubDetector::P2OTB ||
2012                pixdet->subDetector() == GeomDetEnumerators::SubDetector::P2OTEC) {
2013       // If phase 2 outer tracker, hardcoded values as they have been so far
2014       pixelEfficiency = 0.999;
2015       columnEfficiency = 0.999;
2016       chipEfficiency = 0.999;
2017     }  // if barrel/forward
2018   } else {  // Load precomputed factors from Database
2019     pixelEfficiency = eff.PixelGeomFactors.at(detID);
2020     columnEfficiency = eff.ColGeomFactors.at(detID) * eff.pu_scale[eff.iPU.at(detID)];
2021     chipEfficiency = eff.ChipGeomFactors.at(detID);
2022     if (isPhase1) {
2023       for (unsigned int i_roc = 0; i_roc < eff.PixelGeomFactorsROCStdPixels.at(detID).size(); ++i_roc) {
2024         pixelEfficiencyROCStdPixels[i_roc] = eff.PixelGeomFactorsROCStdPixels.at(detID).at(i_roc);
2025         pixelEfficiencyROCBigPixels[i_roc] = eff.PixelGeomFactorsROCBigPixels.at(detID).at(i_roc);
2026       }
2027     }  // is Phase 1
2028   }
2029 
2030 #ifdef TP_DEBUG
2031   LogDebug("Pixel Digitizer") << " enter pixel_inefficiency " << pixelEfficiency << " " << columnEfficiency << " "
2032                               << chipEfficiency;
2033 #endif
2034 
2035   // Initilize the index converter
2036   //PixelIndices indexConverter(numColumns,numRows);
2037 
2038   int chipIndex = 0;
2039   int rowROC = 0;
2040   int colROC = 0;
2041   std::map<int, int, std::less<int> > chips, columns, pixelStd, pixelBig;
2042   std::map<int, int, std::less<int> >::iterator iter;
2043 
2044   // Find out the number of columns and rocs hits
2045   // Loop over hit pixels, amplitude in electrons, channel = coded row,col
2046   for (signal_map_const_iterator i = theSignal.begin(); i != theSignal.end(); ++i) {
2047     int chan = i->first;
2048     std::pair<int, int> ip = PixelDigi::channelToPixel(chan);
2049     int row = ip.first;   // X in row
2050     int col = ip.second;  // Y is in col
2051     //transform to ROC index coordinates
2052     pIndexConverter.transformToROC(col, row, chipIndex, colROC, rowROC);
2053     int dColInChip = pIndexConverter.DColumn(colROC);  // get ROC dcol from ROC col
2054     //dcol in mod
2055     int dColInDet = pIndexConverter.DColumnInModule(dColInChip, chipIndex);
2056 
2057     chips[chipIndex]++;
2058     columns[dColInDet]++;
2059     if (isPhase1) {
2060       if (topol->isItBigPixelInX(row) || topol->isItBigPixelInY(col))
2061         pixelBig[chipIndex]++;
2062       else
2063         pixelStd[chipIndex]++;
2064     }
2065   }
2066 
2067   // Delete some ROC hits.
2068   for (iter = chips.begin(); iter != chips.end(); iter++) {
2069     //float rand  = RandFlat::shoot();
2070     float rand = CLHEP::RandFlat::shoot(engine);
2071     if (rand > chipEfficiency)
2072       chips[iter->first] = 0;
2073   }
2074 
2075   // Delete some Dcol hits.
2076   for (iter = columns.begin(); iter != columns.end(); iter++) {
2077     //float rand  = RandFlat::shoot();
2078     float rand = CLHEP::RandFlat::shoot(engine);
2079     if (rand > columnEfficiency)
2080       columns[iter->first] = 0;
2081   }
2082 
2083   // Delete some pixel hits based on DCDC issue damage.
2084   if (isPhase1) {
2085     for (iter = pixelStd.begin(); iter != pixelStd.end(); iter++) {
2086       float rand = CLHEP::RandFlat::shoot(engine);
2087       if (rand > pixelEfficiencyROCStdPixels[iter->first])
2088         pixelStd[iter->first] = 0;
2089     }
2090 
2091     for (iter = pixelBig.begin(); iter != pixelBig.end(); iter++) {
2092       float rand = CLHEP::RandFlat::shoot(engine);
2093       if (rand > pixelEfficiencyROCBigPixels[iter->first])
2094         pixelBig[iter->first] = 0;
2095     }
2096   }
2097 
2098   // Now loop again over pixels to kill some of them.
2099   // Loop over hit pixels, amplitude in electrons, channel = coded row,col
2100   for (signal_map_iterator i = theSignal.begin(); i != theSignal.end(); ++i) {
2101     //    int chan = i->first;
2102     std::pair<int, int> ip = PixelDigi::channelToPixel(i->first);  //get pixel pos
2103     int row = ip.first;                                            // X in row
2104     int col = ip.second;                                           // Y is in col
2105     //transform to ROC index coordinates
2106     pIndexConverter.transformToROC(col, row, chipIndex, colROC, rowROC);
2107     int dColInChip = pIndexConverter.DColumn(colROC);  //get ROC dcol from ROC col
2108     //dcol in mod
2109     int dColInDet = pIndexConverter.DColumnInModule(dColInChip, chipIndex);
2110 
2111     //float rand  = RandFlat::shoot();
2112     float rand = CLHEP::RandFlat::shoot(engine);
2113     if (chips[chipIndex] == 0 || columns[dColInDet] == 0 || rand > pixelEfficiency ||
2114         (pixelStd.count(chipIndex) && pixelStd[chipIndex] == 0) ||
2115         (pixelBig.count(chipIndex) && pixelBig[chipIndex] == 0)) {
2116       // make pixel amplitude =0, pixel will be lost at clusterization
2117       i->second.set(0.);  // reset amplitude,
2118     }  // end if
2119     if (isPhase1) {
2120       if ((pixelStd.count(chipIndex) && pixelStd[chipIndex] == 0) ||
2121           (pixelBig.count(chipIndex) && pixelBig[chipIndex] == 0) || (badRocsFromFEDChannels.at(chipIndex) == 1)) {
2122         //============================================================
2123         // make pixel amplitude =0, pixel will be lost at clusterization
2124         i->second.set(0.);  // reset amplitude,
2125       }  // end if
2126     }  // is Phase 1
2127     if (KillBadFEDChannels && badRocsFromFEDChannels.at(chipIndex) == 1) {
2128       i->second.set(0.);
2129     }
2130   }  // end pixel loop
2131 }  // end pixel_indefficiency
2132 
2133 //***************************************************************************************
2134 // Simulate pixel aging with an exponential function
2135 //**************************************************************************************
2136 
2137 float SiPixelDigitizerAlgorithm::pixel_aging(const PixelAging& aging,
2138                                              const PixelGeomDetUnit* pixdet,
2139                                              const TrackerTopology* tTopo) const {
2140   uint32_t detID = pixdet->geographicalId().rawId();
2141 
2142   // Predefined damage parameter (no aging)
2143   float pseudoRadDamage = 0.0f;
2144 
2145   // setup the chip indices conversion
2146   if (pixdet->subDetector() == GeomDetEnumerators::SubDetector::PixelBarrel ||
2147       pixdet->subDetector() == GeomDetEnumerators::SubDetector::P1PXB) {  // barrel layers
2148     int layerIndex = tTopo->layer(detID);
2149 
2150     pseudoRadDamage = aging.thePixelPseudoRadDamage[layerIndex - 1];
2151 
2152     LogDebug("Pixel Digitizer") << "pixel_aging: "
2153                                 << "\n";
2154     LogDebug("Pixel Digitizer") << "Subid " << pixdet->subDetector() << " layerIndex " << layerIndex << " ladder "
2155                                 << tTopo->pxbLadder(detID) << " module  " << tTopo->pxbModule(detID) << "\n";
2156 
2157   } else if (pixdet->subDetector() == GeomDetEnumerators::SubDetector::PixelEndcap ||
2158              pixdet->subDetector() == GeomDetEnumerators::SubDetector::P1PXEC ||
2159              pixdet->subDetector() == GeomDetEnumerators::SubDetector::P2PXEC) {  // forward disks
2160     unsigned int diskIndex =
2161         tTopo->layer(detID) + aging.FPixIndex;  // Use diskIndex-1 later to stay consistent with BPix
2162 
2163     pseudoRadDamage = aging.thePixelPseudoRadDamage[diskIndex - 1];
2164 
2165     LogDebug("Pixel Digitizer") << "pixel_aging: "
2166                                 << "\n";
2167     LogDebug("Pixel Digitizer") << "Subid " << pixdet->subDetector() << " diskIndex " << diskIndex << "\n";
2168   } else if (pixdet->subDetector() == GeomDetEnumerators::SubDetector::P2OTB ||
2169              pixdet->subDetector() == GeomDetEnumerators::SubDetector::P2OTEC) {
2170     // if phase 2 OT hardcoded value as it has always been
2171     pseudoRadDamage = 0.f;
2172   }  // if barrel/forward
2173 
2174   LogDebug("Pixel Digitizer") << " pseudoRadDamage " << pseudoRadDamage << "\n";
2175   LogDebug("Pixel Digitizer") << " end pixel_aging "
2176                               << "\n";
2177 
2178   return pseudoRadDamage;
2179 #ifdef TP_DEBUG
2180   LogDebug("Pixel Digitizer") << " enter pixel_aging " << pseudoRadDamage;
2181 #endif
2182 }
2183 
2184 //***********************************************************************
2185 
2186 // Fluctuate the gain and offset for the amplitude calibration
2187 // Use gaussian smearing.
2188 //float SiPixelDigitizerAlgorithm::missCalibrate(const float amp) const {
2189 //float gain  = RandGaussQ::shoot(1.,theGainSmearing);
2190 //float offset  = RandGaussQ::shoot(0.,theOffsetSmearing);
2191 //float newAmp = amp * gain + offset;
2192 // More complex misscalibration
2193 float SiPixelDigitizerAlgorithm::missCalibrate(uint32_t detID,
2194                                                const TrackerTopology* tTopo,
2195                                                const PixelGeomDetUnit* pixdet,
2196                                                int col,
2197                                                int row,
2198                                                const float signalInElectrons) const {
2199   // Central values
2200   //const float p0=0.00352, p1=0.868, p2=112., p3=113.; // pix(0,0,0)
2201   //  const float p0=0.00382, p1=0.886, p2=112.7, p3=113.0; // average roc=0
2202   //const float p0=0.00492, p1=1.998, p2=90.6, p3=134.1; // average roc=6
2203   // Smeared (rms)
2204   //const float s0=0.00020, s1=0.051, s2=5.4, s3=4.4; // average roc=0
2205   //const float s0=0.00015, s1=0.043, s2=3.2, s3=3.1; // col average roc=0
2206 
2207   // Make 2 sets of parameters for Fpix and BPIx:
2208 
2209   float p0 = 0.0f;
2210   float p1 = 0.0f;
2211   float p2 = 0.0f;
2212   float p3 = 0.0f;
2213 
2214   if (pixdet->type().isTrackerPixel() && pixdet->type().isBarrel()) {  // barrel layers
2215     p0 = BPix_p0;
2216     p1 = BPix_p1;
2217     p2 = BPix_p2;
2218     p3 = BPix_p3;
2219   } else if (pixdet->type().isTrackerPixel()) {  // forward disks
2220     p0 = FPix_p0;
2221     p1 = FPix_p1;
2222     p2 = FPix_p2;
2223     p3 = FPix_p3;
2224   } else {
2225     throw cms::Exception("NotAPixelGeomDetUnit") << "Not a pixel geomdet unit" << detID;
2226   }
2227 
2228   float newAmp = 0.f;  //Modified signal
2229 
2230   // Convert electrons to VCAL units
2231   float signal = (signalInElectrons - electronsPerVCAL_Offset) / electronsPerVCAL;
2232 
2233   // New gains/offsets are needed for phase1 L1
2234   int layer = 0;
2235   if (DetId(detID).subdetId() == 1)
2236     layer = tTopo->pxbLayer(detID);
2237   if (layer == 1)
2238     signal = (signalInElectrons - electronsPerVCAL_L1_Offset) / electronsPerVCAL_L1;
2239 
2240   // Simulate the analog response with fixed parametrization
2241   newAmp = p3 + p2 * tanh(p0 * signal - p1);
2242 
2243   // Use the pixel-by-pixel calibrations
2244   //transform to ROC index coordinates
2245   //int chipIndex=0, colROC=0, rowROC=0;
2246   //std::unique_ptr<PixelIndices> pIndexConverter(new PixelIndices(numColumns,numRows));
2247   //pIndexConverter->transformToROC(col,row,chipIndex,colROC,rowROC);
2248 
2249   // Use calibration from a file
2250   //int chanROC = PixelIndices::pixelToChannelROC(rowROC,colROC); // use ROC coordinates
2251   //float pp0=0, pp1=0,pp2=0,pp3=0;
2252   //map<int,CalParameters,std::less<int> >::const_iterator it=calmap.find(chanROC);
2253   //CalParameters y  = (*it).second;
2254   //pp0 = y.p0;
2255   //pp1 = y.p1;
2256   //pp2 = y.p2;
2257   //pp3 = y.p3;
2258 
2259   //
2260   // Use random smearing
2261   // Randomize the pixel response
2262   //float pp0  = RandGaussQ::shoot(p0,s0);
2263   //float pp1  = RandGaussQ::shoot(p1,s1);
2264   //float pp2  = RandGaussQ::shoot(p2,s2);
2265   //float pp3  = RandGaussQ::shoot(p3,s3);
2266 
2267   //newAmp = pp3 + pp2 * tanh(pp0*signal - pp1); // Final signal
2268 
2269   LogDebug("Pixel Digitizer") << " misscalibrate " << col << " " << row
2270                               << " "
2271                               // <<chipIndex<<" " <<colROC<<" " <<rowROC<<" "
2272                               << signalInElectrons << " " << signal << " " << newAmp << " "
2273                               << (signalInElectrons / theElectronPerADC) << "\n";
2274 
2275   return newAmp;
2276 }
2277 //******************************************************************************
2278 
2279 // Set the drift direction accoring to the Bfield in local det-unit frame
2280 // Works for both barrel and forward pixels.
2281 // Replace the sign convention to fit M.Swartz's formulaes.
2282 // Configurations for barrel and foward pixels possess different tanLorentzAngleperTesla
2283 // parameter value
2284 
2285 LocalVector SiPixelDigitizerAlgorithm::DriftDirection(const PixelGeomDetUnit* pixdet,
2286                                                       const GlobalVector& bfield,
2287                                                       const DetId& detId) const {
2288   Frame detFrame(pixdet->surface().position(), pixdet->surface().rotation());
2289   LocalVector Bfield = detFrame.toLocal(bfield);
2290 
2291   float alpha2_FPix;
2292   float alpha2_BPix;
2293   float alpha2;
2294 
2295   //float dir_x = -tanLorentzAnglePerTesla * Bfield.y();
2296   //float dir_y = +tanLorentzAnglePerTesla * Bfield.x();
2297   //float dir_z = -1.; // E field always in z direction, so electrons go to -z
2298   // The dir_z has to be +/- 1. !
2299   // LocalVector theDriftDirection = LocalVector(dir_x,dir_y,dir_z);
2300 
2301   float dir_x = 0.0f;
2302   float dir_y = 0.0f;
2303   float dir_z = 0.0f;
2304   float scale = 0.0f;
2305 
2306   uint32_t detID = pixdet->geographicalId().rawId();
2307 
2308   // Read Lorentz angle from cfg file:**************************************************************
2309 
2310   if (!use_LorentzAngle_DB_) {
2311     if (alpha2Order) {
2312       alpha2_FPix = tanLorentzAnglePerTesla_FPix * tanLorentzAnglePerTesla_FPix;
2313       alpha2_BPix = tanLorentzAnglePerTesla_BPix * tanLorentzAnglePerTesla_BPix;
2314     } else {
2315       alpha2_FPix = 0.0f;
2316       alpha2_BPix = 0.0f;
2317     }
2318 
2319     if (pixdet->type().isTrackerPixel() && pixdet->type().isBarrel()) {  // barrel layers
2320       dir_x = -(tanLorentzAnglePerTesla_BPix * Bfield.y() + alpha2_BPix * Bfield.z() * Bfield.x());
2321       dir_y = +(tanLorentzAnglePerTesla_BPix * Bfield.x() - alpha2_BPix * Bfield.z() * Bfield.y());
2322       dir_z = -(1 + alpha2_BPix * Bfield.z() * Bfield.z());
2323       scale = -dir_z;
2324     } else if (pixdet->type().isTrackerPixel()) {  // forward disks
2325       dir_x = -(tanLorentzAnglePerTesla_FPix * Bfield.y() + alpha2_FPix * Bfield.z() * Bfield.x());
2326       dir_y = +(tanLorentzAnglePerTesla_FPix * Bfield.x() - alpha2_FPix * Bfield.z() * Bfield.y());
2327       dir_z = -(1 + alpha2_FPix * Bfield.z() * Bfield.z());
2328       scale = -dir_z;
2329     } else {
2330       throw cms::Exception("NotAPixelGeomDetUnit") << "Not a pixel geomdet unit" << detID;
2331     }
2332   }  // end: Read LA from cfg file.
2333 
2334   //Read Lorentz angle from DB:********************************************************************
2335   if (use_LorentzAngle_DB_) {
2336     float lorentzAngle = SiPixelLorentzAngle_->getLorentzAngle(detId);
2337     alpha2 = lorentzAngle * lorentzAngle;
2338     dir_x = -(lorentzAngle * Bfield.y() + alpha2 * Bfield.z() * Bfield.x());
2339     dir_y = +(lorentzAngle * Bfield.x() - alpha2 * Bfield.z() * Bfield.y());
2340     dir_z = -(1 + alpha2 * Bfield.z() * Bfield.z());
2341     scale = -dir_z;
2342   }  // end: Read LA from DataBase.
2343 
2344   LocalVector theDriftDirection = LocalVector(dir_x / scale, dir_y / scale, dir_z / scale);
2345 
2346 #ifdef TP_DEBUG
2347   LogDebug("Pixel Digitizer") << " The drift direction in local coordinate is " << theDriftDirection;
2348 #endif
2349 
2350   return theDriftDirection;
2351 }
2352 
2353 //****************************************************************************************************
2354 
2355 void SiPixelDigitizerAlgorithm::pixel_inefficiency_db(uint32_t detID) {
2356   signal_map_type& theSignal = _signal[detID];
2357 
2358   // Loop over hit pixels, amplitude in electrons, channel = coded row,col
2359   for (signal_map_iterator i = theSignal.begin(); i != theSignal.end(); ++i) {
2360     //    int chan = i->first;
2361     std::pair<int, int> ip = PixelDigi::channelToPixel(i->first);  //get pixel pos
2362     int row = ip.first;                                            // X in row
2363     int col = ip.second;                                           // Y is in col
2364     //transform to ROC index coordinates
2365     if (theSiPixelGainCalibrationService_->isDead(detID, col, row)) {
2366       LogDebug("Pixel Digitizer") << "now in isdead check, row " << detID << " " << col << "," << row << "\n";
2367       // make pixel amplitude =0, pixel will be lost at clusterization
2368       i->second.set(0.);  // reset amplitude,
2369     }  // end if
2370   }  // end pixel loop
2371 }  // end pixel_indefficiency
2372 
2373 //****************************************************************************************************
2374 
2375 void SiPixelDigitizerAlgorithm::module_killing_conf(uint32_t detID) {
2376   bool isbad = false;
2377 
2378   Parameters::const_iterator itDeadModules = DeadModules.begin();
2379 
2380   int detid = detID;
2381   for (; itDeadModules != DeadModules.end(); ++itDeadModules) {
2382     int Dead_detID = itDeadModules->getParameter<int>("Dead_detID");
2383     if (detid == Dead_detID) {
2384       isbad = true;
2385       break;
2386     }
2387   }
2388 
2389   if (!isbad)
2390     return;
2391 
2392   signal_map_type& theSignal = _signal[detID];
2393 
2394   std::string Module = itDeadModules->getParameter<std::string>("Module");
2395 
2396   if (Module == "whole") {
2397     for (signal_map_iterator i = theSignal.begin(); i != theSignal.end(); ++i) {
2398       i->second.set(0.);  // reset amplitude
2399     }
2400   }
2401 
2402   for (signal_map_iterator i = theSignal.begin(); i != theSignal.end(); ++i) {
2403     std::pair<int, int> ip = PixelDigi::channelToPixel(i->first);  //get pixel pos
2404 
2405     if (Module == "tbmA" && ip.first >= 80 && ip.first <= 159) {
2406       i->second.set(0.);
2407     }
2408 
2409     if (Module == "tbmB" && ip.first <= 79) {
2410       i->second.set(0.);
2411     }
2412   }
2413 }
2414 //****************************************************************************************************
2415 void SiPixelDigitizerAlgorithm::module_killing_DB(uint32_t detID) {
2416   // Not SLHC safe for now
2417 
2418   bool isbad = false;
2419 
2420   std::vector<SiPixelQuality::disabledModuleType> disabledModules = SiPixelBadModule_->getBadComponentList();
2421 
2422   SiPixelQuality::disabledModuleType badmodule;
2423 
2424   for (size_t id = 0; id < disabledModules.size(); id++) {
2425     if (detID == disabledModules[id].DetID) {
2426       isbad = true;
2427       badmodule = disabledModules[id];
2428       break;
2429     }
2430   }
2431 
2432   if (!isbad)
2433     return;
2434 
2435   signal_map_type& theSignal = _signal[detID];
2436 
2437   LogDebug("Pixel Digitizer") << "Hit in: " << detID << " errorType " << badmodule.errorType << " BadRocs=" << std::hex
2438                               << SiPixelBadModule_->getBadRocs(detID) << std::dec << " "
2439                               << "\n";
2440   if (badmodule.errorType == 0) {  // this is a whole dead module.
2441 
2442     for (signal_map_iterator i = theSignal.begin(); i != theSignal.end(); ++i) {
2443       i->second.set(0.);  // reset amplitude
2444     }
2445   } else {  // all other module types: half-modules and single ROCs.
2446     // Get Bad ROC position:
2447     //follow the example of getBadRocPositions in CondFormats/SiPixelObjects/src/SiPixelQuality.cc
2448     std::vector<GlobalPixel> badrocpositions(0);
2449     for (unsigned int j = 0; j < 16; j++) {
2450       if (SiPixelBadModule_->IsRocBad(detID, j) == true) {
2451         std::vector<CablingPathToDetUnit> path = map_->pathToDetUnit(detID);
2452         typedef std::vector<CablingPathToDetUnit>::const_iterator IT;
2453         for (IT it = path.begin(); it != path.end(); ++it) {
2454           const PixelROC* myroc = map_->findItem(*it);
2455           if (myroc->idInDetUnit() == j) {
2456             LocalPixel::RocRowCol local = {39, 25};  //corresponding to center of ROC row, col
2457             GlobalPixel global = myroc->toGlobal(LocalPixel(local));
2458             badrocpositions.push_back(global);
2459             break;
2460           }
2461         }
2462       }
2463     }  // end of getBadRocPositions
2464 
2465     for (signal_map_iterator i = theSignal.begin(); i != theSignal.end(); ++i) {
2466       std::pair<int, int> ip = PixelDigi::channelToPixel(i->first);  //get pixel pos
2467 
2468       for (std::vector<GlobalPixel>::const_iterator it = badrocpositions.begin(); it != badrocpositions.end(); ++it) {
2469         if (it->row >= 80 && ip.first >= 80) {
2470           if ((std::abs(ip.second - it->col) < 26)) {
2471             i->second.set(0.);
2472           } else if (it->row == 120 && ip.second - it->col == 26) {
2473             i->second.set(0.);
2474           } else if (it->row == 119 && it->col - ip.second == 26) {
2475             i->second.set(0.);
2476           }
2477         } else if (it->row < 80 && ip.first < 80) {
2478           if ((std::abs(ip.second - it->col) < 26)) {
2479             i->second.set(0.);
2480           } else if (it->row == 40 && ip.second - it->col == 26) {
2481             i->second.set(0.);
2482           } else if (it->row == 39 && it->col - ip.second == 26) {
2483             i->second.set(0.);
2484           }
2485         }
2486       }
2487     }
2488   }
2489 }
2490 
2491 /******************************************************************/
2492 
2493 void SiPixelDigitizerAlgorithm::lateSignalReweight(const PixelGeomDetUnit* pixdet,
2494                                                    std::vector<PixelDigi>& digis,
2495                                                    std::vector<PixelSimHitExtraInfo>& newClass_Sim_extra,
2496                                                    const TrackerTopology* tTopo,
2497                                                    CLHEP::HepRandomEngine* engine) {
2498   // Function to apply the Charge Reweighting on top of digi in case of PU from mixing library
2499   // for time dependent MC
2500   std::vector<PixelDigi> New_digis;
2501   uint32_t detID = pixdet->geographicalId().rawId();
2502 
2503   if (UseReweighting) {
2504     LogError("PixelDigitizer ") << " ********************************  \n";
2505     LogError("PixelDigitizer ") << " ********************************  \n";
2506     LogError("PixelDigitizer ") << " *****  INCONSISTENCY !!!   *****  \n";
2507     LogError("PixelDigitizer ")
2508         << " applyLateReweighting_ and UseReweighting can not be true at the same time for PU ! \n";
2509     LogError("PixelDigitizer ") << " ---> DO NOT APPLY CHARGE REWEIGHTING TWICE !!! \n";
2510     LogError("PixelDigitizer ") << " ******************************** \n";
2511     LogError("PixelDigitizer ") << " ******************************** \n";
2512     return;
2513   }
2514 
2515   float thePixelThresholdInE = 0.;
2516   if (theNoiseInElectrons > 0.) {
2517     if (pixdet->type().isTrackerPixel() && pixdet->type().isBarrel()) {  // Barrel modules
2518       int lay = tTopo->layer(detID);
2519       if (addThresholdSmearing) {
2520         if (pixdet->subDetector() == GeomDetEnumerators::SubDetector::PixelBarrel ||
2521             pixdet->subDetector() == GeomDetEnumerators::SubDetector::P1PXB) {
2522           if (lay == 1) {
2523             thePixelThresholdInE = CLHEP::RandGaussQ::shoot(
2524                 engine, theThresholdInE_BPix_L1, theThresholdSmearing_BPix_L1);  // gaussian smearing
2525           } else if (lay == 2) {
2526             thePixelThresholdInE = CLHEP::RandGaussQ::shoot(
2527                 engine, theThresholdInE_BPix_L2, theThresholdSmearing_BPix_L2);  // gaussian smearing
2528           } else {
2529             thePixelThresholdInE =
2530                 CLHEP::RandGaussQ::shoot(engine, theThresholdInE_BPix, theThresholdSmearing_BPix);  // gaussian smearing
2531           }
2532         }
2533       } else {
2534         if (pixdet->subDetector() == GeomDetEnumerators::SubDetector::PixelBarrel ||
2535             pixdet->subDetector() == GeomDetEnumerators::SubDetector::P1PXB) {
2536           if (lay == 1) {
2537             thePixelThresholdInE = theThresholdInE_BPix_L1;
2538           } else if (lay == 2) {
2539             thePixelThresholdInE = theThresholdInE_BPix_L2;
2540           } else {
2541             thePixelThresholdInE = theThresholdInE_BPix;  // no smearing
2542           }
2543         }
2544       }
2545 
2546     } else if (pixdet->type().isTrackerPixel()) {  // Forward disks modules
2547 
2548       if (addThresholdSmearing) {
2549         thePixelThresholdInE =
2550             CLHEP::RandGaussQ::shoot(engine, theThresholdInE_FPix, theThresholdSmearing_FPix);  // gaussian smearing
2551       } else {
2552         thePixelThresholdInE = theThresholdInE_FPix;  // no smearing
2553       }
2554 
2555     } else {
2556       throw cms::Exception("NotAPixelGeomDetUnit") << "Not a pixel geomdet unit" << detID;
2557     }
2558   }
2559 
2560   // loop on the SimHit extra info class
2561   // apply the reweighting for that SimHit on a cluster way
2562   bool reweighted = false;
2563   std::vector<PixelSimHitExtraInfo>::iterator loopTempSH;
2564   for (loopTempSH = newClass_Sim_extra.begin(); loopTempSH != newClass_Sim_extra.end(); ++loopTempSH) {
2565     signal_map_type theDigiSignal;
2566     PixelSimHitExtraInfo TheNewInfo = *loopTempSH;
2567     reweighted = TheNewSiPixelChargeReweightingAlgorithmClass->lateSignalReweight<digitizerUtility::Amplitude>(
2568         pixdet, digis, TheNewInfo, theDigiSignal, tTopo, engine);
2569     if (!reweighted) {
2570       // loop on the non-reweighthed digis associated to the considered SimHit
2571       std::vector<PixelDigi>::const_iterator loopDigi;
2572       for (loopDigi = digis.begin(); loopDigi != digis.end(); ++loopDigi) {
2573         unsigned int chan = loopDigi->channel();
2574         // check if that digi is related to the SimHit
2575         if (loopTempSH->isInTheList(chan)) {
2576           float corresponding_charge = loopDigi->adc();
2577           theDigiSignal[chan] += digitizerUtility::Amplitude(corresponding_charge, corresponding_charge);
2578         }
2579       }
2580     }
2581 
2582     // transform theDigiSignal into digis
2583     int Thresh_inADC = int(thePixelThresholdInE / theElectronPerADC);
2584     for (signal_map_const_iterator i = theDigiSignal.begin(); i != theDigiSignal.end(); ++i) {
2585       float signalInADC = (*i).second;  // signal in ADC
2586       if (signalInADC > 0.) {
2587         if (signalInADC >= Thresh_inADC) {
2588           int chan = (*i).first;  // channel number
2589           std::pair<int, int> ip = PixelDigi::channelToPixel(chan);
2590           int adc = int(signalInADC);
2591           // add MissCalibration
2592           if (doMissCalInLateCR) {
2593             int row = ip.first;
2594             int col = ip.second;
2595             adc =
2596                 int(missCalibrate(detID, tTopo, pixdet, col, row, signalInADC * theElectronPerADC));  //full misscalib.
2597           }
2598 
2599           if (adc > theAdcFullScLateCR)
2600             adc = theAdcFullScLateCR;  // Check maximum value
2601 
2602 #ifdef TP_DEBUG
2603           LogDebug("Pixel Digitizer") << (*i).first << " " << (*i).second << " " << signalInADC << " " << adc
2604                                       << ip.first << " " << ip.second;
2605 #endif
2606 
2607           // Load digis
2608           New_digis.emplace_back(ip.first, ip.second, adc);
2609         }
2610       }
2611     }  // end loop on theDigiSignal
2612     theDigiSignal.clear();
2613   }
2614   digis.clear();
2615   digis = New_digis;
2616 }
2617 
2618 /******************************************************************/
2619 
2620 void SiPixelDigitizerAlgorithm::lateSignalReweight(const PixelGeomDetUnit* pixdet,
2621                                                    std::vector<PixelDigi>& digis,
2622                                                    std::vector<PixelSimHitExtraInfoLite>& newClass_Sim_extra,
2623                                                    const TrackerTopology* tTopo,
2624                                                    CLHEP::HepRandomEngine* engine) {
2625   // Function to apply the Charge Reweighting on top of digi in case of PU from mixing library
2626   // for time dependent MC
2627   std::vector<PixelDigi> New_digis;
2628   uint32_t detID = pixdet->geographicalId().rawId();
2629 
2630   if (UseReweighting) {
2631     LogError("PixelDigitizer ") << " ********************************  \n";
2632     LogError("PixelDigitizer ") << " ********************************  \n";
2633     LogError("PixelDigitizer ") << " *****  INCONSISTENCY !!!   *****  \n";
2634     LogError("PixelDigitizer ")
2635         << " applyLateReweighting_ and UseReweighting can not be true at the same time for PU ! \n";
2636     LogError("PixelDigitizer ") << " ---> DO NOT APPLY CHARGE REWEIGHTING TWICE !!! \n";
2637     LogError("PixelDigitizer ") << " ******************************** \n";
2638     LogError("PixelDigitizer ") << " ******************************** \n";
2639     return;
2640   }
2641 
2642   float thePixelThresholdInE = 0.;
2643   if (theNoiseInElectrons > 0.) {
2644     if (pixdet->type().isTrackerPixel() && pixdet->type().isBarrel()) {  // Barrel modules
2645       int lay = tTopo->layer(detID);
2646       if (addThresholdSmearing) {
2647         if (pixdet->subDetector() == GeomDetEnumerators::SubDetector::PixelBarrel ||
2648             pixdet->subDetector() == GeomDetEnumerators::SubDetector::P1PXB) {
2649           if (lay == 1) {
2650             thePixelThresholdInE = CLHEP::RandGaussQ::shoot(
2651                 engine, theThresholdInE_BPix_L1, theThresholdSmearing_BPix_L1);  // gaussian smearing
2652           } else if (lay == 2) {
2653             thePixelThresholdInE = CLHEP::RandGaussQ::shoot(
2654                 engine, theThresholdInE_BPix_L2, theThresholdSmearing_BPix_L2);  // gaussian smearing
2655           } else {
2656             thePixelThresholdInE =
2657                 CLHEP::RandGaussQ::shoot(engine, theThresholdInE_BPix, theThresholdSmearing_BPix);  // gaussian smearing
2658           }
2659         }
2660       } else {
2661         if (pixdet->subDetector() == GeomDetEnumerators::SubDetector::PixelBarrel ||
2662             pixdet->subDetector() == GeomDetEnumerators::SubDetector::P1PXB) {
2663           if (lay == 1) {
2664             thePixelThresholdInE = theThresholdInE_BPix_L1;
2665           } else if (lay == 2) {
2666             thePixelThresholdInE = theThresholdInE_BPix_L2;
2667           } else {
2668             thePixelThresholdInE = theThresholdInE_BPix;  // no smearing
2669           }
2670         }
2671       }
2672 
2673     } else if (pixdet->type().isTrackerPixel()) {  // Forward disks modules
2674 
2675       if (addThresholdSmearing) {
2676         thePixelThresholdInE =
2677             CLHEP::RandGaussQ::shoot(engine, theThresholdInE_FPix, theThresholdSmearing_FPix);  // gaussian smearing
2678       } else {
2679         thePixelThresholdInE = theThresholdInE_FPix;  // no smearing
2680       }
2681 
2682     } else {
2683       throw cms::Exception("NotAPixelGeomDetUnit") << "Not a pixel geomdet unit" << detID;
2684     }
2685   }
2686 
2687   // loop on the SimHit extra info class
2688   // apply the reweighting for that SimHit on a cluster way
2689   bool reweighted = false;
2690   std::vector<PixelSimHitExtraInfoLite>::iterator loopTempSH;
2691   for (loopTempSH = newClass_Sim_extra.begin(); loopTempSH != newClass_Sim_extra.end(); ++loopTempSH) {
2692     signal_map_type theDigiSignal;
2693     PixelSimHitExtraInfoLite TheNewInfo = *loopTempSH;
2694     reweighted = TheNewSiPixelChargeReweightingAlgorithmClass->lateSignalReweight<digitizerUtility::Amplitude>(
2695         pixdet, digis, TheNewInfo, theDigiSignal, tTopo, engine);
2696     if (!reweighted) {
2697       // loop on the non-reweighthed digis associated to the considered SimHit
2698       std::vector<PixelDigi>::const_iterator loopDigi;
2699       for (loopDigi = digis.begin(); loopDigi != digis.end(); ++loopDigi) {
2700         unsigned int chan = loopDigi->channel();
2701         // check if that digi is related to the SimHit
2702         if (loopTempSH->isInTheList(chan)) {
2703           float corresponding_charge = loopDigi->adc();
2704           theDigiSignal[chan] += digitizerUtility::Amplitude(corresponding_charge, corresponding_charge);
2705         }
2706       }
2707     }
2708 
2709     // transform theDigiSignal into digis
2710     int Thresh_inADC = int(thePixelThresholdInE / theElectronPerADC);
2711     for (signal_map_const_iterator i = theDigiSignal.begin(); i != theDigiSignal.end(); ++i) {
2712       float signalInADC = (*i).second;  // signal in ADC
2713       if (signalInADC > 0.) {
2714         if (signalInADC >= Thresh_inADC) {
2715           int chan = (*i).first;  // channel number
2716           std::pair<int, int> ip = PixelDigi::channelToPixel(chan);
2717           int adc = int(signalInADC);
2718           // add MissCalibration
2719           if (doMissCalInLateCR) {
2720             int row = ip.first;
2721             int col = ip.second;
2722             adc =
2723                 int(missCalibrate(detID, tTopo, pixdet, col, row, signalInADC * theElectronPerADC));  //full misscalib.
2724           }
2725 
2726           if (adc > theAdcFullScLateCR)
2727             adc = theAdcFullScLateCR;  // Check maximum value
2728 
2729 #ifdef TP_DEBUG
2730           LogDebug("Pixel Digitizer") << (*i).first << " " << (*i).second << " " << signalInADC << " " << adc
2731                                       << ip.first << " " << ip.second;
2732 #endif
2733 
2734           // Load digis
2735           New_digis.emplace_back(ip.first, ip.second, adc);
2736         }
2737       }
2738     }  // end loop on theDigiSignal
2739     theDigiSignal.clear();
2740   }
2741   digis.clear();
2742   digis = New_digis;
2743 }