CMS 3D CMS Logo
HcalPulseShapes.cc
Go to the documentation of this file.
1 #include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShapes.h"
2 #include "FWCore/Utilities/interface/Exception.h"
3 #include "FWCore/Framework/interface/ESHandle.h"
4 #include "FWCore/Framework/interface/EventSetup.h"
5 #include "FWCore/MessageLogger/interface/MessageLogger.h"
6 #include "CLHEP/Random/RandFlat.h"
7 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.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 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);
119  theShapes[124] = &hpdBV30ShapeMC_v2;
120 
121  // HF and SiPM
122 
123  computeHFShape();
124  computeSiPMShapeHO();
125  computeSiPMShapeData2017();
126  computeSiPMShapeData2018();
127 
128  theShapes[201] = &siPMShapeHO_;
129  theShapes[202] = theShapes[201];
130  theShapes[203] = &(computeSiPMShapeHE203());
131  theShapes[205] = &siPMShapeData2017_;
132  theShapes[206] = &(computeSiPMShapeHE206());
133  theShapes[207] = &siPMShapeData2018_;
134  theShapes[301] = &hfShape_;
135  //theShapes[401] = new CaloCachedShapeIntegrator(&theZDCShape);
136 }
137 
138 HcalPulseShapes::~HcalPulseShapes() {}
139 
140 void HcalPulseShapes::beginRun(edm::EventSetup const& es) {
141  edm::ESHandle<HcalDbService> conditions;
142  es.get<HcalDbRecord>().get(conditions);
143  theDbService = conditions.product();
144 }
145 
146 void HcalPulseShapes::beginRun(const HcalDbService* conditions) { theDbService = conditions; }
147 
148 //void HcalPulseShapes::computeHPDShape()
149 void HcalPulseShapes::computeHPDShape(
150  float ts1, float ts2, float ts3, float thpd, float tpre, float wd1, float wd2, float wd3, Shape& tmphpdShape_) {
151  // pulse shape time constants in ns
152  /*
153  const float ts1 = 8.; // scintillation time constants : 1,2,3
154  const float ts2 = 10.;
155  const float ts3 = 29.3;
156  const float thpd = 4.; // HPD current collection drift time
157  const float tpre = 9.; // preamp time constant (refit on TB04 data)
158 
159  const float wd1 = 2.; // relative weights of decay exponents
160  const float wd2 = 0.7;
161  const float wd3 = 1.;
162 */
163  // pulse shape components over a range of time 0 ns to 255 ns in 1 ns steps
164  unsigned int nbin = 256;
165  tmphpdShape_.setNBin(nbin);
166  std::vector<float> ntmp(nbin, 0.0); // zeroing output pulse shape
167  std::vector<float> nth(nbin, 0.0); // zeroing HPD drift shape
168  std::vector<float> ntp(nbin, 0.0); // zeroing Binkley preamp shape
169  std::vector<float> ntd(nbin, 0.0); // zeroing Scintillator decay shape
170 
171  unsigned int i, j, k;
172  float norm;
173 
174  // HPD starts at I and rises to 2I in thpd of time
175  norm = 0.0;
176  for (j = 0; j < thpd && j < nbin; j++) {
177  nth[j] = 1.0 + ((float)j) / thpd;
178  norm += nth[j];
179  }
180  // normalize integrated current to 1.0
181  for (j = 0; j < thpd && j < nbin; j++) {
182  nth[j] /= norm;
183  }
184 
185  // Binkley shape over 6 time constants
186  norm = 0.0;
187  for (j = 0; j < 6 * tpre && j < nbin; j++) {
188  ntp[j] = ((float)j) * exp(-((float)(j * j)) / (tpre * tpre));
189  norm += ntp[j];
190  }
191  // normalize pulse area to 1.0
192  for (j = 0; j < 6 * tpre && j < nbin; j++) {
193  ntp[j] /= norm;
194  }
195 
196  // ignore stochastic variation of photoelectron emission
197  // <...>
198 
199  // effective tile plus wave-length shifter decay time over 4 time constants
200  unsigned int tmax = 6 * (int)ts3;
201 
202  norm = 0.0;
203  for (j = 0; j < tmax && j < nbin; j++) {
204  ntd[j] = wd1 * exp(-((float)j) / ts1) + wd2 * exp(-((float)j) / ts2) + wd3 * exp(-((float)j) / ts3);
205  norm += ntd[j];
206  }
207  // normalize pulse area to 1.0
208  for (j = 0; j < tmax && j < nbin; j++) {
209  ntd[j] /= norm;
210  }
211 
212  unsigned int t1, t2, t3, t4;
213  for (i = 0; i < tmax && i < nbin; i++) {
214  t1 = i;
215  // t2 = t1 + top*rand;
216  // ignoring jitter from optical path length
217  t2 = t1;
218  for (j = 0; j < thpd && j < nbin; j++) {
219  t3 = t2 + j;
220  for (k = 0; k < 4 * tpre && k < nbin; k++) { // here "4" is set deliberately,
221  t4 = t3 + k; // as in test fortran toy MC ...
222  if (t4 < nbin) {
223  unsigned int ntb = t4;
224  ntmp[ntb] += ntd[i] * nth[j] * ntp[k];
225  }
226  }
227  }
228  }
229 
230  // normalize for 1 GeV pulse height
231  norm = 0.;
232  for (i = 0; i < nbin; i++) {
233  norm += ntmp[i];
234  }
235 
236  for (i = 0; i < nbin; i++) {
237  ntmp[i] /= norm;
238  }
239 
240  for (i = 0; i < nbin; i++) {
241  tmphpdShape_.setShapeBin(i, ntmp[i]);
242  }
243 }
244 
245 void HcalPulseShapes::computeHFShape() {
246  // first create pulse shape over a range of time 0 ns to 255 ns in 1 ns steps
247  unsigned int nbin = 256;
248  hfShape_.setNBin(nbin);
249  std::vector<float> ntmp(nbin, 0.0); //
250 
251  const float k0 = 0.7956; // shape parameters
252  const float p2 = 1.355;
253  const float p4 = 2.327;
254  const float p1 = 4.3; // position parameter
255 
256  float norm = 0.0;
257 
258  for (unsigned int j = 0; j < 25 && j < nbin; ++j) {
259  float r0 = j - p1;
260  float sigma0 = (r0 < 0) ? p2 : p2 * p4;
261  r0 /= sigma0;
262  if (r0 < k0)
263  ntmp[j] = exp(-0.5 * r0 * r0);
264  else
265  ntmp[j] = exp(0.5 * k0 * k0 - k0 * r0);
266  norm += ntmp[j];
267  }
268  // normalize pulse area to 1.0
269  for (unsigned int j = 0; j < 25 && j < nbin; ++j) {
270  ntmp[j] /= norm;
271  hfShape_.setShapeBin(j, ntmp[j]);
272  }
273 }
274 
275 void HcalPulseShapes::computeSiPMShapeData2018() {
276  //Combination of all phase scan data (May,Jul,Oct2017)
277  //runs: 294736-294740, 294929-294950, 298594-298598 and 305744-305758
278 
279  unsigned int nbin = 250;
280 
281  std::array<float, 250> nt{
282  {5.22174e-12, 7.04852e-10, 3.49584e-08, 7.78029e-07, 9.11847e-06, 6.39666e-05, 0.000297587, 0.000996661,
283  0.00256618, 0.00535396, 0.00944073, 0.0145521, 0.020145, 0.0255936, 0.0303632, 0.0341078,
284  0.0366849, 0.0381183, 0.0385392, 0.0381327, 0.0370956, 0.0356113, 0.0338366, 0.0318978,
285  0.029891, 0.0278866, 0.0259336, 0.0240643, 0.0222981, 0.0206453, 0.0191097, 0.0176902,
286  0.0163832, 0.0151829, 0.0140826, 0.0130752, 0.0121533, 0.01131, 0.0105382, 0.00983178,
287  0.00918467, 0.00859143, 0.00804709, 0.0075471, 0.00708733, 0.00666406, 0.00627393, 0.00591389,
288  0.00558122, 0.00527344, 0.00498834, 0.00472392, 0.00447837, 0.00425007, 0.00403754, 0.00383947,
289  0.00365465, 0.00348199, 0.00332052, 0.00316934, 0.00302764, 0.0028947, 0.00276983, 0.00265242,
290  0.00254193, 0.00243785, 0.00233971, 0.00224709, 0.0021596, 0.00207687, 0.0019986, 0.00192447,
291  0.00185421, 0.00178756, 0.0017243, 0.00166419, 0.00160705, 0.00155268, 0.00150093, 0.00145162,
292  0.00140461, 0.00135976, 0.00131696, 0.00127607, 0.00123699, 0.00119962, 0.00116386, 0.00112963,
293  0.00109683, 0.0010654, 0.00103526, 0.00100634, 0.000978578, 0.000951917, 0.000926299, 0.000901672,
294  0.000877987, 0.000855198, 0.00083326, 0.000812133, 0.000791778, 0.000772159, 0.000753242, 0.000734994,
295  0.000717384, 0.000700385, 0.000683967, 0.000668107, 0.000652779, 0.00063796, 0.000623629, 0.000609764,
296  0.000596346, 0.000583356, 0.000570777, 0.000558592, 0.000546785, 0.00053534, 0.000524243, 0.000513481,
297  0.00050304, 0.000492907, 0.000483072, 0.000473523, 0.000464248, 0.000455238, 0.000446483, 0.000437974,
298  0.0004297, 0.000421655, 0.00041383, 0.000406216, 0.000398807, 0.000391595, 0.000384574, 0.000377736,
299  0.000371076, 0.000364588, 0.000358266, 0.000352104, 0.000346097, 0.00034024, 0.000334528, 0.000328956,
300  0.00032352, 0.000318216, 0.000313039, 0.000307986, 0.000303052, 0.000298234, 0.000293528, 0.000288931,
301  0.000284439, 0.00028005, 0.000275761, 0.000271567, 0.000267468, 0.000263459, 0.000259538, 0.000255703,
302  0.000251951, 0.00024828, 0.000244688, 0.000241172, 0.00023773, 0.000234361, 0.000231061, 0.00022783,
303  0.000224666, 0.000221566, 0.000218528, 0.000215553, 0.000212636, 0.000209778, 0.000206977, 0.00020423,
304  0.000201537, 0.000198896, 0.000196307, 0.000193767, 0.000191275, 0.000188831, 0.000186432, 0.000184079,
305  0.000181769, 0.000179502, 0.000177277, 0.000175092, 0.000172947, 0.000170841, 0.000168772, 0.000166741,
306  0.000164745, 0.000162785, 0.000160859, 0.000158967, 0.000157108, 0.00015528, 0.000153484, 0.000151719,
307  0.000149984, 0.000148278, 0.000146601, 0.000144951, 0.000143329, 0.000141734, 0.000140165, 0.000138622,
308  0.000137104, 0.00013561, 0.000134141, 0.000132695, 0.000131272, 0.000129871, 0.000128493, 0.000127136,
309  0.000125801, 0.000124486, 0.000123191, 0.000121917, 0.000120662, 0.000119426, 0.000118209, 0.00011701,
310  0.000115829, 0.000114665, 0.000113519, 0.00011239, 0.000111278, 0.000110182, 0.000109102, 0.000108037,
311  0.000106988, 0.000105954, 0.000104935, 0.00010393, 0.000102939, 0.000101963, 0.000101, 0.000100051,
312  9.91146e-05, 9.81915e-05, 9.7281e-05, 9.63831e-05, 9.54975e-05, 9.46239e-05, 9.37621e-05, 9.2912e-05,
313  9.20733e-05, 9.12458e-05}};
314 
315  siPMShapeData2018_.setNBin(nbin);
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;
324  siPMShapeData2018_.setShapeBin(j, nt[j]);
325  }
326 }
327 
328 void HcalPulseShapes::computeSiPMShapeData2017() {
329  //From Jay Lawhorn: derived from data Edward Laird phase scan may2017
330  //https://indico.cern.ch/event/641978/contributions/2604491/attachments/1468666/2271582/17-05-31-hcal-hep17-pulse-shape.pdf
331  //Run numbers are 294736-294740 and 294929-294950
332 
333  unsigned int nbin = 250;
334 
335  std::array<float, 250> nt{
336  {3.97958e-29, 1.11634e-22, 9.96106e-18, 6.25334e-14, 5.08863e-11, 8.59141e-09, 4.32285e-07, 8.56617e-06,
337  8.28549e-05, 0.000461447, 0.00168052, 0.00441395, 0.00901637, 0.0151806, 0.0220314, 0.028528,
338  0.0338471, 0.0375578, 0.0395985, 0.0401567, 0.0395398, 0.0380776, 0.0360669, 0.0337474,
339  0.0312984, 0.0288457, 0.0264721, 0.0242276, 0.0221393, 0.0202181, 0.0184647, 0.0168731,
340  0.0154335, 0.0141346, 0.0129639, 0.0119094, 0.0109594, 0.0101031, 0.0093305, 0.00863267,
341  0.0080015, 0.00742977, 0.00691107, 0.00643969, 0.00601059, 0.00561931, 0.00526188, 0.00493483,
342  0.00463505, 0.00435981, 0.00410667, 0.00387348, 0.00365832, 0.00345949, 0.00327547, 0.0031049,
343  0.00294656, 0.00279938, 0.00266237, 0.00253467, 0.00241548, 0.0023041, 0.00219989, 0.00210227,
344  0.00201072, 0.00192476, 0.00184397, 0.00176795, 0.00169634, 0.00162884, 0.00156512, 0.00150494,
345  0.00144803, 0.00139418, 0.00134317, 0.00129481, 0.00124894, 0.00120537, 0.00116398, 0.00112461,
346  0.00108715, 0.00105147, 0.00101747, 0.000985042, 0.000954096, 0.000924545, 0.000896308, 0.000869311,
347  0.000843482, 0.000818758, 0.000795077, 0.000772383, 0.000750623, 0.000729747, 0.00070971, 0.000690466,
348  0.000671977, 0.000654204, 0.00063711, 0.000620663, 0.000604831, 0.000589584, 0.000574894, 0.000560735,
349  0.000547081, 0.00053391, 0.0005212, 0.000508929, 0.000497078, 0.000485628, 0.000474561, 0.000463862,
350  0.000453514, 0.000443501, 0.000433811, 0.000424429, 0.000415343, 0.00040654, 0.00039801, 0.000389741,
351  0.000381722, 0.000373944, 0.000366398, 0.000359074, 0.000351964, 0.00034506, 0.000338353, 0.000331838,
352  0.000325505, 0.00031935, 0.000313365, 0.000307544, 0.000301881, 0.000296371, 0.000291009, 0.000285788,
353  0.000280705, 0.000275755, 0.000270932, 0.000266233, 0.000261653, 0.00025719, 0.000252837, 0.000248593,
354  0.000244454, 0.000240416, 0.000236475, 0.00023263, 0.000228876, 0.000225212, 0.000221633, 0.000218138,
355  0.000214724, 0.000211389, 0.00020813, 0.000204945, 0.000201831, 0.000198787, 0.000195811, 0.0001929,
356  0.000190053, 0.000187268, 0.000184543, 0.000181876, 0.000179266, 0.000176711, 0.00017421, 0.000171761,
357  0.000169363, 0.000167014, 0.000164713, 0.000162459, 0.00016025, 0.000158086, 0.000155964, 0.000153885,
358  0.000151847, 0.000149848, 0.000147888, 0.000145966, 0.000144081, 0.000142232, 0.000140418, 0.000138638,
359  0.000136891, 0.000135177, 0.000133494, 0.000131843, 0.000130221, 0.00012863, 0.000127066, 0.000125531,
360  0.000124023, 0.000122543, 0.000121088, 0.000119658, 0.000118254, 0.000116874, 0.000115518, 0.000114185,
361  0.000112875, 0.000111587, 0.000110321, 0.000109076, 0.000107851, 0.000106648, 0.000105464, 0.000104299,
362  0.000103154, 0.000102027, 0.000100918, 9.98271e-05, 9.87537e-05, 9.76974e-05, 9.66578e-05, 9.56346e-05,
363  9.46274e-05, 9.3636e-05, 9.26599e-05, 9.16989e-05, 9.07526e-05, 8.98208e-05, 8.89032e-05, 8.79995e-05,
364  8.71093e-05, 8.62325e-05, 8.53688e-05, 8.45179e-05, 8.36796e-05, 8.28536e-05, 8.20397e-05, 8.12376e-05,
365  8.04471e-05, 7.96681e-05, 7.89002e-05, 7.81433e-05, 7.73972e-05, 7.66616e-05, 7.59364e-05, 7.52213e-05,
366  7.45163e-05, 7.3821e-05, 7.31354e-05, 7.24592e-05, 7.17923e-05, 7.11345e-05, 7.04856e-05, 6.98455e-05,
367  6.9214e-05, 6.8591e-05}};
368 
369  siPMShapeData2017_.setNBin(nbin);
370 
371  double norm = 0.;
372  for (unsigned int j = 0; j < nbin; ++j) {
373  norm += (nt[j] > 0) ? nt[j] : 0.;
374  }
375 
376  for (unsigned int j = 0; j < nbin; ++j) {
377  nt[j] /= norm;
378  siPMShapeData2017_.setShapeBin(j, nt[j]);
379  }
380 }
381 
382 void HcalPulseShapes::computeSiPMShapeHO() {
383  unsigned int nbin = 128;
384 
385  //From Jake Anderson: toy MC convolution of SiPM pulse + WLS fiber shape + SiPM nonlinear response
386  std::array<float, 128> nt{
387  {2.782980485851731e-6, 4.518134885954626e-5, 2.7689305197392056e-4, 9.18328418900969e-4,
388  .002110072599166349, .003867856860331454, .006120046224897771, .008754774090536956,
389  0.0116469503358586, .01467007449455966, .01770489955229477, .02064621450689512,
390  .02340678093764222, .02591874610854916, .02813325527435303, 0.0300189241965647,
391  .03155968107671164, .03275234052577155, .03360415306318798, .03413048377960748,
392  .03435270899678218, .03429637464659661, .03398962975487166, .03346192884394954,
393  .03274298516247742, .03186195009136525, .03084679116113031, 0.0297238406141036,
394  .02851748748929785, .02724998816332392, .02594137274487424, .02460942736731527,
395  .02326973510736116, .02193576080366117, 0.0206189674254987, .01932895378564653,
396  0.0180736052958666, .01685925112650875, 0.0156908225633535, .01457200857138456,
397  .01350540559602467, .01249265947824805, .01153459805300423, .01063135355597282,
398  .009782474412011936, .008987026319784546, 0.00824368281357106, .007550805679909604,
399  .006906515742762193, .006308754629755056, .005755338185695127, .005244002229973356,
400  .004772441359900532, .004338341490928299, .003939406800854143, 0.00357338171220501,
401  0.0032380685079891, .002931341133259233, .002651155690306086, .002395558090237333,
402  .002162689279320922, .001950788415487319, .001758194329648101, .001583345567913682,
403  .001424779275191974, .001281129147671334, 0.00115112265163774, .001033577678808199,
404  9.273987838127585e-4, 8.315731274976846e-4, 7.451662302008696e-4, 6.673176219006913e-4,
405  5.972364609644049e-4, 5.341971801529036e-4, 4.775352065178378e-4, 4.266427928961177e-4,
406  3.8096498904225923e-4, 3.3999577417327287e-4, 3.032743659102713e-4, 2.703817158798329e-4,
407  2.4093719775272793e-4, 2.145954900503894e-4, 1.9104365317752797e-4, 1.6999839784346724e-4,
408  1.5120354022478893e-4, 1.3442763782650755e-4, 1.1946179895521507e-4, 1.0611765796993575e-4,
409  9.422550797617687e-5, 8.363258233342666e-5, 7.420147621931836e-5, 6.580869950304933e-5,
410  5.834335229919868e-5, 5.17059147771959e-5, 4.5807143072062634e-5, 4.0567063461299446e-5,
411  3.591405732740723e-5, 3.178402980354131e-5, 2.811965539165646e-5, 2.4869694240316126e-5,
412  2.1988373166730962e-5, 1.9434825899529382e-5, 1.717258740121378e-5, 1.5169137499243157e-5,
413  1.339548941011129e-5, 1.1825819079078403e-5, 1.0437131581057595e-5, 9.208961130078894e-6,
414  8.12310153137994e-6, 7.163364176588591e-6, 6.315360932244386e-6, 5.566309502463164e-6,
415  4.904859063429651e-6, 4.320934164082596e-6, 3.8055950719111903e-6, 3.350912911083174e-6,
416  2.9498580949517117e-6, 2.596200697612328e-6, 2.2844215378879293e-6, 2.0096328693141094e-6,
417  1.7675076766686654e-6, 1.5542166787225756e-6, 1.366372225473431e-6, 1.200978365778838e-6,
418  1.0553864128982371e-6, 9.272554464808518e-7, 8.145171945902259e-7, 7.153448381918271e-7}};
419 
420  siPMShapeHO_.setNBin(nbin);
421 
422  double norm = 0.;
423  for (unsigned int j = 0; j < nbin; ++j) {
424  norm += (nt[j] > 0) ? nt[j] : 0.;
425  }
426 
427  for (unsigned int j = 0; j < nbin; ++j) {
428  nt[j] /= norm;
429  siPMShapeHO_.setShapeBin(j, nt[j]);
430  }
431 }
432 
433 const HcalPulseShape& HcalPulseShapes::computeSiPMShapeHE203() {
434  //numerical convolution of SiPM pulse + WLS fiber shape
435  static const HcalPulseShape siPMShapeMC2017(
436  normalize(convolve(nBinsSiPM_, analyticPulseShapeSiPMHE, Y11203), nBinsSiPM_), nBinsSiPM_);
437  return siPMShapeMC2017;
438 }
439 
440 const HcalPulseShape& HcalPulseShapes::computeSiPMShapeHE206() {
441  //numerical convolution of SiPM pulse + WLS fiber shape
442  //shift: aligning 206 phase closer to 205 in order to have good reco agreement
443  static const HcalPulseShape siPMShapeMC2018(
444  normalizeShift(convolve(nBinsSiPM_, analyticPulseShapeSiPMHE, Y11206), nBinsSiPM_, -2), nBinsSiPM_);
445  return siPMShapeMC2018;
446 }
447 
448 const HcalPulseShapes::Shape& HcalPulseShapes::getShape(int shapeType) const {
449  ShapeMap::const_iterator shapeMapItr = theShapes.find(shapeType);
450  if (shapeMapItr == theShapes.end()) {
451  throw cms::Exception("HcalPulseShapes") << "unknown shapeType";
452  return hpdShape_; // should not return this, but...
453  } else {
454  return *(shapeMapItr->second);
455  }
456 }
457 
458 const HcalPulseShapes::Shape& HcalPulseShapes::shape(const HcalDetId& detId) const {
459  if (!theDbService) {
460  return defaultShape(detId);
461  }
462  int shapeType = theDbService->getHcalMCParam(detId)->signalShape();
463 
464  ShapeMap::const_iterator shapeMapItr = theShapes.find(shapeType);
465  if (shapeMapItr == theShapes.end()) {
466  return defaultShape(detId);
467  } else {
468  return *(shapeMapItr->second);
469  }
470 }
471 
472 const HcalPulseShapes::Shape& HcalPulseShapes::shapeForReco(const HcalDetId& detId) const {
473  if (!theDbService) {
474  return defaultShape(detId);
475  }
476  int shapeType = theDbService->getHcalRecoParam(detId.rawId())->pulseShapeID();
477 
478  ShapeMap::const_iterator shapeMapItr = theShapes.find(shapeType);
479  if (shapeMapItr == theShapes.end()) {
480  return defaultShape(detId);
481  } else {
482  return *(shapeMapItr->second);
483  }
484 }
485 
486 const HcalPulseShapes::Shape& HcalPulseShapes::defaultShape(const HcalDetId& detId) const {
487  edm::LogWarning("HcalPulseShapes") << "Cannot find HCAL MC Params ";
488  HcalSubdetector subdet = detId.subdet();
489  switch (subdet) {
490  case HcalBarrel:
491  return hbShape();
492  case HcalEndcap:
493  return heShape();
494  case HcalForward:
495  return hfShape();
496  case HcalOuter:
497  //FIXME doesn't look for SiPMs
498  return hoShape(false);
499  default:
500  throw cms::Exception("HcalPulseShapes") << "unknown detId";
501  break;
502  }
503 }
504 
505 //SiPM helpers
506 
507 inline double gexp(double t, double A, double c, double t0, double s) {
508  static double const root2(sqrt(2));
509  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);
510 }
511 
512 inline double onePulse(double t, double A, double sigma, double theta, double m) {
513  return (t < theta) ? 0 : A * TMath::LogNormal(t, sigma, theta, m);
514 }
515 
516 double HcalPulseShapes::analyticPulseShapeSiPMHO(double t) {
517  // HO SiPM pulse shape fit from Jake Anderson ca. 2013
518  double A1(0.08757), c1(-0.5257), t01(2.4013), s1(0.6721);
519  double A2(0.007598), c2(-0.1501), t02(6.9412), s2(0.8710);
520  return gexp(t, A1, c1, t01, s1) + gexp(t, A2, c2, t02, s2);
521 }
522 
523 double HcalPulseShapes::analyticPulseShapeSiPMHE(double t) {
524  // taken from fit to laser measurement taken by Iouri M. in Spring 2016.
525  double A1(5.204 / 6.94419), sigma1_shape(0.5387), theta1_loc(-0.3976), m1_scale(4.428);
526  double A2(1.855 / 6.94419), sigma2_shape(0.8132), theta2_loc(7.025), m2_scale(12.29);
527  return onePulse(t, A1, sigma1_shape, theta1_loc, m1_scale) + onePulse(t, A2, sigma2_shape, theta2_loc, m2_scale);
528 }
529 
530 double HcalPulseShapes::generatePhotonTime(CLHEP::HepRandomEngine* engine, unsigned int signalShape) {
531  if (signalShape == 206)
532  return generatePhotonTime206(engine);
533  else
534  return generatePhotonTime203(engine);
535 }
536 
537 double HcalPulseShapes::generatePhotonTime203(CLHEP::HepRandomEngine* engine) {
538  double result(0.);
539  while (true) {
540  result = CLHEP::RandFlat::shoot(engine, HcalPulseShapes::Y11RANGE_);
541  if (CLHEP::RandFlat::shoot(engine, HcalPulseShapes::Y11MAX203_) < HcalPulseShapes::Y11203(result))
542  return result;
543  }
544 }
545 
546 double HcalPulseShapes::generatePhotonTime206(CLHEP::HepRandomEngine* engine) {
547  double result(0.);
548  while (true) {
549  result = CLHEP::RandFlat::shoot(engine, HcalPulseShapes::Y11RANGE_);
550  if (CLHEP::RandFlat::shoot(engine, HcalPulseShapes::Y11MAX206_) < HcalPulseShapes::Y11206(result))
551  return result;
552  }
553 }
554 
555 //Original scintillator+Y11 fit from Vasken's 2001 measurement
556 double HcalPulseShapes::Y11203(double t) { return exp(-0.0635 - 0.1518 * t + log(t) * 2.528) / 2485.9; }
557 
558 //New scintillator+Y11 model from Vasken's 2017 measurement plus a Landau correction term
559 double HcalPulseShapes::Y11206(double t) {
560  //Shifting phase to have better comparison of digi shape with data
561  //If necessary, further digi phase adjustment can be done here:
562  //SimCalorimetry/HcalSimProducers/python/hcalSimParameters_cfi.py
563  //by changing "timePhase"
564  double shift = 7.2;
565 
566  //Fit From Deconvolved Data
567  double A, n, t0, fit;
568  A = 0.104204;
569  n = 0.44064;
570  t0 = 10.0186;
571  if (t > shift)
572  fit = A * (1 - exp(-(t - shift) / n)) * exp(-(t - shift) / t0);
573  else
574  fit = 0.0;
575 
576  //Correction Term
577  double norm, mpv, sigma, corTerm;
578  norm = 0.0809882;
579  mpv = 0;
580  sigma = 20;
581  if (t > shift)
582  corTerm = norm * TMath::Landau((t - shift), mpv, sigma);
583  else
584  corTerm = 0.0;
585 
586  //Overall Y11
587  double frac = 0.11;
588  double val = (1 - frac) * fit + frac * corTerm;
589 
590  if (val >= 0)
591  return val;
592  else
593  return 0.0;
594 }
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:486
HcalPulseShapes::theShapes
ShapeMap theShapes
Definition: HcalPulseShapes.h:116
HcalPulseShapes::~HcalPulseShapes
~HcalPulseShapes()
Definition: HcalPulseShapes.cc:138
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:112
HcalPulseShapes::nBinsSiPM_
static const int nBinsSiPM_
Definition: HcalPulseShapes.h:42
HcalPulseShapes::shapeForReco
const Shape & shapeForReco(const HcalDetId &detId) const
Definition: HcalPulseShapes.cc:472
L1TowerCalibrationProducer_cfi.fit
fit
Definition: L1TowerCalibrationProducer_cfi.py:36
indexGen.s2
s2
Definition: indexGen.py:107
HcalPulseShapes::computeHFShape
void computeHFShape()
Definition: HcalPulseShapes.cc:245
HcalPulseShapes::computeSiPMShapeData2017
void computeSiPMShapeData2017()
Definition: HcalPulseShapes.cc:328
HcalBarrel
Definition: HcalAssistant.h:33
HcalPulseShapes::computeSiPMShapeHE206
const HcalPulseShape & computeSiPMShapeHE206()
Definition: HcalPulseShapes.cc:440
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:109
visualization-live-secondInstance_cfg.m
m
Definition: visualization-live-secondInstance_cfg.py:72
HcalPulseShapes::generatePhotonTime203
static double generatePhotonTime203(CLHEP::HepRandomEngine *engine)
Definition: HcalPulseShapes.cc:537
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
HcalPulseShapes::computeSiPMShapeHO
void computeSiPMShapeHO()
Definition: HcalPulseShapes.cc:382
HcalPulseShapes::Y11206
static double Y11206(double t)
Definition: HcalPulseShapes.cc:559
HcalPulseShapes::getShape
const Shape & getShape(int shapeType) const
Definition: HcalPulseShapes.cc:448
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:507
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:556
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:110
edm::LogWarning
Definition: MessageLogger.h:141
HcalPulseShapes::hpdBV30Shape_v2
Shape hpdBV30Shape_v2
Definition: HcalPulseShapes.h:113
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:512
HcalPulseShapes::hpdBV30ShapeMC_v2
Shape hpdBV30ShapeMC_v2
Definition: HcalPulseShapes.h:113
A
HcalPulseShapes::analyticPulseShapeSiPMHE
static double analyticPulseShapeSiPMHE(double t)
Definition: HcalPulseShapes.cc:523
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:140
alignmentValidation.c1
c1
do drawing
Definition: alignmentValidation.py:1025
p4
double p4[4]
Definition: TauolaWrapper.h:92
HcalPulseShapes::computeSiPMShapeData2018
void computeSiPMShapeData2018()
Definition: HcalPulseShapes.cc:275
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
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:530
HcalPulseShapes::hfShape
const Shape & hfShape() const
Definition: HcalPulseShapes.h:32
HcalPulseShapes::hpdShape_v2
Shape hpdShape_v2
Definition: HcalPulseShapes.h:111
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:458
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:149
HcalPulseShapes::siPMShapeHO_
Shape siPMShapeHO_
Definition: HcalPulseShapes.h:109
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:111
Exception
Definition: hltDiff.cc:246
HcalPulseShapes::theDbService
const HcalDbService * theDbService
Definition: HcalPulseShapes.h:114
ZMuMuCategoriesSequences_cff.nbin
nbin
Definition: ZMuMuCategoriesSequences_cff.py:25
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:546
HcalPulseShapes::computeSiPMShapeHE203
const HcalPulseShape & computeSiPMShapeHE203()
Definition: HcalPulseShapes.cc:433
HcalPulseShapes::hpdShape_v3
Shape hpdShape_v3
Definition: HcalPulseShapes.h:112
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:110
fit
Definition: CombinedChiSquaredLikelihood.h:6
HcalPulseShapes::analyticPulseShapeSiPMHO
static double analyticPulseShapeSiPMHO(double t)
Definition: HcalPulseShapes.cc:516
HcalPulseShapes::hfShape_
Shape hfShape_
Definition: HcalPulseShapes.h:109
HcalPulseShape::setShapeBin
void setShapeBin(int i, float f)
Definition: HcalPulseShape.cc:13
Generated for CMSSW Reference Manual by doxygen 1.8.16