File indexing completed on 2022-02-16 06:16:01
0001 #include "RecoLocalCalo/EcalRecProducers/plugins/EcalUncalibRecHitWorkerMultiFit.h"
0002
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/Run.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008
0009 #include <FWCore/ParameterSet/interface/ConfigurationDescriptions.h>
0010 #include <FWCore/ParameterSet/interface/ParameterSetDescription.h>
0011 #include <FWCore/ParameterSet/interface/EmptyGroupDescription.h>
0012
0013 EcalUncalibRecHitWorkerMultiFit::EcalUncalibRecHitWorkerMultiFit(const edm::ParameterSet& ps, edm::ConsumesCollector& c)
0014 : EcalUncalibRecHitWorkerBaseClass(ps, c) {
0015
0016 std::vector<int32_t> activeBXs = ps.getParameter<std::vector<int32_t>>("activeBXs");
0017 activeBX.resize(activeBXs.size());
0018 for (unsigned int ibx = 0; ibx < activeBXs.size(); ++ibx) {
0019 activeBX.coeffRef(ibx) = activeBXs[ibx];
0020 }
0021
0022
0023 ampErrorCalculation_ = ps.getParameter<bool>("ampErrorCalculation");
0024 useLumiInfoRunHeader_ = ps.getParameter<bool>("useLumiInfoRunHeader");
0025
0026 if (useLumiInfoRunHeader_) {
0027 bunchSpacing_ = c.consumes<unsigned int>(edm::InputTag("bunchSpacingProducer"));
0028 bunchSpacingManual_ = 0;
0029 } else {
0030 bunchSpacingManual_ = ps.getParameter<int>("bunchSpacing");
0031 }
0032
0033 doPrefitEB_ = ps.getParameter<bool>("doPrefitEB");
0034 doPrefitEE_ = ps.getParameter<bool>("doPrefitEE");
0035
0036 prefitMaxChiSqEB_ = ps.getParameter<double>("prefitMaxChiSqEB");
0037 prefitMaxChiSqEE_ = ps.getParameter<double>("prefitMaxChiSqEE");
0038
0039 dynamicPedestalsEB_ = ps.getParameter<bool>("dynamicPedestalsEB");
0040 dynamicPedestalsEE_ = ps.getParameter<bool>("dynamicPedestalsEE");
0041 mitigateBadSamplesEB_ = ps.getParameter<bool>("mitigateBadSamplesEB");
0042 mitigateBadSamplesEE_ = ps.getParameter<bool>("mitigateBadSamplesEE");
0043 gainSwitchUseMaxSampleEB_ = ps.getParameter<bool>("gainSwitchUseMaxSampleEB");
0044 gainSwitchUseMaxSampleEE_ = ps.getParameter<bool>("gainSwitchUseMaxSampleEE");
0045 selectiveBadSampleCriteriaEB_ = ps.getParameter<bool>("selectiveBadSampleCriteriaEB");
0046 selectiveBadSampleCriteriaEE_ = ps.getParameter<bool>("selectiveBadSampleCriteriaEE");
0047 addPedestalUncertaintyEB_ = ps.getParameter<double>("addPedestalUncertaintyEB");
0048 addPedestalUncertaintyEE_ = ps.getParameter<double>("addPedestalUncertaintyEE");
0049 simplifiedNoiseModelForGainSwitch_ = ps.getParameter<bool>("simplifiedNoiseModelForGainSwitch");
0050 pedsToken_ = c.esConsumes<EcalPedestals, EcalPedestalsRcd>();
0051 gainsToken_ = c.esConsumes<EcalGainRatios, EcalGainRatiosRcd>();
0052 noiseConvariancesToken_ = c.esConsumes<EcalSamplesCorrelation, EcalSamplesCorrelationRcd>();
0053 pulseShapesToken_ = c.esConsumes<EcalPulseShapes, EcalPulseShapesRcd>();
0054 pulseConvariancesToken_ = c.esConsumes<EcalPulseCovariances, EcalPulseCovariancesRcd>();
0055 sampleMaskToken_ = c.esConsumes<EcalSampleMask, EcalSampleMaskRcd>();
0056 grpsToken_ = c.esConsumes<EcalWeightXtalGroups, EcalWeightXtalGroupsRcd>();
0057 wgtsToken_ = c.esConsumes<EcalTBWeights, EcalTBWeightsRcd>();
0058 timeCorrBiasToken_ = c.esConsumes<EcalTimeBiasCorrections, EcalTimeBiasCorrectionsRcd>();
0059 itimeToken_ = c.esConsumes<EcalTimeCalibConstants, EcalTimeCalibConstantsRcd>();
0060 offtimeToken_ = c.esConsumes<EcalTimeOffsetConstant, EcalTimeOffsetConstantRcd>();
0061
0062
0063 auto const& timeAlgoName = ps.getParameter<std::string>("timealgo");
0064 if (timeAlgoName == "RatioMethod")
0065 timealgo_ = ratioMethod;
0066 else if (timeAlgoName == "WeightsMethod")
0067 timealgo_ = weightsMethod;
0068 else if (timeAlgoName == "crossCorrelationMethod") {
0069 timealgo_ = crossCorrelationMethod;
0070 double startTime = ps.getParameter<double>("crossCorrelationStartTime");
0071 double stopTime = ps.getParameter<double>("crossCorrelationStopTime");
0072 double targetTimePrecision = ps.getParameter<double>("crossCorrelationTargetTimePrecision");
0073 computeCC_ = std::make_unique<EcalUncalibRecHitTimingCCAlgo>(startTime, stopTime, targetTimePrecision);
0074 } else if (timeAlgoName != "None")
0075 edm::LogError("EcalUncalibRecHitError") << "No time estimation algorithm defined";
0076
0077
0078 EBtimeFitParameters_ = ps.getParameter<std::vector<double>>("EBtimeFitParameters");
0079 EEtimeFitParameters_ = ps.getParameter<std::vector<double>>("EEtimeFitParameters");
0080 EBamplitudeFitParameters_ = ps.getParameter<std::vector<double>>("EBamplitudeFitParameters");
0081 EEamplitudeFitParameters_ = ps.getParameter<std::vector<double>>("EEamplitudeFitParameters");
0082 EBtimeFitLimits_.first = ps.getParameter<double>("EBtimeFitLimits_Lower");
0083 EBtimeFitLimits_.second = ps.getParameter<double>("EBtimeFitLimits_Upper");
0084 EEtimeFitLimits_.first = ps.getParameter<double>("EEtimeFitLimits_Lower");
0085 EEtimeFitLimits_.second = ps.getParameter<double>("EEtimeFitLimits_Upper");
0086 EBtimeConstantTerm_ = ps.getParameter<double>("EBtimeConstantTerm");
0087 EEtimeConstantTerm_ = ps.getParameter<double>("EEtimeConstantTerm");
0088 EBtimeNconst_ = ps.getParameter<double>("EBtimeNconst");
0089 EEtimeNconst_ = ps.getParameter<double>("EEtimeNconst");
0090 outOfTimeThreshG12pEB_ = ps.getParameter<double>("outOfTimeThresholdGain12pEB");
0091 outOfTimeThreshG12mEB_ = ps.getParameter<double>("outOfTimeThresholdGain12mEB");
0092 outOfTimeThreshG61pEB_ = ps.getParameter<double>("outOfTimeThresholdGain61pEB");
0093 outOfTimeThreshG61mEB_ = ps.getParameter<double>("outOfTimeThresholdGain61mEB");
0094 outOfTimeThreshG12pEE_ = ps.getParameter<double>("outOfTimeThresholdGain12pEE");
0095 outOfTimeThreshG12mEE_ = ps.getParameter<double>("outOfTimeThresholdGain12mEE");
0096 outOfTimeThreshG61pEE_ = ps.getParameter<double>("outOfTimeThresholdGain61pEE");
0097 outOfTimeThreshG61mEE_ = ps.getParameter<double>("outOfTimeThresholdGain61mEE");
0098 amplitudeThreshEB_ = ps.getParameter<double>("amplitudeThresholdEB");
0099 amplitudeThreshEE_ = ps.getParameter<double>("amplitudeThresholdEE");
0100 }
0101
0102 void EcalUncalibRecHitWorkerMultiFit::set(const edm::EventSetup& es) {
0103
0104 gains = es.getHandle(gainsToken_);
0105 peds = es.getHandle(pedsToken_);
0106
0107
0108 if (!ampErrorCalculation_)
0109 multiFitMethod_.disableErrorCalculation();
0110 noisecovariances = es.getHandle(noiseConvariancesToken_);
0111 pulseshapes = es.getHandle(pulseShapesToken_);
0112 pulsecovariances = es.getHandle(pulseConvariancesToken_);
0113
0114
0115 grps = es.getHandle(grpsToken_);
0116 wgts = es.getHandle(wgtsToken_);
0117
0118
0119 sampleMaskHand_ = es.getHandle(sampleMaskToken_);
0120
0121
0122 itime = es.getHandle(itimeToken_);
0123 offtime = es.getHandle(offtimeToken_);
0124
0125
0126 timeCorrBias_ = es.getHandle(timeCorrBiasToken_);
0127
0128 int nnoise = SampleVector::RowsAtCompileTime;
0129 SampleMatrix& noisecorEBg12 = noisecors_[1][0];
0130 SampleMatrix& noisecorEBg6 = noisecors_[1][1];
0131 SampleMatrix& noisecorEBg1 = noisecors_[1][2];
0132 SampleMatrix& noisecorEEg12 = noisecors_[0][0];
0133 SampleMatrix& noisecorEEg6 = noisecors_[0][1];
0134 SampleMatrix& noisecorEEg1 = noisecors_[0][2];
0135
0136 for (int i = 0; i < nnoise; ++i) {
0137 for (int j = 0; j < nnoise; ++j) {
0138 int vidx = std::abs(j - i);
0139 noisecorEBg12(i, j) = noisecovariances->EBG12SamplesCorrelation[vidx];
0140 noisecorEEg12(i, j) = noisecovariances->EEG12SamplesCorrelation[vidx];
0141 noisecorEBg6(i, j) = noisecovariances->EBG6SamplesCorrelation[vidx];
0142 noisecorEEg6(i, j) = noisecovariances->EEG6SamplesCorrelation[vidx];
0143 noisecorEBg1(i, j) = noisecovariances->EBG1SamplesCorrelation[vidx];
0144 noisecorEEg1(i, j) = noisecovariances->EEG1SamplesCorrelation[vidx];
0145 }
0146 }
0147 }
0148
0149 void EcalUncalibRecHitWorkerMultiFit::set(const edm::Event& evt) {
0150 unsigned int bunchspacing = 450;
0151
0152 if (useLumiInfoRunHeader_) {
0153 edm::Handle<unsigned int> bunchSpacingH;
0154 evt.getByToken(bunchSpacing_, bunchSpacingH);
0155 bunchspacing = *bunchSpacingH;
0156 } else {
0157 bunchspacing = bunchSpacingManual_;
0158 }
0159
0160 if (useLumiInfoRunHeader_ || bunchSpacingManual_ > 0) {
0161 if (bunchspacing == 25) {
0162 activeBX.resize(10);
0163 activeBX << -5, -4, -3, -2, -1, 0, 1, 2, 3, 4;
0164 } else {
0165
0166 activeBX.resize(5);
0167 activeBX << -4, -2, 0, 2, 4;
0168 }
0169 }
0170 }
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181 double EcalUncalibRecHitWorkerMultiFit::timeCorrection(float ampli,
0182 const std::vector<float>& amplitudeBins,
0183 const std::vector<float>& shiftBins) {
0184
0185
0186 double theCorrection = 0;
0187
0188
0189 if (amplitudeBins.empty()) {
0190 edm::LogError("EcalRecHitError") << "timeCorrAmplitudeBins is empty, forcing no time bias corrections.";
0191
0192 return 0;
0193 }
0194
0195 if (amplitudeBins.size() != shiftBins.size()) {
0196 edm::LogError("EcalRecHitError") << "Size of timeCorrAmplitudeBins different from "
0197 "timeCorrShiftBins. Forcing no time bias corrections. ";
0198
0199 return 0;
0200 }
0201
0202
0203 int myBin = -1;
0204 for (int bin = 0; bin < (int)amplitudeBins.size(); bin++) {
0205 if (ampli > amplitudeBins[bin]) {
0206 myBin = bin;
0207 } else {
0208 break;
0209 }
0210 }
0211
0212 if (myBin == -1) {
0213 theCorrection = shiftBins[0];
0214 } else if (myBin == ((int)(amplitudeBins.size() - 1))) {
0215 theCorrection = shiftBins[myBin];
0216 } else {
0217
0218 theCorrection = (shiftBins[myBin + 1] - shiftBins[myBin]);
0219 theCorrection *= (((double)ampli) - amplitudeBins[myBin]) / (amplitudeBins[myBin + 1] - amplitudeBins[myBin]);
0220 theCorrection += shiftBins[myBin];
0221 }
0222
0223
0224 constexpr double inv25 = 1. / 25.;
0225 return theCorrection * inv25;
0226 }
0227
0228 void EcalUncalibRecHitWorkerMultiFit::run(const edm::Event& evt,
0229 const EcalDigiCollection& digis,
0230 EcalUncalibratedRecHitCollection& result) {
0231 if (digis.empty())
0232 return;
0233
0234
0235 DetId detid(digis.begin()->id());
0236 bool barrel = (detid.subdetId() == EcalBarrel);
0237
0238 multiFitMethod_.setSimplifiedNoiseModelForGainSwitch(simplifiedNoiseModelForGainSwitch_);
0239 if (barrel) {
0240 multiFitMethod_.setDoPrefit(doPrefitEB_);
0241 multiFitMethod_.setPrefitMaxChiSq(prefitMaxChiSqEB_);
0242 multiFitMethod_.setDynamicPedestals(dynamicPedestalsEB_);
0243 multiFitMethod_.setMitigateBadSamples(mitigateBadSamplesEB_);
0244 multiFitMethod_.setGainSwitchUseMaxSample(gainSwitchUseMaxSampleEB_);
0245 multiFitMethod_.setSelectiveBadSampleCriteria(selectiveBadSampleCriteriaEB_);
0246 multiFitMethod_.setAddPedestalUncertainty(addPedestalUncertaintyEB_);
0247 } else {
0248 multiFitMethod_.setDoPrefit(doPrefitEE_);
0249 multiFitMethod_.setPrefitMaxChiSq(prefitMaxChiSqEE_);
0250 multiFitMethod_.setDynamicPedestals(dynamicPedestalsEE_);
0251 multiFitMethod_.setMitigateBadSamples(mitigateBadSamplesEE_);
0252 multiFitMethod_.setGainSwitchUseMaxSample(gainSwitchUseMaxSampleEE_);
0253 multiFitMethod_.setSelectiveBadSampleCriteria(selectiveBadSampleCriteriaEE_);
0254 multiFitMethod_.setAddPedestalUncertainty(addPedestalUncertaintyEE_);
0255 }
0256
0257 FullSampleVector fullpulse(FullSampleVector::Zero());
0258 FullSampleMatrix fullpulsecov(FullSampleMatrix::Zero());
0259
0260 result.reserve(result.size() + digis.size());
0261 for (auto itdg = digis.begin(); itdg != digis.end(); ++itdg) {
0262 DetId detid(itdg->id());
0263
0264 const EcalSampleMask* sampleMask_ = sampleMaskHand_.product();
0265
0266
0267 float offsetTime = 0;
0268
0269 const EcalPedestals::Item* aped = nullptr;
0270 const EcalMGPAGainRatio* aGain = nullptr;
0271 const EcalXtalGroupId* gid = nullptr;
0272 const EcalPulseShapes::Item* aPulse = nullptr;
0273 const EcalPulseCovariances::Item* aPulseCov = nullptr;
0274
0275 if (barrel) {
0276 unsigned int hashedIndex = EBDetId(detid).hashedIndex();
0277 aped = &peds->barrel(hashedIndex);
0278 aGain = &gains->barrel(hashedIndex);
0279 gid = &grps->barrel(hashedIndex);
0280 aPulse = &pulseshapes->barrel(hashedIndex);
0281 aPulseCov = &pulsecovariances->barrel(hashedIndex);
0282 offsetTime = offtime->getEBValue();
0283 } else {
0284 unsigned int hashedIndex = EEDetId(detid).hashedIndex();
0285 aped = &peds->endcap(hashedIndex);
0286 aGain = &gains->endcap(hashedIndex);
0287 gid = &grps->endcap(hashedIndex);
0288 aPulse = &pulseshapes->endcap(hashedIndex);
0289 aPulseCov = &pulsecovariances->endcap(hashedIndex);
0290 offsetTime = offtime->getEEValue();
0291 }
0292
0293 double pedVec[3] = {aped->mean_x12, aped->mean_x6, aped->mean_x1};
0294 double pedRMSVec[3] = {aped->rms_x12, aped->rms_x6, aped->rms_x1};
0295 double gainRatios[3] = {1., aGain->gain12Over6(), aGain->gain6Over1() * aGain->gain12Over6()};
0296
0297 for (int i = 0; i < EcalPulseShape::TEMPLATESAMPLES; ++i)
0298 fullpulse(i + 7) = aPulse->pdfval[i];
0299
0300 for (int i = 0; i < EcalPulseShape::TEMPLATESAMPLES; i++)
0301 for (int j = 0; j < EcalPulseShape::TEMPLATESAMPLES; j++)
0302 fullpulsecov(i + 7, j + 7) = aPulseCov->covval[i][j];
0303
0304
0305 EcalTimeCalibConstantMap::const_iterator it = itime->find(detid);
0306 EcalTimeCalibConstant itimeconst = 0;
0307 if (it != itime->end()) {
0308 itimeconst = (*it);
0309 } else {
0310 edm::LogError("EcalRecHitError") << "No time intercalib const found for xtal " << detid.rawId()
0311 << "! something wrong with EcalTimeCalibConstants in your DB? ";
0312 }
0313
0314 int lastSampleBeforeSaturation = -2;
0315 for (unsigned int iSample = 0; iSample < EcalDataFrame::MAXSAMPLES; iSample++) {
0316 if (((EcalDataFrame)(*itdg)).sample(iSample).gainId() == 0) {
0317 lastSampleBeforeSaturation = iSample - 1;
0318 break;
0319 }
0320 }
0321
0322
0323
0324 if (lastSampleBeforeSaturation == 4) {
0325 result.emplace_back((*itdg).id(), 4095 * 12, 0, 0, 0);
0326 auto& uncalibRecHit = result.back();
0327 uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kSaturated);
0328
0329 uncalibRecHit.setChi2(0);
0330 } else if (lastSampleBeforeSaturation >=
0331 -1) {
0332 int gainId = ((EcalDataFrame)(*itdg)).sample(5).gainId();
0333 if (gainId == 0)
0334 gainId = 3;
0335 auto pedestal = pedVec[gainId - 1];
0336 auto gainratio = gainRatios[gainId - 1];
0337 double amplitude = ((double)(((EcalDataFrame)(*itdg)).sample(5).adc()) - pedestal) * gainratio;
0338 result.emplace_back((*itdg).id(), amplitude, 0, 0, 0);
0339 auto& uncalibRecHit = result.back();
0340 uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kSaturated);
0341
0342 uncalibRecHit.setChi2(0);
0343 } else {
0344
0345 const SampleMatrixGainArray& noisecors = noisecor(barrel);
0346
0347 result.push_back(multiFitMethod_.makeRecHit(*itdg, aped, aGain, noisecors, fullpulse, fullpulsecov, activeBX));
0348 auto& uncalibRecHit = result.back();
0349
0350
0351 if (timealgo_ == ratioMethod) {
0352
0353 constexpr float clockToNsConstant = 25.;
0354 constexpr float invClockToNs = 1. / clockToNsConstant;
0355 if (not barrel) {
0356 ratioMethod_endcap_.init(*itdg, *sampleMask_, pedVec, pedRMSVec, gainRatios);
0357 ratioMethod_endcap_.computeTime(EEtimeFitParameters_, EEtimeFitLimits_, EEamplitudeFitParameters_);
0358 ratioMethod_endcap_.computeAmplitude(EEamplitudeFitParameters_);
0359 EcalUncalibRecHitRatioMethodAlgo<EEDataFrame>::CalculatedRecHit crh =
0360 ratioMethod_endcap_.getCalculatedRecHit();
0361 double theTimeCorrectionEE = timeCorrection(
0362 uncalibRecHit.amplitude(), timeCorrBias_->EETimeCorrAmplitudeBins, timeCorrBias_->EETimeCorrShiftBins);
0363
0364 uncalibRecHit.setJitter(crh.timeMax - 5 + theTimeCorrectionEE);
0365 uncalibRecHit.setJitterError(
0366 std::sqrt(std::pow(crh.timeError, 2) + std::pow(EEtimeConstantTerm_ * invClockToNs, 2)));
0367
0368
0369 if (uncalibRecHit.amplitude() > pedRMSVec[0] * amplitudeThreshEE_) {
0370 float outOfTimeThreshP = outOfTimeThreshG12pEE_;
0371 float outOfTimeThreshM = outOfTimeThreshG12mEE_;
0372
0373
0374
0375 if (uncalibRecHit.amplitude() > 3000.) {
0376 for (int iSample = 0; iSample < EEDataFrame::MAXSAMPLES; iSample++) {
0377 int GainId = ((EcalDataFrame)(*itdg)).sample(iSample).gainId();
0378 if (GainId != 1) {
0379 outOfTimeThreshP = outOfTimeThreshG61pEE_;
0380 outOfTimeThreshM = outOfTimeThreshG61mEE_;
0381 break;
0382 }
0383 }
0384 }
0385 float correctedTime = (crh.timeMax - 5) * clockToNsConstant + itimeconst + offsetTime;
0386 float cterm = EEtimeConstantTerm_;
0387 float sigmaped = pedRMSVec[0];
0388 float nterm = EEtimeNconst_ * sigmaped / uncalibRecHit.amplitude();
0389 float sigmat = std::sqrt(nterm * nterm + cterm * cterm);
0390 if ((correctedTime > sigmat * outOfTimeThreshP) || (correctedTime < -sigmat * outOfTimeThreshM)) {
0391 uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kOutOfTime);
0392 }
0393 }
0394
0395 } else {
0396 ratioMethod_barrel_.init(*itdg, *sampleMask_, pedVec, pedRMSVec, gainRatios);
0397 ratioMethod_barrel_.fixMGPAslew(*itdg);
0398 ratioMethod_barrel_.computeTime(EBtimeFitParameters_, EBtimeFitLimits_, EBamplitudeFitParameters_);
0399 ratioMethod_barrel_.computeAmplitude(EBamplitudeFitParameters_);
0400 EcalUncalibRecHitRatioMethodAlgo<EBDataFrame>::CalculatedRecHit crh =
0401 ratioMethod_barrel_.getCalculatedRecHit();
0402
0403 double theTimeCorrectionEB = timeCorrection(
0404 uncalibRecHit.amplitude(), timeCorrBias_->EBTimeCorrAmplitudeBins, timeCorrBias_->EBTimeCorrShiftBins);
0405
0406 uncalibRecHit.setJitter(crh.timeMax - 5 + theTimeCorrectionEB);
0407 uncalibRecHit.setJitterError(std::hypot(crh.timeError, EBtimeConstantTerm_ / clockToNsConstant));
0408
0409
0410 if (uncalibRecHit.amplitude() > pedRMSVec[0] * amplitudeThreshEB_) {
0411 float outOfTimeThreshP = outOfTimeThreshG12pEB_;
0412 float outOfTimeThreshM = outOfTimeThreshG12mEB_;
0413
0414
0415
0416 if (uncalibRecHit.amplitude() > 3000.) {
0417 for (int iSample = 0; iSample < EBDataFrame::MAXSAMPLES; iSample++) {
0418 int GainId = ((EcalDataFrame)(*itdg)).sample(iSample).gainId();
0419 if (GainId != 1) {
0420 outOfTimeThreshP = outOfTimeThreshG61pEB_;
0421 outOfTimeThreshM = outOfTimeThreshG61mEB_;
0422 break;
0423 }
0424 }
0425 }
0426 float correctedTime = (crh.timeMax - 5) * clockToNsConstant + itimeconst + offsetTime;
0427 float cterm = EBtimeConstantTerm_;
0428 float sigmaped = pedRMSVec[0];
0429 float nterm = EBtimeNconst_ * sigmaped / uncalibRecHit.amplitude();
0430 float sigmat = std::sqrt(nterm * nterm + cterm * cterm);
0431 if ((correctedTime > sigmat * outOfTimeThreshP) || (correctedTime < -sigmat * outOfTimeThreshM)) {
0432 uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kOutOfTime);
0433 }
0434 }
0435 }
0436 } else if (timealgo_ == weightsMethod) {
0437
0438 std::vector<double> amplitudes;
0439 for (unsigned int ibx = 0; ibx < activeBX.size(); ++ibx)
0440 amplitudes.push_back(uncalibRecHit.outOfTimeAmplitude(ibx));
0441
0442 EcalTBWeights::EcalTDCId tdcid(1);
0443 EcalTBWeights::EcalTBWeightMap const& wgtsMap = wgts->getMap();
0444 EcalTBWeights::EcalTBWeightMap::const_iterator wit;
0445 wit = wgtsMap.find(std::make_pair(*gid, tdcid));
0446 if (wit == wgtsMap.end()) {
0447 edm::LogError("EcalUncalibRecHitError")
0448 << "No weights found for EcalGroupId: " << gid->id() << " and EcalTDCId: " << tdcid
0449 << "\n skipping digi with id: " << detid.rawId();
0450 result.pop_back();
0451 continue;
0452 }
0453 const EcalWeightSet& wset = wit->second;
0454
0455 const EcalWeightSet::EcalWeightMatrix& mat1 = wset.getWeightsBeforeGainSwitch();
0456 const EcalWeightSet::EcalWeightMatrix& mat2 = wset.getWeightsAfterGainSwitch();
0457
0458 weights[0] = &mat1;
0459 weights[1] = &mat2;
0460
0461 double timerh;
0462 if (detid.subdetId() == EcalEndcap) {
0463 timerh = weightsMethod_endcap_.time(*itdg, amplitudes, aped, aGain, fullpulse, weights);
0464 } else {
0465 timerh = weightsMethod_barrel_.time(*itdg, amplitudes, aped, aGain, fullpulse, weights);
0466 }
0467 uncalibRecHit.setJitter(timerh);
0468 uncalibRecHit.setJitterError(0.);
0469
0470 } else if (timealgo_ == crossCorrelationMethod) {
0471 std::vector<double> amplitudes(activeBX.size());
0472 for (unsigned int ibx = 0; ibx < activeBX.size(); ++ibx)
0473 amplitudes[ibx] = uncalibRecHit.outOfTimeAmplitude(ibx);
0474
0475 float jitterError = 0.;
0476 float jitter = computeCC_->computeTimeCC(*itdg, amplitudes, aped, aGain, fullpulse, uncalibRecHit, jitterError);
0477
0478 uncalibRecHit.setJitter(jitter);
0479 uncalibRecHit.setJitterError(jitterError);
0480
0481 } else {
0482 uncalibRecHit.setJitter(0.);
0483 uncalibRecHit.setJitterError(0.);
0484 }
0485 }
0486
0487
0488 auto& uncalibRecHit = result.back();
0489 if (((EcalDataFrame)(*itdg)).hasSwitchToGain6())
0490 uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kHasSwitchToGain6);
0491 if (((EcalDataFrame)(*itdg)).hasSwitchToGain1())
0492 uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kHasSwitchToGain1);
0493 }
0494 }
0495
0496 edm::ParameterSetDescription EcalUncalibRecHitWorkerMultiFit::getAlgoDescription() {
0497 edm::ParameterSetDescription psd;
0498 psd.addNode(edm::ParameterDescription<std::vector<int>>("activeBXs", {-5, -4, -3, -2, -1, 0, 1, 2, 3, 4}, true) and
0499 edm::ParameterDescription<bool>("ampErrorCalculation", true, true) and
0500 edm::ParameterDescription<bool>("useLumiInfoRunHeader", true, true) and
0501 edm::ParameterDescription<int>("bunchSpacing", 0, true) and
0502 edm::ParameterDescription<bool>("doPrefitEB", false, true) and
0503 edm::ParameterDescription<bool>("doPrefitEE", false, true) and
0504 edm::ParameterDescription<double>("prefitMaxChiSqEB", 25., true) and
0505 edm::ParameterDescription<double>("prefitMaxChiSqEE", 10., true) and
0506 edm::ParameterDescription<bool>("dynamicPedestalsEB", false, true) and
0507 edm::ParameterDescription<bool>("dynamicPedestalsEE", false, true) and
0508 edm::ParameterDescription<bool>("mitigateBadSamplesEB", false, true) and
0509 edm::ParameterDescription<bool>("mitigateBadSamplesEE", false, true) and
0510 edm::ParameterDescription<bool>("gainSwitchUseMaxSampleEB", false, true) and
0511 edm::ParameterDescription<bool>("gainSwitchUseMaxSampleEE", false, true) and
0512 edm::ParameterDescription<bool>("selectiveBadSampleCriteriaEB", false, true) and
0513 edm::ParameterDescription<bool>("selectiveBadSampleCriteriaEE", false, true) and
0514 edm::ParameterDescription<double>("addPedestalUncertaintyEB", 0., true) and
0515 edm::ParameterDescription<double>("addPedestalUncertaintyEE", 0., true) and
0516 edm::ParameterDescription<bool>("simplifiedNoiseModelForGainSwitch", true, true) and
0517 edm::ParameterDescription<std::string>("timealgo", "RatioMethod", true) and
0518 edm::ParameterDescription<std::vector<double>>("EBtimeFitParameters",
0519 {-2.015452e+00,
0520 3.130702e+00,
0521 -1.234730e+01,
0522 4.188921e+01,
0523 -8.283944e+01,
0524 9.101147e+01,
0525 -5.035761e+01,
0526 1.105621e+01},
0527 true) and
0528 edm::ParameterDescription<std::vector<double>>("EEtimeFitParameters",
0529 {-2.390548e+00,
0530 3.553628e+00,
0531 -1.762341e+01,
0532 6.767538e+01,
0533 -1.332130e+02,
0534 1.407432e+02,
0535 -7.541106e+01,
0536 1.620277e+01},
0537 true) and
0538 edm::ParameterDescription<std::vector<double>>("EBamplitudeFitParameters", {1.138, 1.652}, true) and
0539 edm::ParameterDescription<std::vector<double>>("EEamplitudeFitParameters", {1.890, 1.400}, true) and
0540 edm::ParameterDescription<double>("EBtimeFitLimits_Lower", 0.2, true) and
0541 edm::ParameterDescription<double>("EBtimeFitLimits_Upper", 1.4, true) and
0542 edm::ParameterDescription<double>("EEtimeFitLimits_Lower", 0.2, true) and
0543 edm::ParameterDescription<double>("EEtimeFitLimits_Upper", 1.4, true) and
0544 edm::ParameterDescription<double>("EBtimeConstantTerm", .6, true) and
0545 edm::ParameterDescription<double>("EEtimeConstantTerm", 1.0, true) and
0546 edm::ParameterDescription<double>("EBtimeNconst", 28.5, true) and
0547 edm::ParameterDescription<double>("EEtimeNconst", 31.8, true) and
0548 edm::ParameterDescription<double>("outOfTimeThresholdGain12pEB", 5, true) and
0549 edm::ParameterDescription<double>("outOfTimeThresholdGain12mEB", 5, true) and
0550 edm::ParameterDescription<double>("outOfTimeThresholdGain61pEB", 5, true) and
0551 edm::ParameterDescription<double>("outOfTimeThresholdGain61mEB", 5, true) and
0552 edm::ParameterDescription<double>("outOfTimeThresholdGain12pEE", 1000, true) and
0553 edm::ParameterDescription<double>("outOfTimeThresholdGain12mEE", 1000, true) and
0554 edm::ParameterDescription<double>("outOfTimeThresholdGain61pEE", 1000, true) and
0555 edm::ParameterDescription<double>("outOfTimeThresholdGain61mEE", 1000, true) and
0556 edm::ParameterDescription<double>("amplitudeThresholdEB", 10, true) and
0557 edm::ParameterDescription<double>("amplitudeThresholdEE", 10, true) and
0558 edm::ParameterDescription<double>("crossCorrelationStartTime", -25.0, true) and
0559 edm::ParameterDescription<double>("crossCorrelationStopTime", 25.0, true) and
0560 edm::ParameterDescription<double>("crossCorrelationTargetTimePrecision", 0.01, true));
0561
0562 return psd;
0563 }
0564
0565 #include "FWCore/Framework/interface/MakerMacros.h"
0566 #include "RecoLocalCalo/EcalRecProducers/interface/EcalUncalibRecHitWorkerFactory.h"
0567 DEFINE_EDM_PLUGIN(EcalUncalibRecHitWorkerFactory, EcalUncalibRecHitWorkerMultiFit, "EcalUncalibRecHitWorkerMultiFit");
0568 #include "RecoLocalCalo/EcalRecProducers/interface/EcalUncalibRecHitFillDescriptionWorkerFactory.h"
0569 DEFINE_EDM_PLUGIN(EcalUncalibRecHitFillDescriptionWorkerFactory,
0570 EcalUncalibRecHitWorkerMultiFit,
0571 "EcalUncalibRecHitWorkerMultiFit");