File indexing completed on 2024-04-06 11:58:08
0001 #include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShapes.h"
0002 #include "FWCore/Utilities/interface/Exception.h"
0003 #include "FWCore/Framework/interface/ConsumesCollector.h"
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "CLHEP/Random/RandFlat.h"
0007
0008
0009 #include <cmath>
0010 #include <iostream>
0011 #include <fstream>
0012 #include "TMath.h"
0013
0014 HcalPulseShapes::HcalPulseShapes() : theDbService(nullptr), theShapes() {
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036 float ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3;
0037
0038
0039 ts1 = 8.;
0040 ts2 = 10.;
0041 ts3 = 29.3;
0042 thpd = 4.0;
0043 tpre = 9.0;
0044 wd1 = 2.0;
0045 wd2 = 0.7;
0046 wd3 = 1.0;
0047 computeHPDShape(ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3, hpdShape_);
0048 theShapes[101] = &hpdShape_;
0049 theShapes[102] = theShapes[101];
0050
0051
0052 ts1 = 8.;
0053 ts2 = 10.;
0054 ts3 = 25.0;
0055 thpd = 4.0;
0056 tpre = 9.0;
0057 wd1 = 2.0;
0058 wd2 = 0.7;
0059 wd3 = 1.0;
0060 computeHPDShape(ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3, hpdShape_v2);
0061 theShapes[103] = &hpdShape_v2;
0062
0063 ts1 = 8.;
0064 ts2 = 10.;
0065 ts3 = 29.3;
0066 thpd = 4.0;
0067 tpre = 7.0;
0068 wd1 = 2.0;
0069 wd2 = 0.7;
0070 wd3 = 1.0;
0071 computeHPDShape(ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3, hpdShapeMC_v2);
0072 theShapes[123] = &hpdShapeMC_v2;
0073
0074
0075 ts1 = 8.;
0076 ts2 = 19.;
0077 ts3 = 29.3;
0078 thpd = 4.0;
0079 tpre = 9.0;
0080 wd1 = 2.0;
0081 wd2 = 0.7;
0082 wd3 = 0.32;
0083 computeHPDShape(ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3, hpdShape_v3);
0084 theShapes[105] = &hpdShape_v3;
0085
0086 ts1 = 8.;
0087 ts2 = 10.;
0088 ts3 = 22.3;
0089 thpd = 4.0;
0090 tpre = 7.0;
0091 wd1 = 2.0;
0092 wd2 = 0.7;
0093 wd3 = 1.0;
0094 computeHPDShape(ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3, hpdShapeMC_v3);
0095 theShapes[125] = &hpdShapeMC_v3;
0096
0097
0098
0099 ts1 = 8.;
0100 ts2 = 12.;
0101 ts3 = 31.7;
0102 thpd = 9.0;
0103 tpre = 9.0;
0104 wd1 = 2.0;
0105 wd2 = 0.7;
0106 wd3 = 1.0;
0107 computeHPDShape(ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3, hpdBV30Shape_v2);
0108 theShapes[104] = &hpdBV30Shape_v2;
0109
0110 ts1 = 8.;
0111 ts2 = 12.;
0112 ts3 = 31.7;
0113 thpd = 9.0;
0114 tpre = 9.0;
0115 wd1 = 2.0;
0116 wd2 = 0.7;
0117 wd3 = 1.0;
0118 computeHPDShape(ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3, hpdBV30ShapeMC_v2);
0119 theShapes[124] = &hpdBV30ShapeMC_v2;
0120
0121
0122
0123 computeHFShape();
0124 computeSiPMShapeHO();
0125 computeSiPMShapeData2017();
0126 computeSiPMShapeData2018();
0127 computeSiPMShapeMCRecoRun3();
0128
0129 theShapes[201] = &siPMShapeHO_;
0130 theShapes[202] = theShapes[201];
0131 theShapes[203] = &(computeSiPMShapeHE203());
0132 theShapes[205] = &siPMShapeData2017_;
0133 theShapes[206] = &(computeSiPMShapeHE206());
0134 theShapes[207] = &siPMShapeData2018_;
0135 theShapes[208] = &siPMShapeMCRecoRun3_;
0136 theShapes[301] = &hfShape_;
0137
0138 }
0139
0140 HcalPulseShapes::HcalPulseShapes(edm::ConsumesCollector iC) : HcalPulseShapes() {
0141 theDbServiceToken = iC.esConsumes<edm::Transition::BeginRun>();
0142 }
0143
0144 HcalPulseShapes::~HcalPulseShapes() {}
0145
0146 void HcalPulseShapes::beginRun(edm::EventSetup const& es) { theDbService = &es.getData(theDbServiceToken); }
0147
0148 void HcalPulseShapes::beginRun(const HcalDbService* conditions) { theDbService = conditions; }
0149
0150
0151 void HcalPulseShapes::computeHPDShape(
0152 float ts1, float ts2, float ts3, float thpd, float tpre, float wd1, float wd2, float wd3, Shape& tmphpdShape_) {
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166 unsigned int nbin = 256;
0167 tmphpdShape_.setNBin(nbin);
0168 std::vector<float> ntmp(nbin, 0.0);
0169 std::vector<float> nth(nbin, 0.0);
0170 std::vector<float> ntp(nbin, 0.0);
0171 std::vector<float> ntd(nbin, 0.0);
0172
0173 unsigned int i, j, k;
0174 float norm;
0175
0176
0177 norm = 0.0;
0178 for (j = 0; j < thpd && j < nbin; j++) {
0179 nth[j] = 1.0 + ((float)j) / thpd;
0180 norm += nth[j];
0181 }
0182
0183 for (j = 0; j < thpd && j < nbin; j++) {
0184 nth[j] /= norm;
0185 }
0186
0187
0188 norm = 0.0;
0189 for (j = 0; j < 6 * tpre && j < nbin; j++) {
0190 ntp[j] = ((float)j) * exp(-((float)(j * j)) / (tpre * tpre));
0191 norm += ntp[j];
0192 }
0193
0194 for (j = 0; j < 6 * tpre && j < nbin; j++) {
0195 ntp[j] /= norm;
0196 }
0197
0198
0199
0200
0201
0202 unsigned int tmax = 6 * (int)ts3;
0203
0204 norm = 0.0;
0205 for (j = 0; j < tmax && j < nbin; j++) {
0206 ntd[j] = wd1 * exp(-((float)j) / ts1) + wd2 * exp(-((float)j) / ts2) + wd3 * exp(-((float)j) / ts3);
0207 norm += ntd[j];
0208 }
0209
0210 for (j = 0; j < tmax && j < nbin; j++) {
0211 ntd[j] /= norm;
0212 }
0213
0214 unsigned int t1, t2, t3, t4;
0215 for (i = 0; i < tmax && i < nbin; i++) {
0216 t1 = i;
0217
0218
0219 t2 = t1;
0220 for (j = 0; j < thpd && j < nbin; j++) {
0221 t3 = t2 + j;
0222 for (k = 0; k < 4 * tpre && k < nbin; k++) {
0223 t4 = t3 + k;
0224 if (t4 < nbin) {
0225 unsigned int ntb = t4;
0226 ntmp[ntb] += ntd[i] * nth[j] * ntp[k];
0227 }
0228 }
0229 }
0230 }
0231
0232
0233 norm = 0.;
0234 for (i = 0; i < nbin; i++) {
0235 norm += ntmp[i];
0236 }
0237
0238 for (i = 0; i < nbin; i++) {
0239 ntmp[i] /= norm;
0240 }
0241
0242 for (i = 0; i < nbin; i++) {
0243 tmphpdShape_.setShapeBin(i, ntmp[i]);
0244 }
0245 }
0246
0247 void HcalPulseShapes::computeHFShape() {
0248
0249 unsigned int nbin = 256;
0250 hfShape_.setNBin(nbin);
0251 std::vector<float> ntmp(nbin, 0.0);
0252
0253 const float k0 = 0.7956;
0254 const float p2 = 1.355;
0255 const float p4 = 2.327;
0256 const float p1 = 4.3;
0257
0258 float norm = 0.0;
0259
0260 for (unsigned int j = 0; j < 25 && j < nbin; ++j) {
0261 float r0 = j - p1;
0262 float sigma0 = (r0 < 0) ? p2 : p2 * p4;
0263 r0 /= sigma0;
0264 if (r0 < k0)
0265 ntmp[j] = exp(-0.5 * r0 * r0);
0266 else
0267 ntmp[j] = exp(0.5 * k0 * k0 - k0 * r0);
0268 norm += ntmp[j];
0269 }
0270
0271 for (unsigned int j = 0; j < 25 && j < nbin; ++j) {
0272 ntmp[j] /= norm;
0273 hfShape_.setShapeBin(j, ntmp[j]);
0274 }
0275 }
0276 void HcalPulseShapes::computeSiPMShapeMCRecoRun3() {
0277
0278
0279 unsigned int nbin = 250;
0280 std::array<float, 250> nt{
0281 {0, 0, 0, 0, 0, 0, 0, 0.000117468,
0282 0.0031549, 0.0117368, 0.0219974, 0.0305776, 0.0365429, 0.0400524, 0.0415915, 0.0416765,
0283 0.0408111, 0.0394627, 0.0379353, 0.0363688, 0.0348152, 0.0332891, 0.0317923, 0.0303237,
0284 0.028883, 0.0274714, 0.0260914, 0.0247462, 0.0234392, 0.0221738, 0.0209531, 0.0197793,
0285 0.0186544, 0.0175796, 0.0165556, 0.0155823, 0.0146596, 0.0137866, 0.0129623, 0.0121853,
0286 0.0114539, 0.0107665, 0.0101213, 0.0095162, 0.00894934, 0.00841873, 0.0079224, 0.00745841,
0287 0.00702487, 0.00661995, 0.00624189, 0.00588898, 0.00555961, 0.00525223, 0.00496539, 0.0046977,
0288 0.00444786, 0.00421464, 0.00399689, 0.00379353, 0.00360355, 0.00342602, 0.00326004, 0.0031048,
0289 0.00295954, 0.00282355, 0.00269616, 0.00257676, 0.00246479, 0.00235972, 0.00226106, 0.00216834,
0290 0.00208117, 0.00199914, 0.00192189, 0.0018491, 0.00178044, 0.00171565, 0.00165445, 0.00159659,
0291 0.00154186, 0.00149003, 0.00144092, 0.00139435, 0.00135015, 0.00130816, 0.00126825, 0.00123027,
0292 0.00119412, 0.00115966, 0.0011268, 0.00109544, 0.00106548, 0.00103685, 0.00100946, 0.000983242,
0293 0.000958125, 0.000934047, 0.000910949, 0.000888775, 0.000867475, 0.000847, 0.000827306, 0.000808352,
0294 0.000790097, 0.000772506, 0.000755545, 0.000739182, 0.000723387, 0.000708132, 0.00069339, 0.000679138,
0295 0.000665352, 0.00065201, 0.000639091, 0.000626577, 0.00061445, 0.000602692, 0.000591287, 0.00058022,
0296 0.000569477, 0.000559044, 0.000548908, 0.000539058, 0.000529481, 0.000520167, 0.000511106, 0.000502288,
0297 0.000493704, 0.000485344, 0.000477201, 0.000469266, 0.000459912, 0.000448544, 0.000437961, 0.000428079,
0298 0.000418825, 0.000410133, 0.000401945, 0.00039421, 0.000386883, 0.000379924, 0.000373298, 0.000366973,
0299 0.000360922, 0.00035512, 0.000349545, 0.000344179, 0.000339003, 0.000334002, 0.000329163, 0.000324475,
0300 0.000319925, 0.000315504, 0.000311204, 0.000307017, 0.000302935, 0.000298954, 0.000295066, 0.000291267,
0301 0.000287553, 0.000283919, 0.000280361, 0.000276877, 0.000273462, 0.000270114, 0.000266831, 0.000263609,
0302 0.000260447, 0.000257343, 0.000254295, 0.0002513, 0.000248358, 0.000245467, 0.000242625, 0.000239831,
0303 0.000237083, 0.000234381, 0.000231723, 0.000229109, 0.000226536, 0.000224005, 0.000221514, 0.000219062,
0304 0.000216648, 0.000214272, 0.000211933, 0.00020963, 0.000207362, 0.000205129, 0.000202929, 0.000200763,
0305 0.000198629, 0.000196526, 0.000194455, 0.000192415, 0.000190405, 0.000188424, 0.000186472, 0.000184548,
0306 0.000182653, 0.000180784, 0.000178943, 0.000177127, 0.000175338, 0.000173574, 0.000171835, 0.00017012,
0307 0.000168429, 0.000166762, 0.000165119, 0.000163498, 0.000161899, 0.000160322, 0.000158767, 0.000157233,
0308 0.000155721, 0.000154228, 0.000152756, 0.000151304, 0.000149871, 0.000148457, 0.000147062, 0.000145686,
0309 0.000144327, 0.000142987, 0.000141664, 0.000140359, 0.000139071, 0.000137799, 0.000136544, 0.000135305,
0310 0.000134082, 0.000132874, 0.000131682, 0.000130505, 0.000129344, 0.000128196, 0.000127064, 0.000125945,
0311 0.00012484, 0.00012375, 0.000122672, 0.000121608, 0.000120558, 0.00011952, 0.000118495, 0.000117482,
0312 0.000116482, 0.000115493}};
0313
0314 siPMShapeMCRecoRun3_.setNBin(nbin);
0315
0316 double norm = 0.;
0317 for (unsigned int j = 0; j < nbin; ++j) {
0318 norm += (nt[j] > 0) ? nt[j] : 0.;
0319 }
0320
0321 for (unsigned int j = 0; j < nbin; ++j) {
0322 nt[j] /= norm;
0323 siPMShapeMCRecoRun3_.setShapeBin(j, nt[j]);
0324 }
0325 }
0326
0327 void HcalPulseShapes::computeSiPMShapeData2018() {
0328
0329
0330
0331 unsigned int nbin = 250;
0332
0333 std::array<float, 250> nt{
0334 {5.22174e-12, 7.04852e-10, 3.49584e-08, 7.78029e-07, 9.11847e-06, 6.39666e-05, 0.000297587, 0.000996661,
0335 0.00256618, 0.00535396, 0.00944073, 0.0145521, 0.020145, 0.0255936, 0.0303632, 0.0341078,
0336 0.0366849, 0.0381183, 0.0385392, 0.0381327, 0.0370956, 0.0356113, 0.0338366, 0.0318978,
0337 0.029891, 0.0278866, 0.0259336, 0.0240643, 0.0222981, 0.0206453, 0.0191097, 0.0176902,
0338 0.0163832, 0.0151829, 0.0140826, 0.0130752, 0.0121533, 0.01131, 0.0105382, 0.00983178,
0339 0.00918467, 0.00859143, 0.00804709, 0.0075471, 0.00708733, 0.00666406, 0.00627393, 0.00591389,
0340 0.00558122, 0.00527344, 0.00498834, 0.00472392, 0.00447837, 0.00425007, 0.00403754, 0.00383947,
0341 0.00365465, 0.00348199, 0.00332052, 0.00316934, 0.00302764, 0.0028947, 0.00276983, 0.00265242,
0342 0.00254193, 0.00243785, 0.00233971, 0.00224709, 0.0021596, 0.00207687, 0.0019986, 0.00192447,
0343 0.00185421, 0.00178756, 0.0017243, 0.00166419, 0.00160705, 0.00155268, 0.00150093, 0.00145162,
0344 0.00140461, 0.00135976, 0.00131696, 0.00127607, 0.00123699, 0.00119962, 0.00116386, 0.00112963,
0345 0.00109683, 0.0010654, 0.00103526, 0.00100634, 0.000978578, 0.000951917, 0.000926299, 0.000901672,
0346 0.000877987, 0.000855198, 0.00083326, 0.000812133, 0.000791778, 0.000772159, 0.000753242, 0.000734994,
0347 0.000717384, 0.000700385, 0.000683967, 0.000668107, 0.000652779, 0.00063796, 0.000623629, 0.000609764,
0348 0.000596346, 0.000583356, 0.000570777, 0.000558592, 0.000546785, 0.00053534, 0.000524243, 0.000513481,
0349 0.00050304, 0.000492907, 0.000483072, 0.000473523, 0.000464248, 0.000455238, 0.000446483, 0.000437974,
0350 0.0004297, 0.000421655, 0.00041383, 0.000406216, 0.000398807, 0.000391595, 0.000384574, 0.000377736,
0351 0.000371076, 0.000364588, 0.000358266, 0.000352104, 0.000346097, 0.00034024, 0.000334528, 0.000328956,
0352 0.00032352, 0.000318216, 0.000313039, 0.000307986, 0.000303052, 0.000298234, 0.000293528, 0.000288931,
0353 0.000284439, 0.00028005, 0.000275761, 0.000271567, 0.000267468, 0.000263459, 0.000259538, 0.000255703,
0354 0.000251951, 0.00024828, 0.000244688, 0.000241172, 0.00023773, 0.000234361, 0.000231061, 0.00022783,
0355 0.000224666, 0.000221566, 0.000218528, 0.000215553, 0.000212636, 0.000209778, 0.000206977, 0.00020423,
0356 0.000201537, 0.000198896, 0.000196307, 0.000193767, 0.000191275, 0.000188831, 0.000186432, 0.000184079,
0357 0.000181769, 0.000179502, 0.000177277, 0.000175092, 0.000172947, 0.000170841, 0.000168772, 0.000166741,
0358 0.000164745, 0.000162785, 0.000160859, 0.000158967, 0.000157108, 0.00015528, 0.000153484, 0.000151719,
0359 0.000149984, 0.000148278, 0.000146601, 0.000144951, 0.000143329, 0.000141734, 0.000140165, 0.000138622,
0360 0.000137104, 0.00013561, 0.000134141, 0.000132695, 0.000131272, 0.000129871, 0.000128493, 0.000127136,
0361 0.000125801, 0.000124486, 0.000123191, 0.000121917, 0.000120662, 0.000119426, 0.000118209, 0.00011701,
0362 0.000115829, 0.000114665, 0.000113519, 0.00011239, 0.000111278, 0.000110182, 0.000109102, 0.000108037,
0363 0.000106988, 0.000105954, 0.000104935, 0.00010393, 0.000102939, 0.000101963, 0.000101, 0.000100051,
0364 9.91146e-05, 9.81915e-05, 9.7281e-05, 9.63831e-05, 9.54975e-05, 9.46239e-05, 9.37621e-05, 9.2912e-05,
0365 9.20733e-05, 9.12458e-05}};
0366
0367 siPMShapeData2018_.setNBin(nbin);
0368
0369 double norm = 0.;
0370 for (unsigned int j = 0; j < nbin; ++j) {
0371 norm += (nt[j] > 0) ? nt[j] : 0.;
0372 }
0373
0374 for (unsigned int j = 0; j < nbin; ++j) {
0375 nt[j] /= norm;
0376 siPMShapeData2018_.setShapeBin(j, nt[j]);
0377 }
0378 }
0379
0380 void HcalPulseShapes::computeSiPMShapeData2017() {
0381
0382
0383
0384
0385 unsigned int nbin = 250;
0386
0387 std::array<float, 250> nt{
0388 {3.97958e-29, 1.11634e-22, 9.96106e-18, 6.25334e-14, 5.08863e-11, 8.59141e-09, 4.32285e-07, 8.56617e-06,
0389 8.28549e-05, 0.000461447, 0.00168052, 0.00441395, 0.00901637, 0.0151806, 0.0220314, 0.028528,
0390 0.0338471, 0.0375578, 0.0395985, 0.0401567, 0.0395398, 0.0380776, 0.0360669, 0.0337474,
0391 0.0312984, 0.0288457, 0.0264721, 0.0242276, 0.0221393, 0.0202181, 0.0184647, 0.0168731,
0392 0.0154335, 0.0141346, 0.0129639, 0.0119094, 0.0109594, 0.0101031, 0.0093305, 0.00863267,
0393 0.0080015, 0.00742977, 0.00691107, 0.00643969, 0.00601059, 0.00561931, 0.00526188, 0.00493483,
0394 0.00463505, 0.00435981, 0.00410667, 0.00387348, 0.00365832, 0.00345949, 0.00327547, 0.0031049,
0395 0.00294656, 0.00279938, 0.00266237, 0.00253467, 0.00241548, 0.0023041, 0.00219989, 0.00210227,
0396 0.00201072, 0.00192476, 0.00184397, 0.00176795, 0.00169634, 0.00162884, 0.00156512, 0.00150494,
0397 0.00144803, 0.00139418, 0.00134317, 0.00129481, 0.00124894, 0.00120537, 0.00116398, 0.00112461,
0398 0.00108715, 0.00105147, 0.00101747, 0.000985042, 0.000954096, 0.000924545, 0.000896308, 0.000869311,
0399 0.000843482, 0.000818758, 0.000795077, 0.000772383, 0.000750623, 0.000729747, 0.00070971, 0.000690466,
0400 0.000671977, 0.000654204, 0.00063711, 0.000620663, 0.000604831, 0.000589584, 0.000574894, 0.000560735,
0401 0.000547081, 0.00053391, 0.0005212, 0.000508929, 0.000497078, 0.000485628, 0.000474561, 0.000463862,
0402 0.000453514, 0.000443501, 0.000433811, 0.000424429, 0.000415343, 0.00040654, 0.00039801, 0.000389741,
0403 0.000381722, 0.000373944, 0.000366398, 0.000359074, 0.000351964, 0.00034506, 0.000338353, 0.000331838,
0404 0.000325505, 0.00031935, 0.000313365, 0.000307544, 0.000301881, 0.000296371, 0.000291009, 0.000285788,
0405 0.000280705, 0.000275755, 0.000270932, 0.000266233, 0.000261653, 0.00025719, 0.000252837, 0.000248593,
0406 0.000244454, 0.000240416, 0.000236475, 0.00023263, 0.000228876, 0.000225212, 0.000221633, 0.000218138,
0407 0.000214724, 0.000211389, 0.00020813, 0.000204945, 0.000201831, 0.000198787, 0.000195811, 0.0001929,
0408 0.000190053, 0.000187268, 0.000184543, 0.000181876, 0.000179266, 0.000176711, 0.00017421, 0.000171761,
0409 0.000169363, 0.000167014, 0.000164713, 0.000162459, 0.00016025, 0.000158086, 0.000155964, 0.000153885,
0410 0.000151847, 0.000149848, 0.000147888, 0.000145966, 0.000144081, 0.000142232, 0.000140418, 0.000138638,
0411 0.000136891, 0.000135177, 0.000133494, 0.000131843, 0.000130221, 0.00012863, 0.000127066, 0.000125531,
0412 0.000124023, 0.000122543, 0.000121088, 0.000119658, 0.000118254, 0.000116874, 0.000115518, 0.000114185,
0413 0.000112875, 0.000111587, 0.000110321, 0.000109076, 0.000107851, 0.000106648, 0.000105464, 0.000104299,
0414 0.000103154, 0.000102027, 0.000100918, 9.98271e-05, 9.87537e-05, 9.76974e-05, 9.66578e-05, 9.56346e-05,
0415 9.46274e-05, 9.3636e-05, 9.26599e-05, 9.16989e-05, 9.07526e-05, 8.98208e-05, 8.89032e-05, 8.79995e-05,
0416 8.71093e-05, 8.62325e-05, 8.53688e-05, 8.45179e-05, 8.36796e-05, 8.28536e-05, 8.20397e-05, 8.12376e-05,
0417 8.04471e-05, 7.96681e-05, 7.89002e-05, 7.81433e-05, 7.73972e-05, 7.66616e-05, 7.59364e-05, 7.52213e-05,
0418 7.45163e-05, 7.3821e-05, 7.31354e-05, 7.24592e-05, 7.17923e-05, 7.11345e-05, 7.04856e-05, 6.98455e-05,
0419 6.9214e-05, 6.8591e-05}};
0420
0421 siPMShapeData2017_.setNBin(nbin);
0422
0423 double norm = 0.;
0424 for (unsigned int j = 0; j < nbin; ++j) {
0425 norm += (nt[j] > 0) ? nt[j] : 0.;
0426 }
0427
0428 for (unsigned int j = 0; j < nbin; ++j) {
0429 nt[j] /= norm;
0430 siPMShapeData2017_.setShapeBin(j, nt[j]);
0431 }
0432 }
0433
0434 void HcalPulseShapes::computeSiPMShapeHO() {
0435 unsigned int nbin = 128;
0436
0437
0438 std::array<float, 128> nt{
0439 {2.782980485851731e-6, 4.518134885954626e-5, 2.7689305197392056e-4, 9.18328418900969e-4,
0440 .002110072599166349, .003867856860331454, .006120046224897771, .008754774090536956,
0441 0.0116469503358586, .01467007449455966, .01770489955229477, .02064621450689512,
0442 .02340678093764222, .02591874610854916, .02813325527435303, 0.0300189241965647,
0443 .03155968107671164, .03275234052577155, .03360415306318798, .03413048377960748,
0444 .03435270899678218, .03429637464659661, .03398962975487166, .03346192884394954,
0445 .03274298516247742, .03186195009136525, .03084679116113031, 0.0297238406141036,
0446 .02851748748929785, .02724998816332392, .02594137274487424, .02460942736731527,
0447 .02326973510736116, .02193576080366117, 0.0206189674254987, .01932895378564653,
0448 0.0180736052958666, .01685925112650875, 0.0156908225633535, .01457200857138456,
0449 .01350540559602467, .01249265947824805, .01153459805300423, .01063135355597282,
0450 .009782474412011936, .008987026319784546, 0.00824368281357106, .007550805679909604,
0451 .006906515742762193, .006308754629755056, .005755338185695127, .005244002229973356,
0452 .004772441359900532, .004338341490928299, .003939406800854143, 0.00357338171220501,
0453 0.0032380685079891, .002931341133259233, .002651155690306086, .002395558090237333,
0454 .002162689279320922, .001950788415487319, .001758194329648101, .001583345567913682,
0455 .001424779275191974, .001281129147671334, 0.00115112265163774, .001033577678808199,
0456 9.273987838127585e-4, 8.315731274976846e-4, 7.451662302008696e-4, 6.673176219006913e-4,
0457 5.972364609644049e-4, 5.341971801529036e-4, 4.775352065178378e-4, 4.266427928961177e-4,
0458 3.8096498904225923e-4, 3.3999577417327287e-4, 3.032743659102713e-4, 2.703817158798329e-4,
0459 2.4093719775272793e-4, 2.145954900503894e-4, 1.9104365317752797e-4, 1.6999839784346724e-4,
0460 1.5120354022478893e-4, 1.3442763782650755e-4, 1.1946179895521507e-4, 1.0611765796993575e-4,
0461 9.422550797617687e-5, 8.363258233342666e-5, 7.420147621931836e-5, 6.580869950304933e-5,
0462 5.834335229919868e-5, 5.17059147771959e-5, 4.5807143072062634e-5, 4.0567063461299446e-5,
0463 3.591405732740723e-5, 3.178402980354131e-5, 2.811965539165646e-5, 2.4869694240316126e-5,
0464 2.1988373166730962e-5, 1.9434825899529382e-5, 1.717258740121378e-5, 1.5169137499243157e-5,
0465 1.339548941011129e-5, 1.1825819079078403e-5, 1.0437131581057595e-5, 9.208961130078894e-6,
0466 8.12310153137994e-6, 7.163364176588591e-6, 6.315360932244386e-6, 5.566309502463164e-6,
0467 4.904859063429651e-6, 4.320934164082596e-6, 3.8055950719111903e-6, 3.350912911083174e-6,
0468 2.9498580949517117e-6, 2.596200697612328e-6, 2.2844215378879293e-6, 2.0096328693141094e-6,
0469 1.7675076766686654e-6, 1.5542166787225756e-6, 1.366372225473431e-6, 1.200978365778838e-6,
0470 1.0553864128982371e-6, 9.272554464808518e-7, 8.145171945902259e-7, 7.153448381918271e-7}};
0471
0472 siPMShapeHO_.setNBin(nbin);
0473
0474 double norm = 0.;
0475 for (unsigned int j = 0; j < nbin; ++j) {
0476 norm += (nt[j] > 0) ? nt[j] : 0.;
0477 }
0478
0479 for (unsigned int j = 0; j < nbin; ++j) {
0480 nt[j] /= norm;
0481 siPMShapeHO_.setShapeBin(j, nt[j]);
0482 }
0483 }
0484
0485 const HcalPulseShape& HcalPulseShapes::computeSiPMShapeHE203() {
0486
0487 static const HcalPulseShape siPMShapeMC2017(
0488 normalize(convolve(nBinsSiPM_, analyticPulseShapeSiPMHE, Y11203), nBinsSiPM_), nBinsSiPM_);
0489 return siPMShapeMC2017;
0490 }
0491
0492 const HcalPulseShape& HcalPulseShapes::computeSiPMShapeHE206() {
0493
0494
0495 static const HcalPulseShape siPMShapeMC2018(
0496 normalizeShift(convolve(nBinsSiPM_, analyticPulseShapeSiPMHE, Y11206), nBinsSiPM_, -2), nBinsSiPM_);
0497 return siPMShapeMC2018;
0498 }
0499
0500 const HcalPulseShapes::Shape& HcalPulseShapes::getShape(int shapeType) const {
0501 ShapeMap::const_iterator shapeMapItr = theShapes.find(shapeType);
0502 if (shapeMapItr == theShapes.end()) {
0503 throw cms::Exception("HcalPulseShapes") << "unknown shapeType";
0504 return hpdShape_;
0505 } else {
0506 return *(shapeMapItr->second);
0507 }
0508 }
0509
0510 const HcalPulseShapes::Shape& HcalPulseShapes::shape(const HcalDetId& detId) const {
0511 if (!theDbService) {
0512 return defaultShape(detId);
0513 }
0514 int shapeType = theDbService->getHcalMCParam(detId)->signalShape();
0515
0516 ShapeMap::const_iterator shapeMapItr = theShapes.find(shapeType);
0517 if (shapeMapItr == theShapes.end()) {
0518 return defaultShape(detId);
0519 } else {
0520 return *(shapeMapItr->second);
0521 }
0522 }
0523
0524 const HcalPulseShapes::Shape& HcalPulseShapes::shapeForReco(const HcalDetId& detId) const {
0525 if (!theDbService) {
0526 return defaultShape(detId);
0527 }
0528 int shapeType = theDbService->getHcalRecoParam(detId.rawId())->pulseShapeID();
0529
0530 ShapeMap::const_iterator shapeMapItr = theShapes.find(shapeType);
0531 if (shapeMapItr == theShapes.end()) {
0532 return defaultShape(detId);
0533 } else {
0534 return *(shapeMapItr->second);
0535 }
0536 }
0537
0538 const HcalPulseShapes::Shape& HcalPulseShapes::defaultShape(const HcalDetId& detId) const {
0539 edm::LogWarning("HcalPulseShapes") << "Cannot find HCAL MC Params ";
0540 HcalSubdetector subdet = detId.subdet();
0541 switch (subdet) {
0542 case HcalBarrel:
0543 return hbShape();
0544 case HcalEndcap:
0545 return heShape();
0546 case HcalForward:
0547 return hfShape();
0548 case HcalOuter:
0549
0550 return hoShape(false);
0551 default:
0552 throw cms::Exception("HcalPulseShapes") << "unknown detId";
0553 break;
0554 }
0555 }
0556
0557
0558
0559 inline double gexp(double t, double A, double c, double t0, double s) {
0560 static double const root2(sqrt(2));
0561 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);
0562 }
0563
0564 inline double onePulse(double t, double A, double sigma, double theta, double m) {
0565 return (t < theta) ? 0 : A * TMath::LogNormal(t, sigma, theta, m);
0566 }
0567
0568 double HcalPulseShapes::analyticPulseShapeSiPMHO(double t) {
0569
0570 double A1(0.08757), c1(-0.5257), t01(2.4013), s1(0.6721);
0571 double A2(0.007598), c2(-0.1501), t02(6.9412), s2(0.8710);
0572 return gexp(t, A1, c1, t01, s1) + gexp(t, A2, c2, t02, s2);
0573 }
0574
0575 double HcalPulseShapes::analyticPulseShapeSiPMHE(double t) {
0576
0577 double A1(5.204 / 6.94419), sigma1_shape(0.5387), theta1_loc(-0.3976), m1_scale(4.428);
0578 double A2(1.855 / 6.94419), sigma2_shape(0.8132), theta2_loc(7.025), m2_scale(12.29);
0579 return onePulse(t, A1, sigma1_shape, theta1_loc, m1_scale) + onePulse(t, A2, sigma2_shape, theta2_loc, m2_scale);
0580 }
0581
0582 double HcalPulseShapes::generatePhotonTime(CLHEP::HepRandomEngine* engine, unsigned int signalShape) {
0583 if (signalShape == 206)
0584 return generatePhotonTime206(engine);
0585 else
0586 return generatePhotonTime203(engine);
0587 }
0588
0589 double HcalPulseShapes::generatePhotonTime203(CLHEP::HepRandomEngine* engine) {
0590 double result(0.);
0591 while (true) {
0592 result = CLHEP::RandFlat::shoot(engine, HcalPulseShapes::Y11RANGE_);
0593 if (CLHEP::RandFlat::shoot(engine, HcalPulseShapes::Y11MAX203_) < HcalPulseShapes::Y11203(result))
0594 return result;
0595 }
0596 }
0597
0598 double HcalPulseShapes::generatePhotonTime206(CLHEP::HepRandomEngine* engine) {
0599 double result(0.);
0600 while (true) {
0601 result = CLHEP::RandFlat::shoot(engine, HcalPulseShapes::Y11RANGE_);
0602 if (CLHEP::RandFlat::shoot(engine, HcalPulseShapes::Y11MAX206_) < HcalPulseShapes::Y11206(result))
0603 return result;
0604 }
0605 }
0606
0607
0608 double HcalPulseShapes::Y11203(double t) { return exp(-0.0635 - 0.1518 * t + log(t) * 2.528) / 2485.9; }
0609
0610
0611 double HcalPulseShapes::Y11206(double t) {
0612
0613
0614
0615
0616 double shift = 7.2;
0617
0618
0619 double A, n, t0, fit;
0620 A = 0.104204;
0621 n = 0.44064;
0622 t0 = 10.0186;
0623 if (t > shift)
0624 fit = A * (1 - exp(-(t - shift) / n)) * exp(-(t - shift) / t0);
0625 else
0626 fit = 0.0;
0627
0628
0629 double norm, mpv, sigma, corTerm;
0630 norm = 0.0809882;
0631 mpv = 0;
0632 sigma = 20;
0633 if (t > shift)
0634 corTerm = norm * TMath::Landau((t - shift), mpv, sigma);
0635 else
0636 corTerm = 0.0;
0637
0638
0639 double frac = 0.11;
0640 double val = (1 - frac) * fit + frac * corTerm;
0641
0642 if (val >= 0)
0643 return val;
0644 else
0645 return 0.0;
0646 }