Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-09-22 22:37:05

0001 #ifndef PFAnalyzer_H
0002 #define PFAnalyzer_H
0003 
0004 /** \class JetMETAnalyzer
0005  *
0006  *  DQM PF candidate analysis monitoring
0007  *
0008  *  \author J. Roloff - Brown University
0009  *
0010  */
0011 
0012 #include <memory>
0013 #include <fstream>
0014 #include <utility>
0015 #include <string>
0016 #include <cmath>
0017 #include <map>
0018 
0019 #include "FWCore/Framework/interface/MakerMacros.h"
0020 #include "FWCore/Utilities/interface/EDGetToken.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/EventSetup.h"
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0026 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0027 #include "FWCore/Common/interface/TriggerNames.h"
0028 
0029 #include "DataFormats/ParticleFlowReco/interface/PFBlockFwd.h"
0030 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
0031 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0032 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0033 
0034 #include "DataFormats/JetReco/interface/Jet.h"
0035 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0036 #include "DataFormats/JetReco/interface/PFJet.h"
0037 
0038 #include "DataFormats/VertexReco/interface/Vertex.h"
0039 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0040 #include "DataFormats/TrackReco/interface/Track.h"
0041 
0042 #include "DataFormats/Common/interface/TriggerResults.h"
0043 #include "DataFormats/Math/interface/deltaR.h"
0044 
0045 #include "DQMServices/Core/interface/DQMStore.h"
0046 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0047 
0048 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0049 class PFAnalyzer;
0050 
0051 class PFAnalyzer : public DQMEDAnalyzer {
0052 public:
0053   /// Constructor
0054   PFAnalyzer(const edm::ParameterSet&);
0055 
0056   /// Destructor
0057   ~PFAnalyzer() override;
0058 
0059   /// Inizialize parameters for histo binning
0060   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0061 
0062   /// Get the analysis
0063   void analyze(const edm::Event&, const edm::EventSetup&) override;
0064 
0065   /// Initialize run-based parameters
0066   void dqmBeginRun(const edm::Run&, const edm::EventSetup&) override;
0067 
0068 private:
0069   struct binInfo;
0070   // A map between an observable name and a function that obtains that observable from a  PFCandidate.
0071   // This allows us to construct more complicated observables easily, and have it more configurable
0072   // in the config file.
0073   std::map<std::string, std::function<double(const reco::PFCandidate)>> m_funcMap;
0074   std::map<std::string, std::function<double(const reco::PFCandidateCollection, reco::PFCandidate::ParticleType pfType)>>
0075       m_eventFuncMap;
0076   std::map<std::string,
0077            std::function<double(const std::vector<reco::PFCandidatePtr> pfCands, reco::PFCandidate::ParticleType pfType)>>
0078       m_jetWideFuncMap;
0079   std::map<std::string, std::function<double(const reco::PFCandidate, const reco::PFJet)>> m_pfInJetFuncMap;
0080   std::map<std::string, std::function<double(const reco::PFJet)>> m_jetFuncMap;
0081 
0082   binInfo getBinInfo(std::string);
0083 
0084   // Book MonitorElements
0085   void bookMESetSelection(std::string, DQMStore::IBooker&);
0086 
0087   bool passesEventSelection(const edm::Event& iEvent);
0088 
0089   int getPFBin(const reco::PFCandidate pfCand, int i);
0090   int getJetBin(const reco::PFJet jetCand, int i);
0091 
0092   int getBinNumber(double binVal, std::vector<double> bins);
0093   int getBinNumbers(std::vector<double> binVal, std::vector<std::vector<double>> bins);
0094   std::vector<double> getBinList(std::string binString);
0095 
0096   std::vector<std::string> getAllSuffixes(std::vector<std::string> observables,
0097                                           std::vector<std::vector<double>> binnings);
0098   std::string stringWithDecimals(int bin, std::vector<double> bins);
0099 
0100   std::string getSuffix(std::vector<int> binList,
0101                         std::vector<std::string> observables,
0102                         std::vector<std::vector<double>> binnings);
0103 
0104   static double getEnergySpectrum(const reco::PFCandidate pfCand, const reco::PFJet jet) {
0105     if (!jet.pt())
0106       return -1;
0107     return pfCand.pt() / jet.pt();
0108   }
0109 
0110   static double getNPFC(const reco::PFCandidateCollection pfCands, reco::PFCandidate::ParticleType pfType) {
0111     int nPF = 0;
0112     for (auto pfCand : pfCands) {
0113       // We use X to indicate all
0114       if (pfCand.particleId() == pfType || pfType == reco::PFCandidate::ParticleType::X) {
0115         nPF++;
0116       }
0117     }
0118     return nPF;
0119   }
0120 
0121   static double getNPFCinJet(const std::vector<reco::PFCandidatePtr> pfCands, reco::PFCandidate::ParticleType pfType) {
0122     int nPF = 0;
0123     for (auto pfCand : pfCands) {
0124       if (!pfCand)
0125         continue;
0126       // We use X to indicate all
0127       if (pfCand->particleId() == pfType || pfType == reco::PFCandidate::ParticleType::X)
0128         nPF++;
0129     }
0130     return nPF;
0131   }
0132 
0133   // Various functions designed to get information from a PF canddidate
0134   static double getPt(const reco::PFCandidate pfCand) { return pfCand.pt(); }
0135   static double getEnergy(const reco::PFCandidate pfCand) { return pfCand.energy(); }
0136   static double getEta(const reco::PFCandidate pfCand) { return pfCand.eta(); }
0137   static double getAbsEta(const reco::PFCandidate pfCand) { return std::abs(pfCand.eta()); }
0138   static double getPhi(const reco::PFCandidate pfCand) { return pfCand.phi(); }
0139 
0140   static double getHadCalibration(const reco::PFCandidate pfCand) {
0141     if (pfCand.rawHcalEnergy() == 0)
0142       return -1;
0143     return pfCand.hcalEnergy() / pfCand.rawHcalEnergy();
0144   }
0145 
0146   static double getTime(const reco::PFCandidate pfCand) { return pfCand.time(); }
0147 
0148   static double getHcalEnergy_depth1(const reco::PFCandidate pfCand) { return pfCand.hcalDepthEnergyFraction(1); }
0149   static double getHcalEnergy_depth2(const reco::PFCandidate pfCand) { return pfCand.hcalDepthEnergyFraction(2); }
0150   static double getHcalEnergy_depth3(const reco::PFCandidate pfCand) { return pfCand.hcalDepthEnergyFraction(3); }
0151   static double getHcalEnergy_depth4(const reco::PFCandidate pfCand) { return pfCand.hcalDepthEnergyFraction(4); }
0152   static double getHcalEnergy_depth5(const reco::PFCandidate pfCand) { return pfCand.hcalDepthEnergyFraction(5); }
0153   static double getHcalEnergy_depth6(const reco::PFCandidate pfCand) { return pfCand.hcalDepthEnergyFraction(6); }
0154   static double getHcalEnergy_depth7(const reco::PFCandidate pfCand) { return pfCand.hcalDepthEnergyFraction(7); }
0155 
0156   static double getEcalEnergy(const reco::PFCandidate pfCand) { return pfCand.ecalEnergy(); }
0157   static double getRawEcalEnergy(const reco::PFCandidate pfCand) { return pfCand.rawEcalEnergy(); }
0158   static double getHcalEnergy(const reco::PFCandidate pfCand) { return pfCand.hcalEnergy(); }
0159   static double getRawHcalEnergy(const reco::PFCandidate pfCand) { return pfCand.rawHcalEnergy(); }
0160   static double getHOEnergy(const reco::PFCandidate pfCand) { return pfCand.hoEnergy(); }
0161   static double getRawHOEnergy(const reco::PFCandidate pfCand) { return pfCand.rawHoEnergy(); }
0162 
0163   static double getMVAIsolated(const reco::PFCandidate pfCand) { return pfCand.mva_Isolated(); }
0164   static double getMVAEPi(const reco::PFCandidate pfCand) { return pfCand.mva_e_pi(); }
0165   static double getMVAEMu(const reco::PFCandidate pfCand) { return pfCand.mva_e_mu(); }
0166   static double getMVAPiMu(const reco::PFCandidate pfCand) { return pfCand.mva_pi_mu(); }
0167   static double getMVANothingGamma(const reco::PFCandidate pfCand) { return pfCand.mva_nothing_gamma(); }
0168   static double getMVANothingNH(const reco::PFCandidate pfCand) { return pfCand.mva_nothing_nh(); }
0169   static double getMVAGammaNH(const reco::PFCandidate pfCand) { return pfCand.mva_gamma_nh(); }
0170 
0171   static double getDNNESigIsolated(const reco::PFCandidate pfCand) { return pfCand.dnn_e_sigIsolated(); }
0172   static double getDNNESigNonIsolated(const reco::PFCandidate pfCand) { return pfCand.dnn_e_sigNonIsolated(); }
0173   static double getDNNEBkgNonIsolated(const reco::PFCandidate pfCand) { return pfCand.dnn_e_bkgNonIsolated(); }
0174   static double getDNNEBkgTauIsolated(const reco::PFCandidate pfCand) { return pfCand.dnn_e_bkgTau(); }
0175   static double getDNNEBkgPhotonIsolated(const reco::PFCandidate pfCand) { return pfCand.dnn_e_bkgPhoton(); }
0176 
0177   static double getECalEFrac(const reco::PFCandidate pfCand) { return pfCand.ecalEnergy() / pfCand.energy(); }
0178   static double getHCalEFrac(const reco::PFCandidate pfCand) { return pfCand.hcalEnergy() / pfCand.energy(); }
0179   static double getTrackPt(const reco::PFCandidate pfCand) {
0180     if (pfCand.trackRef().isNonnull())
0181       return (pfCand.trackRef())->pt();
0182     return 0;
0183   }
0184 
0185   static double getEoverP(const reco::PFCandidate pfCand) {
0186     double energy = 0;
0187     int maxElement = pfCand.elementsInBlocks().size();
0188     for (int e = 0; e < maxElement; ++e) {
0189       // Get elements from block
0190       reco::PFBlockRef blockRef = pfCand.elementsInBlocks()[e].first;
0191       const edm::OwnVector<reco::PFBlockElement>& elements = blockRef->elements();
0192       for (unsigned iEle = 0; iEle < elements.size(); iEle++) {
0193         if (elements[iEle].index() == pfCand.elementsInBlocks()[e].second) {
0194           if (elements[iEle].type() == reco::PFBlockElement::HCAL ||
0195               elements[iEle].type() == reco::PFBlockElement::ECAL) {  // Element is HB or HE
0196             reco::PFClusterRef clusterref = elements[iEle].clusterRef();
0197             reco::PFCluster cluster = *clusterref;
0198             energy += cluster.energy();
0199           }
0200         }
0201       }
0202     }
0203     return energy / pfCand.p();
0204   }
0205 
0206   static double getHCalEnergy(const reco::PFCandidate pfCand) {
0207     double energy = 0;
0208     int maxElement = pfCand.elementsInBlocks().size();
0209     for (int e = 0; e < maxElement; ++e) {
0210       // Get elements from block
0211       reco::PFBlockRef blockRef = pfCand.elementsInBlocks()[e].first;
0212       const edm::OwnVector<reco::PFBlockElement>& elements = blockRef->elements();
0213       for (unsigned iEle = 0; iEle < elements.size(); iEle++) {
0214         if (elements[iEle].index() == pfCand.elementsInBlocks()[e].second) {
0215           if (elements[iEle].type() == reco::PFBlockElement::HCAL) {  // Element is HB or HE
0216             // Get cluster and hits
0217             reco::PFClusterRef clusterref = elements[iEle].clusterRef();
0218             reco::PFCluster cluster = *clusterref;
0219             //std::vector<std::pair<DetId, float>> hitsAndFracs = cluster.hitsAndFractions();
0220             energy += cluster.energy();
0221           }
0222         }
0223       }
0224     }
0225     return energy;
0226   }
0227 
0228   static double getECalEnergy(const reco::PFCandidate pfCand) {
0229     double energy = 0;
0230     int maxElement = pfCand.elementsInBlocks().size();
0231     for (int e = 0; e < maxElement; ++e) {
0232       // Get elements from block
0233       reco::PFBlockRef blockRef = pfCand.elementsInBlocks()[e].first;
0234       const edm::OwnVector<reco::PFBlockElement>& elements = blockRef->elements();
0235       for (unsigned iEle = 0; iEle < elements.size(); iEle++) {
0236         if (elements[iEle].index() == pfCand.elementsInBlocks()[e].second) {
0237           if (elements[iEle].type() == reco::PFBlockElement::ECAL) {  // Element is HB or HE
0238             // Get cluster and hits
0239             reco::PFClusterRef clusterref = elements[iEle].clusterRef();
0240             // When we don't have isolated tracks, this will be a bit useless, since the energy is shared across multiple tracks
0241             reco::PFCluster cluster = *clusterref;
0242             energy += cluster.energy();
0243           }
0244         }
0245       }
0246     }
0247     return energy;
0248   }
0249 
0250   static double getNTracksInBlock(const reco::PFCandidate pfCand) {
0251     // We need this function to return a double, even though this is an integer value
0252     double nTrack = 0;
0253     int maxElement = pfCand.elementsInBlocks().size();
0254     for (int e = 0; e < maxElement; ++e) {
0255       // Get elements from block
0256       reco::PFBlockRef blockRef = pfCand.elementsInBlocks()[e].first;
0257       const edm::OwnVector<reco::PFBlockElement>& elements = blockRef->elements();
0258       for (unsigned iEle = 0; iEle < elements.size(); iEle++) {
0259         if (elements[iEle].index() == pfCand.elementsInBlocks()[e].second) {
0260           if (elements[iEle].type() == reco::PFBlockElement::TRACK) {  // Element is HB or HE
0261             nTrack += 1;
0262           }
0263         }
0264       }
0265     }
0266     return nTrack;
0267   }
0268 
0269   static double getJetPt(const reco::PFJet jet) { return jet.pt(); }
0270 
0271   edm::EDGetTokenT<reco::PFCandidateCollection> thePfCandidateCollection_;
0272   edm::EDGetTokenT<std::vector<reco::Vertex>> vertexToken_;
0273   edm::EDGetTokenT<reco::PFJetCollection> pfJetsToken_;
0274   edm::InputTag srcWeights;
0275 
0276   edm::EDGetTokenT<edm::ValueMap<float>> weightsToken_;
0277   edm::ValueMap<float> const* weights_;
0278 
0279   edm::EDGetTokenT<GenEventInfoProduct> tok_ew_;
0280 
0281   edm::InputTag theTriggerResultsLabel_;
0282   edm::InputTag vertexTag_;
0283   edm::InputTag highPtJetExpr_;
0284   edm::EDGetTokenT<edm::TriggerResults> triggerResultsToken_;
0285 
0286   std::vector<std::vector<std::string>> m_allSuffixes;
0287   std::vector<std::vector<std::string>> m_allJetSuffixes;
0288 
0289   // The directory where the output is stored
0290   std::string m_directory;
0291 
0292   // All of the histograms, stored as a map between the histogram name and the histogram
0293   std::map<std::string, MonitorElement*> map_of_MEs;
0294 
0295   //check later if we need only one set of parameters
0296   edm::ParameterSet parameters_;
0297 
0298   typedef std::vector<std::string> vstring;
0299   typedef std::vector<double> vDouble;
0300   // Information on which observables to make histograms for.
0301   // In the config file, this should come as a comma-separated list of
0302   // the observable name, the number of bins for the histogram, and
0303   // the lowest and highest values for the histogram.
0304   // The observable name should have an entry in m_funcMap to define how
0305   // it can be retrieved from a PFCandidate.
0306   vstring m_observables;
0307   vstring m_eventObservables;
0308   vstring m_pfInJetObservables;
0309 
0310   vstring m_observableNames;
0311   vstring m_eventObservableNames;
0312   vstring m_pfInJetObservableNames;
0313 
0314   // Information on what cuts should be applied to PFCandidates that are
0315   // being monitored. In the config file, this should come as a comma-separated list of
0316   // the observable name, and the lowest and highest values for the histogram.
0317   // The observable name should have an entry in m_funcMap to define how
0318   // it can be retrieved from a PFCandidate.
0319   vstring m_cutList;
0320   std::vector<std::vector<std::string>> m_fullCutList;
0321   std::vector<std::vector<std::vector<double>>> m_binList;
0322 
0323   // Information on what cuts should be applied to PFJets, in the case that we
0324   // match PFCs to jets.In the config file, this should come as a comma-separated list of
0325   // the observable name, and the lowest and highest values for the histogram.
0326   // The observable name should have an entry in m_jetFuncMap to define how
0327   // it can be retrieved from a PFJet.
0328   vstring m_jetCutList;
0329   std::vector<std::vector<std::string>> m_fullJetCutList;
0330   std::vector<std::vector<std::vector<double>>> m_jetBinList;
0331 
0332   vDouble m_npvBins;
0333 
0334   // The dR radius used to match PFCs to jets.
0335   // Making this configurable is useful in case you want to look at the core of a jet.
0336   double m_matchingRadius;
0337 
0338   std::vector<std::string> m_pfNames;
0339 };
0340 
0341 struct PFAnalyzer::binInfo {
0342   std::string observable;
0343   std::string axisName;
0344   int nBins;
0345   double binMin;
0346   double binMax;
0347 };
0348 
0349 DEFINE_FWK_MODULE(PFAnalyzer);
0350 #endif