Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:18:45 

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