File indexing completed on 2024-09-07 04:37:36
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "CondFormats/DataRecord/interface/EcalGainRatiosRcd.h"
0010 #include "CondFormats/DataRecord/interface/EcalPedestalsRcd.h"
0011 #include "CondFormats/DataRecord/interface/EcalPulseCovariancesRcd.h"
0012 #include "CondFormats/DataRecord/interface/EcalPulseShapesRcd.h"
0013 #include "CondFormats/DataRecord/interface/EcalSampleMaskRcd.h"
0014 #include "CondFormats/DataRecord/interface/EcalSamplesCorrelationRcd.h"
0015 #include "CondFormats/DataRecord/interface/EcalTBWeightsRcd.h"
0016 #include "CondFormats/DataRecord/interface/EcalTimeBiasCorrectionsRcd.h"
0017 #include "CondFormats/DataRecord/interface/EcalTimeCalibConstantsRcd.h"
0018 #include "CondFormats/DataRecord/interface/EcalTimeOffsetConstantRcd.h"
0019 #include "CondFormats/DataRecord/interface/EcalWeightXtalGroupsRcd.h"
0020 #include "CondFormats/EcalObjects/interface/EcalGainRatios.h"
0021 #include "CondFormats/EcalObjects/interface/EcalPedestals.h"
0022 #include "CondFormats/EcalObjects/interface/EcalPulseCovariances.h"
0023 #include "CondFormats/EcalObjects/interface/EcalPulseShapes.h"
0024 #include "CondFormats/EcalObjects/interface/EcalSampleMask.h"
0025 #include "CondFormats/EcalObjects/interface/EcalSamplesCorrelation.h"
0026 #include "CondFormats/EcalObjects/interface/EcalTBWeights.h"
0027 #include "CondFormats/EcalObjects/interface/EcalTimeBiasCorrections.h"
0028 #include "CondFormats/EcalObjects/interface/EcalTimeCalibConstants.h"
0029 #include "CondFormats/EcalObjects/interface/EcalTimeOffsetConstant.h"
0030 #include "CondFormats/EcalObjects/interface/EcalWeightXtalGroups.h"
0031 #include "FWCore/Framework/interface/ESHandle.h"
0032 #include "FWCore/Framework/interface/Event.h"
0033 #include "FWCore/Framework/interface/EventSetup.h"
0034 #include "FWCore/Framework/interface/Run.h"
0035 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0036 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0037 #include "FWCore/ParameterSet/interface/EmptyGroupDescription.h"
0038 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0039 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0040 #include "FWCore/Utilities/interface/ESGetToken.h"
0041 #include "FWCore/Utilities/interface/ESInputTag.h"
0042 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitMultiFitAlgo.h"
0043 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRatioMethodAlgo.h"
0044 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRecChi2Algo.h"
0045 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitTimeWeightsAlgo.h"
0046 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitTimingCCAlgo.h"
0047 #include "RecoLocalCalo/EcalRecAlgos/interface/EigenMatrixTypes.h"
0048 #include "RecoLocalCalo/EcalRecProducers/interface/EcalUncalibRecHitWorkerBaseClass.h"
0049
0050 class EcalUncalibRecHitWorkerMultiFit final : public EcalUncalibRecHitWorkerBaseClass {
0051 public:
0052 EcalUncalibRecHitWorkerMultiFit(const edm::ParameterSet&, edm::ConsumesCollector& c);
0053 EcalUncalibRecHitWorkerMultiFit() {}
0054
0055 private:
0056 void set(const edm::EventSetup& es) override;
0057 void set(const edm::Event& evt) override;
0058 void run(const edm::Event& evt, const EcalDigiCollection& digis, EcalUncalibratedRecHitCollection& result) override;
0059
0060 public:
0061 edm::ParameterSetDescription getAlgoDescription() override;
0062
0063 private:
0064 edm::ESHandle<EcalPedestals> peds;
0065 edm::ESGetToken<EcalPedestals, EcalPedestalsRcd> pedsToken_;
0066 edm::ESHandle<EcalGainRatios> gains;
0067 edm::ESGetToken<EcalGainRatios, EcalGainRatiosRcd> gainsToken_;
0068 edm::ESHandle<EcalSamplesCorrelation> noisecovariances;
0069 edm::ESGetToken<EcalSamplesCorrelation, EcalSamplesCorrelationRcd> noiseConvariancesToken_;
0070 edm::ESHandle<EcalPulseShapes> pulseshapes;
0071 edm::ESGetToken<EcalPulseShapes, EcalPulseShapesRcd> pulseShapesToken_;
0072 edm::ESHandle<EcalPulseCovariances> pulsecovariances;
0073 edm::ESGetToken<EcalPulseCovariances, EcalPulseCovariancesRcd> pulseConvariancesToken_;
0074
0075 double timeCorrection(float ampli, const std::vector<float>& amplitudeBins, const std::vector<float>& shiftBins);
0076
0077 const SampleMatrix& noisecor(bool barrel, int gain) const { return noisecors_[barrel ? 1 : 0][gain]; }
0078 const SampleMatrixGainArray& noisecor(bool barrel) const { return noisecors_[barrel ? 1 : 0]; }
0079
0080
0081 std::array<SampleMatrixGainArray, 2> noisecors_;
0082 BXVector activeBX;
0083 bool ampErrorCalculation_;
0084 bool useLumiInfoRunHeader_;
0085 EcalUncalibRecHitMultiFitAlgo multiFitMethod_;
0086
0087 int bunchSpacingManual_;
0088 edm::EDGetTokenT<unsigned int> bunchSpacing_;
0089
0090
0091 edm::ESHandle<EcalSampleMask> sampleMaskHand_;
0092 edm::ESGetToken<EcalSampleMask, EcalSampleMaskRcd> sampleMaskToken_;
0093
0094
0095 enum TimeAlgo { noMethod, ratioMethod, weightsMethod, crossCorrelationMethod };
0096 TimeAlgo timealgo_ = noMethod;
0097
0098
0099 edm::ESHandle<EcalWeightXtalGroups> grps;
0100 edm::ESGetToken<EcalWeightXtalGroups, EcalWeightXtalGroupsRcd> grpsToken_;
0101 edm::ESHandle<EcalTBWeights> wgts;
0102 edm::ESGetToken<EcalTBWeights, EcalTBWeightsRcd> wgtsToken_;
0103 const EcalWeightSet::EcalWeightMatrix* weights[2];
0104 EcalUncalibRecHitTimeWeightsAlgo<EBDataFrame> weightsMethod_barrel_;
0105 EcalUncalibRecHitTimeWeightsAlgo<EEDataFrame> weightsMethod_endcap_;
0106 bool doPrefitEB_;
0107 bool doPrefitEE_;
0108 double prefitMaxChiSqEB_;
0109 double prefitMaxChiSqEE_;
0110 bool dynamicPedestalsEB_;
0111 bool dynamicPedestalsEE_;
0112 bool mitigateBadSamplesEB_;
0113 bool mitigateBadSamplesEE_;
0114 bool gainSwitchUseMaxSampleEB_;
0115 bool gainSwitchUseMaxSampleEE_;
0116 bool selectiveBadSampleCriteriaEB_;
0117 bool selectiveBadSampleCriteriaEE_;
0118 double addPedestalUncertaintyEB_;
0119 double addPedestalUncertaintyEE_;
0120 bool simplifiedNoiseModelForGainSwitch_;
0121
0122
0123 std::vector<double> EBtimeFitParameters_;
0124 std::vector<double> EEtimeFitParameters_;
0125 std::vector<double> EBamplitudeFitParameters_;
0126 std::vector<double> EEamplitudeFitParameters_;
0127 std::pair<double, double> EBtimeFitLimits_;
0128 std::pair<double, double> EEtimeFitLimits_;
0129
0130 EcalUncalibRecHitRatioMethodAlgo<EBDataFrame> ratioMethod_barrel_;
0131 EcalUncalibRecHitRatioMethodAlgo<EEDataFrame> ratioMethod_endcap_;
0132
0133 double EBtimeConstantTerm_;
0134 double EEtimeConstantTerm_;
0135 double EBtimeNconst_;
0136 double EEtimeNconst_;
0137 double outOfTimeThreshG12pEB_;
0138 double outOfTimeThreshG12mEB_;
0139 double outOfTimeThreshG61pEB_;
0140 double outOfTimeThreshG61mEB_;
0141 double outOfTimeThreshG12pEE_;
0142 double outOfTimeThreshG12mEE_;
0143 double outOfTimeThreshG61pEE_;
0144 double outOfTimeThreshG61mEE_;
0145 double amplitudeThreshEB_;
0146 double amplitudeThreshEE_;
0147 double ebSpikeThresh_;
0148
0149 edm::ESHandle<EcalTimeBiasCorrections> timeCorrBias_;
0150 edm::ESGetToken<EcalTimeBiasCorrections, EcalTimeBiasCorrectionsRcd> timeCorrBiasToken_;
0151
0152 edm::ESHandle<EcalTimeCalibConstants> itime;
0153 edm::ESGetToken<EcalTimeCalibConstants, EcalTimeCalibConstantsRcd> itimeToken_;
0154 edm::ESHandle<EcalTimeOffsetConstant> offtime;
0155 edm::ESGetToken<EcalTimeOffsetConstant, EcalTimeOffsetConstantRcd> offtimeToken_;
0156 std::vector<double> ebPulseShape_;
0157 std::vector<double> eePulseShape_;
0158
0159
0160 bool kPoorRecoFlagEB_;
0161 bool kPoorRecoFlagEE_;
0162 double chi2ThreshEB_;
0163 double chi2ThreshEE_;
0164
0165
0166 std::unique_ptr<EcalUncalibRecHitTimingCCAlgo> computeCC_;
0167 double CCminTimeToBeLateMin_;
0168 double CCminTimeToBeLateMax_;
0169 double CCTimeShiftWrtRations_;
0170 double CCtargetTimePrecision_;
0171 double CCtargetTimePrecisionForDelayedPulses_;
0172 bool crossCorrelationUseSlewCorrectionEB_;
0173 bool crossCorrelationUseSlewCorrectionEE_;
0174 };
0175
0176 EcalUncalibRecHitWorkerMultiFit::EcalUncalibRecHitWorkerMultiFit(const edm::ParameterSet& ps, edm::ConsumesCollector& c)
0177 : EcalUncalibRecHitWorkerBaseClass(ps, c) {
0178
0179 std::vector<int32_t> activeBXs = ps.getParameter<std::vector<int32_t>>("activeBXs");
0180 activeBX.resize(activeBXs.size());
0181 for (unsigned int ibx = 0; ibx < activeBXs.size(); ++ibx) {
0182 activeBX.coeffRef(ibx) = activeBXs[ibx];
0183 }
0184
0185
0186 ampErrorCalculation_ = ps.getParameter<bool>("ampErrorCalculation");
0187 useLumiInfoRunHeader_ = ps.getParameter<bool>("useLumiInfoRunHeader");
0188
0189 if (useLumiInfoRunHeader_) {
0190 bunchSpacing_ = c.consumes<unsigned int>(edm::InputTag("bunchSpacingProducer"));
0191 bunchSpacingManual_ = 0;
0192 } else {
0193 bunchSpacingManual_ = ps.getParameter<int>("bunchSpacing");
0194 }
0195
0196 doPrefitEB_ = ps.getParameter<bool>("doPrefitEB");
0197 doPrefitEE_ = ps.getParameter<bool>("doPrefitEE");
0198
0199 prefitMaxChiSqEB_ = ps.getParameter<double>("prefitMaxChiSqEB");
0200 prefitMaxChiSqEE_ = ps.getParameter<double>("prefitMaxChiSqEE");
0201
0202 dynamicPedestalsEB_ = ps.getParameter<bool>("dynamicPedestalsEB");
0203 dynamicPedestalsEE_ = ps.getParameter<bool>("dynamicPedestalsEE");
0204 mitigateBadSamplesEB_ = ps.getParameter<bool>("mitigateBadSamplesEB");
0205 mitigateBadSamplesEE_ = ps.getParameter<bool>("mitigateBadSamplesEE");
0206 gainSwitchUseMaxSampleEB_ = ps.getParameter<bool>("gainSwitchUseMaxSampleEB");
0207 gainSwitchUseMaxSampleEE_ = ps.getParameter<bool>("gainSwitchUseMaxSampleEE");
0208 selectiveBadSampleCriteriaEB_ = ps.getParameter<bool>("selectiveBadSampleCriteriaEB");
0209 selectiveBadSampleCriteriaEE_ = ps.getParameter<bool>("selectiveBadSampleCriteriaEE");
0210 addPedestalUncertaintyEB_ = ps.getParameter<double>("addPedestalUncertaintyEB");
0211 addPedestalUncertaintyEE_ = ps.getParameter<double>("addPedestalUncertaintyEE");
0212 simplifiedNoiseModelForGainSwitch_ = ps.getParameter<bool>("simplifiedNoiseModelForGainSwitch");
0213 pedsToken_ = c.esConsumes<EcalPedestals, EcalPedestalsRcd>();
0214 gainsToken_ = c.esConsumes<EcalGainRatios, EcalGainRatiosRcd>();
0215 noiseConvariancesToken_ = c.esConsumes<EcalSamplesCorrelation, EcalSamplesCorrelationRcd>();
0216 pulseShapesToken_ = c.esConsumes<EcalPulseShapes, EcalPulseShapesRcd>();
0217 pulseConvariancesToken_ = c.esConsumes<EcalPulseCovariances, EcalPulseCovariancesRcd>();
0218 sampleMaskToken_ = c.esConsumes<EcalSampleMask, EcalSampleMaskRcd>();
0219 grpsToken_ = c.esConsumes<EcalWeightXtalGroups, EcalWeightXtalGroupsRcd>();
0220 wgtsToken_ = c.esConsumes<EcalTBWeights, EcalTBWeightsRcd>();
0221 timeCorrBiasToken_ = c.esConsumes<EcalTimeBiasCorrections, EcalTimeBiasCorrectionsRcd>();
0222 itimeToken_ =
0223 c.esConsumes<EcalTimeCalibConstants, EcalTimeCalibConstantsRcd>(ps.getParameter<edm::ESInputTag>("timeCalibTag"));
0224 offtimeToken_ = c.esConsumes<EcalTimeOffsetConstant, EcalTimeOffsetConstantRcd>(
0225 ps.getParameter<edm::ESInputTag>("timeOffsetTag"));
0226
0227
0228 auto const& timeAlgoName = ps.getParameter<std::string>("timealgo");
0229 if (timeAlgoName == "RatioMethod")
0230 timealgo_ = ratioMethod;
0231 else if (timeAlgoName == "WeightsMethod")
0232 timealgo_ = weightsMethod;
0233 else if (timeAlgoName == "crossCorrelationMethod") {
0234 timealgo_ = crossCorrelationMethod;
0235 double startTime = ps.getParameter<double>("crossCorrelationStartTime");
0236 double stopTime = ps.getParameter<double>("crossCorrelationStopTime");
0237 CCtargetTimePrecision_ = ps.getParameter<double>("crossCorrelationTargetTimePrecision");
0238 CCtargetTimePrecisionForDelayedPulses_ =
0239 ps.getParameter<double>("crossCorrelationTargetTimePrecisionForDelayedPulses");
0240 CCminTimeToBeLateMin_ = ps.getParameter<double>("crossCorrelationMinTimeToBeLateMin") / ecalcctiming::clockToNS;
0241 CCminTimeToBeLateMax_ = ps.getParameter<double>("crossCorrelationMinTimeToBeLateMax") / ecalcctiming::clockToNS;
0242 CCTimeShiftWrtRations_ = ps.getParameter<double>("crossCorrelationTimeShiftWrtRations");
0243 crossCorrelationUseSlewCorrectionEB_ = ps.getParameter<bool>("crossCorrelationUseSlewCorrectionEB");
0244 crossCorrelationUseSlewCorrectionEE_ = ps.getParameter<bool>("crossCorrelationUseSlewCorrectionEE");
0245 computeCC_ = std::make_unique<EcalUncalibRecHitTimingCCAlgo>(startTime, stopTime);
0246 } else if (timeAlgoName != "None")
0247 edm::LogError("EcalUncalibRecHitError") << "No time estimation algorithm defined";
0248
0249
0250 EBtimeFitParameters_ = ps.getParameter<std::vector<double>>("EBtimeFitParameters");
0251 EEtimeFitParameters_ = ps.getParameter<std::vector<double>>("EEtimeFitParameters");
0252 EBamplitudeFitParameters_ = ps.getParameter<std::vector<double>>("EBamplitudeFitParameters");
0253 EEamplitudeFitParameters_ = ps.getParameter<std::vector<double>>("EEamplitudeFitParameters");
0254 EBtimeFitLimits_.first = ps.getParameter<double>("EBtimeFitLimits_Lower");
0255 EBtimeFitLimits_.second = ps.getParameter<double>("EBtimeFitLimits_Upper");
0256 EEtimeFitLimits_.first = ps.getParameter<double>("EEtimeFitLimits_Lower");
0257 EEtimeFitLimits_.second = ps.getParameter<double>("EEtimeFitLimits_Upper");
0258 EBtimeConstantTerm_ = ps.getParameter<double>("EBtimeConstantTerm");
0259 EEtimeConstantTerm_ = ps.getParameter<double>("EEtimeConstantTerm");
0260 EBtimeNconst_ = ps.getParameter<double>("EBtimeNconst");
0261 EEtimeNconst_ = ps.getParameter<double>("EEtimeNconst");
0262 outOfTimeThreshG12pEB_ = ps.getParameter<double>("outOfTimeThresholdGain12pEB");
0263 outOfTimeThreshG12mEB_ = ps.getParameter<double>("outOfTimeThresholdGain12mEB");
0264 outOfTimeThreshG61pEB_ = ps.getParameter<double>("outOfTimeThresholdGain61pEB");
0265 outOfTimeThreshG61mEB_ = ps.getParameter<double>("outOfTimeThresholdGain61mEB");
0266 outOfTimeThreshG12pEE_ = ps.getParameter<double>("outOfTimeThresholdGain12pEE");
0267 outOfTimeThreshG12mEE_ = ps.getParameter<double>("outOfTimeThresholdGain12mEE");
0268 outOfTimeThreshG61pEE_ = ps.getParameter<double>("outOfTimeThresholdGain61pEE");
0269 outOfTimeThreshG61mEE_ = ps.getParameter<double>("outOfTimeThresholdGain61mEE");
0270 amplitudeThreshEB_ = ps.getParameter<double>("amplitudeThresholdEB");
0271 amplitudeThreshEE_ = ps.getParameter<double>("amplitudeThresholdEE");
0272 }
0273
0274 void EcalUncalibRecHitWorkerMultiFit::set(const edm::EventSetup& es) {
0275
0276 gains = es.getHandle(gainsToken_);
0277 peds = es.getHandle(pedsToken_);
0278
0279
0280 if (!ampErrorCalculation_)
0281 multiFitMethod_.disableErrorCalculation();
0282 noisecovariances = es.getHandle(noiseConvariancesToken_);
0283 pulseshapes = es.getHandle(pulseShapesToken_);
0284 pulsecovariances = es.getHandle(pulseConvariancesToken_);
0285
0286
0287 grps = es.getHandle(grpsToken_);
0288 wgts = es.getHandle(wgtsToken_);
0289
0290
0291 sampleMaskHand_ = es.getHandle(sampleMaskToken_);
0292
0293
0294 itime = es.getHandle(itimeToken_);
0295 offtime = es.getHandle(offtimeToken_);
0296
0297
0298 timeCorrBias_ = es.getHandle(timeCorrBiasToken_);
0299
0300 int nnoise = SampleVector::RowsAtCompileTime;
0301 SampleMatrix& noisecorEBg12 = noisecors_[1][0];
0302 SampleMatrix& noisecorEBg6 = noisecors_[1][1];
0303 SampleMatrix& noisecorEBg1 = noisecors_[1][2];
0304 SampleMatrix& noisecorEEg12 = noisecors_[0][0];
0305 SampleMatrix& noisecorEEg6 = noisecors_[0][1];
0306 SampleMatrix& noisecorEEg1 = noisecors_[0][2];
0307
0308 for (int i = 0; i < nnoise; ++i) {
0309 for (int j = 0; j < nnoise; ++j) {
0310 int vidx = std::abs(j - i);
0311 noisecorEBg12(i, j) = noisecovariances->EBG12SamplesCorrelation[vidx];
0312 noisecorEEg12(i, j) = noisecovariances->EEG12SamplesCorrelation[vidx];
0313 noisecorEBg6(i, j) = noisecovariances->EBG6SamplesCorrelation[vidx];
0314 noisecorEEg6(i, j) = noisecovariances->EEG6SamplesCorrelation[vidx];
0315 noisecorEBg1(i, j) = noisecovariances->EBG1SamplesCorrelation[vidx];
0316 noisecorEEg1(i, j) = noisecovariances->EEG1SamplesCorrelation[vidx];
0317 }
0318 }
0319 }
0320
0321 void EcalUncalibRecHitWorkerMultiFit::set(const edm::Event& evt) {
0322 unsigned int bunchspacing = 450;
0323
0324 if (useLumiInfoRunHeader_) {
0325 edm::Handle<unsigned int> bunchSpacingH;
0326 evt.getByToken(bunchSpacing_, bunchSpacingH);
0327 bunchspacing = *bunchSpacingH;
0328 } else {
0329 bunchspacing = bunchSpacingManual_;
0330 }
0331
0332 if (useLumiInfoRunHeader_ || bunchSpacingManual_ > 0) {
0333 if (bunchspacing == 25) {
0334 activeBX.resize(10);
0335 activeBX << -5, -4, -3, -2, -1, 0, 1, 2, 3, 4;
0336 } else {
0337
0338 activeBX.resize(5);
0339 activeBX << -4, -2, 0, 2, 4;
0340 }
0341 }
0342 }
0343
0344
0345
0346
0347
0348
0349
0350
0351
0352
0353 double EcalUncalibRecHitWorkerMultiFit::timeCorrection(float ampli,
0354 const std::vector<float>& amplitudeBins,
0355 const std::vector<float>& shiftBins) {
0356
0357
0358 double theCorrection = 0;
0359
0360
0361 if (amplitudeBins.empty()) {
0362 edm::LogError("EcalRecHitError") << "timeCorrAmplitudeBins is empty, forcing no time bias corrections.";
0363
0364 return 0;
0365 }
0366
0367 if (amplitudeBins.size() != shiftBins.size()) {
0368 edm::LogError("EcalRecHitError") << "Size of timeCorrAmplitudeBins different from "
0369 "timeCorrShiftBins. Forcing no time bias corrections. ";
0370
0371 return 0;
0372 }
0373
0374
0375 int myBin = -1;
0376 for (int bin = 0; bin < (int)amplitudeBins.size(); bin++) {
0377 if (ampli > amplitudeBins[bin]) {
0378 myBin = bin;
0379 } else {
0380 break;
0381 }
0382 }
0383
0384 if (myBin == -1) {
0385 theCorrection = shiftBins[0];
0386 } else if (myBin == ((int)(amplitudeBins.size() - 1))) {
0387 theCorrection = shiftBins[myBin];
0388 } else {
0389
0390 theCorrection = (shiftBins[myBin + 1] - shiftBins[myBin]);
0391 theCorrection *= (((double)ampli) - amplitudeBins[myBin]) / (amplitudeBins[myBin + 1] - amplitudeBins[myBin]);
0392 theCorrection += shiftBins[myBin];
0393 }
0394
0395
0396 constexpr double inv25 = 1. / 25.;
0397 return theCorrection * inv25;
0398 }
0399
0400 void EcalUncalibRecHitWorkerMultiFit::run(const edm::Event& evt,
0401 const EcalDigiCollection& digis,
0402 EcalUncalibratedRecHitCollection& result) {
0403 if (digis.empty())
0404 return;
0405
0406
0407 DetId detid(digis.begin()->id());
0408 bool barrel = (detid.subdetId() == EcalBarrel);
0409
0410 multiFitMethod_.setSimplifiedNoiseModelForGainSwitch(simplifiedNoiseModelForGainSwitch_);
0411 if (barrel) {
0412 multiFitMethod_.setDoPrefit(doPrefitEB_);
0413 multiFitMethod_.setPrefitMaxChiSq(prefitMaxChiSqEB_);
0414 multiFitMethod_.setDynamicPedestals(dynamicPedestalsEB_);
0415 multiFitMethod_.setMitigateBadSamples(mitigateBadSamplesEB_);
0416 multiFitMethod_.setGainSwitchUseMaxSample(gainSwitchUseMaxSampleEB_);
0417 multiFitMethod_.setSelectiveBadSampleCriteria(selectiveBadSampleCriteriaEB_);
0418 multiFitMethod_.setAddPedestalUncertainty(addPedestalUncertaintyEB_);
0419 } else {
0420 multiFitMethod_.setDoPrefit(doPrefitEE_);
0421 multiFitMethod_.setPrefitMaxChiSq(prefitMaxChiSqEE_);
0422 multiFitMethod_.setDynamicPedestals(dynamicPedestalsEE_);
0423 multiFitMethod_.setMitigateBadSamples(mitigateBadSamplesEE_);
0424 multiFitMethod_.setGainSwitchUseMaxSample(gainSwitchUseMaxSampleEE_);
0425 multiFitMethod_.setSelectiveBadSampleCriteria(selectiveBadSampleCriteriaEE_);
0426 multiFitMethod_.setAddPedestalUncertainty(addPedestalUncertaintyEE_);
0427 }
0428
0429 FullSampleVector fullpulse(FullSampleVector::Zero());
0430 FullSampleMatrix fullpulsecov(FullSampleMatrix::Zero());
0431
0432 result.reserve(result.size() + digis.size());
0433 for (auto itdg = digis.begin(); itdg != digis.end(); ++itdg) {
0434 DetId detid(itdg->id());
0435
0436 const EcalSampleMask* sampleMask_ = sampleMaskHand_.product();
0437
0438
0439 float offsetTime = 0;
0440
0441 const EcalPedestals::Item* aped = nullptr;
0442 const EcalMGPAGainRatio* aGain = nullptr;
0443 const EcalXtalGroupId* gid = nullptr;
0444 const EcalPulseShapes::Item* aPulse = nullptr;
0445 const EcalPulseCovariances::Item* aPulseCov = nullptr;
0446
0447 if (barrel) {
0448 unsigned int hashedIndex = EBDetId(detid).hashedIndex();
0449 aped = &peds->barrel(hashedIndex);
0450 aGain = &gains->barrel(hashedIndex);
0451 gid = &grps->barrel(hashedIndex);
0452 aPulse = &pulseshapes->barrel(hashedIndex);
0453 aPulseCov = &pulsecovariances->barrel(hashedIndex);
0454 offsetTime = offtime->getEBValue();
0455 } else {
0456 unsigned int hashedIndex = EEDetId(detid).hashedIndex();
0457 aped = &peds->endcap(hashedIndex);
0458 aGain = &gains->endcap(hashedIndex);
0459 gid = &grps->endcap(hashedIndex);
0460 aPulse = &pulseshapes->endcap(hashedIndex);
0461 aPulseCov = &pulsecovariances->endcap(hashedIndex);
0462 offsetTime = offtime->getEEValue();
0463 }
0464
0465 double pedVec[3] = {aped->mean_x12, aped->mean_x6, aped->mean_x1};
0466 double pedRMSVec[3] = {aped->rms_x12, aped->rms_x6, aped->rms_x1};
0467 double gainRatios[3] = {1., aGain->gain12Over6(), aGain->gain6Over1() * aGain->gain12Over6()};
0468
0469 for (int i = 0; i < EcalPulseShape::TEMPLATESAMPLES; ++i)
0470 fullpulse(i + 7) = aPulse->pdfval[i];
0471
0472 for (int i = 0; i < EcalPulseShape::TEMPLATESAMPLES; i++)
0473 for (int j = 0; j < EcalPulseShape::TEMPLATESAMPLES; j++)
0474 fullpulsecov(i + 7, j + 7) = aPulseCov->covval[i][j];
0475
0476
0477 EcalTimeCalibConstantMap::const_iterator it = itime->find(detid);
0478 EcalTimeCalibConstant itimeconst = 0;
0479 if (it != itime->end()) {
0480 itimeconst = (*it);
0481 } else {
0482 edm::LogError("EcalRecHitError") << "No time intercalib const found for xtal " << detid.rawId()
0483 << "! something wrong with EcalTimeCalibConstants in your DB? ";
0484 }
0485
0486 int lastSampleBeforeSaturation = -2;
0487 for (unsigned int iSample = 0; iSample < EcalDataFrame::MAXSAMPLES; iSample++) {
0488 if (((EcalDataFrame)(*itdg)).sample(iSample).gainId() == 0) {
0489 lastSampleBeforeSaturation = iSample - 1;
0490 break;
0491 }
0492 }
0493
0494
0495
0496 if (lastSampleBeforeSaturation == 4) {
0497 result.emplace_back((*itdg).id(), 4095 * 12, 0, 0, 0);
0498 auto& uncalibRecHit = result.back();
0499 uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kSaturated);
0500
0501 uncalibRecHit.setChi2(0);
0502 } else if (lastSampleBeforeSaturation >=
0503 -1) {
0504 int gainId = ((EcalDataFrame)(*itdg)).sample(5).gainId();
0505 if (gainId == 0)
0506 gainId = 3;
0507 auto pedestal = pedVec[gainId - 1];
0508 auto gainratio = gainRatios[gainId - 1];
0509 double amplitude = ((double)(((EcalDataFrame)(*itdg)).sample(5).adc()) - pedestal) * gainratio;
0510 result.emplace_back((*itdg).id(), amplitude, 0, 0, 0);
0511 auto& uncalibRecHit = result.back();
0512 uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kSaturated);
0513
0514 uncalibRecHit.setChi2(0);
0515 } else {
0516
0517 const SampleMatrixGainArray& noisecors = noisecor(barrel);
0518
0519 result.push_back(multiFitMethod_.makeRecHit(*itdg, aped, aGain, noisecors, fullpulse, fullpulsecov, activeBX));
0520 auto& uncalibRecHit = result.back();
0521
0522
0523 if (timealgo_ == ratioMethod) {
0524
0525 constexpr float clockToNsConstant = 25.;
0526 constexpr float invClockToNs = 1. / clockToNsConstant;
0527 if (not barrel) {
0528 ratioMethod_endcap_.init(*itdg, *sampleMask_, pedVec, pedRMSVec, gainRatios);
0529 ratioMethod_endcap_.computeTime(EEtimeFitParameters_, EEtimeFitLimits_, EEamplitudeFitParameters_);
0530 ratioMethod_endcap_.computeAmplitude(EEamplitudeFitParameters_);
0531 EcalUncalibRecHitRatioMethodAlgo<EEDataFrame>::CalculatedRecHit crh =
0532 ratioMethod_endcap_.getCalculatedRecHit();
0533 double theTimeCorrectionEE = timeCorrection(
0534 uncalibRecHit.amplitude(), timeCorrBias_->EETimeCorrAmplitudeBins, timeCorrBias_->EETimeCorrShiftBins);
0535
0536 uncalibRecHit.setJitter(crh.timeMax - 5 + theTimeCorrectionEE);
0537 uncalibRecHit.setJitterError(
0538 std::sqrt(std::pow(crh.timeError, 2) + std::pow(EEtimeConstantTerm_ * invClockToNs, 2)));
0539
0540
0541 if (uncalibRecHit.amplitude() > pedRMSVec[0] * amplitudeThreshEE_) {
0542 float outOfTimeThreshP = outOfTimeThreshG12pEE_;
0543 float outOfTimeThreshM = outOfTimeThreshG12mEE_;
0544
0545
0546
0547 if (uncalibRecHit.amplitude() > 3000.) {
0548 for (int iSample = 0; iSample < EEDataFrame::MAXSAMPLES; iSample++) {
0549 int GainId = ((EcalDataFrame)(*itdg)).sample(iSample).gainId();
0550 if (GainId != 1) {
0551 outOfTimeThreshP = outOfTimeThreshG61pEE_;
0552 outOfTimeThreshM = outOfTimeThreshG61mEE_;
0553 break;
0554 }
0555 }
0556 }
0557 float correctedTime = (crh.timeMax - 5) * clockToNsConstant + itimeconst + offsetTime;
0558 float cterm = EEtimeConstantTerm_;
0559 float sigmaped = pedRMSVec[0];
0560 float nterm = EEtimeNconst_ * sigmaped / uncalibRecHit.amplitude();
0561 float sigmat = std::sqrt(nterm * nterm + cterm * cterm);
0562 if ((correctedTime > sigmat * outOfTimeThreshP) || (correctedTime < -sigmat * outOfTimeThreshM)) {
0563 uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kOutOfTime);
0564 }
0565 }
0566
0567 } else {
0568 ratioMethod_barrel_.init(*itdg, *sampleMask_, pedVec, pedRMSVec, gainRatios);
0569 ratioMethod_barrel_.fixMGPAslew(*itdg);
0570 ratioMethod_barrel_.computeTime(EBtimeFitParameters_, EBtimeFitLimits_, EBamplitudeFitParameters_);
0571 ratioMethod_barrel_.computeAmplitude(EBamplitudeFitParameters_);
0572 EcalUncalibRecHitRatioMethodAlgo<EBDataFrame>::CalculatedRecHit crh =
0573 ratioMethod_barrel_.getCalculatedRecHit();
0574
0575 double theTimeCorrectionEB = timeCorrection(
0576 uncalibRecHit.amplitude(), timeCorrBias_->EBTimeCorrAmplitudeBins, timeCorrBias_->EBTimeCorrShiftBins);
0577
0578 uncalibRecHit.setJitter(crh.timeMax - 5 + theTimeCorrectionEB);
0579 uncalibRecHit.setJitterError(std::hypot(crh.timeError, EBtimeConstantTerm_ / clockToNsConstant));
0580
0581
0582 if (uncalibRecHit.amplitude() > pedRMSVec[0] * amplitudeThreshEB_) {
0583 float outOfTimeThreshP = outOfTimeThreshG12pEB_;
0584 float outOfTimeThreshM = outOfTimeThreshG12mEB_;
0585
0586
0587
0588 if (uncalibRecHit.amplitude() > 3000.) {
0589 for (int iSample = 0; iSample < EBDataFrame::MAXSAMPLES; iSample++) {
0590 int GainId = ((EcalDataFrame)(*itdg)).sample(iSample).gainId();
0591 if (GainId != 1) {
0592 outOfTimeThreshP = outOfTimeThreshG61pEB_;
0593 outOfTimeThreshM = outOfTimeThreshG61mEB_;
0594 break;
0595 }
0596 }
0597 }
0598 float correctedTime = (crh.timeMax - 5) * clockToNsConstant + itimeconst + offsetTime;
0599 float cterm = EBtimeConstantTerm_;
0600 float sigmaped = pedRMSVec[0];
0601 float nterm = EBtimeNconst_ * sigmaped / uncalibRecHit.amplitude();
0602 float sigmat = std::sqrt(nterm * nterm + cterm * cterm);
0603 if ((correctedTime > sigmat * outOfTimeThreshP) || (correctedTime < -sigmat * outOfTimeThreshM)) {
0604 uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kOutOfTime);
0605 }
0606 }
0607 }
0608 } else if (timealgo_ == weightsMethod) {
0609
0610 std::vector<double> amplitudes;
0611 for (unsigned int ibx = 0; ibx < activeBX.size(); ++ibx)
0612 amplitudes.push_back(uncalibRecHit.outOfTimeAmplitude(ibx));
0613
0614 EcalTBWeights::EcalTDCId tdcid(1);
0615 EcalTBWeights::EcalTBWeightMap const& wgtsMap = wgts->getMap();
0616 EcalTBWeights::EcalTBWeightMap::const_iterator wit;
0617 wit = wgtsMap.find(std::make_pair(*gid, tdcid));
0618 if (wit == wgtsMap.end()) {
0619 edm::LogError("EcalUncalibRecHitError")
0620 << "No weights found for EcalGroupId: " << gid->id() << " and EcalTDCId: " << tdcid
0621 << "\n skipping digi with id: " << detid.rawId();
0622 result.pop_back();
0623 continue;
0624 }
0625 const EcalWeightSet& wset = wit->second;
0626
0627 const EcalWeightSet::EcalWeightMatrix& mat1 = wset.getWeightsBeforeGainSwitch();
0628 const EcalWeightSet::EcalWeightMatrix& mat2 = wset.getWeightsAfterGainSwitch();
0629
0630 weights[0] = &mat1;
0631 weights[1] = &mat2;
0632
0633 double timerh;
0634 if (detid.subdetId() == EcalEndcap) {
0635 timerh = weightsMethod_endcap_.time(*itdg, amplitudes, aped, aGain, fullpulse, weights);
0636 } else {
0637 timerh = weightsMethod_barrel_.time(*itdg, amplitudes, aped, aGain, fullpulse, weights);
0638 }
0639 uncalibRecHit.setJitter(timerh);
0640 uncalibRecHit.setJitterError(0.);
0641
0642 } else if (timealgo_ == crossCorrelationMethod) {
0643 std::vector<double> amplitudes(activeBX.size());
0644 for (unsigned int ibx = 0; ibx < activeBX.size(); ++ibx)
0645 amplitudes[ibx] = uncalibRecHit.outOfTimeAmplitude(ibx);
0646
0647 bool const doSlewCorrection =
0648 barrel ? crossCorrelationUseSlewCorrectionEB_ : crossCorrelationUseSlewCorrectionEE_;
0649
0650 float jitter = computeCC_->computeTimeCC(
0651 *itdg, amplitudes, aped, aGain, fullpulse, CCtargetTimePrecision_, true, doSlewCorrection) +
0652 CCTimeShiftWrtRations_ / ecalcctiming::clockToNS;
0653 float noCorrectedJitter = computeCC_->computeTimeCC(*itdg,
0654 amplitudes,
0655 aped,
0656 aGain,
0657 fullpulse,
0658 CCtargetTimePrecisionForDelayedPulses_,
0659 false,
0660 doSlewCorrection) +
0661 CCTimeShiftWrtRations_ / ecalcctiming::clockToNS;
0662
0663 uncalibRecHit.setJitter(jitter);
0664 uncalibRecHit.setNonCorrectedTime(jitter, noCorrectedJitter);
0665
0666 float retreivedNonCorrectedTime = uncalibRecHit.nonCorrectedTime();
0667 float noCorrectedTime = ecalcctiming::clockToNS * noCorrectedJitter;
0668 if (retreivedNonCorrectedTime > -29.0 && std::abs(retreivedNonCorrectedTime - noCorrectedTime) > 0.05) {
0669 edm::LogError("EcalUncalibRecHitError") << "Problem with noCorrectedJitter: true value:" << noCorrectedTime
0670 << "\t received: " << retreivedNonCorrectedTime << std::endl;
0671 }
0672
0673
0674 float threshold, cterm, timeNconst;
0675 float timeThrP = 0.;
0676 float timeThrM = 0.;
0677 if (barrel) {
0678 threshold = pedRMSVec[0] * amplitudeThreshEB_;
0679 cterm = EBtimeConstantTerm_;
0680 timeNconst = EBtimeNconst_;
0681 timeThrP = outOfTimeThreshG12pEB_;
0682 timeThrM = outOfTimeThreshG12mEB_;
0683 if (uncalibRecHit.amplitude() > 3000.) {
0684 for (int iSample = 0; iSample < EBDataFrame::MAXSAMPLES; iSample++) {
0685 int GainId = ((EcalDataFrame)(*itdg)).sample(iSample).gainId();
0686 if (GainId != 1) {
0687 timeThrP = outOfTimeThreshG61pEB_;
0688 timeThrM = outOfTimeThreshG61mEB_;
0689 break;
0690 }
0691 }
0692 }
0693 } else {
0694 threshold = pedRMSVec[0] * amplitudeThreshEE_;
0695 cterm = EEtimeConstantTerm_;
0696 timeNconst = EEtimeNconst_;
0697 timeThrP = outOfTimeThreshG12pEE_;
0698 timeThrM = outOfTimeThreshG12mEE_;
0699 if (uncalibRecHit.amplitude() > 3000.) {
0700 for (int iSample = 0; iSample < EEDataFrame::MAXSAMPLES; iSample++) {
0701 int GainId = ((EcalDataFrame)(*itdg)).sample(iSample).gainId();
0702 if (GainId != 1) {
0703 timeThrP = outOfTimeThreshG61pEE_;
0704 timeThrM = outOfTimeThreshG61mEE_;
0705 break;
0706 }
0707 }
0708 }
0709 }
0710 if (uncalibRecHit.amplitude() > threshold) {
0711 float correctedTime = noCorrectedJitter * ecalcctiming::clockToNS + itimeconst + offsetTime;
0712 float sigmaped = pedRMSVec[0];
0713 float nterm = timeNconst * sigmaped / uncalibRecHit.amplitude();
0714 float sigmat = std::sqrt(nterm * nterm + cterm * cterm);
0715 if ((correctedTime > sigmat * timeThrP) || (correctedTime < -sigmat * timeThrM))
0716 uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kOutOfTime);
0717 }
0718
0719 } else {
0720 uncalibRecHit.setJitter(0.);
0721 uncalibRecHit.setJitterError(0.);
0722 }
0723 }
0724
0725
0726 auto& uncalibRecHit = result.back();
0727 if (((EcalDataFrame)(*itdg)).hasSwitchToGain6())
0728 uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kHasSwitchToGain6);
0729 if (((EcalDataFrame)(*itdg)).hasSwitchToGain1())
0730 uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kHasSwitchToGain1);
0731 }
0732 }
0733
0734 edm::ParameterSetDescription EcalUncalibRecHitWorkerMultiFit::getAlgoDescription() {
0735 edm::ParameterSetDescription psd;
0736 psd.addNode(edm::ParameterDescription<std::vector<int>>("activeBXs", {-5, -4, -3, -2, -1, 0, 1, 2, 3, 4}, true) and
0737 edm::ParameterDescription<bool>("ampErrorCalculation", true, true) and
0738 edm::ParameterDescription<bool>("useLumiInfoRunHeader", true, true) and
0739 edm::ParameterDescription<int>("bunchSpacing", 0, true) and
0740 edm::ParameterDescription<bool>("doPrefitEB", false, true) and
0741 edm::ParameterDescription<bool>("doPrefitEE", false, true) and
0742 edm::ParameterDescription<double>("prefitMaxChiSqEB", 25., true) and
0743 edm::ParameterDescription<double>("prefitMaxChiSqEE", 10., true) and
0744 edm::ParameterDescription<bool>("dynamicPedestalsEB", false, true) and
0745 edm::ParameterDescription<bool>("dynamicPedestalsEE", false, true) and
0746 edm::ParameterDescription<bool>("mitigateBadSamplesEB", false, true) and
0747 edm::ParameterDescription<bool>("mitigateBadSamplesEE", false, true) and
0748 edm::ParameterDescription<bool>("gainSwitchUseMaxSampleEB", true, true) and
0749 edm::ParameterDescription<bool>("gainSwitchUseMaxSampleEE", false, true) and
0750 edm::ParameterDescription<bool>("selectiveBadSampleCriteriaEB", false, true) and
0751 edm::ParameterDescription<bool>("selectiveBadSampleCriteriaEE", false, true) and
0752 edm::ParameterDescription<double>("addPedestalUncertaintyEB", 0., true) and
0753 edm::ParameterDescription<double>("addPedestalUncertaintyEE", 0., true) and
0754 edm::ParameterDescription<bool>("simplifiedNoiseModelForGainSwitch", true, true) and
0755 edm::ParameterDescription<std::string>("timealgo", "RatioMethod", true) and
0756 edm::ParameterDescription<std::vector<double>>("EBtimeFitParameters",
0757 {-2.015452e+00,
0758 3.130702e+00,
0759 -1.234730e+01,
0760 4.188921e+01,
0761 -8.283944e+01,
0762 9.101147e+01,
0763 -5.035761e+01,
0764 1.105621e+01},
0765 true) and
0766 edm::ParameterDescription<std::vector<double>>("EEtimeFitParameters",
0767 {-2.390548e+00,
0768 3.553628e+00,
0769 -1.762341e+01,
0770 6.767538e+01,
0771 -1.332130e+02,
0772 1.407432e+02,
0773 -7.541106e+01,
0774 1.620277e+01},
0775 true) and
0776 edm::ParameterDescription<std::vector<double>>("EBamplitudeFitParameters", {1.138, 1.652}, true) and
0777 edm::ParameterDescription<std::vector<double>>("EEamplitudeFitParameters", {1.890, 1.400}, true) and
0778 edm::ParameterDescription<edm::ESInputTag>("timeCalibTag", edm::ESInputTag(), true) and
0779 edm::ParameterDescription<edm::ESInputTag>("timeOffsetTag", edm::ESInputTag(), true) and
0780 edm::ParameterDescription<double>("EBtimeFitLimits_Lower", 0.2, true) and
0781 edm::ParameterDescription<double>("EBtimeFitLimits_Upper", 1.4, true) and
0782 edm::ParameterDescription<double>("EEtimeFitLimits_Lower", 0.2, true) and
0783 edm::ParameterDescription<double>("EEtimeFitLimits_Upper", 1.4, true) and
0784 edm::ParameterDescription<double>("EBtimeConstantTerm", .6, true) and
0785 edm::ParameterDescription<double>("EEtimeConstantTerm", 1.0, true) and
0786 edm::ParameterDescription<double>("EBtimeNconst", 28.5, true) and
0787 edm::ParameterDescription<double>("EEtimeNconst", 31.8, true) and
0788 edm::ParameterDescription<double>("outOfTimeThresholdGain12pEB", 5., true) and
0789 edm::ParameterDescription<double>("outOfTimeThresholdGain12mEB", 5., true) and
0790 edm::ParameterDescription<double>("outOfTimeThresholdGain61pEB", 5., true) and
0791 edm::ParameterDescription<double>("outOfTimeThresholdGain61mEB", 5., true) and
0792 edm::ParameterDescription<double>("outOfTimeThresholdGain12pEE", 1000, true) and
0793 edm::ParameterDescription<double>("outOfTimeThresholdGain12mEE", 1000, true) and
0794 edm::ParameterDescription<double>("outOfTimeThresholdGain61pEE", 1000, true) and
0795 edm::ParameterDescription<double>("outOfTimeThresholdGain61mEE", 1000, true) and
0796 edm::ParameterDescription<double>("amplitudeThresholdEB", 10, true) and
0797 edm::ParameterDescription<double>("amplitudeThresholdEE", 10, true) and
0798 edm::ParameterDescription<bool>("crossCorrelationUseSlewCorrectionEB", true, true) and
0799 edm::ParameterDescription<bool>("crossCorrelationUseSlewCorrectionEE", false, true) and
0800 edm::ParameterDescription<double>("crossCorrelationStartTime", -25.0, true) and
0801 edm::ParameterDescription<double>("crossCorrelationStopTime", 25.0, true) and
0802 edm::ParameterDescription<double>("crossCorrelationTargetTimePrecision", 0.01, true) and
0803 edm::ParameterDescription<double>("crossCorrelationTargetTimePrecisionForDelayedPulses", 0.05, true) and
0804 edm::ParameterDescription<double>("crossCorrelationTimeShiftWrtRations", 0., true) and
0805 edm::ParameterDescription<double>("crossCorrelationMinTimeToBeLateMin", 2., true) and
0806 edm::ParameterDescription<double>("crossCorrelationMinTimeToBeLateMax", 5., true));
0807
0808 return psd;
0809 }
0810
0811 #include "FWCore/Framework/interface/MakerMacros.h"
0812 #include "RecoLocalCalo/EcalRecProducers/interface/EcalUncalibRecHitWorkerFactory.h"
0813 DEFINE_EDM_PLUGIN(EcalUncalibRecHitWorkerFactory, EcalUncalibRecHitWorkerMultiFit, "EcalUncalibRecHitWorkerMultiFit");
0814 #include "RecoLocalCalo/EcalRecProducers/interface/EcalUncalibRecHitFillDescriptionWorkerFactory.h"
0815 DEFINE_EDM_PLUGIN(EcalUncalibRecHitFillDescriptionWorkerFactory,
0816 EcalUncalibRecHitWorkerMultiFit,
0817 "EcalUncalibRecHitWorkerMultiFit");