Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-12-21 00:09:15

0001 #ifndef SimTracker_SiPixelDigitizer_SiPixelDigitizerAlgorithm_h
0002 #define SimTracker_SiPixelDigitizer_SiPixelDigitizerAlgorithm_h
0003 
0004 #include <map>
0005 #include <memory>
0006 #include <vector>
0007 #include <iostream>
0008 #include "DataFormats/GeometrySurface/interface/GloballyPositioned.h"
0009 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0010 #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
0011 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0012 #include "SimTracker/Common/interface/DigitizerUtility.h"
0013 #include "SimTracker/SiPixelDigitizer/plugins/PixelDigiAddTempInfo.h"
0014 #include "SimDataFormats/TrackerDigiSimLink/interface/PixelSimHitExtraInfo.h"
0015 #include "DataFormats/Math/interface/approx_exp.h"
0016 #include "SimDataFormats/PileupSummaryInfo/interface/PileupMixingContent.h"
0017 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
0018 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0019 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0020 #include "CondFormats/SiPixelTransient/interface/SiPixelTemplate2D.h"
0021 #include "DataFormats/SiPixelDetId/interface/PixelFEDChannel.h"
0022 #include "CalibTracker/Records/interface/SiPixelFEDChannelContainerESProducerRcd.h"
0023 #include "CondFormats/DataRecord/interface/SiPixelQualityRcd.h"
0024 #include "CondFormats/DataRecord/interface/SiPixelFedCablingMapRcd.h"
0025 #include "CondFormats/DataRecord/interface/SiPixelLorentzAngleSimRcd.h"
0026 #include "CondFormats/DataRecord/interface/SiPixelDynamicInefficiencyRcd.h"
0027 #include "CondFormats/DataRecord/interface/SiPixelStatusScenarioProbabilityRcd.h"
0028 #include "CondFormats/DataRecord/interface/SiPixelStatusScenariosRcd.h"
0029 #include "FWCore/Framework/interface/FrameworkfwdMostUsed.h"
0030 #include "boost/multi_array.hpp"
0031 
0032 typedef boost::multi_array<float, 2> array_2d;
0033 
0034 // forward declarations
0035 
0036 // For the random numbers
0037 namespace CLHEP {
0038   class HepRandomEngine;
0039 }
0040 
0041 class DetId;
0042 class GaussianTailNoiseGenerator;
0043 class PixelDigi;
0044 class PixelDigiSimLink;
0045 class PixelGeomDetUnit;
0046 class SiG4UniversalFluctuation;
0047 class SiPixelFedCablingMap;
0048 class SiPixelGainCalibrationOfflineSimService;
0049 class SiPixelLorentzAngle;
0050 class SiPixelQuality;
0051 class SiPixelDynamicInefficiency;
0052 class TrackerGeometry;
0053 class TrackerTopology;
0054 class SiPixelFEDChannelContainer;
0055 class SiPixelQualityProbabilities;
0056 class SiPixelChargeReweightingAlgorithm;
0057 
0058 class SiPixelDigitizerAlgorithm {
0059 public:
0060   SiPixelDigitizerAlgorithm(const edm::ParameterSet& conf, edm::ConsumesCollector iC);
0061   ~SiPixelDigitizerAlgorithm();
0062 
0063   // initialization that cannot be done in the constructor
0064   void init(const edm::EventSetup& es);
0065 
0066   void initializeEvent() { _signal.clear(); }
0067 
0068   //run the algorithm to digitize a single det
0069   void accumulateSimHits(const std::vector<PSimHit>::const_iterator inputBegin,
0070                          const std::vector<PSimHit>::const_iterator inputEnd,
0071                          const size_t inputBeginGlobalIndex,
0072                          const unsigned int tofBin,
0073                          const PixelGeomDetUnit* pixdet,
0074                          const GlobalVector& bfield,
0075                          const TrackerTopology* tTopo,
0076                          CLHEP::HepRandomEngine*);
0077   void digitize(const PixelGeomDetUnit* pixdet,
0078                 std::vector<PixelDigi>& digis,
0079                 std::vector<PixelDigiSimLink>& simlinks,
0080                 std::vector<PixelDigiAddTempInfo>& newClass_Digi_extra,
0081                 const TrackerTopology* tTopo,
0082                 CLHEP::HepRandomEngine*);
0083   void calculateInstlumiFactor(PileupMixingContent* puInfo);
0084   void fillSimHitMaps(std::vector<PSimHit> simHits, const unsigned int tofBin);
0085   void resetSimHitMaps();
0086   void init_DynIneffDB(const edm::EventSetup&);
0087   std::unique_ptr<PixelFEDChannelCollection> chooseScenario(PileupMixingContent* puInfo, CLHEP::HepRandomEngine*);
0088 
0089   void lateSignalReweight(const PixelGeomDetUnit* pixdet,
0090                           std::vector<PixelDigi>& digis,
0091                           std::vector<PixelSimHitExtraInfo>& newClass_Sim_extra,
0092                           const TrackerTopology* tTopo,
0093                           CLHEP::HepRandomEngine* engine);
0094 
0095   // for premixing
0096   void calculateInstlumiFactor(const std::vector<PileupSummaryInfo>& ps,
0097                                int bunchSpacing);  // TODO: try to remove the duplication of logic...
0098   void setSimAccumulator(const std::map<uint32_t, std::map<int, int> >& signalMap);
0099   std::unique_ptr<PixelFEDChannelCollection> chooseScenario(const std::vector<PileupSummaryInfo>& ps,
0100                                                             CLHEP::HepRandomEngine* engine);
0101 
0102   bool killBadFEDChannels() const;
0103   typedef std::unordered_map<std::string, PixelFEDChannelCollection> PixelFEDChannelCollectionMap;
0104   const PixelFEDChannelCollectionMap* quality_map;
0105 
0106 private:
0107   //Accessing Lorentz angle from DB:
0108   edm::ESGetToken<SiPixelLorentzAngle, SiPixelLorentzAngleSimRcd> SiPixelLorentzAngleToken_;
0109   const SiPixelLorentzAngle* SiPixelLorentzAngle_ = nullptr;
0110 
0111   //Accessing Dead pixel modules from DB:
0112   edm::ESGetToken<SiPixelQuality, SiPixelQualityRcd> SiPixelBadModuleToken_;
0113   const SiPixelQuality* SiPixelBadModule_ = nullptr;
0114 
0115   //Accessing Map and Geom:
0116   const edm::ESGetToken<SiPixelFedCablingMap, SiPixelFedCablingMapRcd> mapToken_;
0117   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
0118   const SiPixelFedCablingMap* map_ = nullptr;
0119   const TrackerGeometry* geom_ = nullptr;
0120 
0121   // Get Dynamic Inefficiency scale factors from DB
0122   edm::ESGetToken<SiPixelDynamicInefficiency, SiPixelDynamicInefficiencyRcd> SiPixelDynamicInefficiencyToken_;
0123   const SiPixelDynamicInefficiency* SiPixelDynamicInefficiency_ = nullptr;
0124 
0125   // For BadFEDChannel simulation
0126   edm::ESGetToken<SiPixelQualityProbabilities, SiPixelStatusScenarioProbabilityRcd> scenarioProbabilityToken_;
0127   edm::ESGetToken<PixelFEDChannelCollectionMap, SiPixelFEDChannelContainerESProducerRcd>
0128       PixelFEDChannelCollectionMapToken_;
0129   const SiPixelQualityProbabilities* scenarioProbability_ = nullptr;
0130   // Define internal classes
0131 
0132   // definition class
0133   //
0134 
0135   // Define a class to hold the calibration parameters per pixel
0136   // Internal
0137   class CalParameters {
0138   public:
0139     float p0;
0140     float p1;
0141     float p2;
0142     float p3;
0143   };
0144   //
0145   // Define a class for 3D ionization points and energy
0146   //
0147   /**
0148    * Internal use only.
0149    */
0150   class EnergyDepositUnit {
0151   public:
0152     EnergyDepositUnit() : _energy(0), _position(0, 0, 0) {}
0153     EnergyDepositUnit(float energy, float x, float y, float z) : _energy(energy), _position(x, y, z) {}
0154     EnergyDepositUnit(float energy, Local3DPoint position) : _energy(energy), _position(position) {}
0155     float x() const { return _position.x(); }
0156     float y() const { return _position.y(); }
0157     float z() const { return _position.z(); }
0158     float energy() const { return _energy; }
0159 
0160   private:
0161     float _energy;
0162     Local3DPoint _position;
0163   };
0164 
0165   //
0166   // define class to store signals on the collection surface
0167   //
0168   /**
0169    * Internal use only.
0170    */
0171 
0172   class SignalPoint {
0173   public:
0174     SignalPoint() : _pos(0, 0), _time(0), _amplitude(0), _sigma_x(1.), _sigma_y(1.), _hitp(nullptr) {}
0175 
0176     SignalPoint(float x, float y, float sigma_x, float sigma_y, float t, float a = 1.0)
0177         : _pos(x, y), _time(t), _amplitude(a), _sigma_x(sigma_x), _sigma_y(sigma_y), _hitp(nullptr) {}
0178 
0179     SignalPoint(float x, float y, float sigma_x, float sigma_y, float t, const PSimHit& hit, float a = 1.0)
0180         : _pos(x, y), _time(t), _amplitude(a), _sigma_x(sigma_x), _sigma_y(sigma_y), _hitp(&hit) {}
0181 
0182     const LocalPoint& position() const { return _pos; }
0183     float x() const { return _pos.x(); }
0184     float y() const { return _pos.y(); }
0185     float sigma_x() const { return _sigma_x; }
0186     float sigma_y() const { return _sigma_y; }
0187     float time() const { return _time; }
0188     float amplitude() const { return _amplitude; }
0189     const PSimHit& hit() { return *_hitp; }
0190     SignalPoint& set_amplitude(float amp) {
0191       _amplitude = amp;
0192       return *this;
0193     }
0194 
0195   private:
0196     LocalPoint _pos;
0197     float _time;
0198     float _amplitude;
0199     float _sigma_x;  // gaussian sigma in the x direction (cm)
0200     float _sigma_y;  //    "       "          y direction (cm) */
0201     const PSimHit* _hitp;
0202   };
0203 
0204   //
0205   // PixelEfficiencies struct
0206   //
0207   /**
0208    * Internal use only.
0209    */
0210   struct PixelEfficiencies {
0211     PixelEfficiencies(const edm::ParameterSet& conf,
0212                       bool AddPixelInefficiency,
0213                       int NumberOfBarrelLayers,
0214                       int NumberOfEndcapDisks);
0215     bool FromConfig;  // If true read from Config, otherwise use Database
0216 
0217     double theInstLumiScaleFactor;
0218     std::vector<double> pu_scale;                       // in config: 0-3 BPix, 4-5 FPix (inner, outer)
0219     std::vector<std::vector<double> > thePUEfficiency;  // Instlumi dependent efficiency
0220 
0221     // Read factors from Configuration
0222     double thePixelEfficiency[20];                     // Single pixel effciency
0223     double thePixelColEfficiency[20];                  // Column effciency
0224     double thePixelChipEfficiency[20];                 // ROC efficiency
0225     std::vector<double> theLadderEfficiency_BPix[20];  // Ladder efficiency
0226     std::vector<double> theModuleEfficiency_BPix[20];  // Module efficiency
0227     double theInnerEfficiency_FPix[20];                // Fpix inner module efficiency
0228     double theOuterEfficiency_FPix[20];                // Fpix outer module efficiency
0229     unsigned int FPixIndex;                            // The Efficiency index for FPix Disks
0230 
0231     // Read factors from DB and fill containers
0232     std::map<uint32_t, double> PixelGeomFactors;
0233     std::map<uint32_t, std::vector<double> > PixelGeomFactorsROCStdPixels;
0234     std::map<uint32_t, std::vector<double> > PixelGeomFactorsROCBigPixels;
0235     std::map<uint32_t, double> ColGeomFactors;
0236     std::map<uint32_t, double> ChipGeomFactors;
0237     std::map<uint32_t, size_t> iPU;
0238 
0239     // constants for ROC level simulation for Phase1
0240     enum shiftEnumerator { FPixRocIdShift = 3, BPixRocIdShift = 6 };
0241     static const int rocIdMaskBits = 0x1F;
0242     void init_from_db(const TrackerGeometry*, const SiPixelDynamicInefficiency*);
0243     bool matches(const DetId&, const DetId&, const std::vector<uint32_t>&);
0244     std::unique_ptr<PixelFEDChannelCollection> PixelFEDChannelCollection_;
0245   };
0246 
0247   //
0248   // PixelAging struct
0249   //
0250   /**
0251     * Internal use only.
0252     */
0253   struct PixelAging {
0254     PixelAging(const edm::ParameterSet& conf, bool AddPixelAging, int NumberOfBarrelLayers, int NumberOfEndcapDisks);
0255     float thePixelPseudoRadDamage[20];  // PseudoRadiation Damage Values for aging studies
0256     unsigned int FPixIndex;             // The Efficiency index for FPix Disks
0257   };
0258 
0259 private:
0260   // Internal typedefs
0261   typedef std::map<int, digitizerUtility::Amplitude, std::less<int> > signal_map_type;  // from Digi.Skel.
0262   typedef signal_map_type::iterator signal_map_iterator;                                // from Digi.Skel.
0263   typedef signal_map_type::const_iterator signal_map_const_iterator;                    // from Digi.Skel.
0264   typedef std::map<uint32_t, signal_map_type> signalMaps;
0265   typedef GloballyPositioned<double> Frame;
0266   typedef std::vector<edm::ParameterSet> Parameters;
0267   typedef boost::multi_array<float, 2> array_2d;
0268 
0269   typedef std::pair<unsigned int, unsigned int> subDetTofBin;
0270   typedef std::map<unsigned int, std::vector<PSimHit> > simhit_map;
0271   simhit_map SimHitMap;
0272   typedef std::map<subDetTofBin, unsigned int> simhit_collectionMap;
0273   simhit_collectionMap SimHitCollMap;
0274 
0275   // Contains the accumulated hit info.
0276   signalMaps _signal;
0277 
0278   const bool makeDigiSimLinks_;
0279   const bool store_SimHitEntryExitPoints_;
0280 
0281   const bool use_ineff_from_db_;
0282   const bool use_module_killing_;   // remove or not the dead pixel modules
0283   const bool use_deadmodule_DB_;    // if we want to get dead pixel modules from the DataBase.
0284   const bool use_LorentzAngle_DB_;  // if we want to get Lorentz angle from the DataBase.
0285 
0286   const Parameters DeadModules;
0287 
0288   std::unique_ptr<SiPixelChargeReweightingAlgorithm> TheNewSiPixelChargeReweightingAlgorithmClass;
0289 
0290 private:
0291   // Variables
0292   //external parameters
0293   // go from Geant energy GeV to number of electrons
0294   const float GeVperElectron;  // 3.7E-09
0295 
0296   //-- drift
0297   const float Sigma0;      //=0.0007  // Charge diffusion in microns for 300 micron Si
0298   const float Dist300;     //=0.0300  // Define 300microns for normalization
0299   const bool alpha2Order;  // Switch on/off of E.B effect
0300 
0301   //-- induce_signal
0302   const float ClusterWidth;  // Gaussian charge cutoff width in sigma units
0303   //-- Allow for upgrades
0304   const int NumberOfBarrelLayers;  // Default = 3
0305   const int NumberOfEndcapDisks;   // Default = 2
0306 
0307   //-- make_digis
0308   const float theElectronPerADC;    // Gain, number of electrons per adc count.
0309   const int theAdcFullScale;        // Saturation count, 255=8bit.
0310   const int theAdcFullScLateCR;     // Saturation count, 255=8bit.
0311   const float theNoiseInElectrons;  // Noise (RMS) in units of electrons.
0312   const float theReadoutNoise;      // Noise of the readount chain in elec,
0313                                     //inludes DCOL-Amp,TBM-Amp, Alt, AOH,OptRec.
0314 
0315   const float theThresholdInE_FPix;     // Pixel threshold in electrons FPix.
0316   const float theThresholdInE_BPix;     // Pixel threshold in electrons BPix.
0317   const float theThresholdInE_BPix_L1;  // In case the BPix layer1 gets a different threshold
0318   const float theThresholdInE_BPix_L2;  // In case the BPix layer2 gets a different threshold
0319 
0320   const double theThresholdSmearing_FPix;
0321   const double theThresholdSmearing_BPix;
0322   const double theThresholdSmearing_BPix_L1;
0323   const double theThresholdSmearing_BPix_L2;
0324 
0325   const float electronsPerVCAL;            // for electrons - VCAL conversion
0326   const float electronsPerVCAL_Offset;     // in misscalibrate()
0327   const float electronsPerVCAL_L1;         // same for Layer 1
0328   const float electronsPerVCAL_L1_Offset;  // same for Layer 1
0329 
0330   const float theTofLowerCut;                // Cut on the particle TOF
0331   const float theTofUpperCut;                // Cut on the particle TOF
0332   const float tanLorentzAnglePerTesla_FPix;  //FPix Lorentz angle tangent per Tesla
0333   const float tanLorentzAnglePerTesla_BPix;  //BPix Lorentz angle tangent per Tesla
0334 
0335   const float FPix_p0;
0336   const float FPix_p1;
0337   const float FPix_p2;
0338   const float FPix_p3;
0339   const float BPix_p0;
0340   const float BPix_p1;
0341   const float BPix_p2;
0342   const float BPix_p3;
0343 
0344   //-- add_noise
0345   const bool addNoise;
0346   const bool addChargeVCALSmearing;
0347   const bool addNoisyPixels;
0348   const bool fluctuateCharge;
0349 
0350   //-- pixel efficiency
0351   const bool AddPixelInefficiency;  // bool to read in inefficiencies
0352   const bool KillBadFEDChannels;
0353   const bool addThresholdSmearing;
0354 
0355   //-- calibration smearing
0356   const bool doMissCalibrate;     // Switch on the calibration smearing
0357   const bool doMissCalInLateCR;   // Switch on the calibration smearing
0358   const float theGainSmearing;    // The sigma of the gain fluctuation (around 1)
0359   const float theOffsetSmearing;  // The sigma of the offset fluct. (around 0)
0360 
0361   // pixel aging
0362   const bool AddPixelAging;
0363   const bool UseReweighting;
0364 
0365   // The PDTable
0366   //HepPDTable *particleTable;
0367   //ParticleDataTable *particleTable;
0368 
0369   //-- charge fluctuation
0370   const double tMax;  // The delta production cut, should be as in OSCAR = 30keV
0371   //                                           cmsim = 100keV
0372 
0373   // The eloss fluctuation class from G4. Is the right place?
0374   const std::unique_ptr<SiG4UniversalFluctuation> fluctuate;  // make a pointer
0375   const std::unique_ptr<GaussianTailNoiseGenerator> theNoiser;
0376 
0377   // To store calibration constants
0378   const std::map<int, CalParameters, std::less<int> > calmap;
0379 
0380   //-- additional member functions
0381   // Private methods
0382   std::map<int, CalParameters, std::less<int> > initCal() const;
0383   void primary_ionization(const PSimHit& hit,
0384                           std::vector<EnergyDepositUnit>& ionization_points,
0385                           CLHEP::HepRandomEngine*) const;
0386   void drift(const PSimHit& hit,
0387              const PixelGeomDetUnit* pixdet,
0388              const GlobalVector& bfield,
0389              const TrackerTopology* tTopo,
0390              const std::vector<EnergyDepositUnit>& ionization_points,
0391              std::vector<SignalPoint>& collection_points) const;
0392   void induce_signal(std::vector<PSimHit>::const_iterator inputBegin,
0393                      std::vector<PSimHit>::const_iterator inputEnd,
0394                      const PSimHit& hit,
0395                      const size_t hitIndex,
0396                      const size_t FirstHitIndex,
0397                      const unsigned int tofBin,
0398                      const PixelGeomDetUnit* pixdet,
0399                      const std::vector<SignalPoint>& collection_points);
0400   void fluctuateEloss(int particleId,
0401                       float momentum,
0402                       float eloss,
0403                       float length,
0404                       int NumberOfSegments,
0405                       float elossVector[],
0406                       CLHEP::HepRandomEngine*) const;
0407   void add_noise(const PixelGeomDetUnit* pixdet, float thePixelThreshold, CLHEP::HepRandomEngine*);
0408   void make_digis(float thePixelThresholdInE,
0409                   uint32_t detID,
0410                   const PixelGeomDetUnit* pixdet,
0411                   std::vector<PixelDigi>& digis,
0412                   std::vector<PixelDigiSimLink>& simlinks,
0413                   std::vector<PixelDigiAddTempInfo>& newClass_Digi_extra,
0414                   const TrackerTopology* tTopo) const;
0415   void pixel_inefficiency(const PixelEfficiencies& eff,
0416                           const PixelGeomDetUnit* pixdet,
0417                           const TrackerTopology* tTopo,
0418                           CLHEP::HepRandomEngine*);
0419 
0420   void pixel_inefficiency_db(uint32_t detID);
0421 
0422   float pixel_aging(const PixelAging& aging, const PixelGeomDetUnit* pixdet, const TrackerTopology* tTopo) const;
0423 
0424   // access to the gain calibration payloads in the db. Only gets initialized if check_dead_pixels_ is set to true.
0425   const std::unique_ptr<SiPixelGainCalibrationOfflineSimService> theSiPixelGainCalibrationService_;
0426   float missCalibrate(
0427       uint32_t detID, const TrackerTopology* tTopo, const PixelGeomDetUnit* pixdet, int col, int row, float amp) const;
0428   LocalVector DriftDirection(const PixelGeomDetUnit* pixdet, const GlobalVector& bfield, const DetId& detId) const;
0429 
0430   void module_killing_conf(
0431       uint32_t detID);  // remove dead modules using the list in the configuration file PixelDigi_cfi.py
0432   void module_killing_DB(uint32_t detID);  // remove dead modules uisng the list in the DB
0433 
0434   PixelEfficiencies pixelEfficiencies_;
0435   const PixelAging pixelAging_;
0436 
0437   double calcQ(float x) const {
0438     // need erf(x/sqrt2)
0439     //float x2=0.5*x*x;
0440     //float a=0.147;
0441     //double erf=sqrt(1.0f-exp( -1.0f*x2*( (4/M_PI)+a*x2)/(1.0+a*x2)));
0442     //if (x<0.) erf*=-1.0;
0443     //return 0.5*(1.0-erf);
0444 
0445     auto xx = std::min(0.5f * x * x, 12.5f);
0446     return 0.5 * (1.0 - std::copysign(std::sqrt(1.f - unsafe_expf<4>(-xx * (1.f + 0.2733f / (1.f + 0.147f * xx)))), x));
0447   }
0448 };
0449 
0450 #endif