Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:21:05 

0001 #include "RecoParticleFlow/PFClusterProducer/interface/PFClusterEMEnergyCorrector.h"
0002 
0003 #include "vdt/vdtMath.h"
0004 #include <array>
0005 
0006 namespace {
0007   typedef reco::PFCluster::EEtoPSAssociation::value_type EEPSPair;
0008   bool sortByKey(const EEPSPair &a, const EEPSPair &b) { return a.first < b.first; }
0009 
0010   double getOffset(const double lo, const double hi) { return lo + 0.5 * (hi - lo); }
0011   double getScale(const double lo, const double hi) { return 0.5 * (hi - lo); }
0012 }  // namespace
0013 
0014 PFClusterEMEnergyCorrector::PFClusterEMEnergyCorrector(const edm::ParameterSet &conf, edm::ConsumesCollector &&cc)
0015     : ecalClusterToolsESGetTokens_{std::move(cc)},
0016       ecalReadoutToolsESGetTokens_{conf, std::move(cc)},
0017       calibrator_(new PFEnergyCalibration) {
0018   applyCrackCorrections_ = conf.getParameter<bool>("applyCrackCorrections");
0019   applyMVACorrections_ = conf.getParameter<bool>("applyMVACorrections");
0020   srfAwareCorrection_ = conf.getParameter<bool>("srfAwareCorrection");
0021   setEnergyUncertainty_ = conf.getParameter<bool>("setEnergyUncertainty");
0022   maxPtForMVAEvaluation_ = conf.getParameter<double>("maxPtForMVAEvaluation");
0023 
0024   if (applyMVACorrections_) {
0025     meanlimlowEB_ = -0.336;
0026     meanlimhighEB_ = 0.916;
0027     meanoffsetEB_ = getOffset(meanlimlowEB_, meanlimhighEB_);
0028     meanscaleEB_ = getScale(meanlimlowEB_, meanlimhighEB_);
0029 
0030     meanlimlowEE_ = -0.336;
0031     meanlimhighEE_ = 0.916;
0032     meanoffsetEE_ = getOffset(meanlimlowEE_, meanlimhighEE_);
0033     meanscaleEE_ = getScale(meanlimlowEE_, meanlimhighEE_);
0034 
0035     sigmalimlowEB_ = 0.001;
0036     sigmalimhighEB_ = 0.4;
0037     sigmaoffsetEB_ = getOffset(sigmalimlowEB_, sigmalimhighEB_);
0038     sigmascaleEB_ = getScale(sigmalimlowEB_, sigmalimhighEB_);
0039 
0040     sigmalimlowEE_ = 0.001;
0041     sigmalimhighEE_ = 0.4;
0042     sigmaoffsetEE_ = getOffset(sigmalimlowEE_, sigmalimhighEE_);
0043     sigmascaleEE_ = getScale(sigmalimlowEE_, sigmalimhighEE_);
0044 
0045     recHitsEB_ = cc.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("recHitsEBLabel"));
0046     recHitsEE_ = cc.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("recHitsEELabel"));
0047     autoDetectBunchSpacing_ = conf.getParameter<bool>("autoDetectBunchSpacing");
0048 
0049     if (autoDetectBunchSpacing_) {
0050       bunchSpacing_ = cc.consumes<unsigned int>(edm::InputTag("bunchSpacingProducer"));
0051       bunchSpacingManual_ = 0;
0052     } else {
0053       bunchSpacingManual_ = conf.getParameter<int>("bunchSpacing");
0054     }
0055 
0056     condnames_mean_25ns_.insert(condnames_mean_25ns_.end(),
0057                                 {"ecalPFClusterCorV2_EB_pfSize1_mean_25ns",
0058                                  "ecalPFClusterCorV2_EB_pfSize2_mean_25ns",
0059                                  "ecalPFClusterCorV2_EB_pfSize3_ptbin1_mean_25ns",
0060                                  "ecalPFClusterCorV2_EB_pfSize3_ptbin2_mean_25ns",
0061                                  "ecalPFClusterCorV2_EB_pfSize3_ptbin3_mean_25ns",
0062                                  "ecalPFClusterCorV2_EE_pfSize1_mean_25ns",
0063                                  "ecalPFClusterCorV2_EE_pfSize2_mean_25ns",
0064                                  "ecalPFClusterCorV2_EE_pfSize3_ptbin1_mean_25ns",
0065                                  "ecalPFClusterCorV2_EE_pfSize3_ptbin2_mean_25ns",
0066                                  "ecalPFClusterCorV2_EE_pfSize3_ptbin3_mean_25ns"});
0067     condnames_sigma_25ns_.insert(condnames_sigma_25ns_.end(),
0068                                  {"ecalPFClusterCorV2_EB_pfSize1_sigma_25ns",
0069                                   "ecalPFClusterCorV2_EB_pfSize2_sigma_25ns",
0070                                   "ecalPFClusterCorV2_EB_pfSize3_ptbin1_sigma_25ns",
0071                                   "ecalPFClusterCorV2_EB_pfSize3_ptbin2_sigma_25ns",
0072                                   "ecalPFClusterCorV2_EB_pfSize3_ptbin3_sigma_25ns",
0073                                   "ecalPFClusterCorV2_EE_pfSize1_sigma_25ns",
0074                                   "ecalPFClusterCorV2_EE_pfSize2_sigma_25ns",
0075                                   "ecalPFClusterCorV2_EE_pfSize3_ptbin1_sigma_25ns",
0076                                   "ecalPFClusterCorV2_EE_pfSize3_ptbin2_sigma_25ns",
0077                                   "ecalPFClusterCorV2_EE_pfSize3_ptbin3_sigma_25ns"});
0078     condnames_mean_50ns_.insert(condnames_mean_50ns_.end(),
0079                                 {"ecalPFClusterCorV2_EB_pfSize1_mean_50ns",
0080                                  "ecalPFClusterCorV2_EB_pfSize2_mean_50ns",
0081                                  "ecalPFClusterCorV2_EB_pfSize3_ptbin1_mean_50ns",
0082                                  "ecalPFClusterCorV2_EB_pfSize3_ptbin2_mean_50ns",
0083                                  "ecalPFClusterCorV2_EB_pfSize3_ptbin3_mean_50ns",
0084                                  "ecalPFClusterCorV2_EE_pfSize1_mean_50ns",
0085                                  "ecalPFClusterCorV2_EE_pfSize2_mean_50ns",
0086                                  "ecalPFClusterCorV2_EE_pfSize3_ptbin1_mean_50ns",
0087                                  "ecalPFClusterCorV2_EE_pfSize3_ptbin2_mean_50ns",
0088                                  "ecalPFClusterCorV2_EE_pfSize3_ptbin3_mean_50ns"});
0089     condnames_sigma_50ns_.insert(condnames_sigma_50ns_.end(),
0090                                  {"ecalPFClusterCorV2_EB_pfSize1_sigma_50ns",
0091                                   "ecalPFClusterCorV2_EB_pfSize2_sigma_50ns",
0092                                   "ecalPFClusterCorV2_EB_pfSize3_ptbin1_sigma_50ns",
0093                                   "ecalPFClusterCorV2_EB_pfSize3_ptbin2_sigma_50ns",
0094                                   "ecalPFClusterCorV2_EB_pfSize3_ptbin3_sigma_50ns",
0095                                   "ecalPFClusterCorV2_EE_pfSize1_sigma_50ns",
0096                                   "ecalPFClusterCorV2_EE_pfSize2_sigma_50ns",
0097                                   "ecalPFClusterCorV2_EE_pfSize3_ptbin1_sigma_50ns",
0098                                   "ecalPFClusterCorV2_EE_pfSize3_ptbin2_sigma_50ns",
0099                                   "ecalPFClusterCorV2_EE_pfSize3_ptbin3_sigma_50ns"});
0100 
0101     if (srfAwareCorrection_) {
0102       sigmalimlowEE_ = 0.001;
0103       sigmalimhighEE_ = 0.1;
0104       sigmaoffsetEE_ = getOffset(sigmalimlowEE_, sigmalimhighEE_);
0105       sigmascaleEE_ = getScale(sigmalimlowEE_, sigmalimhighEE_);
0106 
0107       ebSrFlagToken_ = cc.consumes<EBSrFlagCollection>(conf.getParameter<edm::InputTag>("ebSrFlagLabel"));
0108       eeSrFlagToken_ = cc.consumes<EESrFlagCollection>(conf.getParameter<edm::InputTag>("eeSrFlagLabel"));
0109 
0110       condnames_mean_.insert(condnames_mean_.end(),
0111                              {"ecalPFClusterCor2017V2_EB_ZS_mean_25ns",
0112                               "ecalPFClusterCor2017V2_EB_Full_ptbin1_mean_25ns",
0113                               "ecalPFClusterCor2017V2_EB_Full_ptbin2_mean_25ns",
0114                               "ecalPFClusterCor2017V2_EB_Full_ptbin3_mean_25ns",
0115                               "ecalPFClusterCor2017V2_EE_ZS_mean_25ns",
0116                               "ecalPFClusterCor2017V2_EE_Full_ptbin1_mean_25ns",
0117                               "ecalPFClusterCor2017V2_EE_Full_ptbin2_mean_25ns",
0118                               "ecalPFClusterCor2017V2_EE_Full_ptbin3_mean_25ns"});
0119 
0120       condnames_sigma_.insert(condnames_sigma_.end(),
0121                               {"ecalPFClusterCor2017V2_EB_ZS_sigma_25ns",
0122                                "ecalPFClusterCor2017V2_EB_Full_ptbin1_sigma_25ns",
0123                                "ecalPFClusterCor2017V2_EB_Full_ptbin2_sigma_25ns",
0124                                "ecalPFClusterCor2017V2_EB_Full_ptbin3_sigma_25ns",
0125                                "ecalPFClusterCor2017V2_EE_ZS_sigma_25ns",
0126                                "ecalPFClusterCor2017V2_EE_Full_ptbin1_sigma_25ns",
0127                                "ecalPFClusterCor2017V2_EE_Full_ptbin2_sigma_25ns",
0128                                "ecalPFClusterCor2017V2_EE_Full_ptbin3_sigma_25ns"});
0129 
0130       for (short i = 0; i < (short)condnames_mean_.size(); i++) {
0131         forestMeanTokens_25ns_.emplace_back(cc.esConsumes(edm::ESInputTag("", condnames_mean_[i])));
0132         forestSigmaTokens_25ns_.emplace_back(cc.esConsumes(edm::ESInputTag("", condnames_sigma_[i])));
0133       }
0134     } else {
0135       for (short i = 0; i < (short)condnames_mean_25ns_.size(); i++) {
0136         forestMeanTokens_25ns_.emplace_back(cc.esConsumes(edm::ESInputTag("", condnames_mean_25ns_[i])));
0137         forestSigmaTokens_25ns_.emplace_back(cc.esConsumes(edm::ESInputTag("", condnames_sigma_25ns_[i])));
0138       }
0139       for (short i = 0; i < (short)condnames_mean_50ns_.size(); i++) {
0140         forestMeanTokens_50ns_.emplace_back(cc.esConsumes(edm::ESInputTag("", condnames_mean_50ns_[i])));
0141         forestSigmaTokens_50ns_.emplace_back(cc.esConsumes(edm::ESInputTag("", condnames_sigma_50ns_[i])));
0142       }
0143     }
0144   }
0145 }
0146 
0147 void PFClusterEMEnergyCorrector::correctEnergies(const edm::Event &evt,
0148                                                  const edm::EventSetup &es,
0149                                                  const reco::PFCluster::EEtoPSAssociation &assoc,
0150                                                  reco::PFClusterCollection &cs) {
0151   // First deal with pre-MVA corrections
0152   // Kept for backward compatibility (and for HLT)
0153   if (!applyMVACorrections_) {
0154     for (unsigned int idx = 0; idx < cs.size(); ++idx) {
0155       reco::PFCluster &cluster = cs[idx];
0156       bool iseb = cluster.layer() == PFLayer::ECAL_BARREL;
0157       float ePS1 = 0., ePS2 = 0.;
0158       if (!iseb)
0159         getAssociatedPSEnergy(idx, assoc, ePS1, ePS2);
0160       double correctedEnergy = calibrator_->energyEm(cluster, ePS1, ePS2, applyCrackCorrections_);
0161       cluster.setCorrectedEnergy(correctedEnergy);
0162     }
0163     return;
0164   }
0165 
0166   // Common objects for SRF-aware and old style corrections
0167   EcalClusterLazyTools lazyTool(evt, ecalClusterToolsESGetTokens_.get(es), recHitsEB_, recHitsEE_);
0168   EcalReadoutTools readoutTool(evt, es, ecalReadoutToolsESGetTokens_);
0169 
0170   if (!srfAwareCorrection_) {
0171     int bunchspacing = 450;
0172     if (autoDetectBunchSpacing_) {
0173       edm::Handle<unsigned int> bunchSpacingH;
0174       evt.getByToken(bunchSpacing_, bunchSpacingH);
0175       bunchspacing = *bunchSpacingH;
0176     } else {
0177       bunchspacing = bunchSpacingManual_;
0178     }
0179 
0180     const unsigned int ncor = (bunchspacing == 25) ? condnames_mean_25ns_.size() : condnames_mean_50ns_.size();
0181 
0182     std::vector<edm::ESHandle<GBRForestD> > forestH_mean(ncor);
0183     std::vector<edm::ESHandle<GBRForestD> > forestH_sigma(ncor);
0184 
0185     if (bunchspacing == 25) {
0186       for (unsigned int icor = 0; icor < ncor; ++icor) {
0187         forestH_mean[icor] = es.getHandle(forestMeanTokens_25ns_[icor]);
0188         forestH_sigma[icor] = es.getHandle(forestSigmaTokens_25ns_[icor]);
0189       }
0190     } else {
0191       for (unsigned int icor = 0; icor < ncor; ++icor) {
0192         forestH_mean[icor] = es.getHandle(forestMeanTokens_50ns_[icor]);
0193         forestH_sigma[icor] = es.getHandle(forestSigmaTokens_50ns_[icor]);
0194       }
0195     }
0196 
0197     std::array<float, 5> eval;
0198     for (unsigned int idx = 0; idx < cs.size(); ++idx) {
0199       reco::PFCluster &cluster = cs[idx];
0200       bool iseb = cluster.layer() == PFLayer::ECAL_BARREL;
0201       float ePS1 = 0., ePS2 = 0.;
0202       if (!iseb)
0203         getAssociatedPSEnergy(idx, assoc, ePS1, ePS2);
0204 
0205       double e = cluster.energy();
0206       double pt = cluster.pt();
0207       double evale = e;
0208       if (maxPtForMVAEvaluation_ > 0. && pt > maxPtForMVAEvaluation_) {
0209         evale *= maxPtForMVAEvaluation_ / pt;
0210       }
0211       double invE = (e == 0.) ? 0. : 1. / e;  //guard against dividing by 0.
0212       int size = lazyTool.n5x5(cluster);
0213 
0214       int ietaix = 0;
0215       int iphiiy = 0;
0216       if (iseb) {
0217         EBDetId ebseed(cluster.seed());
0218         ietaix = ebseed.ieta();
0219         iphiiy = ebseed.iphi();
0220       } else {
0221         EEDetId eeseed(cluster.seed());
0222         ietaix = eeseed.ix();
0223         iphiiy = eeseed.iy();
0224       }
0225 
0226       //find index of corrections (0-4 for EB, 5-9 for EE, depending on cluster size and raw pt)
0227       int coridx = std::min(size, 3) - 1;
0228       if (coridx == 2) {
0229         if (pt > 4.5) {
0230           coridx += 1;
0231         }
0232         if (pt > 18.) {
0233           coridx += 1;
0234         }
0235       }
0236       if (!iseb) {
0237         coridx += 5;
0238       }
0239 
0240       const GBRForestD &meanforest = *forestH_mean[coridx].product();
0241       const GBRForestD &sigmaforest = *forestH_sigma[coridx].product();
0242 
0243       //fill array for forest evaluation
0244       eval[0] = evale;
0245       eval[1] = ietaix;
0246       eval[2] = iphiiy;
0247       if (!iseb) {
0248         eval[3] = ePS1 * invE;
0249         eval[4] = ePS2 * invE;
0250       }
0251 
0252       //these are the actual BDT responses
0253       double rawmean = meanforest.GetResponse(eval.data());
0254       double rawsigma = sigmaforest.GetResponse(eval.data());
0255 
0256       //apply transformation to limited output range (matching the training)
0257       double mean = iseb ? meanoffsetEB_ + meanscaleEB_ * vdt::fast_sin(rawmean)
0258                          : meanoffsetEE_ + meanscaleEE_ * vdt::fast_sin(rawmean);
0259       double sigma = iseb ? sigmaoffsetEB_ + sigmascaleEB_ * vdt::fast_sin(rawsigma)
0260                           : sigmaoffsetEE_ + sigmascaleEE_ * vdt::fast_sin(rawsigma);
0261 
0262       //regression target is ln(Etrue/Eraw)
0263       //so corrected energy is ecor=exp(mean)*e, uncertainty is exp(mean)*eraw*sigma=ecor*sigma
0264       double ecor = vdt::fast_exp(mean) * e;
0265       double sigmacor = sigma * ecor;
0266 
0267       cluster.setCorrectedEnergy(ecor);
0268       if (setEnergyUncertainty_)
0269         cluster.setCorrectedEnergyUncertainty(sigmacor);
0270       else
0271         cluster.setCorrectedEnergyUncertainty(0.);
0272     }
0273     return;
0274   }
0275 
0276   // Selective Readout Flags
0277   edm::Handle<EBSrFlagCollection> ebSrFlags;
0278   edm::Handle<EESrFlagCollection> eeSrFlags;
0279   evt.getByToken(ebSrFlagToken_, ebSrFlags);
0280   evt.getByToken(eeSrFlagToken_, eeSrFlags);
0281   if (not ebSrFlags.isValid() or not eeSrFlags.isValid())
0282     edm::LogInfo("PFClusterEMEnergyCorrector") << "SrFlagCollection information is not available. The ECAL PFCluster "
0283                                                   "corrections will assume \"full readout\" for all hits.";
0284 
0285   const unsigned int ncor = forestMeanTokens_25ns_.size();
0286   std::vector<edm::ESHandle<GBRForestD> > forestH_mean(ncor);
0287   std::vector<edm::ESHandle<GBRForestD> > forestH_sigma(ncor);
0288 
0289   for (unsigned int icor = 0; icor < ncor; ++icor) {
0290     forestH_mean[icor] = es.getHandle(forestMeanTokens_25ns_[icor]);
0291     forestH_sigma[icor] = es.getHandle(forestSigmaTokens_25ns_[icor]);
0292   }
0293 
0294   std::array<float, 6> evalEB;
0295   std::array<float, 5> evalEE;
0296 
0297   for (unsigned int idx = 0; idx < cs.size(); ++idx) {
0298     reco::PFCluster &cluster = cs[idx];
0299     bool iseb = cluster.layer() == PFLayer::ECAL_BARREL;
0300     float ePS1 = 0., ePS2 = 0.;
0301     if (!iseb)
0302       getAssociatedPSEnergy(idx, assoc, ePS1, ePS2);
0303 
0304     double e = cluster.energy();
0305     double pt = cluster.pt();
0306     double evale = e;
0307     if (maxPtForMVAEvaluation_ > 0. && pt > maxPtForMVAEvaluation_) {
0308       evale *= maxPtForMVAEvaluation_ / pt;
0309     }
0310     double invE = (e == 0.) ? 0. : 1. / e;  //guard against dividing by 0.
0311     int size = lazyTool.n5x5(cluster);
0312     int reducedHits = size;
0313     if (size >= 3)
0314       reducedHits = 3;
0315 
0316     int ietaix = 0;
0317     int iphiiy = 0;
0318     if (iseb) {
0319       EBDetId ebseed(cluster.seed());
0320       ietaix = ebseed.ieta();
0321       iphiiy = ebseed.iphi();
0322     } else {
0323       EEDetId eeseed(cluster.seed());
0324       ietaix = eeseed.ix();
0325       iphiiy = eeseed.iy();
0326     }
0327 
0328     // Hardcoded number are positions of modules boundaries of ECAL
0329     int signeta = (ietaix > 0) ? 1 : -1;
0330     int ietamod20 = (std::abs(ietaix) < 26) ? ietaix - signeta : (ietaix - 26 * signeta) % 20;
0331     int iphimod20 = (iphiiy - 1) % 20;
0332 
0333     // Assume that hits for which no information is avaiable have a Full Readout (binary 0011)
0334     int clusFlag = 3;
0335     if (iseb) {
0336       if (ebSrFlags.isValid()) {
0337         auto ecalUnit = readoutTool.readOutUnitOf(static_cast<EBDetId>(cluster.seed()));
0338         EBSrFlagCollection::const_iterator srf = ebSrFlags->find(ecalUnit);
0339         if (srf != ebSrFlags->end())
0340           clusFlag = srf->value();
0341       }
0342     } else {
0343       if (eeSrFlags.isValid()) {
0344         auto ecalUnit = readoutTool.readOutUnitOf(static_cast<EEDetId>(cluster.seed()));
0345         EESrFlagCollection::const_iterator srf = eeSrFlags->find(ecalUnit);
0346         if (srf != eeSrFlags->end())
0347           clusFlag = srf->value();
0348       }
0349     }
0350 
0351     // Find index of corrections (0-3 for EB, 4-7 for EE, depending on cluster size and raw pt)
0352     int coridx = 0;
0353     int regind = 0;
0354     if (!iseb)
0355       regind = 4;
0356 
0357     ///////////////////////////////////////////////////////////////////////////////////
0358     ///a hit can be ZS or forced ZS. A hit can be in Full readout or Forced to be FULL readout
0359     ///if it is ZS, then clusFlag (in binary) = 0001
0360     ///if it is forced ZS, then clusFlag (in binary) = 0101
0361     ///if it is FR, then clusFlag (in binary) = 0011
0362     ///if it is forced FR, then clusFlag (in binary) = 0111
0363     ///i.e 3rd bit is set.
0364     ///Even if it is forced, we should mark it is as ZS or FR. To take care of it, just check the LSB and second LSB(SLSB)
0365     ///////////////////////////////////////////////////////////////////////////////////
0366     int ZS_bit = clusFlag >> 0 & 1;
0367     int FR_bit = clusFlag >> 1 & 1;
0368 
0369     if (ZS_bit != 0 && FR_bit == 0)  ///it is clusFlag==1, 5
0370       coridx = 0 + regind;
0371     else {
0372       if (pt < 2.5)
0373         coridx = 1 + regind;
0374       else if (pt >= 2.5 && pt < 6.)
0375         coridx = 2 + regind;
0376       else if (pt >= 6.)
0377         coridx = 3 + regind;
0378     }
0379     if (ZS_bit == 0 || clusFlag > 7) {
0380       edm::LogWarning("PFClusterEMEnergyCorrector")
0381           << "We can only correct regions readout in ZS (flag 1,5) or FULL readout (flag 3,7). Flag " << clusFlag
0382           << " is not recognized."
0383           << "\n"
0384           << "Assuming FULL readout and continuing";
0385     }
0386 
0387     const GBRForestD &meanforest = *forestH_mean[coridx].product();
0388     const GBRForestD &sigmaforest = *forestH_sigma[coridx].product();
0389 
0390     //fill array for forest evaluation
0391     if (iseb) {
0392       evalEB[0] = evale;
0393       evalEB[1] = ietaix;
0394       evalEB[2] = iphiiy;
0395       evalEB[3] = ietamod20;
0396       evalEB[4] = iphimod20;
0397       evalEB[5] = reducedHits;
0398     } else {
0399       evalEE[0] = evale;
0400       evalEE[1] = ietaix;
0401       evalEE[2] = iphiiy;
0402       evalEE[3] = (ePS1 + ePS2) * invE;
0403       evalEE[4] = reducedHits;
0404     }
0405 
0406     //these are the actual BDT responses
0407     double rawmean = 1;
0408     double rawsigma = 0;
0409 
0410     if (iseb) {
0411       rawmean = meanforest.GetResponse(evalEB.data());
0412       rawsigma = sigmaforest.GetResponse(evalEB.data());
0413     } else {
0414       rawmean = meanforest.GetResponse(evalEE.data());
0415       rawsigma = sigmaforest.GetResponse(evalEE.data());
0416     }
0417 
0418     //apply transformation to limited output range (matching the training)
0419     //the training was done with different transformations for EB and EE (width only)
0420     //makes a the code a bit more cumbersome, but it is not a problem per se
0421     double mean = iseb ? meanoffsetEB_ + meanscaleEB_ * vdt::fast_sin(rawmean)
0422                        : meanoffsetEE_ + meanscaleEE_ * vdt::fast_sin(rawmean);
0423     double sigma = iseb ? sigmaoffsetEB_ + sigmascaleEB_ * vdt::fast_sin(rawsigma)
0424                         : sigmaoffsetEE_ + sigmascaleEE_ * vdt::fast_sin(rawsigma);
0425 
0426     //regression target is ln(Etrue/Eraw)
0427     //so corrected energy is ecor=exp(mean)*e, uncertainty is exp(mean)*eraw*sigma=ecor*sigma
0428     double ecor = iseb ? vdt::fast_exp(mean) * e : vdt::fast_exp(mean) * (e + ePS1 + ePS2);
0429     double sigmacor = sigma * ecor;
0430 
0431     LogDebug("PFClusterEMEnergyCorrector")
0432         << "ieta : iphi : ietamod20 : iphimod20 : size : reducedHits = " << ietaix << " " << iphiiy << " " << ietamod20
0433         << " " << iphimod20 << " " << size << " " << reducedHits << "\n"
0434         << "isEB : eraw : ePS1 : ePS2 : (eps1+eps2)/raw : Flag = " << iseb << " " << evale << " " << ePS1 << " " << ePS2
0435         << " " << (ePS1 + ePS2) / evale << " " << clusFlag << "\n"
0436         << "response : correction = " << exp(mean) << " " << ecor;
0437 
0438     cluster.setCorrectedEnergy(ecor);
0439     if (setEnergyUncertainty_)
0440       cluster.setCorrectedEnergyUncertainty(sigmacor);
0441     else
0442       cluster.setCorrectedEnergyUncertainty(0.);
0443   }
0444 }
0445 
0446 void PFClusterEMEnergyCorrector::getAssociatedPSEnergy(const size_t clusIdx,
0447                                                        const reco::PFCluster::EEtoPSAssociation &assoc,
0448                                                        float &e1,
0449                                                        float &e2) {
0450   e1 = 0;
0451   e2 = 0;
0452   auto ee_key_val = std::make_pair(clusIdx, edm::Ptr<reco::PFCluster>());
0453   const auto clustops = std::equal_range(assoc.begin(), assoc.end(), ee_key_val, sortByKey);
0454   for (auto i_ps = clustops.first; i_ps != clustops.second; ++i_ps) {
0455     edm::Ptr<reco::PFCluster> psclus(i_ps->second);
0456     switch (psclus->layer()) {
0457       case PFLayer::PS1:
0458         e1 += psclus->energy();
0459         break;
0460       case PFLayer::PS2:
0461         e2 += psclus->energy();
0462         break;
0463       default:
0464         break;
0465     }
0466   }
0467 }
[ 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!