MuonIsolation/plugins/CutsIsolatorWithCorrection.cc

0001 #include "CutsIsolatorWithCorrection.h"
0002
0003 #include "DataFormats/Common/interface/Handle.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 using namespace muonisolation;
0007
0008 CutsIsolatorWithCorrection::CutsIsolatorWithCorrection(const edm::ParameterSet& par, edm::ConsumesCollector&& iC)
0009     : theCuts(par.getParameter<std::vector<double> >("EtaBounds"),
0010               par.getParameter<std::vector<double> >("ConeSizes"),
0011               par.getParameter<std::vector<double> >("Thresholds")),
0012       theCutsRel(par.getParameter<std::vector<double> >("EtaBoundsRel"),
0013                  par.getParameter<std::vector<double> >("ConeSizesRel"),
0014                  par.getParameter<std::vector<double> >("ThresholdsRel")),
0015       theCutAbsIso(par.getParameter<bool>("CutAbsoluteIso")),
0016       theCutRelativeIso(par.getParameter<bool>("CutRelativeIso")),
0017       theUseRhoCorrection(par.getParameter<bool>("UseRhoCorrection")),
0018       theRhoToken(iC.consumes<double>(par.getParameter<edm::InputTag>("RhoSrc"))),
0019       theRhoMax(par.getParameter<double>("RhoMax")),
0020       theRhoScaleBarrel(par.getParameter<double>("RhoScaleBarrel")),
0021       theRhoScaleEndcap(par.getParameter<double>("RhoScaleEndcap")),
0022       theEffAreaSFBarrel(par.getParameter<double>("EffAreaSFBarrel")),
0023       theEffAreaSFEndcap(par.getParameter<double>("EffAreaSFEndcap")),
0024       theReturnAbsoluteSum(par.getParameter<bool>("ReturnAbsoluteSum")),
0025       theReturnRelativeSum(par.getParameter<bool>("ReturnRelativeSum")),
0026       theAndOrCuts(par.getParameter<bool>("AndOrCuts")) {
0027   if (!(theCutAbsIso || theCutRelativeIso))
0028     throw cms::Exception("BadConfiguration")
0029         << "Something has to be cut: set either CutAbsoluteIso or CutRelativeIso to true";
0030 }
0031
0032 double CutsIsolatorWithCorrection::depSum(const DepositContainer& deposits, double dr, double corr) const {
0033   double dephlt = -corr;
0034   unsigned int nDeps = deposits.size();
0035   //  edm::LogWarning("CutsIsolatorWithCorrection::depSumIn")
0036   //    << "add nDeposit "<< nDeps<< " \t dr "<<dr<<" \t corr "<<corr;
0037   for (unsigned int iDep = 0; iDep < nDeps; ++iDep) {
0038     double lDep = deposits[iDep].dep->depositWithin(dr);
0039     dephlt += lDep;
0040     //    edm::LogWarning("CutsIsolatorWithCorrection::depSumIDep")
0041     //      <<"dep "<<iDep<<" \t added "<<lDep<<" \t sumnow "<<dephlt;
0042   }
0043
0044   return dephlt;
0045 }
0046
0047 MuIsoBaseIsolator::Result CutsIsolatorWithCorrection::result(const DepositContainer& deposits,
0048                                                              const reco::Track& tk,
0049                                                              const edm::Event* ev) const {
0050   Result answer(ISOL_BOOL_TYPE);
0051
0052   bool absDecision = false;
0053   bool relDecision = false;
0054
0055   double rho = 0.0;
0056   double effAreaSF = 1.0;
0057
0058   static const double pi = 3.14159265358979323846;
0059
0060   //  edm::LogWarning("CutsIsolatorWithCorrection::resultIn")
0061   //    <<"Start tk.pt "<<tk.pt()<<" \t tk.eta "<<tk.eta()<<" \t tk.phi "<<tk.phi();
0062
0063   if (theUseRhoCorrection) {
0064     edm::Handle<double> rhoHandle;
0065     ev->getByToken(theRhoToken, rhoHandle);
0066     rho = *(rhoHandle.product());
0067     if (rho < 0.0)
0068       rho = 0.0;
0069     double rhoScale = fabs(tk.eta()) > 1.442 ? theRhoScaleEndcap : theRhoScaleBarrel;
0070     effAreaSF = fabs(tk.eta()) > 1.442 ? theEffAreaSFEndcap : theEffAreaSFBarrel;
0071     //    edm::LogWarning("CutsIsolatorWithCorrection::resultInRho")
0072     //      << "got rho "<<rho<<" vs max "<<theRhoMax<<" will scale by "<<rhoScale;
0073     if (rho > theRhoMax) {
0074       rho = theRhoMax;
0075     }
0076     rho = rho * rhoScale;
0077     //    edm::LogWarning("CutsIsolatorWithCorrection::resultOutRho")<<" final rho "<<rho;
0078   }
0079
0080   if (theCutAbsIso) {
0081     muonisolation::Cuts::CutSpec cuts_here = theCuts(tk.eta());
0082     double conesize = cuts_here.conesize;
0083     double dephlt = depSum(deposits, conesize, rho * conesize * conesize * pi * effAreaSF);
0084     if (theReturnAbsoluteSum)
0085       answer.valFloat = (float)dephlt;
0086     if (dephlt < cuts_here.threshold) {
0087       absDecision = true;
0088     } else {
0089       absDecision = false;
0090     }
0091     //    edm::LogWarning("CutsIsolatorWithCorrection::resultOutAbsIso")
0092     //      <<"compared dephlt "<<dephlt<<" \t with "<<cuts_here.threshold;
0093   } else
0094     absDecision = true;
0095
0096   if (theCutRelativeIso) {
0097     muonisolation::Cuts::CutSpec cuts_here = theCutsRel(tk.eta());
0098     double conesize = cuts_here.conesize;
0099     double dephlt = depSum(deposits, conesize, rho * conesize * conesize * pi * effAreaSF) / tk.pt();
0100     if (theReturnRelativeSum)
0101       answer.valFloat = (float)dephlt;
0102     if (dephlt < cuts_here.threshold) {
0103       relDecision = true;
0104     } else {
0105       relDecision = false;
0106     }
0107     //    edm::LogWarning("CutsIsolatorWithCorrection::resultOutRelIso")
0108     //      <<"compared dephlt "<<dephlt<<" \t with "<<cuts_here.threshold;
0109   } else
0110     relDecision = true;
0111
0112   if (theAndOrCuts) {
0113     answer.valBool = absDecision && relDecision;
0114   } else {
0115     answer.valBool = absDecision || relDecision;
0116   }
0117
0118   //  edm::LogWarning("CutsIsolatorWithCorrection::result")
0119   //    <<"isAbsIsolated "<<absDecision<<" \t isRelIsolated  "<<relDecision
0120   //    <<" \t combined with AND "<<theAndOrCuts
0121   //    <<" = "  << answer.valBool;
0122
0123   return answer;
0124 }