CMS 3D CMS Logo

SiPixelDigitizerAlgorithm.h
Go to the documentation of this file.
1 #ifndef SiPixelDigitizerAlgorithm_h
2 #define SiPixelDigitizerAlgorithm_h
3 
4 #include <map>
5 #include <memory>
6 #include <vector>
7 #include <iostream>
20 #include "boost/multi_array.hpp"
21 
22 typedef boost::multi_array<float, 2> array_2d;
23 
24 // forward declarations
25 
26 
27 // For the random numbers
28 namespace CLHEP {
29  class HepRandomEngine;
30 }
31 
32 namespace edm {
33  class EventSetup;
34  class ParameterSet;
35 }
36 
37 class DetId;
39 class PixelDigi;
40 class PixelDigiSimLink;
41 class PixelGeomDetUnit;
46 class SiPixelQuality;
48 class TrackerGeometry;
49 class TrackerTopology;
50 
52  public:
55 
56  // initialization that cannot be done in the constructor
57  void init(const edm::EventSetup& es);
58 
59  void initializeEvent() {
60  _signal.clear();
61  }
62 
63  //run the algorithm to digitize a single det
64  void accumulateSimHits(const std::vector<PSimHit>::const_iterator inputBegin,
65  const std::vector<PSimHit>::const_iterator inputEnd,
66  const size_t inputBeginGlobalIndex,
67  const unsigned int tofBin,
68  const PixelGeomDetUnit *pixdet,
69  const GlobalVector& bfield,
70  const TrackerTopology *tTopo,
71  CLHEP::HepRandomEngine*);
72  void digitize(const PixelGeomDetUnit *pixdet,
73  std::vector<PixelDigi>& digis,
74  std::vector<PixelDigiSimLink>& simlinks,
75  const TrackerTopology *tTopo,
76  CLHEP::HepRandomEngine*);
77  void calculateInstlumiFactor(PileupMixingContent* puInfo);
78  void init_DynIneffDB(const edm::EventSetup&, const unsigned int&);
79 
80  // for premixing
81  void calculateInstlumiFactor(const std::vector<PileupSummaryInfo> &ps, int bunchSpacing); // TODO: try to remove the duplication of logic...
82  void setSimAccumulator(const std::map<uint32_t, std::map<int, int> >& signalMap);
83 
84  private:
85 
86  //Accessing Lorentz angle from DB:
88 
89  //Accessing Dead pixel modules from DB:
91 
92  //Accessing Map and Geom:
95 
96  // Get Dynamic Inefficiency scale factors from DB
98 
99  // Define internal classes
100 
101  // definition class
102  //
103  class Amplitude {
104  public:
105  Amplitude() : _amp(0.0) {}
106  Amplitude( float amp, float frac) :
107  _amp(amp), _frac(1, frac) {
108  //in case of digi from noisypixels
109  //the MC information are removed
110  if (_frac[0]<-0.5) {
111  _frac.pop_back();
112  }
113  }
114 
115  Amplitude( float amp, const PSimHit* hitp, size_t hitIndex, unsigned int tofBin, float frac) :
116  _amp(amp), _frac(1, frac) {
117 
118  //in case of digi from noisypixels
119  //the MC information are removed
120  if (_frac[0]<-0.5) {
121  _frac.pop_back();
122  }
123  else {
124  _hitInfos.emplace_back(hitp, hitIndex, tofBin);
125  }
126  }
127 
128  // can be used as a float by convers.
129  operator float() const { return _amp;}
130  float ampl() const {return _amp;}
131  const std::vector<float>& individualampl() const {return _frac;}
132  const std::vector<SimHitInfoForLinks>& hitInfos() const { return _hitInfos; }
133 
134  void operator+=( const Amplitude& other) {
135  _amp += other._amp;
136  //in case of contribution of noise to the digi
137  //the MC information are removed
138  if (other._frac[0]>-0.5){
139  if(!other._hitInfos.empty()) {
140  _hitInfos.insert(_hitInfos.end(), other._hitInfos.begin(), other._hitInfos.end());
141  }
142  _frac.insert(_frac.end(), other._frac.begin(), other._frac.end());
143  }
144  }
145  void operator+=( const float& amp) {
146  _amp += amp;
147  }
148 
149  void set (const float amplitude) { // Used to reset the amplitude
150  _amp = amplitude;
151  }
152 /* void setind (const float indamplitude) { // Used to reset the amplitude */
153 /* _frac = idamplitude; */
154 /* } */
155  private:
156  float _amp;
157  std::vector<float> _frac;
158  std::vector<SimHitInfoForLinks> _hitInfos;
159  }; // end class Amplitude
160 
161  // Define a class to hold the calibration parameters per pixel
162  // Internal
164  public:
165  float p0;
166  float p1;
167  float p2;
168  float p3;
169  };
170  //
171  // Define a class for 3D ionization points and energy
172  //
177  public:
178  EnergyDepositUnit(): _energy(0),_position(0,0,0){}
179  EnergyDepositUnit(float energy,float x, float y, float z):
180  _energy(energy),_position(x,y,z){}
182  _energy(energy),_position(position){}
183  float x() const{return _position.x();}
184  float y() const{return _position.y();}
185  float z() const{return _position.z();}
186  float energy() const { return _energy;}
187  private:
188  float _energy;
190  };
191 
192  //
193  // define class to store signals on the collection surface
194  //
199  class SignalPoint {
200  public:
201  SignalPoint() : _pos(0,0), _time(0), _amplitude(0),
202  _sigma_x(1.), _sigma_y(1.), _hitp(nullptr) {}
203 
204  SignalPoint( float x, float y, float sigma_x, float sigma_y,
205  float t, float a=1.0) :
206  _pos(x,y), _time(t), _amplitude(a), _sigma_x(sigma_x),
207  _sigma_y(sigma_y), _hitp(nullptr) {}
208 
209  SignalPoint( float x, float y, float sigma_x, float sigma_y,
210  float t, const PSimHit& hit, float a=1.0) :
211  _pos(x,y), _time(t), _amplitude(a), _sigma_x(sigma_x),
212  _sigma_y(sigma_y),_hitp(&hit) {}
213 
214  const LocalPoint& position() const { return _pos;}
215  float x() const { return _pos.x();}
216  float y() const { return _pos.y();}
217  float sigma_x() const { return _sigma_x;}
218  float sigma_y() const { return _sigma_y;}
219  float time() const { return _time;}
220  float amplitude() const { return _amplitude;}
221  const PSimHit& hit() { return *_hitp;}
222  SignalPoint& set_amplitude( float amp) { _amplitude = amp; return *this;}
223 
224 
225 
226  private:
228  float _time;
229  float _amplitude;
230  float _sigma_x; // gaussian sigma in the x direction (cm)
231  float _sigma_y; // " " y direction (cm) */
232  const PSimHit* _hitp;
233  };
234 
235  //
236  // PixelEfficiencies struct
237  //
242  PixelEfficiencies(const edm::ParameterSet& conf, bool AddPixelInefficiency, int NumberOfBarrelLayers, int NumberOfEndcapDisks);
243  bool FromConfig; // If true read from Config, otherwise use Database
244 
246  std::vector<double> pu_scale; // in config: 0-3 BPix, 4-5 FPix (inner, outer)
247  std::vector<std::vector<double> > thePUEfficiency; // Instlumi dependent efficiency
248 
249  // Read factors from Configuration
250  double thePixelEfficiency[20]; // Single pixel effciency
251  double thePixelColEfficiency[20]; // Column effciency
252  double thePixelChipEfficiency[20]; // ROC efficiency
253  std::vector<double> theLadderEfficiency_BPix[20]; // Ladder efficiency
254  std::vector<double> theModuleEfficiency_BPix[20]; // Module efficiency
255  double theInnerEfficiency_FPix[20]; // Fpix inner module efficiency
256  double theOuterEfficiency_FPix[20]; // Fpix outer module efficiency
257  unsigned int FPixIndex; // The Efficiency index for FPix Disks
258 
259  // Read factors from DB and fill containers
260  std::map<uint32_t, double> PixelGeomFactors;
261  std::map<uint32_t, std::vector<double> > PixelGeomFactorsROCStdPixels;
262  std::map<uint32_t, std::vector<double> > PixelGeomFactorsROCBigPixels;
263  std::map<uint32_t, double> ColGeomFactors;
264  std::map<uint32_t, double> ChipGeomFactors;
265  std::map<uint32_t, size_t > iPU;
266 
267  // constants for ROC level simulation for Phase1
268  enum shiftEnumerator {FPixRocIdShift = 3, BPixRocIdShift = 6};
269  static const int rocIdMaskBits = 0x1F;
271  bool matches(const DetId&, const DetId&, const std::vector<uint32_t >&);
272  };
273 
274  //
275  // PixelAging struct
276  //
280  struct PixelAging {
281  PixelAging(const edm::ParameterSet& conf, bool AddPixelAging, int NumberOfBarrelLayers, int NumberOfEndcapDisks);
282  float thePixelPseudoRadDamage[20]; // PseudoRadiation Damage Values for aging studies
283  unsigned int FPixIndex; // The Efficiency index for FPix Disks
284  };
285 
286  private:
287  // Internal typedefs
288  typedef std::map<int, Amplitude, std::less<int> > signal_map_type; // from Digi.Skel.
289  typedef signal_map_type::iterator signal_map_iterator; // from Digi.Skel.
290  typedef signal_map_type::const_iterator signal_map_const_iterator; // from Digi.Skel.
291  typedef std::map<uint32_t, signal_map_type> signalMaps;
293  typedef std::vector<edm::ParameterSet> Parameters;
294  typedef boost::multi_array<float, 2> array_2d;
295 
296  // Contains the accumulated hit info.
297  signalMaps _signal;
298 
299  const bool makeDigiSimLinks_;
300 
301  const bool use_ineff_from_db_;
302  const bool use_module_killing_; // remove or not the dead pixel modules
303  const bool use_deadmodule_DB_; // if we want to get dead pixel modules from the DataBase.
304  const bool use_LorentzAngle_DB_; // if we want to get Lorentz angle from the DataBase.
305 
306  const Parameters DeadModules;
307 
308  // Variables and objects for the charge reweighting using 2D templates
310  std::vector<bool> xdouble;
311  std::vector<bool> ydouble;
312  std::vector<float> track;
313  int IDnum, IDden;
314 
315  std::vector<SiPixelTemplateStore2D> templateStores_;
316 
319 
320  private:
321 
322 
323  // Variables
324  //external parameters
325  // go from Geant energy GeV to number of electrons
326  const float GeVperElectron; // 3.7E-09
327 
328  //-- drift
329  const float Sigma0; //=0.0007 // Charge diffusion in microns for 300 micron Si
330  const float Dist300; //=0.0300 // Define 300microns for normalization
331  const bool alpha2Order; // Switch on/off of E.B effect
332 
333 
334 
335  //-- induce_signal
336  const float ClusterWidth; // Gaussian charge cutoff width in sigma units
337  //-- Allow for upgrades
338  const int NumberOfBarrelLayers; // Default = 3
339  const int NumberOfEndcapDisks; // Default = 2
340 
341  //-- make_digis
342  const float theElectronPerADC; // Gain, number of electrons per adc count.
343  const int theAdcFullScale; // Saturation count, 255=8bit.
344  const float theNoiseInElectrons; // Noise (RMS) in units of electrons.
345  const float theReadoutNoise; // Noise of the readount chain in elec,
346  //inludes DCOL-Amp,TBM-Amp, Alt, AOH,OptRec.
347 
348  const float theThresholdInE_FPix; // Pixel threshold in electrons FPix.
349  const float theThresholdInE_BPix; // Pixel threshold in electrons BPix.
350  const float theThresholdInE_BPix_L1; // In case the BPix layer1 gets a different threshold
351  const float theThresholdInE_BPix_L2; // In case the BPix layer2 gets a different threshold
352 
357 
358  const float electronsPerVCAL; // for electrons - VCAL conversion
359  const float electronsPerVCAL_Offset; // in misscalibrate()
360  const float electronsPerVCAL_L1; // same for Layer 1
361  const float electronsPerVCAL_L1_Offset;// same for Layer 1
362 
363  const float theTofLowerCut; // Cut on the particle TOF
364  const float theTofUpperCut; // Cut on the particle TOF
365  const float tanLorentzAnglePerTesla_FPix; //FPix Lorentz angle tangent per Tesla
366  const float tanLorentzAnglePerTesla_BPix; //BPix Lorentz angle tangent per Tesla
367 
368  const float FPix_p0;
369  const float FPix_p1;
370  const float FPix_p2;
371  const float FPix_p3;
372  const float BPix_p0;
373  const float BPix_p1;
374  const float BPix_p2;
375  const float BPix_p3;
376 
377 
378  //-- add_noise
379  const bool addNoise;
381  const bool addNoisyPixels;
382  const bool fluctuateCharge;
383  //-- pixel efficiency
384  const bool AddPixelInefficiency; // bool to read in inefficiencies
385 
387 
388  //-- calibration smearing
389  const bool doMissCalibrate; // Switch on the calibration smearing
390  const float theGainSmearing; // The sigma of the gain fluctuation (around 1)
391  const float theOffsetSmearing; // The sigma of the offset fluct. (around 0)
392 
393  // pixel aging
394  const bool AddPixelAging;
395  const bool UseReweighting;
396  const bool PrintClusters;
397  const bool PrintTemplates;
398 
399  // The PDTable
400  //HepPDTable *particleTable;
401  //ParticleDataTable *particleTable;
402 
403  //-- charge fluctuation
404  const double tMax; // The delta production cut, should be as in OSCAR = 30keV
405  // cmsim = 100keV
406 
407  // The eloss fluctuation class from G4. Is the right place?
408  const std::unique_ptr<SiG4UniversalFluctuation> fluctuate; // make a pointer
409  const std::unique_ptr<GaussianTailNoiseGenerator> theNoiser;
410 
411  // To store calibration constants
412  const std::map<int,CalParameters,std::less<int> > calmap;
413 
414 
415  //-- additional member functions
416  // Private methods
417  std::map<int,CalParameters,std::less<int> > initCal() const;
418  void primary_ionization( const PSimHit& hit, std::vector<EnergyDepositUnit>& ionization_points, CLHEP::HepRandomEngine*) const;
419  void drift(const PSimHit& hit,
420  const PixelGeomDetUnit *pixdet,
421  const GlobalVector& bfield,
422  const TrackerTopology *tTopo,
423  const std::vector<EnergyDepositUnit>& ionization_points,
424  std::vector<SignalPoint>& collection_points) const;
425  void induce_signal(std::vector<PSimHit>::const_iterator inputBegin,
426  std::vector<PSimHit>::const_iterator inputEnd,
427  const PSimHit& hit,
428  const size_t hitIndex,
429  const unsigned int tofBin,
430  const PixelGeomDetUnit *pixdet,
431  const std::vector<SignalPoint>& collection_points);
432  void fluctuateEloss(int particleId, float momentum, float eloss,
433  float length, int NumberOfSegments,
434  float elossVector[],
435  CLHEP::HepRandomEngine*) const;
436  void add_noise(const PixelGeomDetUnit *pixdet,
437  float thePixelThreshold,
438  CLHEP::HepRandomEngine*);
439  void make_digis(float thePixelThresholdInE,
440  uint32_t detID,
441  const PixelGeomDetUnit* pixdet,
442  std::vector<PixelDigi>& digis,
443  std::vector<PixelDigiSimLink>& simlinks,
444  const TrackerTopology *tTopo) const;
445  void pixel_inefficiency(const PixelEfficiencies& eff,
446  const PixelGeomDetUnit* pixdet,
447  const TrackerTopology *tTopo,
448  CLHEP::HepRandomEngine*);
449 
450  void pixel_inefficiency_db(uint32_t detID);
451 
452  float pixel_aging(const PixelAging& aging,
453  const PixelGeomDetUnit* pixdet,
454  const TrackerTopology *tTopo) const;
455 
456  // access to the gain calibration payloads in the db. Only gets initialized if check_dead_pixels_ is set to true.
457  const std::unique_ptr<SiPixelGainCalibrationOfflineSimService> theSiPixelGainCalibrationService_;
458  float missCalibrate(uint32_t detID, const TrackerTopology *tTopo, const PixelGeomDetUnit* pixdet, int col, int row, float amp) const;
459  LocalVector DriftDirection(const PixelGeomDetUnit* pixdet,
460  const GlobalVector& bfield,
461  const DetId& detId) const;
462 
463  void module_killing_conf(uint32_t detID); // remove dead modules using the list in the configuration file PixelDigi_cfi.py
464  void module_killing_DB(uint32_t detID); // remove dead modules uisng the list in the DB
465 
466  // methods for charge reweighting in irradiated sensors
467  int PixelTempRewgt2D( int id_gen, int id_rewgt,
468  array_2d& cluster);
469  bool hitSignalReweight(const PSimHit& hit,
470  std::map< int, float, std::less<int> >& hit_signal,
471  const size_t hitIndex,
472  const unsigned int tofBin,
473  const PixelTopology* topol,
474  uint32_t detID,
475  signal_map_type& theSignal,
476  unsigned short int processType);
477  void printCluster(array_2d& cluster);
478  void printCluster(float arr[BXM2][BYM2]);
479  void printCluster(float arr[TXSIZE][TYSIZE]);
480 
481 
484 
485  double calcQ(float x) const {
486  // need erf(x/sqrt2)
487  //float x2=0.5*x*x;
488  //float a=0.147;
489  //double erf=sqrt(1.0f-exp( -1.0f*x2*( (4/M_PI)+a*x2)/(1.0+a*x2)));
490  //if (x<0.) erf*=-1.0;
491  //return 0.5*(1.0-erf);
492 
493  auto xx=std::min(0.5f*x*x,12.5f);
494  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));
495  }
496 
497 
498 };
499 
500 #endif
const SiPixel2DTemplateDBObject * dbobject_den
Amplitude(float amp, const PSimHit *hitp, size_t hitIndex, unsigned int tofBin, float frac)
signal_map_type::const_iterator signal_map_const_iterator
edm::ESHandle< SiPixelFedCablingMap > map_
SignalPoint(float x, float y, float sigma_x, float sigma_y, float t, const PSimHit &hit, float a=1.0)
const std::unique_ptr< SiPixelGainCalibrationOfflineSimService > theSiPixelGainCalibrationService_
const std::vector< float > & individualampl() const
EnergyDepositUnit(float energy, Local3DPoint position)
#define TXSIZE
const std::unique_ptr< SiG4UniversalFluctuation > fluctuate
SignalPoint(float x, float y, float sigma_x, float sigma_y, float t, float a=1.0)
int init
Definition: HydjetWrapper.h:67
LocalVector drift(const StripGeomDetUnit *, const MagneticField &, const SiStripLorentzAngle &)
Definition: ShallowTools.cc:38
std::vector< std::vector< double > > thePUEfficiency
#define nullptr
edm::ESHandle< TrackerGeometry > geom_
const SiPixel2DTemplateDBObject * dbobject_num
#define BXM2
const std::map< int, CalParameters, std::less< int > > calmap
boost::multi_array< float, 2 > array_2d
std::map< uint32_t, std::vector< double > > PixelGeomFactorsROCBigPixels
T sqrt(T t)
Definition: SSEVec.h:18
#define BYM2
double f[11][100]
T min(T a, T b)
Definition: MathUtil.h:58
signal_map_type::iterator signal_map_iterator
EnergyDepositUnit(float energy, float x, float y, float z)
edm::ESHandle< SiPixelLorentzAngle > SiPixelLorentzAngle_
boost::multi_array< float, 2 > array_2d
const std::unique_ptr< GaussianTailNoiseGenerator > theNoiser
Definition: DetId.h:18
std::map< int, Amplitude, std::less< int > > signal_map_type
#define TYSIZE
edm::ESHandle< SiPixelQuality > SiPixelBadModule_
std::map< uint32_t, std::vector< double > > PixelGeomFactorsROCStdPixels
GloballyPositioned< double > Frame
std::vector< SiPixelTemplateStore2D > templateStores_
edm::ESHandle< SiPixelDynamicInefficiency > SiPixelDynamicInefficiency_
std::vector< edm::ParameterSet > Parameters
const std::vector< SimHitInfoForLinks > & hitInfos() const
HLT enums.
std::vector< SimHitInfoForLinks > _hitInfos
double a
Definition: hdecay.h:121
static int position[264][3]
Definition: ReadPGInfo.cc:509
col
Definition: cuy.py:1010
std::map< uint32_t, signal_map_type > signalMaps
Definition: aging.py:1