CMS 3D CMS Logo

HcalPulseShapes.cc
Go to the documentation of this file.
6 #include "CLHEP/Random/RandFlat.h"
8 
9 // #include "CalibCalorimetry/HcalAlgos/interface/HcalDbASCIIIO.h"
10 #include <cmath>
11 #include <iostream>
12 #include <fstream>
13 #include "TMath.h"
14 
15 HcalPulseShapes::HcalPulseShapes() : theDbService(nullptr), theShapes() {
16  /*
17 
18 Reco MC
19 -------------------------------------------------------------------------------------------
20 000 not used (reserved)
21 101 101 hpdShape_ HPD (original version)
22 102 102 =101 HPD BV 30 volts in HBP iphi54
23 103 123 hpdShape_v2,hpdShapeMC_v2 HPD (2011. oct version)
24 104 124 hpdBV30Shape_v2,hpdBV30ShapeMC_v2 HPD bv30 in HBP iph54
25 105 125 hpdShape_v2,hpdShapeMC_v2 HPD (2011.11.12 version)
26 201 201 siPMShapeHO_ SiPMs Zecotec shape (HO)
27 202 202 =201, SiPMs Hamamatsu shape (HO)
28 205 203 siPMShapeData2017_,siPMShapeMC2017_ SiPMs from Data, Hamamatsu shape (HE 2017)
29 207 206 siPMShapeData2018_,siPMShapeMC2018_ SiPMs from Data, Hamamatsu shape (HE 2018)
30 207 208 siPMShapeData2018_,siPMShapeMCRecoRun3_ SiPMs from Data, 2021 MC phase scan
31 301 301 hfShape_ regular HF PMT shape
32 401 401 regular ZDC shape
33 -------------------------------------------------------------------------------------------
34 
35 */
36 
37  float ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3;
38 
39  // HPD Shape Version 1 (used before CMSSW5, until Oct 2011)
40  ts1 = 8.;
41  ts2 = 10.;
42  ts3 = 29.3;
43  thpd = 4.0;
44  tpre = 9.0;
45  wd1 = 2.0;
46  wd2 = 0.7;
47  wd3 = 1.0;
48  computeHPDShape(ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3, hpdShape_);
49  theShapes[101] = &hpdShape_;
50  theShapes[102] = theShapes[101];
51 
52  // HPD Shape Version 2 for CMSSW 5. Nov 2011 (RECO and MC separately)
53  ts1 = 8.;
54  ts2 = 10.;
55  ts3 = 25.0;
56  thpd = 4.0;
57  tpre = 9.0;
58  wd1 = 2.0;
59  wd2 = 0.7;
60  wd3 = 1.0;
61  computeHPDShape(ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3, hpdShape_v2);
62  theShapes[103] = &hpdShape_v2;
63 
64  ts1 = 8.;
65  ts2 = 10.;
66  ts3 = 29.3;
67  thpd = 4.0;
68  tpre = 7.0;
69  wd1 = 2.0;
70  wd2 = 0.7;
71  wd3 = 1.0;
72  computeHPDShape(ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3, hpdShapeMC_v2);
73  theShapes[123] = &hpdShapeMC_v2;
74 
75  // HPD Shape Version 3 for CMSSW 5. Nov 2011 (RECO and MC separately)
76  ts1 = 8.;
77  ts2 = 19.;
78  ts3 = 29.3;
79  thpd = 4.0;
80  tpre = 9.0;
81  wd1 = 2.0;
82  wd2 = 0.7;
83  wd3 = 0.32;
84  computeHPDShape(ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3, hpdShape_v3);
85  theShapes[105] = &hpdShape_v3;
86 
87  ts1 = 8.;
88  ts2 = 10.;
89  ts3 = 22.3;
90  thpd = 4.0;
91  tpre = 7.0;
92  wd1 = 2.0;
93  wd2 = 0.7;
94  wd3 = 1.0;
95  computeHPDShape(ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3, hpdShapeMC_v3);
96  theShapes[125] = &hpdShapeMC_v3;
97 
98  // HPD with Bias Voltage 30 volts, wider pulse. (HBPlus iphi54)
99 
100  ts1 = 8.;
101  ts2 = 12.;
102  ts3 = 31.7;
103  thpd = 9.0;
104  tpre = 9.0;
105  wd1 = 2.0;
106  wd2 = 0.7;
107  wd3 = 1.0;
108  computeHPDShape(ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3, hpdBV30Shape_v2);
109  theShapes[104] = &hpdBV30Shape_v2;
110 
111  ts1 = 8.;
112  ts2 = 12.;
113  ts3 = 31.7;
114  thpd = 9.0;
115  tpre = 9.0;
116  wd1 = 2.0;
117  wd2 = 0.7;
118  wd3 = 1.0;
119  computeHPDShape(ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3, hpdBV30ShapeMC_v2);
121 
122  // HF and SiPM
123 
124  computeHFShape();
129 
130  theShapes[201] = &siPMShapeHO_;
131  theShapes[202] = theShapes[201];
132  theShapes[203] = &(computeSiPMShapeHE203());
134  theShapes[206] = &(computeSiPMShapeHE206());
137  theShapes[301] = &hfShape_;
138  //theShapes[401] = new CaloCachedShapeIntegrator(&theZDCShape);
139 }
140 
142 
144  edm::ESHandle<HcalDbService> conditions;
145  es.get<HcalDbRecord>().get(conditions);
146  theDbService = conditions.product();
147 }
148 
149 void HcalPulseShapes::beginRun(const HcalDbService* conditions) { theDbService = conditions; }
150 
151 //void HcalPulseShapes::computeHPDShape()
153  float ts1, float ts2, float ts3, float thpd, float tpre, float wd1, float wd2, float wd3, Shape& tmphpdShape_) {
154  // pulse shape time constants in ns
155  /*
156  const float ts1 = 8.; // scintillation time constants : 1,2,3
157  const float ts2 = 10.;
158  const float ts3 = 29.3;
159  const float thpd = 4.; // HPD current collection drift time
160  const float tpre = 9.; // preamp time constant (refit on TB04 data)
161 
162  const float wd1 = 2.; // relative weights of decay exponents
163  const float wd2 = 0.7;
164  const float wd3 = 1.;
165 */
166  // pulse shape components over a range of time 0 ns to 255 ns in 1 ns steps
167  unsigned int nbin = 256;
168  tmphpdShape_.setNBin(nbin);
169  std::vector<float> ntmp(nbin, 0.0); // zeroing output pulse shape
170  std::vector<float> nth(nbin, 0.0); // zeroing HPD drift shape
171  std::vector<float> ntp(nbin, 0.0); // zeroing Binkley preamp shape
172  std::vector<float> ntd(nbin, 0.0); // zeroing Scintillator decay shape
173 
174  unsigned int i, j, k;
175  float norm;
176 
177  // HPD starts at I and rises to 2I in thpd of time
178  norm = 0.0;
179  for (j = 0; j < thpd && j < nbin; j++) {
180  nth[j] = 1.0 + ((float)j) / thpd;
181  norm += nth[j];
182  }
183  // normalize integrated current to 1.0
184  for (j = 0; j < thpd && j < nbin; j++) {
185  nth[j] /= norm;
186  }
187 
188  // Binkley shape over 6 time constants
189  norm = 0.0;
190  for (j = 0; j < 6 * tpre && j < nbin; j++) {
191  ntp[j] = ((float)j) * exp(-((float)(j * j)) / (tpre * tpre));
192  norm += ntp[j];
193  }
194  // normalize pulse area to 1.0
195  for (j = 0; j < 6 * tpre && j < nbin; j++) {
196  ntp[j] /= norm;
197  }
198 
199  // ignore stochastic variation of photoelectron emission
200  // <...>
201 
202  // effective tile plus wave-length shifter decay time over 4 time constants
203  unsigned int tmax = 6 * (int)ts3;
204 
205  norm = 0.0;
206  for (j = 0; j < tmax && j < nbin; j++) {
207  ntd[j] = wd1 * exp(-((float)j) / ts1) + wd2 * exp(-((float)j) / ts2) + wd3 * exp(-((float)j) / ts3);
208  norm += ntd[j];
209  }
210  // normalize pulse area to 1.0
211  for (j = 0; j < tmax && j < nbin; j++) {
212  ntd[j] /= norm;
213  }
214 
215  unsigned int t1, t2, t3, t4;
216  for (i = 0; i < tmax && i < nbin; i++) {
217  t1 = i;
218  // t2 = t1 + top*rand;
219  // ignoring jitter from optical path length
220  t2 = t1;
221  for (j = 0; j < thpd && j < nbin; j++) {
222  t3 = t2 + j;
223  for (k = 0; k < 4 * tpre && k < nbin; k++) { // here "4" is set deliberately,
224  t4 = t3 + k; // as in test fortran toy MC ...
225  if (t4 < nbin) {
226  unsigned int ntb = t4;
227  ntmp[ntb] += ntd[i] * nth[j] * ntp[k];
228  }
229  }
230  }
231  }
232 
233  // normalize for 1 GeV pulse height
234  norm = 0.;
235  for (i = 0; i < nbin; i++) {
236  norm += ntmp[i];
237  }
238 
239  for (i = 0; i < nbin; i++) {
240  ntmp[i] /= norm;
241  }
242 
243  for (i = 0; i < nbin; i++) {
244  tmphpdShape_.setShapeBin(i, ntmp[i]);
245  }
246 }
247 
249  // first create pulse shape over a range of time 0 ns to 255 ns in 1 ns steps
250  unsigned int nbin = 256;
252  std::vector<float> ntmp(nbin, 0.0); //
253 
254  const float k0 = 0.7956; // shape parameters
255  const float p2 = 1.355;
256  const float p4 = 2.327;
257  const float p1 = 4.3; // position parameter
258 
259  float norm = 0.0;
260 
261  for (unsigned int j = 0; j < 25 && j < nbin; ++j) {
262  float r0 = j - p1;
263  float sigma0 = (r0 < 0) ? p2 : p2 * p4;
264  r0 /= sigma0;
265  if (r0 < k0)
266  ntmp[j] = exp(-0.5 * r0 * r0);
267  else
268  ntmp[j] = exp(0.5 * k0 * k0 - k0 * r0);
269  norm += ntmp[j];
270  }
271  // normalize pulse area to 1.0
272  for (unsigned int j = 0; j < 25 && j < nbin; ++j) {
273  ntmp[j] /= norm;
274  hfShape_.setShapeBin(j, ntmp[j]);
275  }
276 }
278  //modified shape 206
279  //7.2 ns shift in 206
280  unsigned int nbin = 250;
281  std::array<float, 250> nt{
282  {0, 0, 0, 0, 0, 0, 0, 0.000117468,
283  0.0031549, 0.0117368, 0.0219974, 0.0305776, 0.0365429, 0.0400524, 0.0415915, 0.0416765,
284  0.0408111, 0.0394627, 0.0379353, 0.0363688, 0.0348152, 0.0332891, 0.0317923, 0.0303237,
285  0.028883, 0.0274714, 0.0260914, 0.0247462, 0.0234392, 0.0221738, 0.0209531, 0.0197793,
286  0.0186544, 0.0175796, 0.0165556, 0.0155823, 0.0146596, 0.0137866, 0.0129623, 0.0121853,
287  0.0114539, 0.0107665, 0.0101213, 0.0095162, 0.00894934, 0.00841873, 0.0079224, 0.00745841,
288  0.00702487, 0.00661995, 0.00624189, 0.00588898, 0.00555961, 0.00525223, 0.00496539, 0.0046977,
289  0.00444786, 0.00421464, 0.00399689, 0.00379353, 0.00360355, 0.00342602, 0.00326004, 0.0031048,
290  0.00295954, 0.00282355, 0.00269616, 0.00257676, 0.00246479, 0.00235972, 0.00226106, 0.00216834,
291  0.00208117, 0.00199914, 0.00192189, 0.0018491, 0.00178044, 0.00171565, 0.00165445, 0.00159659,
292  0.00154186, 0.00149003, 0.00144092, 0.00139435, 0.00135015, 0.00130816, 0.00126825, 0.00123027,
293  0.00119412, 0.00115966, 0.0011268, 0.00109544, 0.00106548, 0.00103685, 0.00100946, 0.000983242,
294  0.000958125, 0.000934047, 0.000910949, 0.000888775, 0.000867475, 0.000847, 0.000827306, 0.000808352,
295  0.000790097, 0.000772506, 0.000755545, 0.000739182, 0.000723387, 0.000708132, 0.00069339, 0.000679138,
296  0.000665352, 0.00065201, 0.000639091, 0.000626577, 0.00061445, 0.000602692, 0.000591287, 0.00058022,
297  0.000569477, 0.000559044, 0.000548908, 0.000539058, 0.000529481, 0.000520167, 0.000511106, 0.000502288,
298  0.000493704, 0.000485344, 0.000477201, 0.000469266, 0.000459912, 0.000448544, 0.000437961, 0.000428079,
299  0.000418825, 0.000410133, 0.000401945, 0.00039421, 0.000386883, 0.000379924, 0.000373298, 0.000366973,
300  0.000360922, 0.00035512, 0.000349545, 0.000344179, 0.000339003, 0.000334002, 0.000329163, 0.000324475,
301  0.000319925, 0.000315504, 0.000311204, 0.000307017, 0.000302935, 0.000298954, 0.000295066, 0.000291267,
302  0.000287553, 0.000283919, 0.000280361, 0.000276877, 0.000273462, 0.000270114, 0.000266831, 0.000263609,
303  0.000260447, 0.000257343, 0.000254295, 0.0002513, 0.000248358, 0.000245467, 0.000242625, 0.000239831,
304  0.000237083, 0.000234381, 0.000231723, 0.000229109, 0.000226536, 0.000224005, 0.000221514, 0.000219062,
305  0.000216648, 0.000214272, 0.000211933, 0.00020963, 0.000207362, 0.000205129, 0.000202929, 0.000200763,
306  0.000198629, 0.000196526, 0.000194455, 0.000192415, 0.000190405, 0.000188424, 0.000186472, 0.000184548,
307  0.000182653, 0.000180784, 0.000178943, 0.000177127, 0.000175338, 0.000173574, 0.000171835, 0.00017012,
308  0.000168429, 0.000166762, 0.000165119, 0.000163498, 0.000161899, 0.000160322, 0.000158767, 0.000157233,
309  0.000155721, 0.000154228, 0.000152756, 0.000151304, 0.000149871, 0.000148457, 0.000147062, 0.000145686,
310  0.000144327, 0.000142987, 0.000141664, 0.000140359, 0.000139071, 0.000137799, 0.000136544, 0.000135305,
311  0.000134082, 0.000132874, 0.000131682, 0.000130505, 0.000129344, 0.000128196, 0.000127064, 0.000125945,
312  0.00012484, 0.00012375, 0.000122672, 0.000121608, 0.000120558, 0.00011952, 0.000118495, 0.000117482,
313  0.000116482, 0.000115493}};
314 
316 
317  double norm = 0.;
318  for (unsigned int j = 0; j < nbin; ++j) {
319  norm += (nt[j] > 0) ? nt[j] : 0.;
320  }
321 
322  for (unsigned int j = 0; j < nbin; ++j) {
323  nt[j] /= norm;
325  }
326 }
327 
329  //Combination of all phase scan data (May,Jul,Oct2017)
330  //runs: 294736-294740, 294929-294950, 298594-298598 and 305744-305758
331 
332  unsigned int nbin = 250;
333 
334  std::array<float, 250> nt{
335  {5.22174e-12, 7.04852e-10, 3.49584e-08, 7.78029e-07, 9.11847e-06, 6.39666e-05, 0.000297587, 0.000996661,
336  0.00256618, 0.00535396, 0.00944073, 0.0145521, 0.020145, 0.0255936, 0.0303632, 0.0341078,
337  0.0366849, 0.0381183, 0.0385392, 0.0381327, 0.0370956, 0.0356113, 0.0338366, 0.0318978,
338  0.029891, 0.0278866, 0.0259336, 0.0240643, 0.0222981, 0.0206453, 0.0191097, 0.0176902,
339  0.0163832, 0.0151829, 0.0140826, 0.0130752, 0.0121533, 0.01131, 0.0105382, 0.00983178,
340  0.00918467, 0.00859143, 0.00804709, 0.0075471, 0.00708733, 0.00666406, 0.00627393, 0.00591389,
341  0.00558122, 0.00527344, 0.00498834, 0.00472392, 0.00447837, 0.00425007, 0.00403754, 0.00383947,
342  0.00365465, 0.00348199, 0.00332052, 0.00316934, 0.00302764, 0.0028947, 0.00276983, 0.00265242,
343  0.00254193, 0.00243785, 0.00233971, 0.00224709, 0.0021596, 0.00207687, 0.0019986, 0.00192447,
344  0.00185421, 0.00178756, 0.0017243, 0.00166419, 0.00160705, 0.00155268, 0.00150093, 0.00145162,
345  0.00140461, 0.00135976, 0.00131696, 0.00127607, 0.00123699, 0.00119962, 0.00116386, 0.00112963,
346  0.00109683, 0.0010654, 0.00103526, 0.00100634, 0.000978578, 0.000951917, 0.000926299, 0.000901672,
347  0.000877987, 0.000855198, 0.00083326, 0.000812133, 0.000791778, 0.000772159, 0.000753242, 0.000734994,
348  0.000717384, 0.000700385, 0.000683967, 0.000668107, 0.000652779, 0.00063796, 0.000623629, 0.000609764,
349  0.000596346, 0.000583356, 0.000570777, 0.000558592, 0.000546785, 0.00053534, 0.000524243, 0.000513481,
350  0.00050304, 0.000492907, 0.000483072, 0.000473523, 0.000464248, 0.000455238, 0.000446483, 0.000437974,
351  0.0004297, 0.000421655, 0.00041383, 0.000406216, 0.000398807, 0.000391595, 0.000384574, 0.000377736,
352  0.000371076, 0.000364588, 0.000358266, 0.000352104, 0.000346097, 0.00034024, 0.000334528, 0.000328956,
353  0.00032352, 0.000318216, 0.000313039, 0.000307986, 0.000303052, 0.000298234, 0.000293528, 0.000288931,
354  0.000284439, 0.00028005, 0.000275761, 0.000271567, 0.000267468, 0.000263459, 0.000259538, 0.000255703,
355  0.000251951, 0.00024828, 0.000244688, 0.000241172, 0.00023773, 0.000234361, 0.000231061, 0.00022783,
356  0.000224666, 0.000221566, 0.000218528, 0.000215553, 0.000212636, 0.000209778, 0.000206977, 0.00020423,
357  0.000201537, 0.000198896, 0.000196307, 0.000193767, 0.000191275, 0.000188831, 0.000186432, 0.000184079,
358  0.000181769, 0.000179502, 0.000177277, 0.000175092, 0.000172947, 0.000170841, 0.000168772, 0.000166741,
359  0.000164745, 0.000162785, 0.000160859, 0.000158967, 0.000157108, 0.00015528, 0.000153484, 0.000151719,
360  0.000149984, 0.000148278, 0.000146601, 0.000144951, 0.000143329, 0.000141734, 0.000140165, 0.000138622,
361  0.000137104, 0.00013561, 0.000134141, 0.000132695, 0.000131272, 0.000129871, 0.000128493, 0.000127136,
362  0.000125801, 0.000124486, 0.000123191, 0.000121917, 0.000120662, 0.000119426, 0.000118209, 0.00011701,
363  0.000115829, 0.000114665, 0.000113519, 0.00011239, 0.000111278, 0.000110182, 0.000109102, 0.000108037,
364  0.000106988, 0.000105954, 0.000104935, 0.00010393, 0.000102939, 0.000101963, 0.000101, 0.000100051,
365  9.91146e-05, 9.81915e-05, 9.7281e-05, 9.63831e-05, 9.54975e-05, 9.46239e-05, 9.37621e-05, 9.2912e-05,
366  9.20733e-05, 9.12458e-05}};
367 
369 
370  double norm = 0.;
371  for (unsigned int j = 0; j < nbin; ++j) {
372  norm += (nt[j] > 0) ? nt[j] : 0.;
373  }
374 
375  for (unsigned int j = 0; j < nbin; ++j) {
376  nt[j] /= norm;
378  }
379 }
380 
382  //From Jay Lawhorn: derived from data Edward Laird phase scan may2017
383  //https://indico.cern.ch/event/641978/contributions/2604491/attachments/1468666/2271582/17-05-31-hcal-hep17-pulse-shape.pdf
384  //Run numbers are 294736-294740 and 294929-294950
385 
386  unsigned int nbin = 250;
387 
388  std::array<float, 250> nt{
389  {3.97958e-29, 1.11634e-22, 9.96106e-18, 6.25334e-14, 5.08863e-11, 8.59141e-09, 4.32285e-07, 8.56617e-06,
390  8.28549e-05, 0.000461447, 0.00168052, 0.00441395, 0.00901637, 0.0151806, 0.0220314, 0.028528,
391  0.0338471, 0.0375578, 0.0395985, 0.0401567, 0.0395398, 0.0380776, 0.0360669, 0.0337474,
392  0.0312984, 0.0288457, 0.0264721, 0.0242276, 0.0221393, 0.0202181, 0.0184647, 0.0168731,
393  0.0154335, 0.0141346, 0.0129639, 0.0119094, 0.0109594, 0.0101031, 0.0093305, 0.00863267,
394  0.0080015, 0.00742977, 0.00691107, 0.00643969, 0.00601059, 0.00561931, 0.00526188, 0.00493483,
395  0.00463505, 0.00435981, 0.00410667, 0.00387348, 0.00365832, 0.00345949, 0.00327547, 0.0031049,
396  0.00294656, 0.00279938, 0.00266237, 0.00253467, 0.00241548, 0.0023041, 0.00219989, 0.00210227,
397  0.00201072, 0.00192476, 0.00184397, 0.00176795, 0.00169634, 0.00162884, 0.00156512, 0.00150494,
398  0.00144803, 0.00139418, 0.00134317, 0.00129481, 0.00124894, 0.00120537, 0.00116398, 0.00112461,
399  0.00108715, 0.00105147, 0.00101747, 0.000985042, 0.000954096, 0.000924545, 0.000896308, 0.000869311,
400  0.000843482, 0.000818758, 0.000795077, 0.000772383, 0.000750623, 0.000729747, 0.00070971, 0.000690466,
401  0.000671977, 0.000654204, 0.00063711, 0.000620663, 0.000604831, 0.000589584, 0.000574894, 0.000560735,
402  0.000547081, 0.00053391, 0.0005212, 0.000508929, 0.000497078, 0.000485628, 0.000474561, 0.000463862,
403  0.000453514, 0.000443501, 0.000433811, 0.000424429, 0.000415343, 0.00040654, 0.00039801, 0.000389741,
404  0.000381722, 0.000373944, 0.000366398, 0.000359074, 0.000351964, 0.00034506, 0.000338353, 0.000331838,
405  0.000325505, 0.00031935, 0.000313365, 0.000307544, 0.000301881, 0.000296371, 0.000291009, 0.000285788,
406  0.000280705, 0.000275755, 0.000270932, 0.000266233, 0.000261653, 0.00025719, 0.000252837, 0.000248593,
407  0.000244454, 0.000240416, 0.000236475, 0.00023263, 0.000228876, 0.000225212, 0.000221633, 0.000218138,
408  0.000214724, 0.000211389, 0.00020813, 0.000204945, 0.000201831, 0.000198787, 0.000195811, 0.0001929,
409  0.000190053, 0.000187268, 0.000184543, 0.000181876, 0.000179266, 0.000176711, 0.00017421, 0.000171761,
410  0.000169363, 0.000167014, 0.000164713, 0.000162459, 0.00016025, 0.000158086, 0.000155964, 0.000153885,
411  0.000151847, 0.000149848, 0.000147888, 0.000145966, 0.000144081, 0.000142232, 0.000140418, 0.000138638,
412  0.000136891, 0.000135177, 0.000133494, 0.000131843, 0.000130221, 0.00012863, 0.000127066, 0.000125531,
413  0.000124023, 0.000122543, 0.000121088, 0.000119658, 0.000118254, 0.000116874, 0.000115518, 0.000114185,
414  0.000112875, 0.000111587, 0.000110321, 0.000109076, 0.000107851, 0.000106648, 0.000105464, 0.000104299,
415  0.000103154, 0.000102027, 0.000100918, 9.98271e-05, 9.87537e-05, 9.76974e-05, 9.66578e-05, 9.56346e-05,
416  9.46274e-05, 9.3636e-05, 9.26599e-05, 9.16989e-05, 9.07526e-05, 8.98208e-05, 8.89032e-05, 8.79995e-05,
417  8.71093e-05, 8.62325e-05, 8.53688e-05, 8.45179e-05, 8.36796e-05, 8.28536e-05, 8.20397e-05, 8.12376e-05,
418  8.04471e-05, 7.96681e-05, 7.89002e-05, 7.81433e-05, 7.73972e-05, 7.66616e-05, 7.59364e-05, 7.52213e-05,
419  7.45163e-05, 7.3821e-05, 7.31354e-05, 7.24592e-05, 7.17923e-05, 7.11345e-05, 7.04856e-05, 6.98455e-05,
420  6.9214e-05, 6.8591e-05}};
421 
423 
424  double norm = 0.;
425  for (unsigned int j = 0; j < nbin; ++j) {
426  norm += (nt[j] > 0) ? nt[j] : 0.;
427  }
428 
429  for (unsigned int j = 0; j < nbin; ++j) {
430  nt[j] /= norm;
432  }
433 }
434 
436  unsigned int nbin = 128;
437 
438  //From Jake Anderson: toy MC convolution of SiPM pulse + WLS fiber shape + SiPM nonlinear response
439  std::array<float, 128> nt{
440  {2.782980485851731e-6, 4.518134885954626e-5, 2.7689305197392056e-4, 9.18328418900969e-4,
441  .002110072599166349, .003867856860331454, .006120046224897771, .008754774090536956,
442  0.0116469503358586, .01467007449455966, .01770489955229477, .02064621450689512,
443  .02340678093764222, .02591874610854916, .02813325527435303, 0.0300189241965647,
444  .03155968107671164, .03275234052577155, .03360415306318798, .03413048377960748,
445  .03435270899678218, .03429637464659661, .03398962975487166, .03346192884394954,
446  .03274298516247742, .03186195009136525, .03084679116113031, 0.0297238406141036,
447  .02851748748929785, .02724998816332392, .02594137274487424, .02460942736731527,
448  .02326973510736116, .02193576080366117, 0.0206189674254987, .01932895378564653,
449  0.0180736052958666, .01685925112650875, 0.0156908225633535, .01457200857138456,
450  .01350540559602467, .01249265947824805, .01153459805300423, .01063135355597282,
451  .009782474412011936, .008987026319784546, 0.00824368281357106, .007550805679909604,
452  .006906515742762193, .006308754629755056, .005755338185695127, .005244002229973356,
453  .004772441359900532, .004338341490928299, .003939406800854143, 0.00357338171220501,
454  0.0032380685079891, .002931341133259233, .002651155690306086, .002395558090237333,
455  .002162689279320922, .001950788415487319, .001758194329648101, .001583345567913682,
456  .001424779275191974, .001281129147671334, 0.00115112265163774, .001033577678808199,
457  9.273987838127585e-4, 8.315731274976846e-4, 7.451662302008696e-4, 6.673176219006913e-4,
458  5.972364609644049e-4, 5.341971801529036e-4, 4.775352065178378e-4, 4.266427928961177e-4,
459  3.8096498904225923e-4, 3.3999577417327287e-4, 3.032743659102713e-4, 2.703817158798329e-4,
460  2.4093719775272793e-4, 2.145954900503894e-4, 1.9104365317752797e-4, 1.6999839784346724e-4,
461  1.5120354022478893e-4, 1.3442763782650755e-4, 1.1946179895521507e-4, 1.0611765796993575e-4,
462  9.422550797617687e-5, 8.363258233342666e-5, 7.420147621931836e-5, 6.580869950304933e-5,
463  5.834335229919868e-5, 5.17059147771959e-5, 4.5807143072062634e-5, 4.0567063461299446e-5,
464  3.591405732740723e-5, 3.178402980354131e-5, 2.811965539165646e-5, 2.4869694240316126e-5,
465  2.1988373166730962e-5, 1.9434825899529382e-5, 1.717258740121378e-5, 1.5169137499243157e-5,
466  1.339548941011129e-5, 1.1825819079078403e-5, 1.0437131581057595e-5, 9.208961130078894e-6,
467  8.12310153137994e-6, 7.163364176588591e-6, 6.315360932244386e-6, 5.566309502463164e-6,
468  4.904859063429651e-6, 4.320934164082596e-6, 3.8055950719111903e-6, 3.350912911083174e-6,
469  2.9498580949517117e-6, 2.596200697612328e-6, 2.2844215378879293e-6, 2.0096328693141094e-6,
470  1.7675076766686654e-6, 1.5542166787225756e-6, 1.366372225473431e-6, 1.200978365778838e-6,
471  1.0553864128982371e-6, 9.272554464808518e-7, 8.145171945902259e-7, 7.153448381918271e-7}};
472 
474 
475  double norm = 0.;
476  for (unsigned int j = 0; j < nbin; ++j) {
477  norm += (nt[j] > 0) ? nt[j] : 0.;
478  }
479 
480  for (unsigned int j = 0; j < nbin; ++j) {
481  nt[j] /= norm;
483  }
484 }
485 
487  //numerical convolution of SiPM pulse + WLS fiber shape
488  static const HcalPulseShape siPMShapeMC2017(
490  return siPMShapeMC2017;
491 }
492 
494  //numerical convolution of SiPM pulse + WLS fiber shape
495  //shift: aligning 206 phase closer to 205 in order to have good reco agreement
496  static const HcalPulseShape siPMShapeMC2018(
498  return siPMShapeMC2018;
499 }
500 
502  ShapeMap::const_iterator shapeMapItr = theShapes.find(shapeType);
503  if (shapeMapItr == theShapes.end()) {
504  throw cms::Exception("HcalPulseShapes") << "unknown shapeType";
505  return hpdShape_; // should not return this, but...
506  } else {
507  return *(shapeMapItr->second);
508  }
509 }
510 
512  if (!theDbService) {
513  return defaultShape(detId);
514  }
515  int shapeType = theDbService->getHcalMCParam(detId)->signalShape();
516 
517  ShapeMap::const_iterator shapeMapItr = theShapes.find(shapeType);
518  if (shapeMapItr == theShapes.end()) {
519  return defaultShape(detId);
520  } else {
521  return *(shapeMapItr->second);
522  }
523 }
524 
526  if (!theDbService) {
527  return defaultShape(detId);
528  }
529  int shapeType = theDbService->getHcalRecoParam(detId.rawId())->pulseShapeID();
530 
531  ShapeMap::const_iterator shapeMapItr = theShapes.find(shapeType);
532  if (shapeMapItr == theShapes.end()) {
533  return defaultShape(detId);
534  } else {
535  return *(shapeMapItr->second);
536  }
537 }
538 
540  edm::LogWarning("HcalPulseShapes") << "Cannot find HCAL MC Params ";
541  HcalSubdetector subdet = detId.subdet();
542  switch (subdet) {
543  case HcalBarrel:
544  return hbShape();
545  case HcalEndcap:
546  return heShape();
547  case HcalForward:
548  return hfShape();
549  case HcalOuter:
550  //FIXME doesn't look for SiPMs
551  return hoShape(false);
552  default:
553  throw cms::Exception("HcalPulseShapes") << "unknown detId";
554  break;
555  }
556 }
557 
558 //SiPM helpers
559 
560 inline double gexp(double t, double A, double c, double t0, double s) {
561  static double const root2(sqrt(2));
562  return -A * 0.5 * exp(c * t + 0.5 * c * c * s * s - c * s) * (erf(-0.5 * root2 / s * (t - t0 + c * s * s)) - 1);
563 }
564 
565 inline double onePulse(double t, double A, double sigma, double theta, double m) {
566  return (t < theta) ? 0 : A * TMath::LogNormal(t, sigma, theta, m);
567 }
568 
570  // HO SiPM pulse shape fit from Jake Anderson ca. 2013
571  double A1(0.08757), c1(-0.5257), t01(2.4013), s1(0.6721);
572  double A2(0.007598), c2(-0.1501), t02(6.9412), s2(0.8710);
573  return gexp(t, A1, c1, t01, s1) + gexp(t, A2, c2, t02, s2);
574 }
575 
577  // taken from fit to laser measurement taken by Iouri M. in Spring 2016.
578  double A1(5.204 / 6.94419), sigma1_shape(0.5387), theta1_loc(-0.3976), m1_scale(4.428);
579  double A2(1.855 / 6.94419), sigma2_shape(0.8132), theta2_loc(7.025), m2_scale(12.29);
580  return onePulse(t, A1, sigma1_shape, theta1_loc, m1_scale) + onePulse(t, A2, sigma2_shape, theta2_loc, m2_scale);
581 }
582 
583 double HcalPulseShapes::generatePhotonTime(CLHEP::HepRandomEngine* engine, unsigned int signalShape) {
584  if (signalShape == 206)
585  return generatePhotonTime206(engine);
586  else
587  return generatePhotonTime203(engine);
588 }
589 
590 double HcalPulseShapes::generatePhotonTime203(CLHEP::HepRandomEngine* engine) {
591  double result(0.);
592  while (true) {
593  result = CLHEP::RandFlat::shoot(engine, HcalPulseShapes::Y11RANGE_);
594  if (CLHEP::RandFlat::shoot(engine, HcalPulseShapes::Y11MAX203_) < HcalPulseShapes::Y11203(result))
595  return result;
596  }
597 }
598 
599 double HcalPulseShapes::generatePhotonTime206(CLHEP::HepRandomEngine* engine) {
600  double result(0.);
601  while (true) {
602  result = CLHEP::RandFlat::shoot(engine, HcalPulseShapes::Y11RANGE_);
603  if (CLHEP::RandFlat::shoot(engine, HcalPulseShapes::Y11MAX206_) < HcalPulseShapes::Y11206(result))
604  return result;
605  }
606 }
607 
608 //Original scintillator+Y11 fit from Vasken's 2001 measurement
609 double HcalPulseShapes::Y11203(double t) { return exp(-0.0635 - 0.1518 * t + log(t) * 2.528) / 2485.9; }
610 
611 //New scintillator+Y11 model from Vasken's 2017 measurement plus a Landau correction term
612 double HcalPulseShapes::Y11206(double t) {
613  //Shifting phase to have better comparison of digi shape with data
614  //If necessary, further digi phase adjustment can be done here:
615  //SimCalorimetry/HcalSimProducers/python/hcalSimParameters_cfi.py
616  //by changing "timePhase"
617  double shift = 7.2;
618 
619  //Fit From Deconvolved Data
620  double A, n, t0, fit;
621  A = 0.104204;
622  n = 0.44064;
623  t0 = 10.0186;
624  if (t > shift)
625  fit = A * (1 - exp(-(t - shift) / n)) * exp(-(t - shift) / t0);
626  else
627  fit = 0.0;
628 
629  //Correction Term
630  double norm, mpv, sigma, corTerm;
631  norm = 0.0809882;
632  mpv = 0;
633  sigma = 20;
634  if (t > shift)
635  corTerm = norm * TMath::Landau((t - shift), mpv, sigma);
636  else
637  corTerm = 0.0;
638 
639  //Overall Y11
640  double frac = 0.11;
641  double val = (1 - frac) * fit + frac * corTerm;
642 
643  if (val >= 0)
644  return val;
645  else
646  return 0.0;
647 }
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
HcalPulseShape::setNBin
void setNBin(int n)
Definition: HcalPulseShape.cc:8
RandomServiceHelper.t2
t2
Definition: RandomServiceHelper.py:257
mps_fire.i
i
Definition: mps_fire.py:355
MessageLogger.h
dqmMemoryStats.float
float
Definition: dqmMemoryStats.py:127
dqmiodumpmetadata.n
n
Definition: dqmiodumpmetadata.py:28
reco::ParticleMasses::k0
const double k0
Definition: ParticleMasses.h:7
ESHandle.h
HcalPulseShapes::heShape
const Shape & heShape() const
Definition: HcalPulseShapes.h:31
HcalPulseShapes::defaultShape
const Shape & defaultShape(const HcalDetId &detId) const
in case of conditions problems
Definition: HcalPulseShapes.cc:539
HcalPulseShapes::theShapes
ShapeMap theShapes
Definition: HcalPulseShapes.h:117
HcalPulseShapes::~HcalPulseShapes
~HcalPulseShapes()
Definition: HcalPulseShapes.cc:141
nt
int nt
Definition: AMPTWrapper.h:42
cropTnPTrees.frac
frac
Definition: cropTnPTrees.py:18
HcalPulseShapes.h
HcalPulseShapes::hpdShapeMC_v3
Shape hpdShapeMC_v3
Definition: HcalPulseShapes.h:113
HcalPulseShapes::nBinsSiPM_
static const int nBinsSiPM_
Definition: HcalPulseShapes.h:42
HcalPulseShapes::shapeForReco
const Shape & shapeForReco(const HcalDetId &detId) const
Definition: HcalPulseShapes.cc:525
L1TowerCalibrationProducer_cfi.fit
fit
Definition: L1TowerCalibrationProducer_cfi.py:36
indexGen.s2
s2
Definition: indexGen.py:107
HcalPulseShapes::computeHFShape
void computeHFShape()
Definition: HcalPulseShapes.cc:248
HcalPulseShapes::computeSiPMShapeData2017
void computeSiPMShapeData2017()
Definition: HcalPulseShapes.cc:381
HcalBarrel
Definition: HcalAssistant.h:33
HcalPulseShapes::computeSiPMShapeHE206
const HcalPulseShape & computeSiPMShapeHE206()
Definition: HcalPulseShapes.cc:493
tmax
static const double tmax[3]
Definition: CastorTimeSlew.cc:7
HcalDbService::getHcalMCParam
const HcalMCParam * getHcalMCParam(const HcalGenericDetId &fId) const
Definition: HcalDbService.cc:398
alignCSCRings.s
s
Definition: alignCSCRings.py:92
RandomServiceHelper.t1
t1
Definition: RandomServiceHelper.py:256
edm::EventSetup::get
T get() const
Definition: EventSetup.h:73
FrontierCondition_GT_autoExpress_cfi.t0
t0
Definition: FrontierCondition_GT_autoExpress_cfi.py:148
HcalPulseShapes::hpdShape_
Shape hpdShape_
Definition: HcalPulseShapes.h:110
visualization-live-secondInstance_cfg.m
m
Definition: visualization-live-secondInstance_cfg.py:72
HcalPulseShapes::generatePhotonTime203
static double generatePhotonTime203(CLHEP::HepRandomEngine *engine)
Definition: HcalPulseShapes.cc:590
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
HcalPulseShapes::computeSiPMShapeHO
void computeSiPMShapeHO()
Definition: HcalPulseShapes.cc:435
HcalPulseShapes::Y11206
static double Y11206(double t)
Definition: HcalPulseShapes.cc:612
HcalPulseShapes::getShape
const Shape & getShape(int shapeType) const
Definition: HcalPulseShapes.cc:501
edm::ESHandle< HcalDbService >
p2
double p2[4]
Definition: TauolaWrapper.h:90
HcalPulseShapes::convolve
static std::vector< double > convolve(unsigned nbin, F1 f1, F2 f2)
Definition: HcalPulseShapes.h:57
HcalOuter
Definition: HcalAssistant.h:35
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
gexp
double gexp(double t, double A, double c, double t0, double s)
Definition: HcalPulseShapes.cc:560
HcalPulseShapes::hbShape
const Shape & hbShape() const
Definition: HcalPulseShapes.h:30
dqmdumpme.k
k
Definition: dqmdumpme.py:60
HcalDbRecord.h
HcalPulseShapes::Y11203
static double Y11203(double t)
Definition: HcalPulseShapes.cc:609
OrderedSet.t
t
Definition: OrderedSet.py:90
RandomServiceHelper.t3
t3
Definition: RandomServiceHelper.py:258
HcalPulseShapes::Y11RANGE_
static constexpr float Y11RANGE_
Definition: HcalPulseShapes.h:47
HcalPulseShapes::HcalPulseShapes
HcalPulseShapes()
Definition: HcalPulseShapes.cc:15
HcalDbService::getHcalRecoParam
const HcalRecoParam * getHcalRecoParam(const HcalGenericDetId &fId) const
Definition: HcalDbService.cc:405
HcalPulseShapes::siPMShapeData2018_
Shape siPMShapeData2018_
Definition: HcalPulseShapes.h:111
edm::LogWarning
Definition: MessageLogger.h:141
HcalPulseShapes::hpdBV30Shape_v2
Shape hpdBV30Shape_v2
Definition: HcalPulseShapes.h:114
HcalPulseShapes::normalize
static std::vector< double > normalize(std::vector< double > nt, unsigned nbin)
Definition: HcalPulseShapes.h:69
onePulse
double onePulse(double t, double A, double sigma, double theta, double m)
Definition: HcalPulseShapes.cc:565
HcalPulseShapes::hpdBV30ShapeMC_v2
Shape hpdBV30ShapeMC_v2
Definition: HcalPulseShapes.h:114
A
HcalPulseShapes::analyticPulseShapeSiPMHE
static double analyticPulseShapeSiPMHE(double t)
Definition: HcalPulseShapes.cc:576
HcalDetId::subdet
constexpr HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:138
HcalDetId
Definition: HcalDetId.h:12
createfilelist.int
int
Definition: createfilelist.py:10
HcalPulseShapes::beginRun
void beginRun(edm::EventSetup const &es)
Definition: HcalPulseShapes.cc:143
alignmentValidation.c1
c1
do drawing
Definition: alignmentValidation.py:1025
p4
double p4[4]
Definition: TauolaWrapper.h:92
HcalPulseShapes::computeSiPMShapeData2018
void computeSiPMShapeData2018()
Definition: HcalPulseShapes.cc:328
p1
double p1[4]
Definition: TauolaWrapper.h:89
edm::EventSetup
Definition: EventSetup.h:57
HltBtagPostValidation_cff.c
c
Definition: HltBtagPostValidation_cff.py:31
get
#define get
HcalPulseShapes::Y11MAX203_
static constexpr float Y11MAX203_
Definition: HcalPulseShapes.h:48
HcalPulseShapes::siPMShapeMCRecoRun3_
Shape siPMShapeMCRecoRun3_
Definition: HcalPulseShapes.h:111
HcalPulseShape
Definition: HcalPulseShape.h:6
HcalSubdetector
HcalSubdetector
Definition: HcalAssistant.h:31
HcalPulseShapes::hoShape
const Shape & hoShape(bool sipm=false) const
Definition: HcalPulseShapes.h:33
HcalForward
Definition: HcalAssistant.h:36
HcalPulseShapes::normalizeShift
static std::vector< double > normalizeShift(std::vector< double > nt, unsigned nbin, int shift)
Definition: HcalPulseShapes.h:83
HcalPulseShapes::generatePhotonTime
static double generatePhotonTime(CLHEP::HepRandomEngine *engine, unsigned int signalShape)
Definition: HcalPulseShapes.cc:583
HcalPulseShapes::hfShape
const Shape & hfShape() const
Definition: HcalPulseShapes.h:32
HcalPulseShapes::hpdShape_v2
Shape hpdShape_v2
Definition: HcalPulseShapes.h:112
HcalMCParam::signalShape
unsigned int signalShape() const
Definition: HcalMCParam.h:38
MaterialEffects_cfi.A
A
Definition: MaterialEffects_cfi.py:11
HcalDbService
Definition: HcalDbService.h:26
HcalPulseShapes::shape
const Shape & shape(const HcalDetId &detId) const
automatically figures out which shape to return
Definition: HcalPulseShapes.cc:511
heppy_batch.val
val
Definition: heppy_batch.py:351
DetId::rawId
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
HcalPulseShapes::computeHPDShape
void computeHPDShape(float, float, float, float, float, float, float, float, Shape &)
Definition: HcalPulseShapes.cc:152
HcalPulseShapes::siPMShapeHO_
Shape siPMShapeHO_
Definition: HcalPulseShapes.h:110
edm::shift
static unsigned const int shift
Definition: LuminosityBlockID.cc:7
HcalEndcap
Definition: HcalAssistant.h:34
HcalPulseShapes::hpdShapeMC_v2
Shape hpdShapeMC_v2
Definition: HcalPulseShapes.h:112
Exception
Definition: hltDiff.cc:246
HcalPulseShapes::theDbService
const HcalDbService * theDbService
Definition: HcalPulseShapes.h:115
ZMuMuCategoriesSequences_cff.nbin
nbin
Definition: ZMuMuCategoriesSequences_cff.py:25
HcalPulseShapes::computeSiPMShapeMCRecoRun3
void computeSiPMShapeMCRecoRun3()
Definition: HcalPulseShapes.cc:277
EventSetup.h
Exception.h
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
HcalPulseShapes::Y11MAX206_
static constexpr float Y11MAX206_
Definition: HcalPulseShapes.h:49
HcalPulseShapes::generatePhotonTime206
static double generatePhotonTime206(CLHEP::HepRandomEngine *engine)
Definition: HcalPulseShapes.cc:599
HcalPulseShapes::computeSiPMShapeHE203
const HcalPulseShape & computeSiPMShapeHE203()
Definition: HcalPulseShapes.cc:486
HcalPulseShapes::hpdShape_v3
Shape hpdShape_v3
Definition: HcalPulseShapes.h:113
mps_fire.result
result
Definition: mps_fire.py:303
HcalDbRecord
Definition: HcalDbRecord.h:30
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
HcalPulseShapes::siPMShapeData2017_
Shape siPMShapeData2017_
Definition: HcalPulseShapes.h:111
fit
Definition: CombinedChiSquaredLikelihood.h:6
HcalPulseShapes::analyticPulseShapeSiPMHO
static double analyticPulseShapeSiPMHO(double t)
Definition: HcalPulseShapes.cc:569
HcalPulseShapes::hfShape_
Shape hfShape_
Definition: HcalPulseShapes.h:110
HcalPulseShape::setShapeBin
void setShapeBin(int i, float f)
Definition: HcalPulseShape.cc:13