CMS 3D CMS Logo

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