CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HPDIonFeedbackSim.cc
Go to the documentation of this file.
1 
2 // --------------------------------------------------------
3 // A class to simulated HPD ion feedback noise.
4 // The deliverable of the class is the ion feedback noise
5 // for an HcalDetId units of fC or GeV
6 //
7 // Project: HPD ion feedback
8 // Author: T.Yetkin University of Iowa, Feb. 16, 2010
9 // --------------------------------------------------------
10 
21 
22 #include "CLHEP/Random/RandBinomial.h"
23 #include "CLHEP/Random/RandGaussQ.h"
24 #include "CLHEP/Random/RandPoissonQ.h"
25 
26 using namespace edm;
27 using namespace std;
28 
29 // constants for simulation/parameterization
30 double pe2Charge = 0.333333; // fC/p.e.
31 
33 : theDbService(0), theShapes(shapes)
34 {
35 }
36 
38 {
39 }
40 
41 double HPDIonFeedbackSim::getIonFeedback(DetId detId, double signal, double pedWidth, bool doThermal, bool isInGeV, CLHEP::HepRandomEngine* engine)
42 {
43 
44  // HcalDetId id = detId;
45 
46  double GeVperfC = 1.;
47  if(isInGeV) GeVperfC = 1./fCtoGeV(detId);
48 
49  double charge = signal / GeVperfC;
50 
51  double noise = 0.; // fC
52  if (charge > 3. * pedWidth) { // 3 sigma away from pedestal mean
53  int npe = int (charge / pe2Charge);
54  if(doThermal) {
55  double electronEmission = 0.08;
56  CLHEP::RandPoissonQ theRandPoissonQ(*engine, electronEmission);
57  npe += theRandPoissonQ.fire();
58  }
59 
60  noise = correctPE(detId, npe, engine) - npe;
61  }
62  return (noise * GeVperfC);
63 
64 }
65 
66 
67 double HPDIonFeedbackSim::correctPE(const DetId & detId, double npe, CLHEP::HepRandomEngine* engine) const
68 {
69  double rateInTail = 0.000211988;//read this from XML file
70  double rateInSecondTail = 4.61579e-06;//read this from XML file
71 
72  // three gauss fit is applied to data to get ion feedback distribution
73  // parameters (in fC)
74  // first gaussian
75  // double p0 = 9.53192e+05;
76  // double p1 = -3.13653e-01;
77  // double p2 = 2.78350e+00;
78 
79  // second gaussian
80  // double p3 = 2.41611e+03;
81  double p4 = 2.06117e+01;
82  double p5 = 1.09239e+01;
83 
84  // third gaussian
85  // double p6 = 3.42793e+01;
86  double p7 = 5.45548e+01;
87  double p8 = 1.59696e+01;
88 
89  double noise = 0.; // fC
90  int nFirst = (int)(CLHEP::RandBinomial::shoot(engine, npe, rateInTail));
91  int nSecond = (int)(CLHEP::RandBinomial::shoot(engine, npe, rateInSecondTail));
92 
93  for (int j = 0; j < nFirst; ++j) {
94  noise += CLHEP::RandGaussQ::shoot(engine, p4, p5);
95  }
96  for (int j = 0; j < nSecond; ++j) {
97  noise += CLHEP::RandGaussQ::shoot(engine, p7, p8);
98  }
99 
100  return npe + std::max(noise/pe2Charge, 0.);
101 }
102 
103 
104 void HPDIonFeedbackSim::addThermalNoise(CaloSamples & samples, CLHEP::HepRandomEngine* engine)
105 {
106  // make some chance to add a PE (with a chance of feedback)
107  // for each time sample
108  double meanPE = 0.02;
109  DetId detId(samples.id());
110  int nSamples = samples.size();
111  const CaloVShape * shape = theShapes->shape(detId);
112  for(int i = 0; i < nSamples; ++i)
113  {
114  CLHEP::RandPoissonQ theRandPoissonQ(*engine, meanPE);
115  double npe = theRandPoissonQ.fire();
116  // TODOprobably should time-smear these
117  if(npe > 0.)
118  {
119  // chance of feedback
120  npe = correctPE(detId, npe, engine);
121  for(int j = i; j < nSamples; ++j)
122  {
123  double timeFromPE = (j-i) * 25.;
124  samples[j] += (*shape)(timeFromPE) * npe;
125  }
126  }
127  }
128 }
129 
130 
131 double HPDIonFeedbackSim::fCtoGeV(const DetId & detId) const
132 {
133  assert(theDbService != 0);
134  HcalGenericDetId hcalGenDetId(detId);
135  const HcalGain* gains = theDbService->getGain(hcalGenDetId);
136  const HcalGainWidth* gwidths = theDbService->getGainWidth(hcalGenDetId);
137  if (!gains || !gwidths )
138  {
139  edm::LogError("HcalAmplifier") << "Could not fetch HCAL conditions for channel " << hcalGenDetId;
140  }
141  // only one gain will be recorded per channel, so just use capID 0 for now
142  double result = gains->getValue(0);
143  return result;
144 }
HPDIonFeedbackSim(const edm::ParameterSet &, const CaloShapes *shapes)
need a shaper in order to set thermal noise
const HcalGainWidth * getGainWidth(const HcalGenericDetId &fId) const
int i
Definition: DBlmapReader.cc:9
double getIonFeedback(DetId detId, double signal, double pedWidth, bool doThermal, bool isInGeV, CLHEP::HepRandomEngine *)
double fCtoGeV(const DetId &detId) const
Electronic response of the preamp.
Definition: CaloVShape.h:11
virtual double correctPE(const DetId &detId, double npe, CLHEP::HepRandomEngine *) const
float getValue(int fCapId) const
get value for capId = 0..3
Definition: HcalGain.h:22
double pe2Charge
const HcalDbService * theDbService
void addThermalNoise(CaloSamples &samples, CLHEP::HepRandomEngine *)
double p4[4]
Definition: TauolaWrapper.h:92
tuple result
Definition: query.py:137
int j
Definition: DBlmapReader.cc:9
Definition: DetId.h:18
int size() const
get the size
Definition: CaloSamples.h:24
const HcalGain * getGain(const HcalGenericDetId &fId) const
const CaloShapes * theShapes
DetId id() const
get the (generic) id
Definition: CaloSamples.h:21
virtual const CaloVShape * shape(const DetId &detId) const
Definition: CaloShapes.h:15