Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:30:54 

0001 #include <cmath>
0002 #include <iostream>
0003 #include <typeinfo>
0004 
0005 #include "CLHEP/Random/RandGaussQ.h"
0006 #include "CalibTracker/SiPixelESProducers/interface/SiPixelGainCalibrationOfflineSimService.h"
0007 #include "CondFormats/SiPhase2TrackerObjects/interface/SiPhase2OuterTrackerLorentzAngle.h"
0008 #include "CondFormats/SiPixelObjects/interface/GlobalPixel.h"
0009 #include "CondFormats/SiPixelObjects/interface/LocalPixel.h"
0010 #include "CondFormats/SiPixelObjects/interface/SiPixelLorentzAngle.h"
0011 #include "DataFormats/DetId/interface/DetId.h"
0012 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0013 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0014 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0015 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/ServiceRegistry/interface/Service.h"
0019 #include "FWCore/Utilities/interface/Exception.h"
0020 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0021 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0022 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0023 #include "SimGeneral/NoiseGenerators/interface/GaussianTailNoiseGenerator.h"
0024 #include "SimTracker/Common/interface/SiG4UniversalFluctuation.h"
0025 #include "SimTracker/SiPhase2Digitizer/plugins/Phase2TrackerDigitizerAlgorithm.h"
0026 
0027 using namespace edm;
0028 
0029 namespace ph2tkdigialgo {
0030   // Mass in MeV
0031   constexpr double m_pion = 139.571;
0032   constexpr double m_kaon = 493.677;
0033   constexpr double m_electron = 0.511;
0034   constexpr double m_muon = 105.658;
0035   constexpr double m_proton = 938.272;
0036 }  // namespace ph2tkdigialgo
0037 
0038 Phase2TrackerDigitizerAlgorithm::Phase2TrackerDigitizerAlgorithm(const edm::ParameterSet& conf_common,
0039                                                                  const edm::ParameterSet& conf_specific,
0040                                                                  edm::ConsumesCollector iC)
0041     : _signal(),
0042       makeDigiSimLinks_(conf_common.getUntrackedParameter<bool>("makeDigiSimLinks", true)),
0043       use_ineff_from_db_(conf_specific.getParameter<bool>("Inefficiency_DB")),
0044       use_module_killing_(conf_specific.getParameter<bool>("KillModules")),    // boolean to kill or not modules
0045       use_deadmodule_DB_(conf_specific.getParameter<bool>("DeadModules_DB")),  // boolean to access dead modules from DB
0046       // boolean to access Lorentz angle from DB
0047       use_LorentzAngle_DB_(conf_specific.getParameter<bool>("LorentzAngle_DB")),
0048 
0049       // get dead module from cfg file
0050       deadModules_(use_deadmodule_DB_ ? Parameters() : conf_specific.getParameter<Parameters>("DeadModules")),
0051 
0052       // Common pixel parameters
0053       // These are parameters which are not likely to be changed
0054       GeVperElectron_(3.61E-09),                                      // 1 electron(3.61eV, 1keV(277e, mod 9/06 d.k.
0055       alpha2Order_(conf_specific.getParameter<bool>("Alpha2Order")),  // switch on/off of E.B effect
0056       addXtalk_(conf_specific.getParameter<bool>("AddXTalk")),
0057       // Interstrip Coupling - Not used in PixelDigitizerAlgorithm
0058       interstripCoupling_(conf_specific.getParameter<double>("InterstripCoupling")),
0059 
0060       Sigma0_(conf_specific.getParameter<double>("SigmaZero")),       // Charge diffusion constant 7->3.7
0061       SigmaCoeff_(conf_specific.getParameter<double>("SigmaCoeff")),  // delta in the diffusion across the strip pitch
0062       // (set between 0 to 0.9,  0-->flat Sigma0, 1-->Sigma_min=0 & Sigma_max=2*Sigma0
0063       // D.B.: Dist300 replaced by moduleThickness, may not work with partially depleted sensors but works otherwise
0064       // Dist300(0.0300),                                          //   normalized to 300micron Silicon
0065 
0066       // Charge integration spread on the collection plane
0067       clusterWidth_(conf_specific.getParameter<double>("ClusterWidth")),
0068 
0069       // Allowed modes of readout which has following values :
0070       // 0          ---> Digital or binary readout
0071       // -1         ---> Analog readout, current digitizer (Inner Pixel) (TDR version) with no threshold subtraction
0072       // Analog readout with dual slope with the "second" slope being 1/2^(n-1) and threshold subtraction (n = 1, 2, 3,4)
0073       thePhase2ReadoutMode_(conf_specific.getParameter<int>("Phase2ReadoutMode")),
0074 
0075       // ADC calibration 1adc count(135e.
0076       // Corresponds to 2adc/kev, 270[e/kev]/135[e/adc](2[adc/kev]
0077       // Be careful, this parameter is also used in SiPixelDet.cc to
0078       // calculate the noise in adc counts from noise in electrons.
0079       // Both defaults should be the same.
0080       theElectronPerADC_(conf_specific.getParameter<double>("ElectronPerAdc")),
0081 
0082       // ADC saturation value, 255(8bit adc.
0083       theAdcFullScale_(conf_specific.getParameter<int>("AdcFullScale")),
0084 
0085       // Noise in electrons:
0086       // Pixel cell noise, relevant for generating noisy pixels
0087       theNoiseInElectrons_(conf_specific.getParameter<double>("NoiseInElectrons")),
0088 
0089       // Fill readout noise, including all readout chain, relevant for smearing
0090       theReadoutNoise_(conf_specific.getParameter<double>("ReadoutNoiseInElec")),
0091 
0092       // Threshold in units of noise:
0093       // thePixelThreshold(conf.getParameter<double>("ThresholdInNoiseUnits")),
0094       // Pixel threshold in electron units.
0095       theThresholdInE_Endcap_(conf_specific.getParameter<double>("ThresholdInElectrons_Endcap")),
0096       theThresholdInE_Barrel_(conf_specific.getParameter<double>("ThresholdInElectrons_Barrel")),
0097 
0098       // Add threshold gaussian smearing:
0099       theThresholdSmearing_Endcap_(conf_specific.getParameter<double>("ThresholdSmearing_Endcap")),
0100       theThresholdSmearing_Barrel_(conf_specific.getParameter<double>("ThresholdSmearing_Barrel")),
0101 
0102       // Add HIP Threshold in electron units.
0103       theHIPThresholdInE_Endcap_(conf_specific.getParameter<double>("HIPThresholdInElectrons_Endcap")),
0104       theHIPThresholdInE_Barrel_(conf_specific.getParameter<double>("HIPThresholdInElectrons_Barrel")),
0105 
0106       // theTofCut 12.5, cut in particle TOD +/- 12.5ns
0107       theTofLowerCut_(conf_specific.getParameter<double>("TofLowerCut")),
0108       theTofUpperCut_(conf_specific.getParameter<double>("TofUpperCut")),
0109 
0110       // Get the Lorentz angle from the cfg file:
0111       tanLorentzAnglePerTesla_Endcap_(
0112           use_LorentzAngle_DB_ ? 0.0 : conf_specific.getParameter<double>("TanLorentzAnglePerTesla_Endcap")),
0113       tanLorentzAnglePerTesla_Barrel_(
0114           use_LorentzAngle_DB_ ? 0.0 : conf_specific.getParameter<double>("TanLorentzAnglePerTesla_Barrel")),
0115 
0116       // Add noise
0117       addNoise_(conf_specific.getParameter<bool>("AddNoise")),
0118 
0119       // Add noisy pixels
0120       addNoisyPixels_(conf_specific.getParameter<bool>("AddNoisyPixels")),
0121 
0122       // Fluctuate charge in track subsegments
0123       fluctuateCharge_(conf_specific.getUntrackedParameter<bool>("FluctuateCharge", true)),
0124 
0125       // Control the pixel inefficiency
0126       addPixelInefficiency_(conf_specific.getParameter<bool>("AddInefficiency")),
0127 
0128       // Add threshold gaussian smearing:
0129       addThresholdSmearing_(conf_specific.getParameter<bool>("AddThresholdSmearing")),
0130 
0131       // Add some pseudo-red damage
0132       pseudoRadDamage_(conf_specific.exists("PseudoRadDamage") ? conf_specific.getParameter<double>("PseudoRadDamage")
0133                                                                : double(0.0)),
0134       pseudoRadDamageRadius_(conf_specific.exists("PseudoRadDamageRadius")
0135                                  ? conf_specific.getParameter<double>("PseudoRadDamageRadius")
0136                                  : double(0.0)),
0137       // Add pixel radiation damage
0138       useChargeReweighting_(conf_specific.getParameter<bool>("UseReweighting")),
0139       theSiPixelChargeReweightingAlgorithm_(
0140           useChargeReweighting_ ? std::make_unique<SiPixelChargeReweightingAlgorithm>(conf_specific, iC) : nullptr),
0141 
0142       // delta cutoff in MeV, has to be same as in OSCAR(0.030/cmsim=1.0 MeV
0143       // tMax(0.030), // In MeV.
0144       // tMax(conf.getUntrackedParameter<double>("DeltaProductionCut",0.030)),
0145       tMax_(conf_common.getParameter<double>("DeltaProductionCut")),
0146 
0147       badPixels_(conf_specific.getParameter<Parameters>("CellsToKill")),
0148 
0149       fluctuate_(fluctuateCharge_ ? std::make_unique<SiG4UniversalFluctuation>() : nullptr),
0150       theNoiser_(addNoise_ ? std::make_unique<GaussianTailNoiseGenerator>() : nullptr),
0151       theSiPixelGainCalibrationService_(
0152           use_ineff_from_db_ ? std::make_unique<SiPixelGainCalibrationOfflineSimService>(conf_specific, iC) : nullptr),
0153       subdetEfficiencies_(conf_specific) {
0154   LogDebug("Phase2TrackerDigitizerAlgorithm")
0155       << "Phase2TrackerDigitizerAlgorithm constructed\n"
0156       << "Configuration parameters:\n"
0157       << "Threshold/Gain = "
0158       << "threshold in electron Endcap = " << theThresholdInE_Endcap_
0159       << "\nthreshold in electron Barrel = " << theThresholdInE_Barrel_ << " ElectronPerADC " << theElectronPerADC_
0160       << " ADC Scale (in bits) " << theAdcFullScale_ << " The delta cut-off is set to " << tMax_ << " pix-inefficiency "
0161       << addPixelInefficiency_;
0162 }
0163 
0164 Phase2TrackerDigitizerAlgorithm::~Phase2TrackerDigitizerAlgorithm() {
0165   LogDebug("Phase2TrackerDigitizerAlgorithm") << "Phase2TrackerDigitizerAlgorithm deleted";
0166 }
0167 
0168 Phase2TrackerDigitizerAlgorithm::SubdetEfficiencies::SubdetEfficiencies(const edm::ParameterSet& conf) {
0169   barrel_efficiencies = conf.getParameter<std::vector<double> >("EfficiencyFactors_Barrel");
0170   endcap_efficiencies = conf.getParameter<std::vector<double> >("EfficiencyFactors_Endcap");
0171 }
0172 
0173 void Phase2TrackerDigitizerAlgorithm::accumulateSimHits(std::vector<PSimHit>::const_iterator inputBegin,
0174                                                         std::vector<PSimHit>::const_iterator inputEnd,
0175                                                         const size_t inputBeginGlobalIndex,
0176                                                         const uint32_t tofBin,
0177                                                         const Phase2TrackerGeomDetUnit* pixdet,
0178                                                         const GlobalVector& bfield) {
0179   // produce SignalPoint's for all SimHit's in detector
0180   // Loop over hits
0181   uint32_t detId = pixdet->geographicalId().rawId();
0182   size_t simHitGlobalIndex = inputBeginGlobalIndex;  // This needs to be stored to create the digi-sim link later
0183 
0184   // find the relevant hits
0185   std::vector<PSimHit> matchedSimHits;
0186   std::copy_if(inputBegin, inputEnd, std::back_inserter(matchedSimHits), [detId](auto const& hit) -> bool {
0187     return hit.detUnitId() == detId;
0188   });
0189   // loop over a much reduced set of SimHits
0190   for (auto const& hit : matchedSimHits) {
0191     LogDebug("Phase2DigitizerAlgorithm") << hit.particleType() << " " << hit.pabs() << " " << hit.energyLoss() << " "
0192                                          << hit.tof() << " " << hit.trackId() << " " << hit.processType() << " "
0193                                          << hit.detUnitId() << hit.entryPoint() << " " << hit.exitPoint();
0194     double signalScale = 1.0;
0195     // fill collection_points for this SimHit, indpendent of topology
0196     if (select_hit(hit, (pixdet->surface().toGlobal(hit.localPosition()).mag() * c_inv), signalScale)) {
0197       const auto& ionization_points = primary_ionization(hit);  // fills ionization_points
0198 
0199       // transforms ionization_points -> collection_points
0200       const auto& collection_points = drift(hit, pixdet, bfield, ionization_points);
0201 
0202       // compute induced signal on readout elements and add to _signal
0203       // hit needed only for SimHit<-->Digi link
0204       induce_signal(inputBegin, hit, simHitGlobalIndex, inputBeginGlobalIndex, tofBin, pixdet, collection_points);
0205       LogDebug("Phase2DigitizerAlgorithm") << "Induced signal was computed";
0206     }
0207     ++simHitGlobalIndex;
0208   }
0209 }
0210 // =================================================================
0211 //
0212 // Generate primary ionization along the track segment.
0213 // Divide the track into small sub-segments
0214 //
0215 // =================================================================
0216 std::vector<digitizerUtility::EnergyDepositUnit> Phase2TrackerDigitizerAlgorithm::primary_ionization(
0217     const PSimHit& hit) const {
0218   // Straight line approximation for trajectory inside active media
0219   constexpr float SegmentLength = 0.0010;  // in cm (10 microns)
0220   // Get the 3D segment direction vector
0221   LocalVector direction = hit.exitPoint() - hit.entryPoint();
0222 
0223   float eLoss = hit.energyLoss();  // Eloss in GeV
0224   float length = direction.mag();  // Track length in Silicon
0225 
0226   int NumberOfSegments = static_cast<int>(length / SegmentLength);  // Number of segments
0227   if (NumberOfSegments < 1)
0228     NumberOfSegments = 1;
0229   LogDebug("Phase2TrackerDigitizerAlgorithm")
0230       << "enter primary_ionzation " << NumberOfSegments << " shift = " << hit.exitPoint().x() - hit.entryPoint().x()
0231       << " " << hit.exitPoint().y() - hit.entryPoint().y() << " " << hit.exitPoint().z() - hit.entryPoint().z() << " "
0232       << hit.particleType() << " " << hit.pabs();
0233 
0234   std::vector<float> elossVector;
0235   elossVector.reserve(NumberOfSegments);
0236   if (fluctuateCharge_) {
0237     // Generate fluctuated charge points
0238     elossVector = fluctuateEloss(hit.particleType(), hit.pabs(), eLoss, length, NumberOfSegments);
0239   } else {
0240     float averageEloss = eLoss / NumberOfSegments;
0241     elossVector.resize(NumberOfSegments, averageEloss);
0242   }
0243 
0244   std::vector<digitizerUtility::EnergyDepositUnit> ionization_points;
0245   ionization_points.reserve(NumberOfSegments);  // set size
0246   // loop over segments
0247   for (size_t i = 0; i < elossVector.size(); ++i) {
0248     // Divide the segment into equal length subsegments
0249     Local3DPoint point = hit.entryPoint() + ((i + 0.5) / NumberOfSegments) * direction;
0250     float energy = elossVector[i] / GeVperElectron_;  // Convert charge to elec.
0251 
0252     digitizerUtility::EnergyDepositUnit edu(energy, point);  // define position,energy point
0253     ionization_points.push_back(edu);                        // save
0254     LogDebug("Phase2TrackerDigitizerAlgorithm")
0255         << "For index = " << i << " EnergyDepositUnit-x = " << edu.x() << " EnergyDepositUnit-y = " << edu.y()
0256         << " EnergyDepositUnit-z = " << edu.z() << " EnergyDepositUnit-energy = " << edu.energy();
0257   }
0258   return ionization_points;
0259 }
0260 //==============================================================================
0261 //
0262 // Fluctuate the charge comming from a small (10um) track segment.
0263 // Use the G4 routine. For mip pions for the moment.
0264 //
0265 //==============================================================================
0266 std::vector<float> Phase2TrackerDigitizerAlgorithm::fluctuateEloss(
0267     int pid, float particleMomentum, float eloss, float length, int NumberOfSegs) const {
0268   double particleMass = ph2tkdigialgo::m_pion;  // Mass in MeV, assume pion
0269   pid = std::abs(pid);
0270   if (pid != 211) {  // Mass in MeV
0271     if (pid == 11)
0272       particleMass = ph2tkdigialgo::m_electron;
0273     else if (pid == 13)
0274       particleMass = ph2tkdigialgo::m_muon;
0275     else if (pid == 321)
0276       particleMass = ph2tkdigialgo::m_kaon;
0277     else if (pid == 2212)
0278       particleMass = ph2tkdigialgo::m_proton;
0279   }
0280   // What is the track segment length.
0281   float segmentLength = length / NumberOfSegs;
0282 
0283   // Generate charge fluctuations.
0284   float sum = 0.;
0285   double segmentEloss = (1000. * eloss) / NumberOfSegs;  //eloss in MeV
0286   std::vector<float> elossVector;
0287   elossVector.reserve(NumberOfSegs);
0288   for (int i = 0; i < NumberOfSegs; ++i) {
0289     // The G4 routine needs momentum in MeV, mass in Mev, delta-cut in MeV,
0290     // track segment length in mm, segment eloss in MeV
0291     // Returns fluctuated eloss in MeV
0292     double deltaCutoff = tMax_;  // the cutoff is sometimes redefined inside, so fix it.
0293     float de = fluctuate_->SampleFluctuations(particleMomentum * 1000.,
0294                                               particleMass,
0295                                               deltaCutoff,
0296                                               segmentLength * 10.,
0297                                               segmentEloss,
0298                                               rengine_) /
0299                1000.;  //convert to GeV
0300     elossVector.push_back(de);
0301     sum += de;
0302   }
0303   if (sum > 0.) {  // if fluctuations give eloss>0.
0304     // Rescale to the same total eloss
0305     float ratio = eloss / sum;
0306     std::transform(
0307         std::begin(elossVector), std::end(elossVector), std::begin(elossVector), [&ratio](auto const& c) -> float {
0308           return c * ratio;
0309         });  // use a simple lambda expression
0310   } else {   // if fluctuations gives 0 eloss
0311     float averageEloss = eloss / NumberOfSegs;
0312     elossVector.resize(NumberOfSegs, averageEloss);
0313   }
0314   return elossVector;
0315 }
0316 
0317 // ======================================================================
0318 //
0319 // Drift the charge segments to the sensor surface (collection plane)
0320 // Include the effect of E-field and B-field
0321 //
0322 // =====================================================================
0323 std::vector<digitizerUtility::SignalPoint> Phase2TrackerDigitizerAlgorithm::drift(
0324     const PSimHit& hit,
0325     const Phase2TrackerGeomDetUnit* pixdet,
0326     const GlobalVector& bfield,
0327     const std::vector<digitizerUtility::EnergyDepositUnit>& ionization_points) const {
0328   LogDebug("Phase2TrackerDigitizerAlgorithm") << "Enter drift ";
0329 
0330   std::vector<digitizerUtility::SignalPoint> collection_points;
0331   collection_points.reserve(ionization_points.size());                     // set size
0332   LocalVector driftDir = driftDirection(pixdet, bfield, hit.detUnitId());  // get the charge drift direction
0333   if (driftDir.z() == 0.) {
0334     LogWarning("Phase2TrackerDigitizerAlgorithm") << " pxlx: drift in z is zero ";
0335     return collection_points;
0336   }
0337 
0338   float TanLorenzAngleX = driftDir.x();                                       // tangent of Lorentz angle
0339   float TanLorenzAngleY = 0.;                                                 // force to 0, driftDir.y()/driftDir.z();
0340   float dir_z = driftDir.z();                                                 // The z drift direction
0341   float CosLorenzAngleX = 1. / std::sqrt(1. + std::pow(TanLorenzAngleX, 2));  // cosine to estimate the path length
0342   float CosLorenzAngleY = 1.;
0343   if (alpha2Order_) {
0344     TanLorenzAngleY = driftDir.y();
0345     CosLorenzAngleY = 1. / std::sqrt(1. + std::pow(TanLorenzAngleY, 2));  // cosine
0346   }
0347 
0348   float moduleThickness = pixdet->specificSurface().bounds().thickness();
0349   float stripPitch = pixdet->specificTopology().pitch().first;
0350 
0351   LogDebug("Phase2TrackerDigitizerAlgorithm")
0352       << " Lorentz Tan-X " << TanLorenzAngleX << "  Lorentz Tan-Y " << TanLorenzAngleY << " Lorentz Cos-X "
0353       << CosLorenzAngleX << " Lorentz Cos-Y " << CosLorenzAngleY
0354       << " ticknes * Lorentz Tan-X = " << moduleThickness * TanLorenzAngleX << " drift direction " << driftDir;
0355 
0356   for (auto const& val : ionization_points) {
0357     // position
0358     float SegX = val.x(), SegY = val.y(), SegZ = val.z();
0359 
0360     // Distance from the collection plane
0361     // DriftDistance = (moduleThickness/2. + SegZ); // Drift to -z
0362     // Include explixitely the E drift direction (for CMS dir_z=-1)
0363 
0364     // Distance between charge generation and collection
0365     float driftDistance = moduleThickness / 2. - (dir_z * SegZ);  // Drift to -z
0366 
0367     if (driftDistance < 0.)
0368       driftDistance = 0.;
0369     else if (driftDistance > moduleThickness)
0370       driftDistance = moduleThickness;
0371 
0372     // Assume full depletion now, partial depletion will come later.
0373     float XDriftDueToMagField = driftDistance * TanLorenzAngleX;
0374     float YDriftDueToMagField = driftDistance * TanLorenzAngleY;
0375 
0376     // Shift cloud center
0377     float CloudCenterX = SegX + XDriftDueToMagField;
0378     float CloudCenterY = SegY + YDriftDueToMagField;
0379 
0380     // Calculate how long is the charge drift path
0381     // Actual Drift Lentgh
0382     float driftLength =
0383         std::sqrt(std::pow(driftDistance, 2) + std::pow(XDriftDueToMagField, 2) + std::pow(YDriftDueToMagField, 2));
0384 
0385     // What is the charge diffusion after this path
0386     // Sigma0=0.00037 is for 300um thickness (make sure moduleThickness is in [cm])
0387     float Sigma = std::sqrt(driftLength / moduleThickness) * Sigma0_ * moduleThickness / 0.0300;
0388     // D.B.: sigmaCoeff=0 means no modulation
0389     if (SigmaCoeff_)
0390       Sigma *= (SigmaCoeff_ * std::pow(cos(SegX * M_PI / stripPitch), 2) + 1);
0391     // NB: divided by 4 to get a periodicity of stripPitch
0392 
0393     // Project the diffusion sigma on the collection plane
0394     float Sigma_x = Sigma / CosLorenzAngleX;
0395     float Sigma_y = Sigma / CosLorenzAngleY;
0396 
0397     // Insert a charge loss due to Rad Damage here
0398     float energyOnCollector = val.energy();  // The energy that reaches the collector
0399 
0400     // pseudoRadDamage
0401     if (pseudoRadDamage_) {
0402       float moduleRadius = pixdet->surface().position().perp();
0403       if (moduleRadius <= pseudoRadDamageRadius_) {
0404         float kValue = pseudoRadDamage_ / std::pow(moduleRadius, 2);
0405         energyOnCollector *= exp(-1 * kValue * driftDistance / moduleThickness);
0406       }
0407     }
0408     LogDebug("Phase2TrackerDigitizerAlgorithm")
0409         << "Dift DistanceZ = " << driftDistance << " module thickness = " << moduleThickness
0410         << " Start Energy = " << val.energy() << " Energy after loss= " << energyOnCollector;
0411     digitizerUtility::SignalPoint sp(CloudCenterX, CloudCenterY, Sigma_x, Sigma_y, hit.tof(), energyOnCollector);
0412 
0413     // Load the Charge distribution parameters
0414     collection_points.push_back(sp);
0415   }
0416   return collection_points;
0417 }
0418 
0419 // ====================================================================
0420 //
0421 // Induce the signal on the collection plane of the active sensor area.
0422 void Phase2TrackerDigitizerAlgorithm::induce_signal(
0423     std::vector<PSimHit>::const_iterator inputBegin,
0424     const PSimHit& hit,
0425     const size_t hitIndex,
0426     const size_t firstHitIndex,
0427     const uint32_t tofBin,
0428     const Phase2TrackerGeomDetUnit* pixdet,
0429     const std::vector<digitizerUtility::SignalPoint>& collection_points) {
0430   // X  - Rows, Left-Right, 160, (1.6cm)   for barrel
0431   // Y  - Columns, Down-Up, 416, (6.4cm)
0432   const Phase2TrackerTopology* topol = &pixdet->specificTopology();
0433   uint32_t detID = pixdet->geographicalId().rawId();
0434   signal_map_type& theSignal = _signal[detID];
0435 
0436   LogDebug("Phase2TrackerDigitizerAlgorithm")
0437       << "Enter induce_signal, Pitch size is " << topol->pitch().first << " cm vs " << topol->pitch().second << " cm";
0438 
0439   // local map to store pixels hit by 1 Hit.
0440   using hit_map_type = std::map<int, float, std::less<int> >;
0441   hit_map_type hit_signal;
0442 
0443   // Assign signals to readout channels and store sorted by channel number
0444   // Iterate over collection points on the collection plane
0445   for (auto const& v : collection_points) {
0446     float CloudCenterX = v.position().x();  // Charge position in x
0447     float CloudCenterY = v.position().y();  //                 in y
0448     float SigmaX = v.sigma_x();             // Charge spread in x
0449     float SigmaY = v.sigma_y();             //               in y
0450     float Charge = v.amplitude();           // Charge amplitude
0451 
0452     LogDebug("Phase2TrackerDigitizerAlgorithm") << " cloud " << v.position().x() << " " << v.position().y() << " "
0453                                                 << v.sigma_x() << " " << v.sigma_y() << " " << v.amplitude();
0454 
0455     // Find the maximum cloud spread in 2D plane , assume 3*sigma
0456     float CloudRight = CloudCenterX + clusterWidth_ * SigmaX;
0457     float CloudLeft = CloudCenterX - clusterWidth_ * SigmaX;
0458     float CloudUp = CloudCenterY + clusterWidth_ * SigmaY;
0459     float CloudDown = CloudCenterY - clusterWidth_ * SigmaY;
0460 
0461     // Define 2D cloud limit points
0462     LocalPoint PointRightUp = LocalPoint(CloudRight, CloudUp);
0463     LocalPoint PointLeftDown = LocalPoint(CloudLeft, CloudDown);
0464 
0465     // This points can be located outside the sensor area.
0466     // The conversion to measurement point does not check for that
0467     // so the returned pixel index might be wrong (outside range).
0468     // We rely on the limits check below to fix this.
0469     // But remember whatever we do here THE CHARGE OUTSIDE THE ACTIVE
0470     // PIXEL ARE IS LOST, it should not be collected.
0471 
0472     // Convert the 2D points to pixel indices
0473     MeasurementPoint mp = topol->measurementPosition(PointRightUp);
0474     int IPixRightUpX = static_cast<int>(std::floor(mp.x()));  // cast reqd.
0475     int IPixRightUpY = static_cast<int>(std::floor(mp.y()));
0476     LogDebug("Phase2TrackerDigitizerAlgorithm")
0477         << " right-up " << PointRightUp << " " << mp.x() << " " << mp.y() << " " << IPixRightUpX << " " << IPixRightUpY;
0478 
0479     mp = topol->measurementPosition(PointLeftDown);
0480     int IPixLeftDownX = static_cast<int>(std::floor(mp.x()));
0481     int IPixLeftDownY = static_cast<int>(std::floor(mp.y()));
0482     LogDebug("Phase2TrackerDigitizerAlgorithm") << " left-down " << PointLeftDown << " " << mp.x() << " " << mp.y()
0483                                                 << " " << IPixLeftDownX << " " << IPixLeftDownY;
0484 
0485     // Check detector limits to correct for pixels outside range.
0486     int numColumns = topol->ncolumns();  // det module number of cols&rows
0487     int numRows = topol->nrows();
0488 
0489     IPixRightUpX = numRows > IPixRightUpX ? IPixRightUpX : numRows - 1;
0490     IPixRightUpY = numColumns > IPixRightUpY ? IPixRightUpY : numColumns - 1;
0491     IPixLeftDownX = 0 < IPixLeftDownX ? IPixLeftDownX : 0;
0492     IPixLeftDownY = 0 < IPixLeftDownY ? IPixLeftDownY : 0;
0493 
0494     // First integrate charge strips in x
0495     hit_map_type x;
0496     for (int ix = IPixLeftDownX; ix <= IPixRightUpX; ++ix) {  // loop over x index
0497       float xLB, LowerBound;
0498       // Why is set to 0 if ix=0, does it meen that we accept charge
0499       // outside the sensor?
0500       if (ix == 0 || SigmaX == 0.) {  // skip for surface segemnts
0501         LowerBound = 0.;
0502       } else {
0503         mp = MeasurementPoint(ix, 0.0);
0504         xLB = topol->localPosition(mp).x();
0505         LowerBound = 1 - calcQ((xLB - CloudCenterX) / SigmaX);
0506       }
0507 
0508       float xUB, UpperBound;
0509       if (ix == numRows - 1 || SigmaX == 0.) {
0510         UpperBound = 1.;
0511       } else {
0512         mp = MeasurementPoint(ix + 1, 0.0);
0513         xUB = topol->localPosition(mp).x();
0514         UpperBound = 1. - calcQ((xUB - CloudCenterX) / SigmaX);
0515       }
0516       float TotalIntegrationRange = UpperBound - LowerBound;  // get strip
0517       x.emplace(ix, TotalIntegrationRange);                   // save strip integral
0518     }
0519 
0520     // Now integrate strips in y
0521     hit_map_type y;
0522     for (int iy = IPixLeftDownY; iy <= IPixRightUpY; ++iy) {  // loop over y index
0523       float yLB, LowerBound;
0524       if (iy == 0 || SigmaY == 0.) {
0525         LowerBound = 0.;
0526       } else {
0527         mp = MeasurementPoint(0.0, iy);
0528         yLB = topol->localPosition(mp).y();
0529         LowerBound = 1. - calcQ((yLB - CloudCenterY) / SigmaY);
0530       }
0531 
0532       float yUB, UpperBound;
0533       if (iy == numColumns - 1 || SigmaY == 0.) {
0534         UpperBound = 1.;
0535       } else {
0536         mp = MeasurementPoint(0.0, iy + 1);
0537         yUB = topol->localPosition(mp).y();
0538         UpperBound = 1. - calcQ((yUB - CloudCenterY) / SigmaY);
0539       }
0540 
0541       float TotalIntegrationRange = UpperBound - LowerBound;
0542       y.emplace(iy, TotalIntegrationRange);  // save strip integral
0543     }
0544 
0545     // Get the 2D charge integrals by folding x and y strips
0546     for (int ix = IPixLeftDownX; ix <= IPixRightUpX; ++ix) {    // loop over x index
0547       for (int iy = IPixLeftDownY; iy <= IPixRightUpY; ++iy) {  // loop over y index
0548         float ChargeFraction = Charge * x[ix] * y[iy];
0549         int chanFired = -1;
0550         if (ChargeFraction > 0.) {
0551           chanFired =
0552               pixelFlag_ ? PixelDigi::pixelToChannel(ix, iy) : Phase2TrackerDigi::pixelToChannel(ix, iy);  // Get index
0553           // Load the amplitude
0554           hit_signal[chanFired] += ChargeFraction;
0555         }
0556 
0557         mp = MeasurementPoint(ix, iy);
0558         LocalPoint lp = topol->localPosition(mp);
0559         int chan = topol->channel(lp);
0560 
0561         LogDebug("Phase2TrackerDigitizerAlgorithm")
0562             << " pixel " << ix << " " << iy << " - "
0563             << " " << chanFired << " " << ChargeFraction << " " << mp.x() << " " << mp.y() << " " << lp.x() << " "
0564             << lp.y() << " "  // givex edge position
0565             << chan;          // edge belongs to previous ?
0566       }
0567     }
0568   }
0569   // Fill the global map with all hit pixels from this event
0570   bool reweighted = false;
0571   size_t referenceIndex4CR = 0;
0572 
0573   if (useChargeReweighting_) {
0574     if (hit.processType() == 0) {
0575       referenceIndex4CR = hitIndex;
0576       reweighted =
0577           theSiPixelChargeReweightingAlgorithm_->hitSignalReweight<digitizerUtility::Ph2Amplitude>(hit,
0578                                                                                                    hit_signal,
0579                                                                                                    hitIndex,
0580                                                                                                    referenceIndex4CR,
0581                                                                                                    tofBin,
0582                                                                                                    topol,
0583                                                                                                    detID,
0584                                                                                                    theSignal,
0585                                                                                                    hit.processType(),
0586                                                                                                    makeDigiSimLinks_);
0587     } else {
0588       // If it's not the primary particle, use the first hit in the collection as SimHit, which should be the corresponding primary.
0589       referenceIndex4CR = firstHitIndex;
0590       reweighted =
0591           theSiPixelChargeReweightingAlgorithm_->hitSignalReweight<digitizerUtility::Ph2Amplitude>((*inputBegin),
0592                                                                                                    hit_signal,
0593                                                                                                    hitIndex,
0594                                                                                                    referenceIndex4CR,
0595                                                                                                    tofBin,
0596                                                                                                    topol,
0597                                                                                                    detID,
0598                                                                                                    theSignal,
0599                                                                                                    hit.processType(),
0600                                                                                                    makeDigiSimLinks_);
0601     }
0602   }
0603 
0604   float corr_time = hit.tof() - pixdet->surface().toGlobal(hit.localPosition()).mag() * c_inv;
0605   if (!reweighted) {
0606     for (auto const& hit_s : hit_signal) {
0607       int chan = hit_s.first;
0608       theSignal[chan] += (makeDigiSimLinks_ ? digitizerUtility::Ph2Amplitude(
0609                                                   hit_s.second, &hit, hit_s.second, corr_time, hitIndex, tofBin)
0610                                             : digitizerUtility::Ph2Amplitude(hit_s.second, nullptr, hit_s.second));
0611     }
0612   }
0613 }  // end of induce_signal function
0614 
0615 // ======================================================================
0616 //
0617 //  Add electronic noise to pixel charge
0618 //
0619 // ======================================================================
0620 void Phase2TrackerDigitizerAlgorithm::add_noise(const Phase2TrackerGeomDetUnit* pixdet) {
0621   uint32_t detID = pixdet->geographicalId().rawId();
0622   signal_map_type& theSignal = _signal[detID];
0623   for (auto& s : theSignal) {
0624     float noise = gaussDistribution_->fire();
0625     if ((s.second.ampl() + noise) < 0.)
0626       s.second.set(0);
0627     else
0628       s.second += noise;
0629   }
0630 }
0631 
0632 // ======================================================================
0633 //
0634 //  Add  Cross-talk contribution
0635 //
0636 // ======================================================================
0637 void Phase2TrackerDigitizerAlgorithm::add_cross_talk(const Phase2TrackerGeomDetUnit* pixdet) {
0638   uint32_t detID = pixdet->geographicalId().rawId();
0639   signal_map_type& theSignal = _signal[detID];
0640   signal_map_type signalNew;
0641   const Phase2TrackerTopology* topol = &pixdet->specificTopology();
0642   int numRows = topol->nrows();
0643 
0644   for (auto& s : theSignal) {
0645     float signalInElectrons = s.second.ampl();  // signal in electrons
0646 
0647     std::pair<int, int> hitChan;
0648     if (pixelFlag_)
0649       hitChan = PixelDigi::channelToPixel(s.first);
0650     else
0651       hitChan = Phase2TrackerDigi::channelToPixel(s.first);
0652 
0653     float signalInElectrons_Xtalk = signalInElectrons * interstripCoupling_;
0654     // subtract the charge which will be shared
0655     s.second.set(signalInElectrons - signalInElectrons_Xtalk);
0656 
0657     if (hitChan.first != 0) {
0658       auto XtalkPrev = std::make_pair(hitChan.first - 1, hitChan.second);
0659       int chanXtalkPrev = pixelFlag_ ? PixelDigi::pixelToChannel(XtalkPrev.first, XtalkPrev.second)
0660                                      : Phase2TrackerDigi::pixelToChannel(XtalkPrev.first, XtalkPrev.second);
0661       signalNew.emplace(chanXtalkPrev, digitizerUtility::Ph2Amplitude(signalInElectrons_Xtalk, nullptr, -1.0));
0662     }
0663     if (hitChan.first < numRows - 1) {
0664       auto XtalkNext = std::make_pair(hitChan.first + 1, hitChan.second);
0665       int chanXtalkNext = pixelFlag_ ? PixelDigi::pixelToChannel(XtalkNext.first, XtalkNext.second)
0666                                      : Phase2TrackerDigi::pixelToChannel(XtalkNext.first, XtalkNext.second);
0667       signalNew.emplace(chanXtalkNext, digitizerUtility::Ph2Amplitude(signalInElectrons_Xtalk, nullptr, -1.0));
0668     }
0669   }
0670   for (auto const& l : signalNew) {
0671     int chan = l.first;
0672     auto iter = theSignal.find(chan);
0673     if (iter != theSignal.end()) {
0674       theSignal[chan] += l.second.ampl();
0675     } else {
0676       theSignal.emplace(chan, digitizerUtility::Ph2Amplitude(l.second.ampl(), nullptr, -1.0));
0677     }
0678   }
0679 }
0680 
0681 // ======================================================================
0682 //
0683 //  Add noise on non-hit cells
0684 //
0685 // ======================================================================
0686 void Phase2TrackerDigitizerAlgorithm::add_noisy_cells(const Phase2TrackerGeomDetUnit* pixdet, float thePixelThreshold) {
0687   uint32_t detID = pixdet->geographicalId().rawId();
0688   signal_map_type& theSignal = _signal[detID];
0689   const Phase2TrackerTopology* topol = &pixdet->specificTopology();
0690 
0691   int numColumns = topol->ncolumns();  // det module number of cols&rows
0692   int numRows = topol->nrows();
0693 
0694   int numberOfPixels = numRows * numColumns;
0695   std::map<int, float, std::less<int> > otherPixels;
0696 
0697   theNoiser_->generate(numberOfPixels,
0698                        thePixelThreshold,     //thr. in un. of nois
0699                        theNoiseInElectrons_,  // noise in elec.
0700                        otherPixels,
0701                        rengine_);
0702 
0703   LogDebug("Phase2TrackerDigitizerAlgorithm")
0704       << " Add noisy pixels " << numRows << " " << numColumns << " " << theNoiseInElectrons_ << " "
0705       << theThresholdInE_Endcap_ << "  " << theThresholdInE_Barrel_ << " " << numberOfPixels << " "
0706       << otherPixels.size();
0707 
0708   // Add noisy pixels
0709   for (auto const& el : otherPixels) {
0710     int iy = el.first / numRows;
0711     if (iy < 0 || iy > numColumns - 1)
0712       LogWarning("Phase2TrackerDigitizerAlgorithm") << " error in iy " << iy;
0713 
0714     int ix = el.first - iy * numRows;
0715     if (ix < 0 || ix > numRows - 1)
0716       LogWarning("Phase2TrackerDigitizerAlgorithm") << " error in ix " << ix;
0717 
0718     int chan = pixelFlag_ ? PixelDigi::pixelToChannel(ix, iy) : Phase2TrackerDigi::pixelToChannel(ix, iy);
0719 
0720     LogDebug("Phase2TrackerDigitizerAlgorithm")
0721         << " Storing noise = " << el.first << " " << el.second << " " << ix << " " << iy << " " << chan;
0722 
0723     if (theSignal[chan] == 0)
0724       theSignal[chan] = digitizerUtility::Ph2Amplitude(el.second, nullptr, -1.);
0725   }
0726 }
0727 // ============================================================================
0728 //
0729 // Simulate the readout inefficiencies.
0730 // Delete a selected number of single pixels, dcols and rocs.
0731 void Phase2TrackerDigitizerAlgorithm::pixel_inefficiency(const SubdetEfficiencies& eff,
0732                                                          const Phase2TrackerGeomDetUnit* pixdet,
0733                                                          const TrackerTopology* tTopo) {
0734   uint32_t detID = pixdet->geographicalId().rawId();
0735   signal_map_type& theSignal = _signal[detID];  // check validity
0736 
0737   // Predefined efficiencies
0738   float subdetEfficiency = 1.0;
0739 
0740   // setup the chip indices conversion
0741   uint32_t Subid = DetId(detID).subdetId();
0742   if (Subid == PixelSubdetector::PixelBarrel || Subid == StripSubdetector::TOB) {  // barrel layers
0743     uint32_t layerIndex = tTopo->pxbLayer(detID);
0744     if (layerIndex - 1 < eff.barrel_efficiencies.size())
0745       subdetEfficiency = eff.barrel_efficiencies[layerIndex - 1];
0746   } else {  // forward disks
0747     uint32_t diskIndex = 2 * tTopo->pxfDisk(detID) - tTopo->pxfSide(detID);
0748     if (diskIndex - 1 < eff.endcap_efficiencies.size())
0749       subdetEfficiency = eff.endcap_efficiencies[diskIndex - 1];
0750   }
0751 
0752   LogDebug("Phase2TrackerDigitizerAlgorithm") << " enter pixel_inefficiency " << subdetEfficiency;
0753 
0754   // Now loop again over pixels to kill some of them.
0755   // Loop over hits, amplitude in electrons, channel = coded row,col
0756   for (auto& s : theSignal) {
0757     float rand = rengine_->flat();
0758     if (rand > subdetEfficiency) {
0759       // make amplitude =0
0760       s.second.set(0.);  // reset amplitude
0761     }
0762   }
0763 }
0764 
0765 void Phase2TrackerDigitizerAlgorithm::initializeEvent(CLHEP::HepRandomEngine& eng) {
0766   if (addNoise_ || addPixelInefficiency_ || fluctuateCharge_ || addThresholdSmearing_) {
0767     gaussDistribution_ = std::make_unique<CLHEP::RandGaussQ>(eng, 0., theNoiseInElectrons_);
0768   }
0769   // Threshold smearing with gaussian distribution:
0770   if (addThresholdSmearing_) {
0771     smearedThreshold_Endcap_ =
0772         std::make_unique<CLHEP::RandGaussQ>(eng, theThresholdInE_Endcap_, theThresholdSmearing_Endcap_);
0773     smearedThreshold_Barrel_ =
0774         std::make_unique<CLHEP::RandGaussQ>(eng, theThresholdInE_Barrel_, theThresholdSmearing_Barrel_);
0775   }
0776   rengine_ = &eng;
0777   _signal.clear();
0778 }
0779 
0780 // =======================================================================================
0781 //
0782 // Set the drift direction accoring to the Bfield in local det-unit frame
0783 // Works for both barrel and forward pixels.
0784 // Replace the sign convention to fit M.Swartz's formulaes.
0785 // Configurations for barrel and foward pixels possess different tanLorentzAngleperTesla
0786 // parameter value
0787 
0788 LocalVector Phase2TrackerDigitizerAlgorithm::driftDirection(const Phase2TrackerGeomDetUnit* pixdet,
0789                                                             const GlobalVector& bfield,
0790                                                             const DetId& detId) const {
0791   Frame detFrame(pixdet->surface().position(), pixdet->surface().rotation());
0792   LocalVector Bfield = detFrame.toLocal(bfield);
0793 
0794   float dir_x = 0.0;
0795   float dir_y = 0.0;
0796   float dir_z = 0.0;
0797   float scale = 1.0;
0798 
0799   uint32_t detID = pixdet->geographicalId().rawId();
0800   uint32_t Sub_detid = DetId(detID).subdetId();
0801 
0802   // Read Lorentz angle from DB:
0803   if (use_LorentzAngle_DB_) {
0804     bool isPixel = (Sub_detid == PixelSubdetector::PixelBarrel || Sub_detid == PixelSubdetector::PixelEndcap);
0805     if (isPixel) {
0806       LogDebug("Phase2TrackerDigitizerAlgorithm") << "Read Lorentz angle from DB for Phase-2 IT";
0807     } else {
0808       LogDebug("Phase2TrackerDigitizerAlgorithm") << "Read Lorentz angle from DB for Phase-2 OT";
0809     }
0810 
0811     float lorentzAngle =
0812         isPixel ? siPixelLorentzAngle_->getLorentzAngle(detId) : siPhase2OTLorentzAngle_->getLorentzAngle(detId);
0813     float alpha2 = std::pow(lorentzAngle, 2);
0814 
0815     dir_x = -(lorentzAngle * Bfield.y() + alpha2 * Bfield.z() * Bfield.x());
0816     dir_y = +(lorentzAngle * Bfield.x() - alpha2 * Bfield.z() * Bfield.y());
0817     dir_z = -(1 + alpha2 * std::pow(Bfield.z(), 2));
0818     scale = (1 + alpha2 * std::pow(Bfield.z(), 2));
0819   } else {
0820     // Read Lorentz angle from cfg file:
0821     LogDebug("Phase2TrackerDigitizerAlgorithm") << "Read Lorentz angle from cfg file";
0822     float alpha2_Endcap = 0.0;
0823     float alpha2_Barrel = 0.0;
0824     if (alpha2Order_) {
0825       alpha2_Endcap = std::pow(tanLorentzAnglePerTesla_Endcap_, 2);
0826       alpha2_Barrel = std::pow(tanLorentzAnglePerTesla_Barrel_, 2);
0827     }
0828 
0829     if (Sub_detid == PixelSubdetector::PixelBarrel || Sub_detid == StripSubdetector::TOB) {  // barrel layers
0830       dir_x = -(tanLorentzAnglePerTesla_Barrel_ * Bfield.y() + alpha2_Barrel * Bfield.z() * Bfield.x());
0831       dir_y = +(tanLorentzAnglePerTesla_Barrel_ * Bfield.x() - alpha2_Barrel * Bfield.z() * Bfield.y());
0832       dir_z = -(1 + alpha2_Barrel * std::pow(Bfield.z(), 2));
0833       scale = (1 + alpha2_Barrel * std::pow(Bfield.z(), 2));
0834 
0835     } else {  // forward disks
0836       dir_x = -(tanLorentzAnglePerTesla_Endcap_ * Bfield.y() + alpha2_Endcap * Bfield.z() * Bfield.x());
0837       dir_y = +(tanLorentzAnglePerTesla_Endcap_ * Bfield.x() - alpha2_Endcap * Bfield.z() * Bfield.y());
0838       dir_z = -(1 + alpha2_Endcap * std::pow(Bfield.z(), 2));
0839       scale = (1 + alpha2_Endcap * std::pow(Bfield.z(), 2));
0840     }
0841   }
0842 
0843   LocalVector theDriftDirection = LocalVector(dir_x / scale, dir_y / scale, dir_z / scale);
0844   LogDebug("Phase2TrackerDigitizerAlgorithm") << " The drift direction in local coordinate is " << theDriftDirection;
0845   return theDriftDirection;
0846 }
0847 
0848 // =============================================================================
0849 
0850 void Phase2TrackerDigitizerAlgorithm::pixel_inefficiency_db(uint32_t detID) {
0851   signal_map_type& theSignal = _signal[detID];  // check validity
0852 
0853   // Loop over hit pixels, amplitude in electrons, channel = coded row,col
0854   for (auto& s : theSignal) {
0855     std::pair<int, int> ip;
0856     if (pixelFlag_)
0857       ip = PixelDigi::channelToPixel(s.first);  //get pixel pos
0858     else
0859       ip = Phase2TrackerDigi::channelToPixel(s.first);  //get pixel pos
0860 
0861     int row = ip.first;   // X in row
0862     int col = ip.second;  // Y is in col
0863     // transform to ROC index coordinates
0864     if (theSiPixelGainCalibrationService_->isDead(detID, col, row))
0865       s.second.set(0.);  // reset amplitude
0866   }
0867 }
0868 
0869 // ==========================================================================
0870 
0871 void Phase2TrackerDigitizerAlgorithm::module_killing_conf(uint32_t detID) {
0872   bool isbad = false;
0873   int detid = detID;
0874   std::string Module;
0875   for (auto const& det_m : deadModules_) {
0876     int Dead_detID = det_m.getParameter<int>("Dead_detID");
0877     Module = det_m.getParameter<std::string>("Module");
0878     if (detid == Dead_detID) {
0879       isbad = true;
0880       break;
0881     }
0882   }
0883 
0884   if (!isbad)
0885     return;
0886 
0887   signal_map_type& theSignal = _signal[detID];  // check validity
0888   for (auto& s : theSignal) {
0889     std::pair<int, int> ip;
0890     if (pixelFlag_)
0891       ip = PixelDigi::channelToPixel(s.first);
0892     else
0893       ip = Phase2TrackerDigi::channelToPixel(s.first);  //get pixel pos
0894 
0895     if (Module == "whole" || (Module == "tbmA" && ip.first >= 80 && ip.first <= 159) ||
0896         (Module == "tbmB" && ip.first <= 79))
0897       s.second.set(0.);
0898   }
0899 }
0900 // For premixing
0901 void Phase2TrackerDigitizerAlgorithm::loadAccumulator(uint32_t detId, const std::map<int, float>& accumulator) {
0902   auto& theSignal = _signal[detId];
0903   // the input channel is always with PixelDigi definition
0904   // if needed, that has to be converted to Phase2TrackerDigi convention
0905   for (const auto& elem : accumulator) {
0906     auto inserted = theSignal.emplace(elem.first, digitizerUtility::Ph2Amplitude(elem.second, nullptr));
0907     if (!inserted.second) {
0908       throw cms::Exception("LogicError") << "Signal was already set for DetId " << detId;
0909     }
0910   }
0911 }
0912 
0913 void Phase2TrackerDigitizerAlgorithm::digitize(const Phase2TrackerGeomDetUnit* pixdet,
0914                                                std::map<int, digitizerUtility::DigiSimInfo>& digi_map,
0915                                                const TrackerTopology* tTopo) {
0916   uint32_t detID = pixdet->geographicalId().rawId();
0917   auto it = _signal.find(detID);
0918   if (it == _signal.end())
0919     return;
0920 
0921   const signal_map_type& theSignal = _signal[detID];
0922 
0923   uint32_t Sub_detid = DetId(detID).subdetId();
0924 
0925   float theThresholdInE = 0.;
0926   float theHIPThresholdInE = 0.;
0927   // Define Threshold
0928   if (Sub_detid == PixelSubdetector::PixelBarrel || Sub_detid == StripSubdetector::TOB) {  // Barrel modules
0929     theThresholdInE = addThresholdSmearing_ ? smearedThreshold_Barrel_->fire()             // gaussian smearing
0930                                             : theThresholdInE_Barrel_;                     // no smearing
0931     theHIPThresholdInE = theHIPThresholdInE_Barrel_;
0932   } else {                                                                      // Forward disks modules
0933     theThresholdInE = addThresholdSmearing_ ? smearedThreshold_Endcap_->fire()  // gaussian smearing
0934                                             : theThresholdInE_Endcap_;          // no smearing
0935     theHIPThresholdInE = theHIPThresholdInE_Endcap_;
0936   }
0937 
0938   //  if (addNoise) add_noise(pixdet, theThresholdInE/theNoiseInElectrons_);  // generate noise
0939   if (addNoise_)
0940     add_noise(pixdet);  // generate noise
0941   if (addXtalk_)
0942     add_cross_talk(pixdet);
0943   if (addNoisyPixels_) {
0944     float thresholdInNoiseUnits = 99.9;
0945     if (theNoiseInElectrons_)
0946       thresholdInNoiseUnits = theThresholdInE / theNoiseInElectrons_;
0947     add_noisy_cells(pixdet, thresholdInNoiseUnits);
0948   }
0949 
0950   // Do only if needed
0951   if (addPixelInefficiency_ && !theSignal.empty()) {
0952     if (use_ineff_from_db_)
0953       pixel_inefficiency_db(detID);
0954     else
0955       pixel_inefficiency(subdetEfficiencies_, pixdet, tTopo);
0956   }
0957   if (use_module_killing_) {
0958     if (use_deadmodule_DB_)  // remove dead modules using DB
0959       module_killing_DB(pixdet);
0960     else  // remove dead modules using the list in cfg file
0961       module_killing_conf(detID);
0962   }
0963 
0964   // Digitize if the signal is greater than threshold
0965   for (auto const& s : theSignal) {
0966     const digitizerUtility::Ph2Amplitude& sig_data = s.second;
0967     float signalInElectrons = sig_data.ampl();
0968 
0969     const auto& info_list = sig_data.simInfoList();
0970     const digitizerUtility::SimHitInfo* hitInfo = nullptr;
0971     if (!info_list.empty())
0972       hitInfo = std::max_element(info_list.begin(), info_list.end())->second.get();
0973 
0974     if (isAboveThreshold(hitInfo, signalInElectrons, theThresholdInE)) {  // check threshold
0975       digitizerUtility::DigiSimInfo info;
0976       info.sig_tot = convertSignalToAdc(detID, signalInElectrons, theThresholdInE);  // adc
0977       info.ot_bit = signalInElectrons > theHIPThresholdInE ? true : false;
0978       if (makeDigiSimLinks_) {
0979         for (auto const& l : sig_data.simInfoList()) {
0980           float charge_frac = l.first / signalInElectrons;
0981           if (l.first > -5.0)
0982             info.simInfoList.push_back({charge_frac, l.second.get()});
0983         }
0984       }
0985       digi_map.insert({s.first, info});
0986     }
0987   }
0988 }
0989 //
0990 // Scale the Signal using Dual Slope option
0991 //
0992 int Phase2TrackerDigitizerAlgorithm::convertSignalToAdc(uint32_t detID, float signal_in_elec, float threshold) {
0993   int signal_in_adc;
0994   int temp_signal;
0995   const int max_limit = 10;
0996   if (thePhase2ReadoutMode_ == 0) {
0997     signal_in_adc = theAdcFullScale_;
0998   } else {
0999     if (thePhase2ReadoutMode_ == -1) {
1000       temp_signal = std::min(static_cast<int>(std::round(signal_in_elec / theElectronPerADC_)), theAdcFullScale_);
1001     } else {
1002       // calculate the kink point and the slope
1003       int dualslope_param = std::min(std::abs(thePhase2ReadoutMode_), max_limit);
1004       int kink_point = static_cast<int>(theAdcFullScale_ / 2) + 1;
1005       // C-ROC: first valid ToT code above threshold is 0b0000
1006       temp_signal = std::floor((signal_in_elec - threshold) / theElectronPerADC_);
1007       if (temp_signal > kink_point)
1008         temp_signal =
1009             static_cast<int>(std::floor((temp_signal - kink_point) / (pow(2, dualslope_param - 1)))) + kink_point;
1010     }
1011     signal_in_adc = std::min(temp_signal, theAdcFullScale_);
1012     LogTrace("Phase2TrackerDigitizerAlgorithm")
1013         << " DetId " << detID << " signal_in_elec " << signal_in_elec << " threshold " << threshold
1014         << " signal_above_thr " << signal_in_elec - threshold << " temp conversion "
1015         << std::floor((signal_in_elec - threshold) / theElectronPerADC_) + 1 << " signal after slope correction "
1016         << temp_signal << " signal_in_adc " << signal_in_adc;
1017   }
1018   return signal_in_adc;
1019 }
1020 
1021 float Phase2TrackerDigitizerAlgorithm::calcQ(float x) {
1022   constexpr float p1 = 12.5f;
1023   constexpr float p2 = 0.2733f;
1024   constexpr float p3 = 0.147f;
1025 
1026   auto xx = std::min(0.5f * x * x, p1);
1027   return 0.5f * (1.f - std::copysign(std::sqrt(1.f - unsafe_expf<4>(-xx * (1.f + p2 / (1.f + p3 * xx)))), x));
1028 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!