File indexing completed on 2024-04-06 12:29:32
0001 #include "SimCalorimetry/HcalSimAlgos/interface/HcalTimeSlewSim.h"
0002 #include "SimCalorimetry/CaloSimAlgos/interface/CaloSimParameters.h"
0003 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSimParameters.h"
0004 #include "CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h"
0005 #include "DataFormats/DetId/interface/DetId.h"
0006 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0007 #include "DataFormats/HcalDetId/interface/HcalZDCDetId.h"
0008 #include "FWCore/Utilities/interface/Exception.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include "CLHEP/Random/RandGaussQ.h"
0012
0013 HcalTimeSlewSim::HcalTimeSlewSim(const CaloVSimParameterMap* parameterMap, double minFCToDelay)
0014 : theParameterMap(parameterMap), minFCToDelay_(minFCToDelay) {}
0015
0016
0017 double HcalTimeSlewSim::charge(const CaloSamples& samples) const {
0018 double totalCharge = 0.;
0019 for (int i = 0; i < 4; ++i) {
0020 int bin = i + samples.presamples();
0021 if (bin < samples.size()) {
0022 totalCharge += samples[bin];
0023 }
0024 }
0025 return totalCharge;
0026 }
0027
0028 void HcalTimeSlewSim::delay(CaloSamples& cs,
0029 CLHEP::HepRandomEngine* engine,
0030 const HcalTimeSlew* hcalTimeSlew_delay) const {
0031
0032
0033
0034 DetId detId(cs.id());
0035 if (detId.det() == DetId::Calo && detId.subdetId() == HcalZDCDetId::SubdetectorId)
0036 return;
0037 HcalDetId hcalDetId(detId);
0038
0039 if (hcalDetId.subdet() == HcalBarrel || hcalDetId.subdet() == HcalEndcap || hcalDetId.subdet() == HcalOuter) {
0040 HcalTimeSlew::BiasSetting biasSetting =
0041 (hcalDetId.subdet() == HcalOuter) ? HcalTimeSlew::Slow : HcalTimeSlew::Medium;
0042
0043
0044
0045 int maxbin = cs.size();
0046 CaloSamples data(detId, maxbin);
0047 data = cs;
0048
0049
0050 int soi = cs.presamples();
0051 double eps = 1.e-6;
0052 double scale_factor = 0.6;
0053 double scale = cs[soi] / scale_factor;
0054 double smearns = 0.;
0055 double cut = minFCToDelay_;
0056
0057 const HcalSimParameters& params = static_cast<const HcalSimParameters&>(theParameterMap->simParameters(detId));
0058 if (params.doTimeSmear()) {
0059 double rms = params.timeSmearRMS(scale);
0060 smearns = CLHEP::RandGaussQ::shoot(engine) * rms;
0061 LogDebug("HcalTimeSlewSim") << "TimeSmear charge " << scale << " rms " << rms << " smearns " << smearns;
0062 }
0063
0064
0065 for (int it = 0; it < maxbin - 1;) {
0066 double datai = cs[it];
0067 int nts = 0;
0068 double tshift = smearns;
0069 double totalCharge = datai / scale_factor;
0070
0071
0072 if (totalCharge <= 0.)
0073 totalCharge = eps;
0074 tshift += hcalTimeSlew_delay->delay(totalCharge, biasSetting);
0075 if (tshift <= 0.)
0076 tshift = eps;
0077
0078 if (cut > -999.) {
0079 if ((datai > cut) && (it == 0 || (datai > cs[it - 1]))) {
0080
0081
0082 nts = 2;
0083 if (datai > cs[it + 1])
0084 nts = 3;
0085
0086
0087
0088 for (int j = it; j < it + nts && j < maxbin; ++j) {
0089
0090
0091 double t = j * 25. - tshift;
0092 int firstbin = floor(t / 25.);
0093 double f = t / 25. - firstbin;
0094 int nextbin = firstbin + 1;
0095 double v1 = (firstbin < 0 || firstbin >= maxbin) ? 0. : cs[firstbin];
0096 double v2 = (nextbin < 0 || nextbin >= maxbin) ? 0. : cs[nextbin];
0097
0098
0099 if (nts == 2) {
0100 if (j == it)
0101 data[j] = v2 * f;
0102 else
0103 data[j] = v1 * (1. - f) + v2;
0104 } else {
0105 if (j == it)
0106 data[j] = v2 * f;
0107 if (j == it + 1)
0108 data[j] = v1 * (1. - f) + v2 * f;
0109 if (j == it + nts - 1)
0110 data[j] = v1 * (1. - f) + v2;
0111 }
0112
0113 }
0114 cs = data;
0115
0116 }
0117 } else {
0118 double t = it * 25. - tshift;
0119 int firstbin = floor(t / 25.);
0120 double f = t / 25. - firstbin;
0121 int nextbin = firstbin + 1;
0122 double v2 = (nextbin < 0 || nextbin >= maxbin) ? 0. : data[nextbin];
0123 data[it] = v2 * f;
0124 data[it + 1] += (v2 - data[it]);
0125 cs = data;
0126 }
0127
0128 if (nts < 3)
0129 it++;
0130 else
0131 it += 2;
0132 }
0133
0134
0135 cs = data;
0136 }
0137 }