File indexing completed on 2023-03-17 11:20:58
0001 #ifndef RecoParticleFlow_Benchmark_GenericBenchmark_h
0002 #define RecoParticleFlow_Benchmark_GenericBenchmark_h
0003
0004
0005 #include "RecoParticleFlow/Benchmark/interface/PFBenchmarkAlgo.h"
0006
0007 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0008 #include "DataFormats/Candidate/interface/Candidate.h"
0009 #include "DataFormats/METReco/interface/MET.h"
0010 #include "DQMServices/Core/interface/DQMStore.h"
0011
0012 #include <string>
0013
0014 #include <TH1F.h>
0015 #include <TH2F.h>
0016 #include <TFile.h>
0017
0018 class BenchmarkTree;
0019
0020
0021 class GenericBenchmark {
0022 public:
0023 typedef dqm::legacy::DQMStore DQMStore;
0024 typedef dqm::legacy::MonitorElement MonitorElement;
0025
0026 GenericBenchmark();
0027 virtual ~GenericBenchmark() noexcept(false);
0028
0029 void setup(DQMStore *DQM = nullptr,
0030 bool PlotAgainstReco_ = true,
0031 float minDeltaEt = -100.,
0032 float maxDeltaEt = 50.,
0033 float minDeltaPhi = -0.5,
0034 float maxDeltaPhi = 0.5,
0035 bool doMetPlots = false);
0036
0037 template <typename C>
0038 void fill(const C *RecoCollection,
0039 const C *GenCollection,
0040 bool startFromGen = false,
0041 bool PlotAgainstReco = true,
0042 bool onlyTwoJets = false,
0043 double recPt_cut = -1.,
0044 double minEta_cut = -1.,
0045 double maxEta_cut = -1.,
0046 double deltaR_cut = -1.);
0047
0048 void write(std::string Filename);
0049
0050 void setfile(TFile *file);
0051
0052 private:
0053 bool accepted(const reco::Candidate *particle, double ptCut, double minEtaCut, double maxEtaCut) const;
0054
0055 void fillHistos(const reco::Candidate *genParticle,
0056 const reco::Candidate *recParticle,
0057 double deltaR_cut,
0058 bool plotAgainstReco);
0059
0060 TFile *file_;
0061
0062 TH1F *hDeltaEt;
0063 TH1F *hDeltaEx;
0064 TH1F *hDeltaEy;
0065 TH2F *hDeltaEtvsEt;
0066 TH2F *hDeltaEtOverEtvsEt;
0067 TH2F *hDeltaEtvsEta;
0068 TH2F *hDeltaEtOverEtvsEta;
0069 TH2F *hDeltaEtvsPhi;
0070 TH2F *hDeltaEtOverEtvsPhi;
0071 TH2F *hDeltaEtvsDeltaR;
0072 TH2F *hDeltaEtOverEtvsDeltaR;
0073
0074 TH2F *hEtRecvsEt;
0075 TH2F *hEtRecOverTrueEtvsTrueEt;
0076
0077 TH1F *hDeltaEta;
0078 TH2F *hDeltaEtavsEt;
0079 TH2F *hDeltaEtavsEta;
0080
0081 TH1F *hDeltaPhi;
0082 TH2F *hDeltaPhivsEt;
0083 TH2F *hDeltaPhivsEta;
0084
0085 TH1F *hDeltaR;
0086 TH2F *hDeltaRvsEt;
0087 TH2F *hDeltaRvsEta;
0088
0089 TH1F *hNRec;
0090
0091 TH1F *hEtGen;
0092 TH1F *hEtaGen;
0093 TH2F *hEtvsEtaGen;
0094 TH2F *hEtvsPhiGen;
0095 TH1F *hPhiGen;
0096
0097 TH1F *hNGen;
0098
0099 TH1F *hEtSeen;
0100 TH1F *hEtaSeen;
0101 TH1F *hPhiSeen;
0102 TH2F *hEtvsEtaSeen;
0103 TH2F *hEtvsPhiSeen;
0104
0105 TH1F *hEtRec;
0106 TH1F *hExRec;
0107 TH1F *hEyRec;
0108 TH1F *hPhiRec;
0109
0110 TH1F *hSumEt;
0111 TH1F *hTrueSumEt;
0112 TH2F *hDeltaSetvsSet;
0113 TH2F *hDeltaMexvsSet;
0114 TH2F *hDeltaSetOverSetvsSet;
0115 TH2F *hRecSetvsTrueSet;
0116 TH2F *hRecSetOverTrueSetvsTrueSet;
0117 TH2F *hTrueMexvsTrueSet;
0118
0119 BenchmarkTree *tree_;
0120
0121 bool doMetPlots_;
0122
0123 protected:
0124 DQMStore *dbe_;
0125 PFBenchmarkAlgo *algo_;
0126 };
0127
0128 template <typename C>
0129 void GenericBenchmark::fill(const C *RecoCollection,
0130 const C *GenCollection,
0131 bool startFromGen,
0132 bool PlotAgainstReco,
0133 bool onlyTwoJets,
0134 double recPt_cut,
0135 double minEta_cut,
0136 double maxEta_cut,
0137 double deltaR_cut) {
0138
0139
0140
0141
0142
0143
0144
0145
0146 if (!startFromGen) {
0147 int nRec = 0;
0148 for (unsigned int i = 0; i < RecoCollection->size(); i++) {
0149
0150 const reco::Candidate *particle = &(*RecoCollection)[i];
0151
0152 assert(particle != nullptr);
0153 if (!accepted(particle, recPt_cut, minEta_cut, maxEta_cut))
0154 continue;
0155
0156
0157 if (onlyTwoJets) {
0158 unsigned highJets = 0;
0159 for (unsigned j = 0; j < RecoCollection->size(); j++) {
0160 const reco::Candidate *otherParticle = &(*RecoCollection)[j];
0161 if (j != i && otherParticle->pt() > particle->pt())
0162 highJets++;
0163 }
0164 if (highJets > 1)
0165 continue;
0166 }
0167 nRec++;
0168
0169 const reco::Candidate *gen_particle = algo_->matchByDeltaR(particle, GenCollection);
0170 if (gen_particle == nullptr)
0171 continue;
0172
0173
0174 fillHistos(gen_particle, particle, deltaR_cut, PlotAgainstReco);
0175 }
0176
0177 hNRec->Fill(nRec);
0178 }
0179
0180
0181
0182
0183
0184
0185 int nGen = 0;
0186 for (unsigned int i = 0; i < GenCollection->size(); i++) {
0187 const reco::Candidate *gen_particle = &(*GenCollection)[i];
0188
0189 if (!accepted(gen_particle, recPt_cut, minEta_cut, maxEta_cut)) {
0190 continue;
0191 }
0192
0193 hEtaGen->Fill(gen_particle->eta());
0194 hPhiGen->Fill(gen_particle->phi());
0195 hEtGen->Fill(gen_particle->et());
0196 hEtvsEtaGen->Fill(gen_particle->eta(), gen_particle->et());
0197 hEtvsPhiGen->Fill(gen_particle->phi(), gen_particle->et());
0198
0199 const reco::Candidate *rec_particle = algo_->matchByDeltaR(gen_particle, RecoCollection);
0200 nGen++;
0201 if (!rec_particle)
0202 continue;
0203
0204
0205 double deltaR = algo_->deltaR(rec_particle, gen_particle);
0206 if (deltaR > deltaR_cut && deltaR_cut != -1.)
0207 continue;
0208
0209 hEtaSeen->Fill(gen_particle->eta());
0210 hPhiSeen->Fill(gen_particle->phi());
0211 hEtSeen->Fill(gen_particle->et());
0212 hEtvsEtaSeen->Fill(gen_particle->eta(), gen_particle->et());
0213 hEtvsPhiSeen->Fill(gen_particle->phi(), gen_particle->et());
0214
0215 hPhiRec->Fill(rec_particle->phi());
0216 hEtRec->Fill(rec_particle->et());
0217 hExRec->Fill(rec_particle->px());
0218 hEyRec->Fill(rec_particle->py());
0219
0220 hEtRecvsEt->Fill(gen_particle->et(), rec_particle->et());
0221 if (gen_particle->et() != 0.0)
0222 hEtRecOverTrueEtvsTrueEt->Fill(gen_particle->et(), rec_particle->et() / gen_particle->et());
0223
0224 if (startFromGen)
0225 fillHistos(gen_particle, rec_particle, deltaR_cut, PlotAgainstReco);
0226 }
0227 hNGen->Fill(nGen);
0228 }
0229
0230 #endif