Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:44 

0001 /** \class EcalUncalibRecHitRecGlobalAlgo
0002  *  Template used to compute amplitude, pedestal using a weights method
0003  *                           time using a ratio method
0004  *                           chi2 using express method
0005  *
0006  *  \author R. Bruneliere - A. Zabi
0007  */
0008 
0009 #include "CondFormats/DataRecord/interface/EcalGainRatiosRcd.h"
0010 #include "CondFormats/DataRecord/interface/EcalPedestalsRcd.h"
0011 #include "CondFormats/DataRecord/interface/EcalSampleMaskRcd.h"
0012 #include "CondFormats/DataRecord/interface/EcalTBWeightsRcd.h"
0013 #include "CondFormats/DataRecord/interface/EcalTimeBiasCorrectionsRcd.h"
0014 #include "CondFormats/DataRecord/interface/EcalTimeCalibConstantsRcd.h"
0015 #include "CondFormats/DataRecord/interface/EcalTimeOffsetConstantRcd.h"
0016 #include "CondFormats/DataRecord/interface/EcalWeightXtalGroupsRcd.h"
0017 #include "CondFormats/EcalObjects/interface/EcalGainRatios.h"
0018 #include "CondFormats/EcalObjects/interface/EcalPedestals.h"
0019 #include "CondFormats/EcalObjects/interface/EcalSampleMask.h"
0020 #include "CondFormats/EcalObjects/interface/EcalTBWeights.h"
0021 #include "CondFormats/EcalObjects/interface/EcalTimeBiasCorrections.h"
0022 #include "CondFormats/EcalObjects/interface/EcalTimeCalibConstants.h"
0023 #include "CondFormats/EcalObjects/interface/EcalTimeOffsetConstant.h"
0024 #include "CondFormats/EcalObjects/interface/EcalWeightXtalGroups.h"
0025 #include "FWCore/Framework/interface/ESHandle.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/EventSetup.h"
0028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0029 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0032 #include "FWCore/Utilities/interface/ESGetToken.h"
0033 #include "FWCore/Utilities/interface/ESInputTag.h"
0034 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRatioMethodAlgo.h"
0035 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRecChi2Algo.h"
0036 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRecWeightsAlgo.h"
0037 #include "RecoLocalCalo/EcalRecProducers/interface/EcalUncalibRecHitWorkerRunOneDigiBase.h"
0038 #include "SimCalorimetry/EcalSimAlgos/interface/EBShape.h"
0039 #include "SimCalorimetry/EcalSimAlgos/interface/EEShape.h"
0040 
0041 class EcalUncalibRecHitWorkerGlobal : public EcalUncalibRecHitWorkerRunOneDigiBase {
0042 public:
0043   EcalUncalibRecHitWorkerGlobal(const edm::ParameterSet&, edm::ConsumesCollector& c);
0044   EcalUncalibRecHitWorkerGlobal() = default;  // for EcalUncalibRecHitFillDescriptionWorkerFactory
0045 
0046   void set(const edm::EventSetup& es) override;
0047   bool run(const edm::Event& evt,
0048            const EcalDigiCollection::const_iterator& digi,
0049            EcalUncalibratedRecHitCollection& result) override;
0050 
0051   edm::ParameterSetDescription getAlgoDescription() override;
0052 
0053 protected:
0054   double pedVec[3];
0055   double pedRMSVec[3];
0056   double gainRatios[3];
0057 
0058   edm::ESGetToken<EcalPedestals, EcalPedestalsRcd> tokenPeds_;
0059   edm::ESGetToken<EcalGainRatios, EcalGainRatiosRcd> tokenGains_;
0060   edm::ESHandle<EcalPedestals> peds_;
0061   edm::ESHandle<EcalGainRatios> gains_;
0062 
0063   template <class C>
0064   int isSaturated(const C& digi);
0065 
0066   double timeCorrection(float ampli, const std::vector<float>& amplitudeBins, const std::vector<float>& shiftBins);
0067 
0068   // weights method
0069   edm::ESGetToken<EcalWeightXtalGroups, EcalWeightXtalGroupsRcd> tokenGrps_;
0070   edm::ESGetToken<EcalTBWeights, EcalTBWeightsRcd> tokenWgts_;
0071   edm::ESHandle<EcalWeightXtalGroups> grps_;
0072   edm::ESHandle<EcalTBWeights> wgts_;
0073   const EcalWeightSet::EcalWeightMatrix* weights[2];
0074   const EcalWeightSet::EcalChi2WeightMatrix* chi2mat[2];
0075   EcalUncalibRecHitRecWeightsAlgo<EBDataFrame> weightsMethod_barrel_;
0076   EcalUncalibRecHitRecWeightsAlgo<EEDataFrame> weightsMethod_endcap_;
0077   EEShape testbeamEEShape;  // used in the chi2
0078   EBShape testbeamEBShape;  // can be replaced by simple shape arrays of float in the future
0079 
0080   // determie which of the samples must actually be used by ECAL local reco
0081   edm::ESGetToken<EcalSampleMask, EcalSampleMaskRcd> tokenSampleMask_;
0082   edm::ESHandle<EcalSampleMask> sampleMaskHand_;
0083 
0084   // ratio method
0085   std::vector<double> EBtimeFitParameters_;
0086   std::vector<double> EEtimeFitParameters_;
0087   std::vector<double> EBamplitudeFitParameters_;
0088   std::vector<double> EEamplitudeFitParameters_;
0089   std::pair<double, double> EBtimeFitLimits_;
0090   std::pair<double, double> EEtimeFitLimits_;
0091 
0092   EcalUncalibRecHitRatioMethodAlgo<EBDataFrame> ratioMethod_barrel_;
0093   EcalUncalibRecHitRatioMethodAlgo<EEDataFrame> ratioMethod_endcap_;
0094 
0095   double EBtimeConstantTerm_;
0096   double EBtimeNconst_;
0097   double EEtimeConstantTerm_;
0098   double EEtimeNconst_;
0099   double outOfTimeThreshG12pEB_;
0100   double outOfTimeThreshG12mEB_;
0101   double outOfTimeThreshG61pEB_;
0102   double outOfTimeThreshG61mEB_;
0103   double outOfTimeThreshG12pEE_;
0104   double outOfTimeThreshG12mEE_;
0105   double outOfTimeThreshG61pEE_;
0106   double outOfTimeThreshG61mEE_;
0107   double amplitudeThreshEB_;
0108   double amplitudeThreshEE_;
0109   double ebSpikeThresh_;
0110 
0111   edm::ESGetToken<EcalTimeBiasCorrections, EcalTimeBiasCorrectionsRcd> tokenTimeCorrBias_;
0112   edm::ESHandle<EcalTimeBiasCorrections> timeCorrBias_;
0113 
0114   edm::ESGetToken<EcalTimeCalibConstants, EcalTimeCalibConstantsRcd> tokenItime_;
0115   edm::ESGetToken<EcalTimeOffsetConstant, EcalTimeOffsetConstantRcd> tokenOfftime_;
0116   edm::ESHandle<EcalTimeCalibConstants> itime_;
0117   edm::ESHandle<EcalTimeOffsetConstant> offtime_;
0118   std::vector<double> ebPulseShape_;
0119   std::vector<double> eePulseShape_;
0120 
0121   // chi2 method
0122   bool kPoorRecoFlagEB_;
0123   bool kPoorRecoFlagEE_;
0124   double chi2ThreshEB_;
0125   double chi2ThreshEE_;
0126   std::vector<double> EBchi2Parameters_;
0127   std::vector<double> EEchi2Parameters_;
0128 };
0129 
0130 EcalUncalibRecHitWorkerGlobal::EcalUncalibRecHitWorkerGlobal(const edm::ParameterSet& ps, edm::ConsumesCollector& c)
0131     : EcalUncalibRecHitWorkerRunOneDigiBase(ps, c),
0132       tokenPeds_(c.esConsumes<EcalPedestals, EcalPedestalsRcd>()),
0133       tokenGains_(c.esConsumes<EcalGainRatios, EcalGainRatiosRcd>()),
0134       tokenGrps_(c.esConsumes<EcalWeightXtalGroups, EcalWeightXtalGroupsRcd>()),
0135       tokenWgts_(c.esConsumes<EcalTBWeights, EcalTBWeightsRcd>()),
0136       testbeamEEShape(c),
0137       testbeamEBShape(c),
0138       tokenSampleMask_(c.esConsumes<EcalSampleMask, EcalSampleMaskRcd>()),
0139       tokenTimeCorrBias_(c.esConsumes<EcalTimeBiasCorrections, EcalTimeBiasCorrectionsRcd>()),
0140       tokenItime_(c.esConsumes<EcalTimeCalibConstants, EcalTimeCalibConstantsRcd>(
0141           ps.getParameter<edm::ESInputTag>("timeCalibTag"))),
0142       tokenOfftime_(c.esConsumes<EcalTimeOffsetConstant, EcalTimeOffsetConstantRcd>(
0143           ps.getParameter<edm::ESInputTag>("timeOffsetTag"))) {
0144   // ratio method parameters
0145   EBtimeFitParameters_ = ps.getParameter<std::vector<double>>("EBtimeFitParameters");
0146   EEtimeFitParameters_ = ps.getParameter<std::vector<double>>("EEtimeFitParameters");
0147   EBamplitudeFitParameters_ = ps.getParameter<std::vector<double>>("EBamplitudeFitParameters");
0148   EEamplitudeFitParameters_ = ps.getParameter<std::vector<double>>("EEamplitudeFitParameters");
0149   EBtimeFitLimits_.first = ps.getParameter<double>("EBtimeFitLimits_Lower");
0150   EBtimeFitLimits_.second = ps.getParameter<double>("EBtimeFitLimits_Upper");
0151   EEtimeFitLimits_.first = ps.getParameter<double>("EEtimeFitLimits_Lower");
0152   EEtimeFitLimits_.second = ps.getParameter<double>("EEtimeFitLimits_Upper");
0153   EBtimeConstantTerm_ = ps.getParameter<double>("EBtimeConstantTerm");
0154   EBtimeNconst_ = ps.getParameter<double>("EBtimeNconst");
0155   EEtimeConstantTerm_ = ps.getParameter<double>("EEtimeConstantTerm");
0156   EEtimeNconst_ = ps.getParameter<double>("EEtimeNconst");
0157   outOfTimeThreshG12pEB_ = ps.getParameter<double>("outOfTimeThresholdGain12pEB");
0158   outOfTimeThreshG12mEB_ = ps.getParameter<double>("outOfTimeThresholdGain12mEB");
0159   outOfTimeThreshG61pEB_ = ps.getParameter<double>("outOfTimeThresholdGain61pEB");
0160   outOfTimeThreshG61mEB_ = ps.getParameter<double>("outOfTimeThresholdGain61mEB");
0161   outOfTimeThreshG12pEE_ = ps.getParameter<double>("outOfTimeThresholdGain12pEE");
0162   outOfTimeThreshG12mEE_ = ps.getParameter<double>("outOfTimeThresholdGain12mEE");
0163   outOfTimeThreshG61pEE_ = ps.getParameter<double>("outOfTimeThresholdGain61pEE");
0164   outOfTimeThreshG61mEE_ = ps.getParameter<double>("outOfTimeThresholdGain61mEE");
0165   amplitudeThreshEB_ = ps.getParameter<double>("amplitudeThresholdEB");
0166   amplitudeThreshEE_ = ps.getParameter<double>("amplitudeThresholdEE");
0167 
0168   // spike threshold
0169   ebSpikeThresh_ = ps.getParameter<double>("ebSpikeThreshold");
0170 
0171   ebPulseShape_ = ps.getParameter<std::vector<double>>("ebPulseShape");
0172   eePulseShape_ = ps.getParameter<std::vector<double>>("eePulseShape");
0173 
0174   // chi2 parameters
0175   kPoorRecoFlagEB_ = ps.getParameter<bool>("kPoorRecoFlagEB");
0176   kPoorRecoFlagEE_ = ps.getParameter<bool>("kPoorRecoFlagEE");
0177 
0178   chi2ThreshEB_ = ps.getParameter<double>("chi2ThreshEB_");
0179   chi2ThreshEE_ = ps.getParameter<double>("chi2ThreshEE_");
0180   EBchi2Parameters_ = ps.getParameter<std::vector<double>>("EBchi2Parameters");
0181   EEchi2Parameters_ = ps.getParameter<std::vector<double>>("EEchi2Parameters");
0182 }
0183 
0184 void EcalUncalibRecHitWorkerGlobal::set(const edm::EventSetup& es) {
0185   // common setup
0186   gains_ = es.getHandle(tokenGains_);
0187   peds_ = es.getHandle(tokenPeds_);
0188 
0189   // for the weights method
0190   grps_ = es.getHandle(tokenGrps_);
0191   wgts_ = es.getHandle(tokenWgts_);
0192 
0193   // which of the samples need be used
0194   sampleMaskHand_ = es.getHandle(tokenSampleMask_);
0195 
0196   // for the ratio method
0197 
0198   itime_ = es.getHandle(tokenItime_);
0199   offtime_ = es.getHandle(tokenOfftime_);
0200 
0201   // for the time correction methods
0202   timeCorrBias_ = es.getHandle(tokenTimeCorrBias_);
0203 
0204   // for the DB Ecal Pulse Sim Shape
0205   testbeamEEShape.setEventSetup(es);
0206   testbeamEBShape.setEventSetup(es);
0207 }
0208 
0209 // check saturation: 5 samples with gainId = 0
0210 template <class C>
0211 int EcalUncalibRecHitWorkerGlobal::isSaturated(const C& dataFrame) {
0212   //bool saturated_ = 0;
0213   int cnt;
0214   for (int j = 0; j < C::MAXSAMPLES - 5; ++j) {
0215     cnt = 0;
0216     for (int i = j; i < (j + 5) && i < C::MAXSAMPLES; ++i) {
0217       if (dataFrame.sample(i).gainId() == 0)
0218         ++cnt;
0219     }
0220     if (cnt == 5)
0221       return j - 1;  // the last unsaturated sample
0222   }
0223   return -1;  // no saturation found
0224 }
0225 
0226 /**
0227  * Amplitude-dependent time corrections; EE and EB have separate corrections:
0228  * EXtimeCorrAmplitudes (ADC) and EXtimeCorrShifts (ns) need to have the same number of elements
0229  * Bins must be ordered in amplitude. First-last bins take care of under-overflows.
0230  *
0231  * The algorithm is the same for EE and EB, only the correction vectors are different.
0232  *
0233  * @return Jitter (in clock cycles) which will be added to UncalibRechit.setJitter(), 0 if no correction is applied.
0234  */
0235 double EcalUncalibRecHitWorkerGlobal::timeCorrection(float ampli,
0236                                                      const std::vector<float>& amplitudeBins,
0237                                                      const std::vector<float>& shiftBins) {
0238   // computed initially in ns. Than turned in the BX's, as
0239   // EcalUncalibratedRecHit need be.
0240   double theCorrection = 0;
0241 
0242   // sanity check for arrays
0243   if (amplitudeBins.empty()) {
0244     edm::LogError("EcalRecHitError") << "timeCorrAmplitudeBins is empty, forcing no time bias corrections.";
0245 
0246     return 0;
0247   }
0248 
0249   if (amplitudeBins.size() != shiftBins.size()) {
0250     edm::LogError("EcalRecHitError") << "Size of timeCorrAmplitudeBins different from "
0251                                         "timeCorrShiftBins. Forcing no time bias corrections. ";
0252 
0253     return 0;
0254   }
0255 
0256   int myBin = -1;
0257   for (int bin = 0; bin < (int)amplitudeBins.size(); bin++) {
0258     if (ampli > amplitudeBins.at(bin)) {
0259       myBin = bin;
0260     } else {
0261       break;
0262     }
0263   }
0264 
0265   if (myBin == -1) {
0266     theCorrection = shiftBins.at(0);
0267   } else if (myBin == ((int)(amplitudeBins.size() - 1))) {
0268     theCorrection = shiftBins.at(myBin);
0269   } else if (-1 < myBin && myBin < ((int)amplitudeBins.size() - 1)) {
0270     // interpolate linearly between two assingned points
0271     theCorrection = (shiftBins.at(myBin + 1) - shiftBins.at(myBin));
0272     theCorrection *=
0273         (((double)ampli) - amplitudeBins.at(myBin)) / (amplitudeBins.at(myBin + 1) - amplitudeBins.at(myBin));
0274     theCorrection += shiftBins.at(myBin);
0275   } else {
0276     edm::LogError("EcalRecHitError") << "Assigning time correction impossible. Setting it to 0 ";
0277     theCorrection = 0.;
0278   }
0279 
0280   // convert ns into clocks
0281   return theCorrection / 25.;
0282 }
0283 
0284 bool EcalUncalibRecHitWorkerGlobal::run(const edm::Event& evt,
0285                                         const EcalDigiCollection::const_iterator& itdg,
0286                                         EcalUncalibratedRecHitCollection& result) {
0287   DetId detid(itdg->id());
0288 
0289   const EcalSampleMask* sampleMask_ = sampleMaskHand_.product();
0290 
0291   // intelligence for recHit computation
0292   EcalUncalibratedRecHit uncalibRecHit;
0293 
0294   const EcalPedestals::Item* aped = nullptr;
0295   const EcalMGPAGainRatio* aGain = nullptr;
0296   const EcalXtalGroupId* gid = nullptr;
0297   float offsetTime = 0;
0298 
0299   if (detid.subdetId() == EcalEndcap) {
0300     unsigned int hashedIndex = EEDetId(detid).hashedIndex();
0301     aped = &peds_->endcap(hashedIndex);
0302     aGain = &gains_->endcap(hashedIndex);
0303     gid = &grps_->endcap(hashedIndex);
0304     offsetTime = offtime_->getEEValue();
0305   } else {
0306     unsigned int hashedIndex = EBDetId(detid).hashedIndex();
0307     aped = &peds_->barrel(hashedIndex);
0308     aGain = &gains_->barrel(hashedIndex);
0309     gid = &grps_->barrel(hashedIndex);
0310     offsetTime = offtime_->getEBValue();
0311   }
0312 
0313   pedVec[0] = aped->mean_x12;
0314   pedVec[1] = aped->mean_x6;
0315   pedVec[2] = aped->mean_x1;
0316   pedRMSVec[0] = aped->rms_x12;
0317   pedRMSVec[1] = aped->rms_x6;
0318   pedRMSVec[2] = aped->rms_x1;
0319   gainRatios[0] = 1.;
0320   gainRatios[1] = aGain->gain12Over6();
0321   gainRatios[2] = aGain->gain6Over1() * aGain->gain12Over6();
0322 
0323   // compute the right bin of the pulse shape using time calibration constants
0324   EcalTimeCalibConstantMap::const_iterator it = itime_->find(detid);
0325   EcalTimeCalibConstant itimeconst = 0;
0326   if (it != itime_->end()) {
0327     itimeconst = (*it);
0328   } else {
0329     edm::LogError("EcalRecHitError") << "No time intercalib const found for xtal " << detid.rawId()
0330                                      << "! something wrong with EcalTimeCalibConstants in your DB? ";
0331   }
0332 
0333   // === amplitude computation ===
0334   int leadingSample = -1;
0335   if (detid.subdetId() == EcalEndcap) {
0336     leadingSample = ((EcalDataFrame)(*itdg)).lastUnsaturatedSample();
0337   } else {
0338     leadingSample = ((EcalDataFrame)(*itdg)).lastUnsaturatedSample();
0339   }
0340 
0341   if (leadingSample == 4) {  // saturation on the expected max sample
0342     uncalibRecHit = EcalUncalibratedRecHit((*itdg).id(), 4095 * 12, 0, 0, 0);
0343     uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kSaturated);
0344     // do not propagate the default chi2 = -1 value to the calib rechit (mapped to 64), set it to 0 when saturation
0345     uncalibRecHit.setChi2(0);
0346   } else if (leadingSample >= 0) {  // saturation on other samples: cannot extrapolate from the fourth one
0347     double pedestal = 0.;
0348     double gainratio = 1.;
0349     int gainId = ((EcalDataFrame)(*itdg)).sample(5).gainId();
0350 
0351     if (gainId == 0 || gainId == 3) {
0352       pedestal = aped->mean_x1;
0353       gainratio = aGain->gain6Over1() * aGain->gain12Over6();
0354     } else if (gainId == 1) {
0355       pedestal = aped->mean_x12;
0356       gainratio = 1.;
0357     } else if (gainId == 2) {
0358       pedestal = aped->mean_x6;
0359       gainratio = aGain->gain12Over6();
0360     }
0361     double amplitude = ((double)(((EcalDataFrame)(*itdg)).sample(5).adc()) - pedestal) * gainratio;
0362     uncalibRecHit = EcalUncalibratedRecHit((*itdg).id(), amplitude, 0, 0, 0);
0363     uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kSaturated);
0364     // do not propagate the default chi2 = -1 value to the calib rechit (mapped to 64), set it to 0 when saturation
0365     uncalibRecHit.setChi2(0);
0366   } else {
0367     // weights method
0368     EcalTBWeights::EcalTDCId tdcid(1);
0369     EcalTBWeights::EcalTBWeightMap const& wgtsMap = wgts_->getMap();
0370     EcalTBWeights::EcalTBWeightMap::const_iterator wit;
0371     wit = wgtsMap.find(std::make_pair(*gid, tdcid));
0372     if (wit == wgtsMap.end()) {
0373       edm::LogError("EcalUncalibRecHitError")
0374           << "No weights found for EcalGroupId: " << gid->id() << " and  EcalTDCId: " << tdcid
0375           << "\n  skipping digi with id: " << detid.rawId();
0376 
0377       return false;
0378     }
0379     const EcalWeightSet& wset = wit->second;  // this is the EcalWeightSet
0380 
0381     const EcalWeightSet::EcalWeightMatrix& mat1 = wset.getWeightsBeforeGainSwitch();
0382     const EcalWeightSet::EcalWeightMatrix& mat2 = wset.getWeightsAfterGainSwitch();
0383 
0384     weights[0] = &mat1;
0385     weights[1] = &mat2;
0386 
0387     // get uncalibrated recHit from weights
0388     if (detid.subdetId() == EcalEndcap) {
0389       uncalibRecHit = weightsMethod_endcap_.makeRecHit(*itdg, pedVec, pedRMSVec, gainRatios, weights, testbeamEEShape);
0390     } else {
0391       uncalibRecHit = weightsMethod_barrel_.makeRecHit(*itdg, pedVec, pedRMSVec, gainRatios, weights, testbeamEBShape);
0392     }
0393 
0394     // === time computation ===
0395     // ratio method
0396     float const clockToNsConstant = 25.;
0397     if (detid.subdetId() == EcalEndcap) {
0398       ratioMethod_endcap_.init(*itdg, *sampleMask_, pedVec, pedRMSVec, gainRatios);
0399       ratioMethod_endcap_.computeTime(EEtimeFitParameters_, EEtimeFitLimits_, EEamplitudeFitParameters_);
0400       ratioMethod_endcap_.computeAmplitude(EEamplitudeFitParameters_);
0401       EcalUncalibRecHitRatioMethodAlgo<EEDataFrame>::CalculatedRecHit crh = ratioMethod_endcap_.getCalculatedRecHit();
0402       double theTimeCorrectionEE = timeCorrection(
0403           uncalibRecHit.amplitude(), timeCorrBias_->EETimeCorrAmplitudeBins, timeCorrBias_->EETimeCorrShiftBins);
0404 
0405       uncalibRecHit.setJitter(crh.timeMax - 5 + theTimeCorrectionEE);
0406       uncalibRecHit.setJitterError(
0407           std::sqrt(pow(crh.timeError, 2) + std::pow(EEtimeConstantTerm_, 2) / std::pow(clockToNsConstant, 2)));
0408       // consider flagging as kOutOfTime only if above noise
0409       if (uncalibRecHit.amplitude() > pedRMSVec[0] * amplitudeThreshEE_) {
0410         float outOfTimeThreshP = outOfTimeThreshG12pEE_;
0411         float outOfTimeThreshM = outOfTimeThreshG12mEE_;
0412         // determine if gain has switched away from gainId==1 (x12 gain)
0413         // and determine cuts (number of 'sigmas') to ose for kOutOfTime
0414         // >3k ADC is necessasry condition for gain switch to occur
0415         if (uncalibRecHit.amplitude() > 3000.) {
0416           for (int iSample = 0; iSample < EEDataFrame::MAXSAMPLES; iSample++) {
0417             int GainId = ((EcalDataFrame)(*itdg)).sample(iSample).gainId();
0418             if (GainId != 1) {
0419               outOfTimeThreshP = outOfTimeThreshG61pEE_;
0420               outOfTimeThreshM = outOfTimeThreshG61mEE_;
0421               break;
0422             }
0423           }
0424         }
0425         float correctedTime = (crh.timeMax - 5) * clockToNsConstant + itimeconst + offsetTime;
0426         float cterm = EEtimeConstantTerm_;
0427         float sigmaped = pedRMSVec[0];  // approx for lower gains
0428         float nterm = EEtimeNconst_ * sigmaped / uncalibRecHit.amplitude();
0429         float sigmat = std::sqrt(nterm * nterm + cterm * cterm);
0430         if ((correctedTime > sigmat * outOfTimeThreshP) || (correctedTime < (-1. * sigmat * outOfTimeThreshM))) {
0431           uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kOutOfTime);
0432         }
0433       }
0434 
0435     } else {
0436       ratioMethod_barrel_.init(*itdg, *sampleMask_, pedVec, pedRMSVec, gainRatios);
0437       ratioMethod_barrel_.fixMGPAslew(*itdg);
0438       ratioMethod_barrel_.computeTime(EBtimeFitParameters_, EBtimeFitLimits_, EBamplitudeFitParameters_);
0439       ratioMethod_barrel_.computeAmplitude(EBamplitudeFitParameters_);
0440       EcalUncalibRecHitRatioMethodAlgo<EBDataFrame>::CalculatedRecHit crh = ratioMethod_barrel_.getCalculatedRecHit();
0441 
0442       double theTimeCorrectionEB = timeCorrection(
0443           uncalibRecHit.amplitude(), timeCorrBias_->EBTimeCorrAmplitudeBins, timeCorrBias_->EBTimeCorrShiftBins);
0444 
0445       uncalibRecHit.setJitter(crh.timeMax - 5 + theTimeCorrectionEB);
0446 
0447       uncalibRecHit.setJitterError(
0448           std::sqrt(std::pow(crh.timeError, 2) + std::pow(EBtimeConstantTerm_, 2) / std::pow(clockToNsConstant, 2)));
0449       // consider flagging as kOutOfTime only if above noise
0450       if (uncalibRecHit.amplitude() > pedRMSVec[0] * amplitudeThreshEB_) {
0451         float outOfTimeThreshP = outOfTimeThreshG12pEB_;
0452         float outOfTimeThreshM = outOfTimeThreshG12mEB_;
0453         // determine if gain has switched away from gainId==1 (x12 gain)
0454         // and determine cuts (number of 'sigmas') to ose for kOutOfTime
0455         // >3k ADC is necessasry condition for gain switch to occur
0456         if (uncalibRecHit.amplitude() > 3000.) {
0457           for (int iSample = 0; iSample < EBDataFrame::MAXSAMPLES; iSample++) {
0458             int GainId = ((EcalDataFrame)(*itdg)).sample(iSample).gainId();
0459             if (GainId != 1) {
0460               outOfTimeThreshP = outOfTimeThreshG61pEB_;
0461               outOfTimeThreshM = outOfTimeThreshG61mEB_;
0462               break;
0463             }
0464           }
0465         }
0466         float correctedTime = (crh.timeMax - 5) * clockToNsConstant + itimeconst + offsetTime;
0467         float cterm = EBtimeConstantTerm_;
0468         float sigmaped = pedRMSVec[0];  // approx for lower gains
0469         float nterm = EBtimeNconst_ * sigmaped / uncalibRecHit.amplitude();
0470         float sigmat = std::sqrt(nterm * nterm + cterm * cterm);
0471         if ((correctedTime > sigmat * outOfTimeThreshP) || (correctedTime < (-1. * sigmat * outOfTimeThreshM))) {
0472           uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kOutOfTime);
0473         }
0474       }
0475     }
0476 
0477     // === chi2express ===
0478     if (detid.subdetId() == EcalEndcap) {
0479       double amplitude = uncalibRecHit.amplitude();
0480       double amplitudeOutOfTime = 0.;
0481       double jitter = uncalibRecHit.jitter();
0482 
0483       EcalUncalibRecHitRecChi2Algo<EEDataFrame> chi2expressEE_(*itdg,
0484                                                                amplitude,
0485                                                                (itimeconst + offsetTime),
0486                                                                amplitudeOutOfTime,
0487                                                                jitter,
0488                                                                pedVec,
0489                                                                pedRMSVec,
0490                                                                gainRatios,
0491                                                                testbeamEEShape,
0492                                                                EEchi2Parameters_);
0493       double chi2 = chi2expressEE_.chi2();
0494       uncalibRecHit.setChi2(chi2);
0495 
0496       if (kPoorRecoFlagEE_) {
0497         if (chi2 > chi2ThreshEE_) {
0498           // first check if all samples are ok, if not don't use chi2 to flag
0499           bool samplesok = true;
0500           for (int sample = 0; sample < EcalDataFrame::MAXSAMPLES; ++sample) {
0501             if (!sampleMask_->useSampleEE(sample)) {
0502               samplesok = false;
0503               break;
0504             }
0505           }
0506           if (samplesok)
0507             uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kPoorReco);
0508         }
0509       }
0510 
0511     } else {
0512       double amplitude = uncalibRecHit.amplitude();
0513       double amplitudeOutOfTime = 0.;
0514       double jitter = uncalibRecHit.jitter();
0515 
0516       EcalUncalibRecHitRecChi2Algo<EBDataFrame> chi2expressEB_(*itdg,
0517                                                                amplitude,
0518                                                                (itimeconst + offsetTime),
0519                                                                amplitudeOutOfTime,
0520                                                                jitter,
0521                                                                pedVec,
0522                                                                pedRMSVec,
0523                                                                gainRatios,
0524                                                                testbeamEBShape,
0525                                                                EBchi2Parameters_);
0526       double chi2 = chi2expressEB_.chi2();
0527       uncalibRecHit.setChi2(chi2);
0528 
0529       if (kPoorRecoFlagEB_) {
0530         if (chi2 > chi2ThreshEB_) {
0531           // first check if all samples are ok, if not don't use chi2 to flag
0532           bool samplesok = true;
0533           for (int sample = 0; sample < EcalDataFrame::MAXSAMPLES; ++sample) {
0534             if (!sampleMask_->useSampleEB(sample)) {
0535               samplesok = false;
0536               break;
0537             }
0538           }
0539           if (samplesok)
0540             uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kPoorReco);
0541         }
0542       }
0543     }
0544   }
0545 
0546   // set flags if gain switch has occurred
0547   if (((EcalDataFrame)(*itdg)).hasSwitchToGain6())
0548     uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kHasSwitchToGain6);
0549   if (((EcalDataFrame)(*itdg)).hasSwitchToGain1())
0550     uncalibRecHit.setFlagBit(EcalUncalibratedRecHit::kHasSwitchToGain1);
0551 
0552   // put the recHit in the collection
0553   if (detid.subdetId() == EcalEndcap) {
0554     result.push_back(uncalibRecHit);
0555   } else {
0556     result.push_back(uncalibRecHit);
0557   }
0558 
0559   return true;
0560 }
0561 
0562 edm::ParameterSetDescription EcalUncalibRecHitWorkerGlobal::getAlgoDescription() {
0563   edm::ParameterSetDescription psd;
0564   psd.addNode(
0565       edm::ParameterDescription<std::vector<double>>(
0566           "eePulseShape", {5.2e-05, -5.26e-05, 6.66e-05, 0.1168, 0.7575, 1.0, 0.8876, 0.6732, 0.4741, 0.3194}, true) and
0567       edm::ParameterDescription<std::vector<double>>(
0568           "EBtimeFitParameters",
0569           {-2.015452, 3.130702, -12.3473, 41.88921, -82.83944, 91.01147, -50.35761, 11.05621},
0570           true) and
0571       edm::ParameterDescription<double>("outOfTimeThresholdGain61pEB", 5, true) and
0572       edm::ParameterDescription<double>("amplitudeThresholdEE", 10, true) and
0573       edm::ParameterDescription<double>("EBtimeConstantTerm", 0.6, true) and
0574       edm::ParameterDescription<double>("outOfTimeThresholdGain61pEE", 1000, true) and
0575       edm::ParameterDescription<double>("ebSpikeThreshold", 1.042, true) and
0576       edm::ParameterDescription<double>("EBtimeNconst", 28.5, true) and
0577       edm::ParameterDescription<bool>("kPoorRecoFlagEB", true, true) and
0578       edm::ParameterDescription<std::vector<double>>(
0579           "ebPulseShape", {5.2e-05, -5.26e-05, 6.66e-05, 0.1168, 0.7575, 1.0, 0.8876, 0.6732, 0.4741, 0.3194}, true) and
0580       edm::ParameterDescription<double>("EBtimeFitLimits_Lower", 0.2, true) and
0581       edm::ParameterDescription<bool>("kPoorRecoFlagEE", false, true) and
0582       edm::ParameterDescription<double>("chi2ThreshEB_", 36.0, true) and
0583       edm::ParameterDescription<edm::ESInputTag>("timeCalibTag", edm::ESInputTag(), true) and
0584       edm::ParameterDescription<edm::ESInputTag>("timeOffsetTag", edm::ESInputTag(), true) and
0585       edm::ParameterDescription<std::vector<double>>(
0586           "EEtimeFitParameters",
0587           {-2.390548, 3.553628, -17.62341, 67.67538, -133.213, 140.7432, -75.41106, 16.20277},
0588           true) and
0589       edm::ParameterDescription<double>("outOfTimeThresholdGain61mEE", 1000, true) and
0590       edm::ParameterDescription<std::vector<double>>("EEchi2Parameters", {2.122, 0.022, 2.122, 0.022}, true) and
0591       edm::ParameterDescription<double>("outOfTimeThresholdGain12mEE", 1000, true) and
0592       edm::ParameterDescription<double>("outOfTimeThresholdGain12mEB", 5, true) and
0593       edm::ParameterDescription<double>("EEtimeFitLimits_Upper", 1.4, true) and
0594       edm::ParameterDescription<double>("EEtimeFitLimits_Lower", 0.2, true) and
0595       edm::ParameterDescription<std::vector<double>>("EEamplitudeFitParameters", {1.89, 1.4}, true) and
0596       edm::ParameterDescription<std::vector<double>>("EBamplitudeFitParameters", {1.138, 1.652}, true) and
0597       edm::ParameterDescription<double>("amplitudeThresholdEB", 10, true) and
0598       edm::ParameterDescription<double>("outOfTimeThresholdGain12pEE", 1000, true) and
0599       edm::ParameterDescription<double>("outOfTimeThresholdGain12pEB", 5, true) and
0600       edm::ParameterDescription<double>("EEtimeNconst", 31.8, true) and
0601       edm::ParameterDescription<double>("outOfTimeThresholdGain61mEB", 5, true) and
0602       edm::ParameterDescription<std::vector<double>>("EBchi2Parameters", {2.122, 0.022, 2.122, 0.022}, true) and
0603       edm::ParameterDescription<double>("EEtimeConstantTerm", 1.0, true) and
0604       edm::ParameterDescription<double>("chi2ThreshEE_", 95.0, true) and
0605       edm::ParameterDescription<double>("EBtimeFitLimits_Upper", 1.4, true));
0606 
0607   return psd;
0608 }
0609 
0610 #include "FWCore/Framework/interface/MakerMacros.h"
0611 #include "RecoLocalCalo/EcalRecProducers/interface/EcalUncalibRecHitWorkerFactory.h"
0612 DEFINE_EDM_PLUGIN(EcalUncalibRecHitWorkerFactory, EcalUncalibRecHitWorkerGlobal, "EcalUncalibRecHitWorkerGlobal");
0613 #include "RecoLocalCalo/EcalRecProducers/interface/EcalUncalibRecHitFillDescriptionWorkerFactory.h"
0614 DEFINE_EDM_PLUGIN(EcalUncalibRecHitFillDescriptionWorkerFactory,
0615                   EcalUncalibRecHitWorkerGlobal,
0616                   "EcalUncalibRecHitWorkerGlobal");
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!