Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-01-16 23:37:29

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