SiPixelDigitizer/plugins/SiPixelDigitizerAlgorithm.cc

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