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