Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-02-16 06:16:01 

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