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