Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:49 

0001 #include <functional>
0002 #include <memory>
0003 
0004 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0005 #include "CommonTools/Utils/interface/StringObjectFunction.h"
0006 #include "DataFormats/Candidate/interface/LeafCandidate.h"
0007 #include "DataFormats/VertexReco/interface/Vertex.h"
0008 #include "RecoTauTag/RecoTau/interface/ConeTools.h"
0009 #include "RecoTauTag/RecoTau/interface/RecoTauQualityCuts.h"
0010 #include "RecoTauTag/RecoTau/interface/RecoTauVertexAssociator.h"
0011 #include "RecoTauTag/RecoTau/interface/TauDiscriminationProducerBase.h"
0012 #include <FWCore/ParameterSet/interface/ConfigurationDescriptions.h>
0013 #include <FWCore/ParameterSet/interface/ParameterSetDescription.h>
0014 
0015 #include "TMath.h"
0016 #include "TFormula.h"
0017 
0018 /* class PFRecoTauDiscriminationByIsolation
0019  * created : Jul 23 2007,
0020  * revised : Thu Aug 13 14:44:40 PDT 2009
0021  * contributors : Ludovic Houchu (IPHC, Strasbourg),
0022  *                Christian Veelken (UC Davis),
0023  *                Evan K. Friis (UC Davis)
0024  *                Michalis Bachtis (UW Madison)
0025  */
0026 
0027 using namespace reco;
0028 using namespace std;
0029 
0030 class PFRecoTauDiscriminationByIsolation : public PFTauDiscriminationProducerBase {
0031 public:
0032   explicit PFRecoTauDiscriminationByIsolation(const edm::ParameterSet& pset)
0033       : PFTauDiscriminationProducerBase(pset),
0034         moduleLabel_(pset.getParameter<std::string>("@module_label")),
0035         qualityCutsPSet_(pset.getParameter<edm::ParameterSet>("qualityCuts")) {
0036     includeTracks_ = pset.getParameter<bool>("ApplyDiscriminationByTrackerIsolation");
0037     includeGammas_ = pset.getParameter<bool>("ApplyDiscriminationByECALIsolation");
0038 
0039     calculateWeights_ = pset.getParameter<bool>("ApplyDiscriminationByWeightedECALIsolation");
0040 
0041     enableHGCalWorkaround_ = pset.getParameter<bool>("enableHGCalWorkaround");
0042 
0043     // RIC: multiply neutral isolation by a flat factor.
0044     //      Useful, for instance, to combine charged and neutral isolations
0045     //      with different relative weights
0046     weightGammas_ = pset.getParameter<double>("WeightECALIsolation");
0047 
0048     // RIC: allow to relax the isolation completely beyond a given tau pt
0049     minPtForNoIso_ = pset.getParameter<double>("minTauPtForNoIso");
0050 
0051     applyOccupancyCut_ = pset.getParameter<bool>("applyOccupancyCut");
0052     maximumOccupancy_ = pset.getParameter<uint32_t>("maximumOccupancy");
0053 
0054     applySumPtCut_ = pset.getParameter<bool>("applySumPtCut");
0055     maximumSumPt_ = pset.getParameter<double>("maximumSumPtCut");
0056 
0057     applyRelativeSumPtCut_ = pset.getParameter<bool>("applyRelativeSumPtCut");
0058     maximumRelativeSumPt_ = pset.getParameter<double>("relativeSumPtCut");
0059     offsetRelativeSumPt_ = pset.getParameter<double>("relativeSumPtOffset");
0060 
0061     storeRawOccupancy_ = pset.getParameter<bool>("storeRawOccupancy");
0062     storeRawSumPt_ = pset.getParameter<bool>("storeRawSumPt");
0063     storeRawPUsumPt_ = pset.getParameter<bool>("storeRawPUsumPt");
0064     storeRawFootprintCorrection_ = pset.getParameter<bool>("storeRawFootprintCorrection");
0065     storeRawPhotonSumPt_outsideSignalCone_ = pset.getParameter<bool>("storeRawPhotonSumPt_outsideSignalCone");
0066 
0067     // Sanity check on requested options.  We can't apply cuts and store the
0068     // raw output at the same time
0069     if (applySumPtCut_ || applyOccupancyCut_ || applyRelativeSumPtCut_) {
0070       if (storeRawSumPt_ || storeRawOccupancy_ || storeRawPUsumPt_) {
0071         throw cms::Exception("BadIsoConfig") << "A 'store raw' and a 'apply cut' option have been set to true "
0072                                              << "simultaneously.  These options are mutually exclusive.";
0073       }
0074     }
0075 
0076     // sanity check2 - can't use weighted and unweighted iso at the same time
0077     if (includeGammas_ && calculateWeights_) {
0078       throw cms::Exception("BasIsoConfig")
0079           << "Both 'ApplyDiscriminationByECALIsolation' and 'ApplyDiscriminationByWeightedECALIsolation' "
0080           << "have been set to true. These options are mutually exclusive.";
0081     }
0082 
0083     // Can only store one type
0084     int numStoreOptions = 0;
0085     if (storeRawSumPt_)
0086       ++numStoreOptions;
0087     if (storeRawOccupancy_)
0088       ++numStoreOptions;
0089     if (storeRawPUsumPt_)
0090       ++numStoreOptions;
0091     if (storeRawFootprintCorrection_)
0092       ++numStoreOptions;
0093     if (storeRawPhotonSumPt_outsideSignalCone_)
0094       ++numStoreOptions;
0095     if (numStoreOptions > 1) {
0096       throw cms::Exception("BadIsoConfig") << "Multiple 'store sum pt' and/or 'store occupancy' options are set."
0097                                            << " These options are mutually exclusive.";
0098     }
0099 
0100     customIsoCone_ = pset.getParameter<double>("customOuterCone");
0101 
0102     applyPhotonPtSumOutsideSignalConeCut_ = pset.getParameter<bool>("applyPhotonPtSumOutsideSignalConeCut");
0103     if (applyPhotonPtSumOutsideSignalConeCut_) {
0104       maxAbsPhotonSumPt_outsideSignalCone_ = pset.getParameter<double>("maxAbsPhotonSumPt_outsideSignalCone");
0105       maxRelPhotonSumPt_outsideSignalCone_ = pset.getParameter<double>("maxRelPhotonSumPt_outsideSignalCone");
0106     }
0107 
0108     applyFootprintCorrection_ = pset.getParameter<bool>("applyFootprintCorrection");
0109     if (applyFootprintCorrection_ || storeRawFootprintCorrection_) {
0110       edm::VParameterSet cfgFootprintCorrections = pset.getParameter<edm::VParameterSet>("footprintCorrections");
0111       for (edm::VParameterSet::const_iterator cfgFootprintCorrection = cfgFootprintCorrections.begin();
0112            cfgFootprintCorrection != cfgFootprintCorrections.end();
0113            ++cfgFootprintCorrection) {
0114         std::string selection = cfgFootprintCorrection->getParameter<std::string>("selection");
0115         std::string offset = cfgFootprintCorrection->getParameter<std::string>("offset");
0116         auto footprintCorrection(std::make_unique<FootprintCorrection>(selection, offset));
0117         footprintCorrections_.push_back(std::move(footprintCorrection));
0118       }
0119     }
0120 
0121     // Get the quality cuts specific to the isolation region
0122     edm::ParameterSet isolationQCuts = qualityCutsPSet_.getParameterSet("isolationQualityCuts");
0123 
0124     qcuts_ = std::make_unique<tau::RecoTauQualityCuts>(isolationQCuts);
0125 
0126     vertexAssociator_ = std::make_unique<tau::RecoTauVertexAssociator>(qualityCutsPSet_, consumesCollector());
0127 
0128     applyDeltaBeta_ = pset.getParameter<bool>("applyDeltaBetaCorrection");
0129 
0130     if (applyDeltaBeta_ || calculateWeights_) {
0131       // Factorize the isolation QCuts into those that are used to
0132       // select PU and those that are not.
0133       std::pair<edm::ParameterSet, edm::ParameterSet> puFactorizedIsoQCuts =
0134           reco::tau::factorizePUQCuts(isolationQCuts);
0135 
0136       // Determine the pt threshold for the PU tracks
0137       // First check if the user specifies explicitly the cut.
0138       // For that the user has to provide a >= 0  value for the PtCutOverride.
0139       bool deltaBetaPUTrackPtCutOverride = pset.getParameter<bool>("deltaBetaPUTrackPtCutOverride");
0140       if (deltaBetaPUTrackPtCutOverride) {
0141         double deltaBetaPUTrackPtCutOverride_val = pset.getParameter<double>("deltaBetaPUTrackPtCutOverride_val");
0142         puFactorizedIsoQCuts.second.addParameter<double>("minTrackPt", deltaBetaPUTrackPtCutOverride_val);
0143       } else {
0144         // Secondly take it from the minGammaEt
0145         puFactorizedIsoQCuts.second.addParameter<double>("minTrackPt",
0146                                                          isolationQCuts.getParameter<double>("minGammaEt"));
0147       }
0148 
0149       pileupQcutsPUTrackSelection_ = std::make_unique<tau::RecoTauQualityCuts>(puFactorizedIsoQCuts.first);
0150 
0151       pileupQcutsGeneralQCuts_ = std::make_unique<tau::RecoTauQualityCuts>(puFactorizedIsoQCuts.second);
0152 
0153       pfCandSrc_ = pset.getParameter<edm::InputTag>("particleFlowSrc");
0154       pfCand_token = consumes<edm::View<reco::Candidate> >(pfCandSrc_);
0155       vertexSrc_ = pset.getParameter<edm::InputTag>("vertexSrc");
0156       vertex_token = consumes<reco::VertexCollection>(vertexSrc_);
0157       deltaBetaCollectionCone_ = pset.getParameter<double>("isoConeSizeForDeltaBeta");
0158       std::string deltaBetaFactorFormula = pset.getParameter<string>("deltaBetaFactor");
0159       deltaBetaFormula_ = std::make_unique<TFormula>("DB_corr", deltaBetaFactorFormula.c_str());
0160     }
0161 
0162     applyRhoCorrection_ = pset.getParameter<bool>("applyRhoCorrection");
0163     if (applyRhoCorrection_) {
0164       rhoProducer_ = pset.getParameter<edm::InputTag>("rhoProducer");
0165       rho_token = consumes<double>(rhoProducer_);
0166       rhoConeSize_ = pset.getParameter<double>("rhoConeSize");
0167       rhoUEOffsetCorrection_ = pset.getParameter<double>("rhoUEOffsetCorrection");
0168     }
0169     useAllPFCands_ = pset.getParameter<bool>("UseAllPFCandsForWeights");
0170 
0171     verbosity_ = pset.getParameter<int>("verbosity");
0172   }
0173 
0174   ~PFRecoTauDiscriminationByIsolation() override {}
0175 
0176   void beginEvent(const edm::Event& evt, const edm::EventSetup& evtSetup) override;
0177   double discriminate(const PFTauRef& pfTau) const override;
0178 
0179   inline double weightedSum(const std::vector<CandidatePtr>& inColl_, double eta, double phi) const {
0180     double out = 1.0;
0181     for (auto const& inObj_ : inColl_) {
0182       double sum = (inObj_->pt() * inObj_->pt()) / (deltaR2(eta, phi, inObj_->eta(), inObj_->phi()));
0183       if (sum > 1.0)
0184         out *= sum;
0185     }
0186     return out;
0187   }
0188 
0189   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0190 
0191 private:
0192   std::string moduleLabel_;
0193 
0194   edm::ParameterSet qualityCutsPSet_;
0195   std::unique_ptr<tau::RecoTauQualityCuts> qcuts_;
0196 
0197   // Inverted QCut which selects tracks with bad DZ/trackWeight
0198   std::unique_ptr<tau::RecoTauQualityCuts> pileupQcutsPUTrackSelection_;
0199   std::unique_ptr<tau::RecoTauQualityCuts> pileupQcutsGeneralQCuts_;
0200 
0201   std::unique_ptr<tau::RecoTauVertexAssociator> vertexAssociator_;
0202 
0203   bool includeTracks_;
0204   bool includeGammas_;
0205   bool calculateWeights_;
0206   bool enableHGCalWorkaround_;
0207   double weightGammas_;
0208   bool applyOccupancyCut_;
0209   uint32_t maximumOccupancy_;
0210   bool applySumPtCut_;
0211   double maximumSumPt_;
0212   bool applyRelativeSumPtCut_;
0213   double maximumRelativeSumPt_;
0214   double offsetRelativeSumPt_;
0215   double customIsoCone_;
0216   // RIC:
0217   double minPtForNoIso_;
0218 
0219   bool applyPhotonPtSumOutsideSignalConeCut_;
0220   double maxAbsPhotonSumPt_outsideSignalCone_;
0221   double maxRelPhotonSumPt_outsideSignalCone_;
0222 
0223   bool applyFootprintCorrection_;
0224   struct FootprintCorrection {
0225     FootprintCorrection(const std::string& selection, const std::string& offset)
0226         : selection_(selection), offset_(offset) {}
0227     ~FootprintCorrection() {}
0228     StringCutObjectSelector<PFTau> selection_;
0229     StringObjectFunction<PFTau> offset_;
0230   };
0231   std::vector<std::unique_ptr<FootprintCorrection> > footprintCorrections_;
0232 
0233   // Options to store the raw value in the discriminator instead of boolean pass/fail flag
0234   bool storeRawOccupancy_;
0235   bool storeRawSumPt_;
0236   bool storeRawPUsumPt_;
0237   bool storeRawFootprintCorrection_;
0238   bool storeRawPhotonSumPt_outsideSignalCone_;
0239 
0240   /* **********************************************************************
0241      **** Pileup Subtraction Parameters ***********************************
0242      **********************************************************************/
0243 
0244   // Delta Beta correction
0245   bool applyDeltaBeta_;
0246   edm::InputTag pfCandSrc_;
0247   edm::EDGetTokenT<edm::View<reco::Candidate> > pfCand_token;
0248   // Keep track of how many vertices are in the event
0249   edm::InputTag vertexSrc_;
0250   edm::EDGetTokenT<reco::VertexCollection> vertex_token;
0251   std::vector<reco::CandidatePtr> chargedPFCandidatesInEvent_;
0252   // Size of cone used to collect PU tracks
0253   double deltaBetaCollectionCone_;
0254   std::unique_ptr<TFormula> deltaBetaFormula_;
0255   double deltaBetaFactorThisEvent_;
0256 
0257   // Rho correction
0258   bool applyRhoCorrection_;
0259   bool useAllPFCands_;
0260   edm::InputTag rhoProducer_;
0261   edm::EDGetTokenT<double> rho_token;
0262   double rhoConeSize_;
0263   double rhoUEOffsetCorrection_;
0264   double rhoCorrectionThisEvent_;
0265   double rhoThisEvent_;
0266 
0267   // Flag to enable/disable debug output
0268   int verbosity_;
0269 };
0270 
0271 void PFRecoTauDiscriminationByIsolation::beginEvent(const edm::Event& event, const edm::EventSetup& eventSetup) {
0272   // NB: The use of the PV in this context is necessitated by its use in
0273   // applying quality cuts to the different objects in the isolation cone
0274   // The vertex associator contains the logic to select the appropriate vertex
0275   // We need to pass it the event so it can load the vertices.
0276   vertexAssociator_->setEvent(event);
0277 
0278   // If we are applying the delta beta correction, we need to get the PF
0279   // candidates from the event so we can find the PU tracks.
0280   if (applyDeltaBeta_ || calculateWeights_) {
0281     // Collect all the PF pile up tracks
0282     edm::Handle<edm::View<reco::Candidate> > pfCandidates;
0283     event.getByToken(pfCand_token, pfCandidates);
0284     chargedPFCandidatesInEvent_.clear();
0285     chargedPFCandidatesInEvent_.reserve(pfCandidates->size());
0286     size_t numPFCandidates = pfCandidates->size();
0287     for (size_t i = 0; i < numPFCandidates; ++i) {
0288       reco::CandidatePtr pfCandidate(pfCandidates, i);
0289       if (pfCandidate->charge() != 0) {
0290         chargedPFCandidatesInEvent_.push_back(pfCandidate);
0291       }
0292     }
0293     // Count all the vertices in the event, to parameterize the DB
0294     // correction factor
0295     edm::Handle<reco::VertexCollection> vertices;
0296     event.getByToken(vertex_token, vertices);
0297     size_t nVtxThisEvent = vertices->size();
0298     deltaBetaFactorThisEvent_ = deltaBetaFormula_->Eval(nVtxThisEvent);
0299   }
0300 
0301   if (applyRhoCorrection_) {
0302     edm::Handle<double> rhoHandle_;
0303     event.getByToken(rho_token, rhoHandle_);
0304     rhoThisEvent_ = (*rhoHandle_ - rhoUEOffsetCorrection_) * (3.14159) * rhoConeSize_ * rhoConeSize_;
0305   }
0306 }
0307 
0308 double PFRecoTauDiscriminationByIsolation::discriminate(const PFTauRef& pfTau) const {
0309   LogDebug("discriminate") << " tau: Pt = " << pfTau->pt() << ", eta = " << pfTau->eta() << ", phi = " << pfTau->phi();
0310   LogDebug("discriminate") << *pfTau;
0311 
0312   // collect the objects we are working with (ie tracks, tracks+gammas, etc)
0313   std::vector<CandidatePtr> isoCharged_;
0314   std::vector<CandidatePtr> isoNeutral_;
0315   std::vector<CandidatePtr> isoPU_;
0316   CandidateCollection isoNeutralWeight_;
0317   std::vector<CandidatePtr> chPV_;
0318   isoCharged_.reserve(pfTau->isolationChargedHadrCands().size());
0319   isoNeutral_.reserve(pfTau->isolationGammaCands().size());
0320   isoPU_.reserve(std::min(100UL, chargedPFCandidatesInEvent_.size()));
0321   isoNeutralWeight_.reserve(pfTau->isolationGammaCands().size());
0322 
0323   chPV_.reserve(std::min(50UL, chargedPFCandidatesInEvent_.size()));
0324 
0325   // Get the primary vertex associated to this tau
0326   reco::VertexRef pv = vertexAssociator_->associatedVertex(*pfTau);
0327   // Let the quality cuts know which the vertex to use when applying selections
0328   // on dz, etc.
0329   if (verbosity_) {
0330     if (pv.isNonnull()) {
0331       LogTrace("discriminate") << "pv: x = " << pv->position().x() << ", y = " << pv->position().y()
0332                                << ", z = " << pv->position().z();
0333     } else {
0334       LogTrace("discriminate") << "pv: N/A";
0335     }
0336     if (pfTau->leadChargedHadrCand().isNonnull()) {
0337       LogTrace("discriminate") << "leadPFChargedHadron:"
0338                                << " Pt = " << pfTau->leadChargedHadrCand()->pt() << ","
0339                                << " eta = " << pfTau->leadChargedHadrCand()->eta() << ","
0340                                << " phi = " << pfTau->leadChargedHadrCand()->phi();
0341     } else {
0342       LogTrace("discriminate") << "leadPFChargedHadron: N/A";
0343     }
0344   }
0345 
0346   // CV: isolation is not well defined in case primary vertex or leading charged hadron do not exist
0347   if (!(pv.isNonnull() && pfTau->leadChargedHadrCand().isNonnull()))
0348     return 0.;
0349 
0350   qcuts_->setPV(pv);
0351   qcuts_->setLeadTrack(*pfTau->leadChargedHadrCand());
0352 
0353   if (applyDeltaBeta_ || calculateWeights_) {
0354     pileupQcutsGeneralQCuts_->setPV(pv);
0355     pileupQcutsGeneralQCuts_->setLeadTrack(*pfTau->leadChargedHadrCand());
0356     pileupQcutsPUTrackSelection_->setPV(pv);
0357     pileupQcutsPUTrackSelection_->setLeadTrack(*pfTau->leadChargedHadrCand());
0358   }
0359 
0360   // Load the tracks if they are being used.
0361   if (includeTracks_) {
0362     for (auto const& cand : pfTau->isolationChargedHadrCands()) {
0363       if (qcuts_->filterCandRef(cand)) {
0364         LogTrace("discriminate") << "adding charged iso cand with pt " << cand->pt();
0365         isoCharged_.push_back(cand);
0366       }
0367     }
0368   }
0369   if (includeGammas_ || calculateWeights_) {
0370     for (auto const& cand : pfTau->isolationGammaCands()) {
0371       if (qcuts_->filterCandRef(cand)) {
0372         LogTrace("discriminate") << "adding neutral iso cand with pt " << cand->pt();
0373         isoNeutral_.push_back(cand);
0374       }
0375     }
0376   }
0377 
0378   typedef reco::tau::cone::DeltaRPtrFilter<CandidatePtr> DRFilter;
0379   typedef reco::tau::cone::DeltaRFilter<Candidate> DRFilter2;
0380 
0381   // If desired, get PU tracks.
0382   if (applyDeltaBeta_ || calculateWeights_) {
0383     // First select by inverted the DZ/track weight cuts. True = invert
0384     if (verbosity_) {
0385       std::cout << "Initial PFCands: " << chargedPFCandidatesInEvent_.size() << std::endl;
0386     }
0387 
0388     std::vector<CandidatePtr> allPU = pileupQcutsPUTrackSelection_->filterCandRefs(chargedPFCandidatesInEvent_, true);
0389 
0390     std::vector<CandidatePtr> allNPU = pileupQcutsPUTrackSelection_->filterCandRefs(chargedPFCandidatesInEvent_);
0391     LogTrace("discriminate") << "After track cuts: " << allPU.size();
0392 
0393     // Now apply the rest of the cuts, like pt, and TIP, tracker hits, etc
0394     if (!useAllPFCands_) {
0395       std::vector<CandidatePtr> cleanPU = pileupQcutsGeneralQCuts_->filterCandRefs(allPU);
0396 
0397       std::vector<CandidatePtr> cleanNPU = pileupQcutsGeneralQCuts_->filterCandRefs(allNPU);
0398 
0399       LogTrace("discriminate") << "After cleaning cuts: " << cleanPU.size();
0400 
0401       // Only select PU tracks inside the isolation cone.
0402       DRFilter deltaBetaFilter(pfTau->p4(), 0, deltaBetaCollectionCone_);
0403       for (auto const& cand : cleanPU) {
0404         if (deltaBetaFilter(cand))
0405           isoPU_.push_back(cand);
0406       }
0407 
0408       for (auto const& cand : cleanNPU) {
0409         if (deltaBetaFilter(cand))
0410           chPV_.push_back(cand);
0411       }
0412       LogTrace("discriminate") << "After cone cuts: " << isoPU_.size() << " " << chPV_.size();
0413     } else {
0414       isoPU_ = std::move(allPU);
0415       chPV_ = std::move(allNPU);
0416     }
0417   }
0418 
0419   if (calculateWeights_) {
0420     for (auto const& isoObject : isoNeutral_) {
0421       if (isoObject->charge() != 0) {
0422         // weight only neutral objects
0423         isoNeutralWeight_.push_back(*isoObject);
0424         continue;
0425       }
0426 
0427       double eta = isoObject->eta();
0428       double phi = isoObject->phi();
0429       double sumNPU = 0.5 * log(weightedSum(chPV_, eta, phi));
0430 
0431       double sumPU = 0.5 * log(weightedSum(isoPU_, eta, phi));
0432       LeafCandidate neutral(*isoObject);
0433       if ((sumNPU + sumPU) > 0)
0434         neutral.setP4(((sumNPU) / (sumNPU + sumPU)) * neutral.p4());
0435 
0436       isoNeutralWeight_.push_back(neutral);
0437     }
0438   }
0439 
0440   // Check if we want a custom iso cone
0441   if (customIsoCone_ >= 0.) {
0442     DRFilter filter(pfTau->p4(), 0, customIsoCone_);
0443     DRFilter2 filter2(pfTau->p4(), 0, customIsoCone_);
0444     std::vector<CandidatePtr> isoCharged_filter;
0445     std::vector<CandidatePtr> isoNeutral_filter;
0446     CandidateCollection isoNeutralWeight_filter;
0447     // Remove all the objects not in our iso cone
0448     for (auto const& isoObject : isoCharged_) {
0449       if (filter(isoObject))
0450         isoCharged_filter.push_back(isoObject);
0451     }
0452     if (!calculateWeights_) {
0453       for (auto const& isoObject : isoNeutral_) {
0454         if (filter(isoObject))
0455           isoNeutral_filter.push_back(isoObject);
0456       }
0457       isoNeutral_ = isoNeutral_filter;
0458     } else {
0459       for (auto const& isoObject : isoNeutralWeight_) {
0460         if (filter2(isoObject))
0461           isoNeutralWeight_filter.push_back(isoObject);
0462       }
0463       isoNeutralWeight_ = isoNeutralWeight_filter;
0464     }
0465     isoCharged_ = isoCharged_filter;
0466   }
0467 
0468   bool failsOccupancyCut = false;
0469   bool failsSumPtCut = false;
0470   bool failsRelativeSumPtCut = false;
0471 
0472   //--- nObjects requirement
0473   int neutrals = isoNeutral_.size();
0474 
0475   if (applyDeltaBeta_) {
0476     neutrals -= TMath::Nint(deltaBetaFactorThisEvent_ * isoPU_.size());
0477   }
0478   if (neutrals < 0) {
0479     neutrals = 0;
0480   }
0481 
0482   size_t nOccupants = isoCharged_.size() + neutrals;
0483 
0484   failsOccupancyCut = (nOccupants > maximumOccupancy_);
0485 
0486   double footprintCorrection_value = 0.;
0487   if (applyFootprintCorrection_ || storeRawFootprintCorrection_) {
0488     for (std::vector<std::unique_ptr<FootprintCorrection> >::const_iterator footprintCorrection =
0489              footprintCorrections_.begin();
0490          footprintCorrection != footprintCorrections_.end();
0491          ++footprintCorrection) {
0492       if ((*footprintCorrection)->selection_(*pfTau)) {
0493         footprintCorrection_value = (*footprintCorrection)->offset_(*pfTau);
0494       }
0495     }
0496   }
0497 
0498   double totalPt = 0.;
0499   double puPt = 0.;
0500   //--- Sum PT requirement
0501   if (applySumPtCut_ || applyRelativeSumPtCut_ || storeRawSumPt_ || storeRawPUsumPt_) {
0502     double chargedPt = 0.;
0503     double neutralPt = 0.;
0504     double weightedNeutralPt = 0.;
0505     for (auto const& isoObject : isoCharged_) {
0506       //-------------------------------------------------------------------------
0507       // CV: fix for Phase-2 HLT tau trigger studies
0508       //    (pT of PFCandidates within HGCal acceptance is significantly higher than track pT !!)
0509       if (enableHGCalWorkaround_) {
0510         double trackPt = (isoObject->bestTrack()) ? isoObject->bestTrack()->pt() : 0.;
0511         double pfCandPt = isoObject->pt();
0512         if (pfCandPt > trackPt) {
0513           chargedPt += trackPt;
0514           neutralPt += pfCandPt - trackPt;
0515         } else {
0516           chargedPt += isoObject->pt();
0517         }
0518       } else {
0519         chargedPt += isoObject->pt();
0520       }
0521       //-------------------------------------------------------------------------
0522     }
0523     if (!calculateWeights_) {
0524       for (auto const& isoObject : isoNeutral_) {
0525         neutralPt += isoObject->pt();
0526       }
0527     } else {
0528       for (auto const& isoObject : isoNeutralWeight_) {
0529         weightedNeutralPt += isoObject.pt();
0530       }
0531     }
0532     for (auto const& isoObject : isoPU_) {
0533       puPt += isoObject->pt();
0534     }
0535     LogTrace("discriminate") << "chargedPt = " << chargedPt;
0536     LogTrace("discriminate") << "neutralPt = " << neutralPt;
0537     LogTrace("discriminate") << "weighted neutral Pt = " << weightedNeutralPt;
0538     LogTrace("discriminate") << "puPt = " << puPt << " (delta-beta corr. = " << (deltaBetaFactorThisEvent_ * puPt)
0539                              << ")";
0540 
0541     if (calculateWeights_) {
0542       neutralPt = weightedNeutralPt;
0543     }
0544 
0545     if (applyDeltaBeta_) {
0546       neutralPt -= (deltaBetaFactorThisEvent_ * puPt);
0547     }
0548 
0549     if (applyFootprintCorrection_) {
0550       neutralPt -= footprintCorrection_value;
0551     }
0552 
0553     if (applyRhoCorrection_) {
0554       neutralPt -= rhoThisEvent_;
0555     }
0556 
0557     if (neutralPt < 0.) {
0558       neutralPt = 0.;
0559     }
0560 
0561     totalPt = chargedPt + weightGammas_ * neutralPt;
0562     LogTrace("discriminate") << "totalPt = " << totalPt << " (cut = " << maximumSumPt_ << ")";
0563 
0564     failsSumPtCut = (totalPt > maximumSumPt_);
0565 
0566     //--- Relative Sum PT requirement
0567     failsRelativeSumPtCut = (totalPt > ((pfTau->pt() - offsetRelativeSumPt_) * maximumRelativeSumPt_));
0568   }
0569 
0570   bool failsPhotonPtSumOutsideSignalConeCut = false;
0571   double photonSumPt_outsideSignalCone = 0.;
0572   if (applyPhotonPtSumOutsideSignalConeCut_ || storeRawPhotonSumPt_outsideSignalCone_) {
0573     const std::vector<reco::CandidatePtr>& signalGammas = pfTau->signalGammaCands();
0574     for (std::vector<reco::CandidatePtr>::const_iterator signalGamma = signalGammas.begin();
0575          signalGamma != signalGammas.end();
0576          ++signalGamma) {
0577       double dR = deltaR(pfTau->eta(), pfTau->phi(), (*signalGamma)->eta(), (*signalGamma)->phi());
0578       if (dR > pfTau->signalConeSize())
0579         photonSumPt_outsideSignalCone += (*signalGamma)->pt();
0580     }
0581     if (photonSumPt_outsideSignalCone > maxAbsPhotonSumPt_outsideSignalCone_ ||
0582         photonSumPt_outsideSignalCone > (maxRelPhotonSumPt_outsideSignalCone_ * pfTau->pt())) {
0583       failsPhotonPtSumOutsideSignalConeCut = true;
0584     }
0585   }
0586 
0587   bool fails = (applyOccupancyCut_ && failsOccupancyCut) || (applySumPtCut_ && failsSumPtCut) ||
0588                (applyRelativeSumPtCut_ && failsRelativeSumPtCut) ||
0589                (applyPhotonPtSumOutsideSignalConeCut_ && failsPhotonPtSumOutsideSignalConeCut);
0590 
0591   if (pfTau->pt() > minPtForNoIso_ && minPtForNoIso_ > 0.) {
0592     return 1.;
0593     LogDebug("discriminate") << "tau pt = " << pfTau->pt() << "\t  min cutoff pt = " << minPtForNoIso_;
0594   }
0595 
0596   // We did error checking in the constructor, so this is safe.
0597   if (storeRawSumPt_) {
0598     return totalPt;
0599   } else if (storeRawPUsumPt_) {
0600     if (applyDeltaBeta_)
0601       return puPt;
0602     else if (applyRhoCorrection_)
0603       return rhoThisEvent_;
0604     else
0605       return 0.;
0606   } else if (storeRawOccupancy_) {
0607     return nOccupants;
0608   } else if (storeRawFootprintCorrection_) {
0609     return footprintCorrection_value;
0610   } else if (storeRawPhotonSumPt_outsideSignalCone_) {
0611     return photonSumPt_outsideSignalCone;
0612   } else {
0613     return (fails ? 0. : 1.);
0614   }
0615 }
0616 
0617 void PFRecoTauDiscriminationByIsolation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0618   // pfRecoTauDiscriminationByIsolation
0619   edm::ParameterSetDescription desc;
0620   desc.add<bool>("storeRawFootprintCorrection", false);
0621   desc.add<edm::InputTag>("PFTauProducer", edm::InputTag("pfRecoTauProducer"));
0622   desc.add<bool>("storeRawOccupancy", false);
0623   desc.add<double>("maximumSumPtCut", 6.0);
0624 
0625   edm::ParameterSetDescription desc_qualityCuts;
0626   reco::tau::RecoTauQualityCuts::fillDescriptions(desc_qualityCuts);
0627   desc.add<edm::ParameterSetDescription>("qualityCuts", desc_qualityCuts);
0628 
0629   desc.add<double>("minTauPtForNoIso", -99.0);
0630   desc.add<double>("maxAbsPhotonSumPt_outsideSignalCone", 1000000000.0);
0631   desc.add<edm::InputTag>("vertexSrc", edm::InputTag("offlinePrimaryVertices"));
0632   desc.add<bool>("applySumPtCut", false);
0633   desc.add<double>("rhoConeSize", 0.5);
0634   desc.add<bool>("ApplyDiscriminationByTrackerIsolation", true);
0635   desc.add<bool>("storeRawPhotonSumPt_outsideSignalCone", false);
0636   desc.add<edm::InputTag>("rhoProducer", edm::InputTag("fixedGridRhoFastjetAll"));
0637   desc.add<bool>("enableHGCalWorkaround", false);
0638 
0639   {
0640     edm::ParameterSetDescription vpsd1;
0641     vpsd1.add<std::string>("selection");
0642     vpsd1.add<std::string>("offset");
0643     desc.addVPSet("footprintCorrections", vpsd1, {});
0644   }
0645 
0646   desc.add<std::string>("deltaBetaFactor", "0.38");
0647   desc.add<bool>("applyFootprintCorrection", false);
0648   desc.add<bool>("UseAllPFCandsForWeights", false);
0649   desc.add<double>("relativeSumPtCut", 0.0);
0650   {
0651     edm::ParameterSetDescription pset_Prediscriminants;
0652     pset_Prediscriminants.add<std::string>("BooleanOperator", "and");
0653     {
0654       edm::ParameterSetDescription psd1;
0655       psd1.add<double>("cut");
0656       psd1.add<edm::InputTag>("Producer");
0657       pset_Prediscriminants.addOptional<edm::ParameterSetDescription>("leadTrack", psd1);
0658     }
0659     {
0660       // encountered this at
0661       // RecoTauTag/Configuration/python/HPSPFTaus_cff.py
0662       // Prediscriminants = requireDecayMode.clone(),
0663       // requireDecayMode = cms.PSet(
0664       //     BooleanOperator = cms.string("and"),
0665       //     decayMode = cms.PSet(
0666       //         Producer = cms.InputTag('hpsPFTauDiscriminationByDecayModeFindingNewDMs'),
0667       //         cut = cms.double(0.5)
0668       //     )
0669       // )
0670       edm::ParameterSetDescription psd1;
0671       psd1.add<double>("cut");
0672       psd1.add<edm::InputTag>("Producer");
0673       pset_Prediscriminants.addOptional<edm::ParameterSetDescription>("decayMode", psd1);
0674     }
0675     {
0676       // encountered this at
0677       // RecoTauTag/Configuration/python/HPSPFTaus_cff.py
0678       // Prediscriminants = requireDecayMode.clone(),
0679       // hpsPFTauDiscriminationByLooseIsolation.Prediscriminants.preIso = cms.PSet(
0680       //     Producer = cms.InputTag("hpsPFTauDiscriminationByLooseChargedIsolation"),
0681       //     cut = cms.double(0.5)
0682       // )
0683       edm::ParameterSetDescription psd1;
0684       psd1.add<double>("cut");
0685       psd1.add<edm::InputTag>("Producer");
0686       pset_Prediscriminants.addOptional<edm::ParameterSetDescription>("preIso", psd1);
0687     }
0688     desc.add<edm::ParameterSetDescription>("Prediscriminants", pset_Prediscriminants);
0689   }
0690 
0691   desc.add<unsigned int>("maximumOccupancy", 0);
0692   desc.add<int>("verbosity", 0);
0693 
0694   desc.add<bool>("applyOccupancyCut", true);
0695   desc.add<bool>("applyDeltaBetaCorrection", false);
0696   desc.add<bool>("applyRelativeSumPtCut", false);
0697   desc.add<bool>("storeRawPUsumPt", false);
0698   desc.add<bool>("applyPhotonPtSumOutsideSignalConeCut", false);
0699   desc.add<bool>("deltaBetaPUTrackPtCutOverride", false);
0700   desc.add<bool>("ApplyDiscriminationByWeightedECALIsolation", false);
0701   desc.add<bool>("storeRawSumPt", false);
0702   desc.add<bool>("ApplyDiscriminationByECALIsolation", true);
0703   desc.add<bool>("applyRhoCorrection", false);
0704 
0705   desc.add<double>("WeightECALIsolation", 1.0);
0706   desc.add<double>("rhoUEOffsetCorrection", 1.0);
0707   desc.add<double>("maxRelPhotonSumPt_outsideSignalCone", 0.1);
0708   desc.add<double>("deltaBetaPUTrackPtCutOverride_val", -1.5);
0709   desc.add<double>("isoConeSizeForDeltaBeta", 0.5);
0710   desc.add<double>("relativeSumPtOffset", 0.0);
0711   desc.add<double>("customOuterCone", -1.0);
0712   desc.add<edm::InputTag>("particleFlowSrc", edm::InputTag("particleFlow"));
0713 
0714   descriptions.add("pfRecoTauDiscriminationByIsolation", desc);
0715 }
0716 
0717 DEFINE_FWK_MODULE(PFRecoTauDiscriminationByIsolation);
[ 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!