CMS 3D CMS Logo
HcalQie.cc
Go to the documentation of this file.
1 // File: HcalQie.cc
3 // Description: Simulation of QIE readout for HCal
5 
6 #include "SimG4CMS/Calo/interface/HcalQie.h"
7 
8 #include "CLHEP/Random/RandGaussQ.h"
9 #include "CLHEP/Random/RandPoissonQ.h"
10 #include "CLHEP/Units/GlobalPhysicalConstants.h"
11 
12 #include <iostream>
13 #include <iomanip>
14 #include <cmath>
15 
16 //#define EDM_ML_DEBUG
17 
18 HcalQie::HcalQie(edm::ParameterSet const& p) {
19  //static SimpleConfigurable<double> p1(4.0, "HcalQie:qToPE");
20  //static SimpleConfigurable<int> p2(6, "HcalQie:BinOfMax");
21  //static SimpleConfigurable<int> p3(2, "HcalQie:SignalBuckets");
22  //static SimpleConfigurable<int> p4(0, "HcalQie:PreSamples");
23  //static SimpleConfigurable<int> p5(10, "HcalQie:NumOfBuckets");
24  //static SimpleConfigurable<double> p6(0.5, "HcalQie:SigmaNoise");
25  //static SimpleConfigurable<double> p7(0.0005,"HcalQie:EDepPerPE");
26  //static SimpleConfigurable<int> p8(4, "HcalQie:BaseLine");
27 
28  edm::ParameterSet m_HQ = p.getParameter<edm::ParameterSet>("HcalQie");
29  qToPE = m_HQ.getParameter<double>("qToPE");
30  binOfMax = m_HQ.getParameter<int>("BinOfMax");
31  signalBuckets = m_HQ.getParameter<int>("SignalBuckets");
32  preSamples = m_HQ.getParameter<int>("PreSamples");
33  numOfBuckets = m_HQ.getParameter<int>("NumOfBuckets");
34  sigma = m_HQ.getParameter<double>("SigmaNoise");
35  eDepPerPE = m_HQ.getParameter<double>("EDepPerPE");
36  int bl = m_HQ.getParameter<int>("BaseLine");
37 
38  shape_ = shape();
39  code_ = code();
40  charge_ = charge();
41  if (signalBuckets == 1) {
42  phase_ = -3;
43  rescale_ = 1.46;
44  } else if (signalBuckets == 3) {
45  phase_ = -1;
46  rescale_ = 1.06;
47  } else if (signalBuckets == 4) {
48  phase_ = 0;
49  rescale_ = 1.03;
50  } else {
51  phase_ = -2;
52  rescale_ = 1.14;
53  signalBuckets = 2;
54  }
55  weight_ = weight(binOfMax, signalBuckets, preSamples, numOfBuckets);
56  baseline = codeToQ(bl);
57  bmin_ = binOfMax - 3;
58  if (bmin_ < 0)
59  bmin_ = 0;
60  if (binOfMax > numOfBuckets)
61  bmax_ = numOfBuckets + 5;
62  else
63  bmax_ = binOfMax + 5;
64 
65  edm::LogVerbatim("HcalSim") << "HcalQie: initialized with binOfMax " << binOfMax << " sample from " << bmin_ << " to "
66  << bmax_ << "; signalBuckets " << signalBuckets << " Baseline/Phase/Scale " << baseline
67  << "/" << phase_ << "/" << rescale_ << "\n Noise " << sigma
68  << "fC fCToPE " << qToPE << " EDepPerPE " << eDepPerPE;
69 }
70 
71 HcalQie::~HcalQie() { edm::LogVerbatim("HcalSim") << "HcalQie:: Deleting Qie"; }
72 
73 std::vector<double> HcalQie::shape() {
74  // pulse shape time constants in ns
75  const float ts1 = 8.; // scintillation time constants : 1,2,3
76  const float ts2 = 10.;
77  const float ts3 = 29.3;
78  const float thpd = 4.; // HPD current collection drift time
79  const float tpre = 5.; // preamp time constant
80 
81  const float wd1 = 2.; // relative weights of decay exponents
82  const float wd2 = 0.7;
83  const float wd3 = 1.;
84 
85  // HPD starts at I and rises to 2I in thpd of time
86  double norm = 0.0;
87  int j, hpd_siz = (int)(thpd);
88  std::vector<double> hpd_drift(hpd_siz);
89  for (j = 0; j < hpd_siz; j++) {
90  double tmp = (double)j + 0.5;
91  hpd_drift[j] = 1.0 + tmp / thpd;
92  norm += hpd_drift[j];
93  }
94  // normalize integrated current to 1.0
95  for (j = 0; j < hpd_siz; j++) {
96  hpd_drift[j] /= norm;
97  }
98 
99  // Binkley shape over 6 time constants
100  int preamp_siz = (int)(6 * tpre);
101  std::vector<double> preamp(preamp_siz);
102  norm = 0;
103  for (j = 0; j < preamp_siz; j++) {
104  double tmp = (double)j + 0.5;
105  preamp[j] = tmp * exp(-(tmp * tmp) / (tpre * tpre));
106  norm += preamp[j];
107  }
108  // normalize pulse area to 1.0
109  for (j = 0; j < preamp_siz; j++) {
110  preamp[j] /= norm;
111  }
112 
113  // ignore stochastic variation of photoelectron emission
114  // <...>
115  // effective tile plus wave-length shifter decay time over 4 time constants
116 
117  int tmax = 6 * (int)ts3;
118  std::vector<double> scnt_decay(tmax);
119  norm = 0;
120  for (j = 0; j < tmax; j++) {
121  double tmp = (double)j + 0.5;
122  scnt_decay[j] = wd1 * exp(-tmp / ts1) + wd2 * exp(-tmp / ts2) + wd3 * exp(-tmp / ts3);
123  norm += scnt_decay[j];
124  }
125  // normalize pulse area to 1.0
126  for (j = 0; j < tmax; j++) {
127  scnt_decay[j] /= norm;
128  }
129 
130  int nsiz = tmax + hpd_siz + preamp_siz + 1;
131  std::vector<double> pulse(nsiz, 0.0); // zeroing output pulse shape
132  norm = 0;
133  int i, k;
134  for (i = 0; i < tmax; i++) {
135  int t1 = i; // and ignore jitter from optical path length
136  for (j = 0; j < hpd_siz; j++) {
137  int t2 = t1 + j;
138  for (k = 0; k < preamp_siz; k++) {
139  int t3 = t2 + k;
140  float tmp = scnt_decay[i] * hpd_drift[j] * preamp[k];
141  pulse[t3] += tmp;
142  norm += tmp;
143  }
144  }
145  }
146 
147  // normalize for 1 GeV pulse height
148  edm::LogVerbatim("HcalSim") << "HcalQie: Convoluted Shape ============== Normalisation " << norm;
149  for (i = 0; i < nsiz; i++) {
150  pulse[i] /= norm;
151 #ifdef EDM_ML_DEBUG
152  edm::LogVerbatim("HcalSim") << "HcalQie: Pulse[" << std::setw(3) << i << "] " << std::setw(8) << pulse[i];
153 #endif
154  }
155 
156  return pulse;
157 }
158 
159 std::vector<int> HcalQie::code() {
160  unsigned int CodeFADCdata[122] = {
161  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
162  21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43,
163  44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 66,
164  67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87,
165  88, 89, 90, 91, 92, 93, 94, 95, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110,
166  111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127};
167 
168  std::vector<int> temp(122);
169  int i;
170  for (i = 0; i < 122; i++)
171  temp[i] = (int)CodeFADCdata[i];
172 
173 #ifdef EDM_ML_DEBUG
174  int siz = temp.size();
175  edm::LogVerbatim("HcalSim") << "HcalQie: Codes in array of size " << siz;
176  for (i = 0; i < siz; i++)
177  edm::LogVerbatim("HcalSim") << "HcalQie: Code[" << std::setw(3) << i << "] " << std::setw(6) << temp[i];
178 #endif
179  return temp;
180 }
181 
182 std::vector<double> HcalQie::charge() {
183  double ChargeFADCdata[122] = {
184  -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5,
185  12.5, 13.5, 14.5, 16.5, 18.5, 20.5, 22.5, 24.5, 26.5, 28.5, 31.5, 34.5, 37.5, 40.5,
186  44.5, 48.5, 52.5, 57.5, 62.5, 67.5, 72.5, 77.5, 82.5, 87.5, 92.5, 97.5, 102.5, 107.5,
187  112.5, 117.5, 122.5, 127.5, 132.5, 142.5, 152.5, 162.5, 172.5, 182.5, 192.5, 202.5, 217.5, 232.5,
188  247.5, 262.5, 282.5, 302.5, 322.5, 347.5, 372.5, 397.5, 422.5, 447.5, 472.5, 497.5, 522.5, 547.5,
189  572.5, 597.5, 622.5, 647.5, 672.5, 697.5, 722.5, 772.5, 822.5, 872.5, 922.5, 972.5, 1022.5, 1072.5,
190  1147.5, 1222.5, 1297.5, 1372.5, 1472.5, 1572.5, 1672.5, 1797.5, 1922.5, 2047.5, 2172.5, 2297.5, 2422.5, 2547.5,
191  2672.5, 2797.5, 2922.5, 3047.5, 3172.5, 3397.5, 3422.5, 3547.5, 3672.5, 3922.5, 4172.5, 4422.5, 4672.5, 4922.5,
192  5172.5, 5422.5, 5797.5, 6172.5, 6547.5, 6922.5, 7422.5, 7922.5, 8422.5, 9047.5};
193 
194  std::vector<double> temp(122);
195  int i;
196  for (i = 0; i < 122; i++)
197  temp[i] = (double)(ChargeFADCdata[i]);
198 
199 #ifdef EDM_ML_DEBUG
200  int siz = temp.size();
201  edm::LogVerbatim("HcalSim") << "HcalQie: Charges in array of size " << siz;
202  for (i = 0; i < siz; i++)
203  edm::LogVerbatim("HcalSim") << "HcalQie: Charge[" << std::setw(3) << i << "] " << std::setw(8) << temp[i];
204 #endif
205  return temp;
206 }
207 
208 std::vector<double> HcalQie::weight(int binOfMax, int mode, int npre, int bucket) {
209  std::vector<double> temp(bucket, 0);
210  int i;
211  for (i = binOfMax - 1; i < binOfMax + mode - 1; i++)
212  temp[i] = 1.;
213  if (npre > 0) {
214  for (i = 0; i < npre; i++) {
215  int j = binOfMax - 2 - i;
216  temp[j] = -(double)mode / (double)npre;
217  }
218  }
219 
220 #ifdef EDM_ML_DEBUG
221  int siz = temp.size();
222  edm::LogVerbatim("HcalSim") << "HcalQie: Weights in array of size " << siz << " and Npre " << npre;
223  for (i = 0; i < siz; i++)
224  edm::LogVerbatim("HcalSim") << "HcalQie: [Weight[" << i << "] = " << temp[i];
225 #endif
226  return temp;
227 }
228 
229 double HcalQie::codeToQ(int ic) {
230  double tmp = 0;
231  for (unsigned int i = 0; i < code_.size(); i++) {
232  if (ic == code_[i]) {
233  double delta;
234  if (i == code_.size() - 1)
235  delta = charge_[i] - charge_[i - 1];
236  else
237  delta = charge_[i + 1] - charge_[i];
238  tmp = charge_[i] + 0.5 * delta;
239  break;
240  }
241  }
242 
243  return tmp;
244 }
245 
246 int HcalQie::getCode(double charge) {
247  int tmp = 0;
248  for (unsigned int i = 0; i < charge_.size(); i++) {
249  if (charge < charge_[i]) {
250  if (i > 0)
251  tmp = code_[i - 1];
252  break;
253  }
254  }
255 
256  return tmp;
257 }
258 
259 double HcalQie::getShape(double time) {
260  double tmp = 0;
261  int k = (int)(time + 0.5);
262  if (k >= 0 && k < ((int)(shape_.size()) - 1))
263  tmp = 0.5 * (shape_[k] + shape_[k + 1]);
264 
265  return tmp;
266 }
267 
268 std::vector<int> HcalQie::getCode(int nht, const std::vector<CaloHit>& hitbuf, CLHEP::HepRandomEngine* engine) {
269  const double bunchSpace = 25.;
270  int nmax = (bmax_ > numOfBuckets ? bmax_ : numOfBuckets);
271  std::vector<double> work(nmax);
272 
273  // Noise in the channel
274  for (int i = 0; i < numOfBuckets; i++)
275  work[i] = CLHEP::RandGaussQ::shoot(engine, baseline, sigma);
276 
277 #ifdef EDM_ML_DEBUG
278  edm::LogVerbatim("HcalSim") << "HcalQie::getCode: Noise with baseline " << baseline << " width " << sigma << " and "
279  << nht << " hits";
280  for (int i = 0; i < numOfBuckets; i++)
281  edm::LogVerbatim("HcalSim") << "HcalQie: Code[" << i << "] = " << work[i];
282 #endif
283  double etot = 0, esum = 0, photons = 0;
284  if (nht > 0) {
285  // Sort the hits
286  std::vector<const CaloHit*> hits(nht);
287  std::vector<const CaloHit*>::iterator k1, k2;
288  int kk;
289  for (kk = 0; kk < nht; kk++) {
290  hits[kk] = &hitbuf[kk];
291  }
292  sort(hits.begin(), hits.end(), CaloHitMore());
293 
294  // Energy deposits
295  for (kk = 0, k1 = hits.begin(); k1 != hits.end(); kk++, k1++) {
296  double ehit = (**k1).e();
297  double jitter = (**k1).t();
298  int jump = 0;
299  for (k2 = k1 + 1; k2 != hits.end() && (jitter - (**k2).t()) < 1. && (jitter - (**k2).t()) > -1.; k2++) {
300  ehit += (**k2).e();
301  jump++;
302  }
303 
304  double avpe = ehit / eDepPerPE;
305  CLHEP::RandPoissonQ randPoissonQ(*engine, avpe);
306  double photo = randPoissonQ.fire();
307  etot += ehit;
308  photons += photo;
309 #ifdef EDM_ML_DEBUG
310  edm::LogVerbatim("HcalSim") << "HcalQie::getCode: Hit " << kk << ":" << kk + jump << " Energy deposit " << ehit
311  << " Time " << jitter << " Average and true no of PE " << avpe << " " << photo;
312 #endif
313  double bintime = jitter - phase_ - bunchSpace * (binOfMax - bmin_);
314 #ifdef EDM_ML_DEBUG
315  edm::LogVerbatim("HcalSim") << "HcalQie::getCode: phase " << phase_ << " binTime " << bintime;
316 #endif
317  std::vector<double> binsum(nmax, 0);
318  double norm = 0, sum = 0.;
319  for (int i = bmin_; i < bmax_; i++) {
320  bintime += bunchSpace;
321  for (int j = 0; j < (int)(bunchSpace); j++) {
322  double tim = bintime + j;
323  double tmp = getShape(tim);
324  binsum[i] += tmp;
325  }
326  sum += binsum[i];
327  }
328 
329  if (sum > 0)
330  norm = (photo / (sum * qToPE));
331 #ifdef EDM_ML_DEBUG
332  edm::LogVerbatim("HcalSim") << "HcalQie::getCode: PE " << photo << " Sum " << sum << " Norm. " << norm;
333 #endif
334  for (int i = bmin_; i < bmax_; i++)
335  work[i] += binsum[i] * norm;
336 
337  kk += jump;
338  k1 += jump;
339  }
340  }
341 
342  std::vector<int> temp(numOfBuckets, 0);
343  for (int i = 0; i < numOfBuckets; i++) {
344  temp[i] = getCode(work[i]);
345  esum += work[i];
346  }
347 #ifdef EDM_ML_DEBUG
348  edm::LogVerbatim("HcalSim") << "HcalQie::getCode: Input " << etot << " GeV; Photons " << photons << "; Output "
349  << esum << " fc";
350 #endif
351  return temp;
352 }
353 
354 double HcalQie::getEnergy(const std::vector<int>& code) {
355  std::vector<double> work(numOfBuckets);
356  double sum = 0;
357  for (int i = 0; i < numOfBuckets; i++) {
358  work[i] = codeToQ(code[i]);
359  sum += work[i] * weight_[i];
360 #ifdef EDM_ML_DEBUG
361  edm::LogVerbatim("HcalSim") << "HcalQie::getEnergy: " << i << " code " << code[i] << " PE " << work[i];
362 #endif
363  }
364 
365  double tmp;
366  if (preSamples == 0) {
367  tmp = (sum - signalBuckets * baseline) * rescale_ * qToPE * eDepPerPE;
368  } else {
369  tmp = sum * rescale_ * qToPE;
370  }
371 #ifdef EDM_ML_DEBUG
372  edm::LogVerbatim("HcalSim") << "HcalQie::getEnergy: PE " << sum * qToPE << " Energy " << tmp << " GeV";
373 #endif
374  return tmp;
375 }
RandomServiceHelper.t2
t2
Definition: RandomServiceHelper.py:257
HcalQie::charge_
std::vector< double > charge_
Definition: HcalQie.h:38
mps_fire.i
i
Definition: mps_fire.py:355
HcalQie::codeToQ
double codeToQ(int ic)
Definition: HcalQie.cc:229
hfClusterShapes_cfi.hits
hits
Definition: hfClusterShapes_cfi.py:5
HcalQie::charge
std::vector< double > charge()
Definition: HcalQie.cc:182
AlCaHLTBitMon_ParallelJobs.p
p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
HcalQie.h
ALCARECOPromptCalibProdSiPixelAli0T_cff.mode
mode
Definition: ALCARECOPromptCalibProdSiPixelAli0T_cff.py:96
pulse
double pulse(double x, double y, double z, double t)
Definition: SiStripPulseShape.cc:49
HcalQie::preSamples
int preSamples
Definition: HcalQie.h:39
HcalQie::shape
std::vector< double > shape()
Definition: HcalQie.cc:73
createJobs.tmp
tmp
align.sh
Definition: createJobs.py:716
groupFilesInBlocks.temp
list temp
Definition: groupFilesInBlocks.py:142
tmax
static const double tmax[3]
Definition: CastorTimeSlew.cc:7
edmIntegrityCheck.work
work
Definition: edmIntegrityCheck.py:75
HcalQie::bmin_
int bmin_
Definition: HcalQie.h:42
RandomServiceHelper.t1
t1
Definition: RandomServiceHelper.py:256
HcalQie::qToPE
double qToPE
Definition: HcalQie.h:41
HcalQie::sigma
double sigma
Definition: HcalQie.h:41
HcalQie::numOfBuckets
int numOfBuckets
Definition: HcalQie.h:39
HcalQie::code_
std::vector< int > code_
Definition: HcalQie.h:37
HcalQie::bmax_
int bmax_
Definition: HcalQie.h:42
HcalQie::code
std::vector< int > code()
Definition: HcalQie.cc:159
HcalQie::HcalQie
HcalQie(edm::ParameterSet const &p)
Definition: HcalQie.cc:18
GetRecoTauVFromDQM_MC_cff.kk
kk
Definition: GetRecoTauVFromDQM_MC_cff.py:84
HcalQie::eDepPerPE
double eDepPerPE
Definition: HcalQie.h:41
dqmdumpme.k
k
Definition: dqmdumpme.py:60
RandomServiceHelper.t3
t3
Definition: RandomServiceHelper.py:258
ALCARECOTkAlJpsiMuMu_cff.charge
charge
Definition: ALCARECOTkAlJpsiMuMu_cff.py:47
HcalQie::rescale_
double rescale_
Definition: HcalQie.h:43
edm::ParameterSet
Definition: ParameterSet.h:36
dumpMFGeometry_cfg.delta
delta
Definition: dumpMFGeometry_cfg.py:25
createfilelist.int
int
Definition: createfilelist.py:10
edm::LogVerbatim
Definition: MessageLogger.h:297
HcalQie::getCode
std::vector< int > getCode(int, const std::vector< CaloHit > &, CLHEP::HepRandomEngine *)
Definition: HcalQie.cc:268
BPHMonitor_cfi.photons
photons
Definition: BPHMonitor_cfi.py:91
HcalQie::phase_
double phase_
Definition: HcalQie.h:43
HcalQie::weight_
std::vector< double > weight_
Definition: HcalQie.h:40
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
HcalQie::baseline
double baseline
Definition: HcalQie.h:41
CaloHitMore
Definition: CaloHit.h:36
HcalQie::getEnergy
double getEnergy(const std::vector< int > &)
Definition: HcalQie.cc:354
HcalQie::shape_
std::vector< double > shape_
Definition: HcalQie.h:36
HcalQie::signalBuckets
int signalBuckets
Definition: HcalQie.h:39
HcalQie::getShape
double getShape(double time)
Definition: HcalQie.cc:259
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
ntuplemaker.time
time
Definition: ntuplemaker.py:310
HcalQie::~HcalQie
virtual ~HcalQie()
Definition: HcalQie.cc:71
HcalQie::binOfMax
int binOfMax
Definition: HcalQie.h:39
HcalQie::weight
std::vector< double > weight(int binofmax, int mode, int npre, int numbucket)
Definition: HcalQie.cc:208
Generated for CMSSW Reference Manual by doxygen 1.8.16