Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-03-29 02:43:21 

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