Back to home page

Project CMSSW displayed by LXR

 
 

    


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 // #include "CalibCalorimetry/HcalAlgos/interface/HcalDbASCIIIO.h"
0009 #include <cmath>
0010 #include <iostream>
0011 #include <fstream>
0012 #include "TMath.h"
0013 
0014 HcalPulseShapes::HcalPulseShapes() : theDbService(nullptr), theShapes() {
0015   /*
0016 
0017 Reco  MC
0018 -------------------------------------------------------------------------------------------
0019 000                                              not used (reserved)
0020 101   101   hpdShape_                            HPD (original version)
0021 102   102   =101                                 HPD BV 30 volts in HBP iphi54
0022 103   123   hpdShape_v2,hpdShapeMC_v2            HPD (2011. oct version)
0023 104   124   hpdBV30Shape_v2,hpdBV30ShapeMC_v2    HPD bv30 in HBP iph54
0024 105   125   hpdShape_v2,hpdShapeMC_v2            HPD (2011.11.12 version)
0025 201   201   siPMShapeHO_                         SiPMs Zecotec shape   (HO)
0026 202   202   =201,                                SiPMs Hamamatsu shape (HO)
0027 205   203   siPMShapeData2017_,siPMShapeMC2017_  SiPMs from Data, Hamamatsu shape (HE 2017)
0028 207   206   siPMShapeData2018_,siPMShapeMC2018_  SiPMs from Data, Hamamatsu shape (HE 2018)
0029 207   208   siPMShapeData2018_,siPMShapeMCRecoRun3_  SiPMs from Data, 2021 MC phase scan
0030 301   301   hfShape_                             regular HF PMT shape
0031 401   401                                        regular ZDC shape
0032 -------------------------------------------------------------------------------------------
0033 
0034 */
0035 
0036   float ts1, ts2, ts3, thpd, tpre, wd1, wd2, wd3;
0037 
0038   //  HPD Shape  Version 1 (used before CMSSW5, until Oct 2011)
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   //  HPD Shape  Version 2 for CMSSW 5. Nov 2011  (RECO and MC separately)
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   //  HPD Shape  Version 3 for CMSSW 5. Nov 2011  (RECO and MC separately)
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   // HPD with Bias Voltage 30 volts, wider pulse.  (HBPlus iphi54)
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   // HF and SiPM
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   //theShapes[401] = new CaloCachedShapeIntegrator(&theZDCShape);
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 //void HcalPulseShapes::computeHPDShape()
0151 void HcalPulseShapes::computeHPDShape(
0152     float ts1, float ts2, float ts3, float thpd, float tpre, float wd1, float wd2, float wd3, Shape& tmphpdShape_) {
0153   // pulse shape time constants in ns
0154   /*
0155   const float ts1  = 8.;          // scintillation time constants : 1,2,3
0156   const float ts2  = 10.;           
0157   const float ts3  = 29.3;         
0158   const float thpd = 4.;          // HPD current collection drift time
0159   const float tpre = 9.;          // preamp time constant (refit on TB04 data)
0160   
0161   const float wd1 = 2.;           // relative weights of decay exponents 
0162   const float wd2 = 0.7;
0163   const float wd3 = 1.;
0164 */
0165   // pulse shape components over a range of time 0 ns to 255 ns in 1 ns steps
0166   unsigned int nbin = 256;
0167   tmphpdShape_.setNBin(nbin);
0168   std::vector<float> ntmp(nbin, 0.0);  // zeroing output pulse shape
0169   std::vector<float> nth(nbin, 0.0);   // zeroing HPD drift shape
0170   std::vector<float> ntp(nbin, 0.0);   // zeroing Binkley preamp shape
0171   std::vector<float> ntd(nbin, 0.0);   // zeroing Scintillator decay shape
0172 
0173   unsigned int i, j, k;
0174   float norm;
0175 
0176   // HPD starts at I and rises to 2I in thpd of time
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   // normalize integrated current to 1.0
0183   for (j = 0; j < thpd && j < nbin; j++) {
0184     nth[j] /= norm;
0185   }
0186 
0187   // Binkley shape over 6 time constants
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   // normalize pulse area to 1.0
0194   for (j = 0; j < 6 * tpre && j < nbin; j++) {
0195     ntp[j] /= norm;
0196   }
0197 
0198   // ignore stochastic variation of photoelectron emission
0199   // <...>
0200 
0201   // effective tile plus wave-length shifter decay time over 4 time constants
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   // normalize pulse area to 1.0
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     //    t2 = t1 + top*rand;
0218     // ignoring jitter from optical path length
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++) {  // here "4" is set deliberately,
0223         t4 = t3 + k;                                // as in test fortran toy MC ...
0224         if (t4 < nbin) {
0225           unsigned int ntb = t4;
0226           ntmp[ntb] += ntd[i] * nth[j] * ntp[k];
0227         }
0228       }
0229     }
0230   }
0231 
0232   // normalize for 1 GeV pulse height
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   // first create pulse shape over a range of time 0 ns to 255 ns in 1 ns steps
0249   unsigned int nbin = 256;
0250   hfShape_.setNBin(nbin);
0251   std::vector<float> ntmp(nbin, 0.0);  //
0252 
0253   const float k0 = 0.7956;  // shape parameters
0254   const float p2 = 1.355;
0255   const float p4 = 2.327;
0256   const float p1 = 4.3;  // position parameter
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   // normalize pulse area to 1.0
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   //modified shape 206
0278   //7.2 ns shift in 206
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   //Combination of all phase scan data (May,Jul,Oct2017)
0329   //runs:  294736-294740, 294929-294950, 298594-298598 and 305744-305758
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   //From Jay Lawhorn: derived from data Edward Laird phase scan may2017
0382   //https://indico.cern.ch/event/641978/contributions/2604491/attachments/1468666/2271582/17-05-31-hcal-hep17-pulse-shape.pdf
0383   //Run numbers are 294736-294740 and 294929-294950
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   //From Jake Anderson: toy MC convolution of SiPM pulse + WLS fiber shape + SiPM nonlinear response
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   //numerical convolution of SiPM pulse + WLS fiber shape
0487   static const HcalPulseShape siPMShapeMC2017(
0488       normalize(convolve(nBinsSiPM_, analyticPulseShapeSiPMHE, Y11203), nBinsSiPM_), nBinsSiPM_);
0489   return siPMShapeMC2017;
0490 }
0491 
0492 const HcalPulseShape& HcalPulseShapes::computeSiPMShapeHE206() {
0493   //numerical convolution of SiPM pulse + WLS fiber shape
0494   //shift: aligning 206 phase closer to 205 in order to have good reco agreement
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_;  // should not return this, but...
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       //FIXME doesn't look for SiPMs
0550       return hoShape(false);
0551     default:
0552       throw cms::Exception("HcalPulseShapes") << "unknown detId";
0553       break;
0554   }
0555 }
0556 
0557 //SiPM helpers
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   // HO SiPM pulse shape fit from Jake Anderson ca. 2013
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   // taken from fit to laser measurement taken by Iouri M. in Spring 2016.
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 //Original scintillator+Y11 fit from Vasken's 2001 measurement
0608 double HcalPulseShapes::Y11203(double t) { return exp(-0.0635 - 0.1518 * t + log(t) * 2.528) / 2485.9; }
0609 
0610 //New scintillator+Y11 model from Vasken's 2017 measurement plus a Landau correction term
0611 double HcalPulseShapes::Y11206(double t) {
0612   //Shifting phase to have better comparison of digi shape with data
0613   //If necessary, further digi phase adjustment can be done here:
0614   //SimCalorimetry/HcalSimProducers/python/hcalSimParameters_cfi.py
0615   //by changing "timePhase"
0616   double shift = 7.2;
0617 
0618   //Fit From Deconvolved Data
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   //Correction Term
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   //Overall Y11
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 }