Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:37:49 

0001 #include "RecoParticleFlow/Benchmark/interface/PFJetBenchmark.h"
0002 #include "RecoParticleFlow/PFTracking/interface/PFTrackAlgoTools.h"
0003 
0004 #include "DataFormats/TrackReco/interface/Track.h"
0005 
0006 // preprocessor macro for booking 1d histos with DQMStore -or- bare Root
0007 #define BOOK1D(name, title, nbinsx, lowx, highx) \
0008   h##name =                                      \
0009       dbe_ ? dbe_->book1D(#name, title, nbinsx, lowx, highx)->getTH1F() : new TH1F(#name, title, nbinsx, lowx, highx)
0010 
0011 // preprocessor macro for booking 2d histos with DQMStore -or- bare Root
0012 #define BOOK2D(name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)                              \
0013   h##name = dbe_ ? dbe_->book2D(#name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)->getTH2F() \
0014                  : new TH2F(#name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)
0015 
0016 //macros for building barrel and endcap histos with one call
0017 #define DBOOK1D(name, title, nbinsx, lowx, highx)         \
0018   BOOK1D(B##name, "Barrel " #title, nbinsx, lowx, highx); \
0019   BOOK1D(E##name, "Endcap " #title, nbinsx, lowx, highx); \
0020   BOOK1D(F##name, "Forward " #title, nbinsx, lowx, highx);
0021 #define DBOOK2D(name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)         \
0022   BOOK2D(B##name, "Barrel " #title, nbinsx, lowx, highx, nbinsy, lowy, highy); \
0023   BOOK2D(E##name, "Endcap " #title, nbinsx, lowx, highx, nbinsy, lowy, highy); \
0024   BOOK2D(F##name, "Forward " #title, nbinsx, lowx, highx, nbinsy, lowy, highy);
0025 
0026 // all versions OK
0027 // preprocesor macro for setting axis titles
0028 #define SETAXES(name, xtitle, ytitle)    \
0029   h##name->GetXaxis()->SetTitle(xtitle); \
0030   h##name->GetYaxis()->SetTitle(ytitle)
0031 
0032 //macro for setting the titles for barrel and endcap together
0033 #define DSETAXES(name, xtitle, ytitle) \
0034   SETAXES(B##name, xtitle, ytitle);    \
0035   SETAXES(E##name, xtitle, ytitle);    \
0036   SETAXES(F##name, xtitle, ytitle);
0037 /*#define SET2AXES(name,xtitle,ytitle) \
0038   hE##name->GetXaxis()->SetTitle(xtitle); hE##name->GetYaxis()->SetTitle(ytitle);  hB##name->GetXaxis()->SetTitle(xtitle); hB##name->GetYaxis()->SetTitle(ytitle)
0039 */
0040 
0041 #define PT (plotAgainstReco_) ? "reconstructed P_{T}" : "generated P_{T}"
0042 #define P (plotAgainstReco_) ? "generated P" : "generated P"
0043 
0044 using namespace reco;
0045 using namespace std;
0046 
0047 PFJetBenchmark::PFJetBenchmark() : file_(nullptr), entry_(0) {}
0048 
0049 PFJetBenchmark::~PFJetBenchmark() {
0050   if (file_)
0051     file_->Close();
0052 }
0053 
0054 void PFJetBenchmark::write() {
0055   // Store the DAQ Histograms
0056   if (!outputFile_.empty()) {
0057     if (dbe_)
0058       dbe_->save(outputFile_);
0059     // use bare Root if no DQM (FWLite applications)
0060     else if (file_) {
0061       file_->Write(outputFile_.c_str());
0062       cout << "Benchmark output written to file " << outputFile_.c_str() << endl;
0063       file_->Close();
0064     }
0065   } else
0066     cout << "No output file specified (" << outputFile_ << "). Results will not be saved!" << endl;
0067 }
0068 
0069 void PFJetBenchmark::setup(string Filename,
0070                            bool debug,
0071                            bool plotAgainstReco,
0072                            bool onlyTwoJets,
0073                            double deltaRMax,
0074                            string benchmarkLabel_,
0075                            double recPt,
0076                            double maxEta,
0077                            DQMStore* dbe_store) {
0078   debug_ = debug;
0079   plotAgainstReco_ = plotAgainstReco;
0080   onlyTwoJets_ = onlyTwoJets;
0081   deltaRMax_ = deltaRMax;
0082   outputFile_ = Filename;
0083   recPt_cut = recPt;
0084   maxEta_cut = maxEta;
0085   file_ = nullptr;
0086   dbe_ = dbe_store;
0087   // print parameters
0088   cout << "PFJetBenchmark Setup parameters ==============================================" << endl;
0089   cout << "Filename to write histograms " << Filename << endl;
0090   cout << "PFJetBenchmark debug " << debug_ << endl;
0091   cout << "plotAgainstReco " << plotAgainstReco_ << endl;
0092   cout << "onlyTwoJets " << onlyTwoJets_ << endl;
0093   cout << "deltaRMax " << deltaRMax << endl;
0094   cout << "benchmarkLabel " << benchmarkLabel_ << endl;
0095   cout << "recPt_cut " << recPt_cut << endl;
0096   cout << "maxEta_cut " << maxEta_cut << endl;
0097 
0098   // Book histogram
0099 
0100   // Establish DQM Store
0101   string path = "PFTask/Benchmarks/" + benchmarkLabel_ + "/";
0102   if (plotAgainstReco)
0103     path += "Reco";
0104   else
0105     path += "Gen";
0106   if (dbe_) {
0107     dbe_->setCurrentFolder(path);
0108   } else {
0109     file_ = new TFile(outputFile_.c_str(), "recreate");
0110     //    TTree * tr = new TTree("PFTast");
0111     //    tr->Branch("Benchmarks/ParticleFlow")
0112     cout << "Info: DQM is not available to provide data storage service. Using TFile to save histograms. " << endl;
0113   }
0114   // Jets inclusive  distributions  (Pt > 20 or specified recPt GeV)
0115   char cutString[35];
0116   sprintf(cutString, "Jet multiplicity P_{T}>%4.1f GeV", recPt_cut);
0117   BOOK1D(Njets, cutString, 50, 0, 50);
0118 
0119   BOOK1D(jetsPt, "Jets P_{T} Distribution", 100, 0, 500);
0120 
0121   sprintf(cutString, "Jets #eta Distribution |#eta|<%4.1f", maxEta_cut);
0122   BOOK1D(jetsEta, cutString, 100, -5, 5);
0123 
0124   BOOK2D(RPtvsEta, "#DeltaP_{T}/P_{T} vs #eta", 200, -5., 5., 200, -1, 1);
0125   BOOK2D(RNeutvsEta, "R_{Neutral} vs #eta", 200, -5., 5., 200, -1, 1);
0126   BOOK2D(RNEUTvsEta, "R_{HCAL+ECAL} vs #eta", 200, -5., 5., 200, -1, 1);
0127   BOOK2D(RNONLvsEta, "R_{HCAL+ECAL - Hcal Only} vs #eta", 200, -5., 5., 200, -1, 1);
0128   BOOK2D(RHCALvsEta, "R_{HCAL} vs #eta", 200, -5., 5., 200, -1, 1);
0129   BOOK2D(RHONLvsEta, "R_{HCAL only} vs #eta", 200, -5., 5., 200, -1, 1);
0130   BOOK2D(RCHEvsEta, "R_{Charged} vs #eta", 200, -5., 5., 200, -1, 1);
0131   BOOK2D(NCHvsEta, "N_{Charged} vs #eta", 200, -5., 5., 200, 0., 200.);
0132   BOOK2D(NCH0vsEta, "N_{Charged} vs #eta, iter 0", 200, -5., 5., 200, 0., 200.);
0133   BOOK2D(NCH1vsEta, "N_{Charged} vs #eta, iter 1", 200, -5., 5., 200, 0., 200.);
0134   BOOK2D(NCH2vsEta, "N_{Charged} vs #eta, iter 2", 200, -5., 5., 200, 0., 200.);
0135   BOOK2D(NCH3vsEta, "N_{Charged} vs #eta, iter 3", 200, -5., 5., 200, 0., 200.);
0136   BOOK2D(NCH4vsEta, "N_{Charged} vs #eta, iter 4", 200, -5., 5., 200, 0., 200.);
0137   BOOK2D(NCH5vsEta, "N_{Charged} vs #eta, iter 5", 200, -5., 5., 200, 0., 200.);
0138   BOOK2D(NCH6vsEta, "N_{Charged} vs #eta, iter 6", 200, -5., 5., 200, 0., 200.);
0139   // delta Pt or E quantities for Barrel
0140   DBOOK1D(RPt, #DeltaP_{T} / P_{T}, 80, -1, 1);
0141   DBOOK1D(RCHE, #DeltaE / E(charged had), 80, -2, 2);
0142   DBOOK1D(RNHE, #DeltaE / E(neutral had), 80, -2, 2);
0143   DBOOK1D(RNEE, #DeltaE / E(neutral em), 80, -2, 2);
0144   DBOOK1D(Rneut, #DeltaE / E(neutral), 80, -2, 2);
0145   DBOOK1D(NCH, #N_{charged}, 200, 0., 200.);
0146   DBOOK2D(RPtvsPt, #DeltaP_{T} / P_{T} vs P_{T}, 250, 0, 500, 200, -1, 1);  //used to be 50 bin for each in x-direction
0147   DBOOK2D(RCHEvsPt, #DeltaE / E(charged had) vs P_{T}, 250, 0, 500, 120, -1, 2);
0148   DBOOK2D(RNHEvsPt, #DeltaE / E(neutral had) vs P_{T}, 250, 0, 500, 120, -1, 2);
0149   DBOOK2D(RNEEvsPt, #DeltaE / E(neutral em) vs P_{T}, 250, 0, 500, 120, -1, 2);
0150   DBOOK2D(RneutvsPt, #DeltaE / E(neutral) vs P_{T}, 250, 0, 500, 120, -1, 2);
0151   DBOOK2D(NCHvsPt, N_{charged} vs P_{T}, 250, 0, 500, 200, 0., 200.);
0152   DBOOK2D(NCH0vsPt, N_{charged} vs P_{T} iter 0, 250, 0, 500, 200, 0., 200.);
0153   DBOOK2D(NCH1vsPt, N_{charged} vs P_{T} iter 1, 250, 0, 500, 200, 0., 200.);
0154   DBOOK2D(NCH2vsPt, N_{charged} vs P_{T} iter 2, 250, 0, 500, 200, 0., 200.);
0155   DBOOK2D(NCH3vsPt, N_{charged} vs P_{T} iter 3, 250, 0, 500, 200, 0., 200.);
0156   DBOOK2D(NCH4vsPt, N_{charged} vs P_{T} iter 4, 250, 0, 500, 200, 0., 200.);
0157   DBOOK2D(NCH5vsPt, N_{charged} vs P_{T} iter 5, 250, 0, 500, 200, 0., 200.);
0158   DBOOK2D(NCH6vsPt, N_{charged} vs P_{T} iter 6, 250, 0, 500, 200, 0., 200.);
0159 
0160   DBOOK2D(RNEUTvsP, #DeltaE / E(ECAL + HCAL) vs P, 250, 0, 1000, 150, -1.5, 1.5);
0161   DBOOK2D(RNONLvsP, #DeltaE / E(ECAL + HCAL - only) vs P, 250, 0, 1000, 150, -1.5, 1.5);
0162   DBOOK2D(RHCALvsP, #DeltaE / E(HCAL) vs P, 250, 0, 1000, 150, -1.5, 1.5);
0163   DBOOK2D(RHONLvsP, #DeltaE / E(HCAL only) vs P, 250, 0, 1000, 150, -1.5, 1.5);
0164   DBOOK1D(RPt20_40, #DeltaP_{T} / P_{T}, 80, -1, 1);
0165   DBOOK1D(RPt40_60, #DeltaP_{T} / P_{T}, 80, -1, 1);
0166   DBOOK1D(RPt60_80, #DeltaP_{T} / P_{T}, 80, -1, 1);
0167   DBOOK1D(RPt80_100, #DeltaP_{T} / P_{T}, 80, -1, 1);
0168   DBOOK1D(RPt100_150, #DeltaP_{T} / P_{T}, 80, -1, 1);
0169   DBOOK1D(RPt150_200, #DeltaP_{T} / P_{T}, 80, -1, 1);
0170   DBOOK1D(RPt200_250, #DeltaP_{T} / P_{T}, 80, -1, 1);
0171   DBOOK1D(RPt250_300, #DeltaP_{T} / P_{T}, 80, -1, 1);
0172   DBOOK1D(RPt300_400, #DeltaP_{T} / P_{T}, 160, -1, 1);
0173   DBOOK1D(RPt400_500, #DeltaP_{T} / P_{T}, 160, -1, 1);
0174   DBOOK1D(RPt500_750, #DeltaP_{T} / P_{T}, 160, -1, 1);
0175   DBOOK1D(RPt750_1250, #DeltaP_{T} / P_{T}, 160, -1, 1);
0176   DBOOK1D(RPt1250_2000, #DeltaP_{T} / P_{T}, 160, -1, 1);
0177   DBOOK1D(RPt2000_5000, #DeltaP_{T} / P_{T}, 160, -1, 1);
0178 
0179   DBOOK2D(DEtavsPt, #Delta #eta vs P_{T}, 1000, 0, 2000, 500, -0.5, 0.5);
0180   DBOOK2D(DPhivsPt, #Delta #phi vs P_{T}, 1000, 0, 2000, 500, -0.5, 0.5);
0181   BOOK2D(DEtavsEta, "#Delta#eta vs P_{T}", 1000, -5., +5., 500, -0.5, 0.5);
0182   BOOK2D(DPhivsEta, "#Delta#phi vs P_{T}", 1000, -5., +5., 500, -0.5, 0.5);
0183 
0184   // Set Axis Titles
0185 
0186   // Jets inclusive  distributions  (Pt > 20 GeV)
0187   SETAXES(Njets, "", "Multiplicity");
0188   SETAXES(jetsPt, PT, "Number of Events");
0189   SETAXES(jetsEta, "#eta", "Number of Events");
0190   SETAXES(RNeutvsEta, "#eta", "#DeltaE/E (Neutral)");
0191   SETAXES(RNEUTvsEta, "#eta", "#DeltaE/E (ECAL+HCAL)");
0192   SETAXES(RNONLvsEta, "#eta", "#DeltaE/E (ECAL+HCAL-only)");
0193   SETAXES(RHCALvsEta, "#eta", "#DeltaE/E (HCAL)");
0194   SETAXES(RHONLvsEta, "#eta", "#DeltaE/E (HCAL Only)");
0195   SETAXES(RCHEvsEta, "#eta", "#DeltaE/E (Charged)");
0196   SETAXES(RPtvsEta, "#eta", "#DeltaP_{T}/P_{T}");
0197   SETAXES(DEtavsEta, "#eta", "#Delta#eta");
0198   SETAXES(DPhivsEta, "#eta", "#Delta#phi");
0199   // delta Pt or E quantities for Barrel and Endcap
0200   DSETAXES(RPt, "#DeltaP_{T}/P_{T}", "Events");
0201   DSETAXES(RPt20_40, "#DeltaP_{T}/P_{T}", "Events");
0202   DSETAXES(RPt40_60, "#DeltaP_{T}/P_{T}", "Events");
0203   DSETAXES(RPt60_80, "#DeltaP_{T}/P_{T}", "Events");
0204   DSETAXES(RPt80_100, "#DeltaP_{T}/P_{T}", "Events");
0205   DSETAXES(RPt100_150, "#DeltaP_{T}/P_{T}", "Events");
0206   DSETAXES(RPt150_200, "#DeltaP_{T}/P_{T}", "Events");
0207   DSETAXES(RPt200_250, "#DeltaP_{T}/P_{T}", "Events");
0208   DSETAXES(RPt250_300, "#DeltaP_{T}/P_{T}", "Events");
0209   DSETAXES(RPt300_400, "#DeltaP_{T}/P_{T}", "Events");
0210   DSETAXES(RPt400_500, "#DeltaP_{T}/P_{T}", "Events");
0211   DSETAXES(RPt500_750, "#DeltaP_{T}/P_{T}", "Events");
0212   DSETAXES(RPt750_1250, "#DeltaP_{T}/P_{T}", "Events");
0213   DSETAXES(RPt1250_2000, "#DeltaP_{T}/P_{T}", "Events");
0214   DSETAXES(RPt2000_5000, "#DeltaP_{T}/P_{T}", "Events");
0215   DSETAXES(RCHE, "#DeltaE/E(charged had)", "Events");
0216   DSETAXES(RNHE, "#DeltaE/E(neutral had)", "Events");
0217   DSETAXES(RNEE, "#DeltaE/E(neutral em)", "Events");
0218   DSETAXES(Rneut, "#DeltaE/E(neutral)", "Events");
0219   DSETAXES(RPtvsPt, PT, "#DeltaP_{T}/P_{T}");
0220   DSETAXES(RCHEvsPt, PT, "#DeltaE/E(charged had)");
0221   DSETAXES(RNHEvsPt, PT, "#DeltaE/E(neutral had)");
0222   DSETAXES(RNEEvsPt, PT, "#DeltaE/E(neutral em)");
0223   DSETAXES(RneutvsPt, PT, "#DeltaE/E(neutral)");
0224   DSETAXES(RHCALvsP, P, "#DeltaE/E(HCAL)");
0225   DSETAXES(RHONLvsP, P, "#DeltaE/E(HCAL-only)");
0226   DSETAXES(RNEUTvsP, P, "#DeltaE/E(ECAL+HCAL)");
0227   DSETAXES(RNONLvsP, P, "#DeltaE/E(ECAL+HCAL-only)");
0228   DSETAXES(DEtavsPt, PT, "#Delta#eta");
0229   DSETAXES(DPhivsPt, PT, "#Delta#phi");
0230 }
0231 
0232 void PFJetBenchmark::process(const reco::PFJetCollection& pfJets, const reco::GenJetCollection& genJets) {
0233   // loop over reco  pf  jets
0234   resPtMax_ = 0.;
0235   resChargedHadEnergyMax_ = 0.;
0236   resNeutralHadEnergyMax_ = 0.;
0237   resNeutralEmEnergyMax_ = 0.;
0238   int NPFJets = 0;
0239 
0240   for (unsigned i = 0; i < pfJets.size(); i++) {
0241     // Count the number of jets with a larger energy
0242     unsigned highJets = 0;
0243     for (unsigned j = 0; j < pfJets.size(); j++) {
0244       if (j != i && pfJets[j].pt() > pfJets[i].pt())
0245         highJets++;
0246     }
0247     if (onlyTwoJets_ && highJets > 1)
0248       continue;
0249 
0250     const reco::PFJet& pfj = pfJets[i];
0251     double rec_pt = pfj.pt();
0252     double rec_eta = pfj.eta();
0253     double rec_phi = pfj.phi();
0254 
0255     // skip PFjets with pt < recPt_cut GeV
0256     if (rec_pt < recPt_cut and recPt_cut != -1.)
0257       continue;
0258     // skip PFjets with eta > maxEta_cut
0259     if (fabs(rec_eta) > maxEta_cut and maxEta_cut != -1.)
0260       continue;
0261 
0262     NPFJets++;
0263 
0264     // fill inclusive PFjet distribution pt > 20 GeV
0265     hNjets->Fill(NPFJets);
0266     hjetsPt->Fill(rec_pt);
0267     hjetsEta->Fill(rec_eta);
0268 
0269     // separate Barrel PFJets from Endcap PFJets
0270     bool Barrel = false;
0271     bool Endcap = false;
0272     bool Forward = false;
0273     if (std::abs(rec_eta) < 1.4)
0274       Barrel = true;
0275     if (std::abs(rec_eta) > 1.6 && std::abs(rec_eta) < 2.4)
0276       Endcap = true;
0277     if (std::abs(rec_eta) > 2.5 && std::abs(rec_eta) < 2.9)
0278       Forward = true;
0279     if (std::abs(rec_eta) > 3.1 && std::abs(rec_eta) < 4.7)
0280       Forward = true;
0281 
0282     // do only barrel for now
0283     //  if(!Barrel) continue;
0284 
0285     // look for the closets gen Jet : truth
0286     const GenJet* truth = algo_->matchByDeltaR(&pfj, &genJets);
0287     if (!truth)
0288       continue;
0289     double deltaR = algo_->deltaR(&pfj, truth);
0290     // check deltaR is small enough
0291     if (deltaR < deltaRMax_ || (abs(rec_eta) > 2.5 && deltaR < 0.2) ||
0292         deltaRMax_ == -1.0) {  //start case deltaR < deltaRMax
0293 
0294       // generate histograms comparing the reco and truth candidate (truth = closest in delta-R)
0295       // get the quantities to place on the denominator and/or divide by
0296       double pt_denom;
0297       double true_E = truth->p();
0298       double true_pt = truth->pt();
0299       double true_eta = truth->eta();
0300       double true_phi = truth->phi();
0301 
0302       if (plotAgainstReco_) {
0303         pt_denom = rec_pt;
0304       } else {
0305         pt_denom = true_pt;
0306       }
0307       // get true specific quantities
0308       double true_ChargedHadEnergy;
0309       double true_NeutralHadEnergy;
0310       double true_NeutralEmEnergy;
0311       gettrue(truth, true_ChargedHadEnergy, true_NeutralHadEnergy, true_NeutralEmEnergy);
0312       double true_NeutralEnergy = true_NeutralHadEnergy + true_NeutralEmEnergy;
0313       double rec_ChargedHadEnergy = pfj.chargedHadronEnergy();
0314       double rec_NeutralHadEnergy = pfj.neutralHadronEnergy();
0315       double rec_NeutralEmEnergy = pfj.neutralEmEnergy();
0316       double rec_NeutralEnergy = rec_NeutralHadEnergy + rec_NeutralEmEnergy;
0317       double rec_ChargedMultiplicity = pfj.chargedMultiplicity();
0318       std::vector<PFCandidatePtr> constituents = pfj.getPFConstituents();
0319       std::vector<unsigned int> chMult(9, static_cast<unsigned int>(0));
0320       for (unsigned ic = 0; ic < constituents.size(); ++ic) {
0321         if (constituents[ic]->particleId() > 3)
0322           continue;
0323         reco::TrackRef trackRef = constituents[ic]->trackRef();
0324         if (trackRef.isNull()) {
0325           //std::cout << "Warning in entry " << entry_
0326           //        << " : a track with Id " << constituents[ic]->particleId()
0327           //        << " has no track ref.." << std::endl;
0328           continue;
0329         }
0330         unsigned int iter = trackRef->algo() > 6 ? 6 : trackRef->algo();
0331         ++(chMult[iter]);
0332       }
0333 
0334       bool plot1 = false;
0335       bool plot2 = false;
0336       bool plot3 = false;
0337       bool plot4 = false;
0338       bool plot5 = false;
0339       bool plot6 = false;
0340       bool plot7 = false;
0341       double cut1 = 0.0001;
0342       double cut2 = 0.0001;
0343       double cut3 = 0.0001;
0344       double cut4 = 0.0001;
0345       double cut5 = 0.0001;
0346       double cut6 = 0.0001;
0347       double cut7 = 0.0001;
0348       double resPt = 0.;
0349       double resChargedHadEnergy = 0.;
0350       double resNeutralHadEnergy = 0.;
0351       double resNeutralEmEnergy = 0.;
0352       double resNeutralEnergy = 0.;
0353 
0354       double resHCALEnergy = 0.;
0355       double resNEUTEnergy = 0.;
0356       if (rec_NeutralHadEnergy > cut6 && rec_ChargedHadEnergy < cut1) {
0357         double true_NEUTEnergy = true_NeutralHadEnergy + true_NeutralEmEnergy;
0358         double true_HCALEnergy = true_NEUTEnergy - rec_NeutralEmEnergy;
0359         double rec_NEUTEnergy = rec_NeutralHadEnergy + rec_NeutralEmEnergy;
0360         double rec_HCALEnergy = rec_NeutralHadEnergy;
0361         resHCALEnergy = (rec_HCALEnergy - true_HCALEnergy) / rec_HCALEnergy;
0362         resNEUTEnergy = (rec_NEUTEnergy - true_NEUTEnergy) / rec_NEUTEnergy;
0363         if (rec_NeutralEmEnergy > cut7) {
0364           plot6 = true;
0365         } else {
0366           plot7 = true;
0367         }
0368       }
0369 
0370       // get relative delta quantities (protect against division by zero!)
0371       if (true_pt > 0.0001) {
0372         resPt = (rec_pt - true_pt) / true_pt;
0373         plot1 = true;
0374       }
0375       if (true_ChargedHadEnergy > cut1) {
0376         resChargedHadEnergy = (rec_ChargedHadEnergy - true_ChargedHadEnergy) / true_ChargedHadEnergy;
0377         plot2 = true;
0378       }
0379       if (true_NeutralHadEnergy > cut2) {
0380         resNeutralHadEnergy = (rec_NeutralHadEnergy - true_NeutralHadEnergy) / true_NeutralHadEnergy;
0381         plot3 = true;
0382       } else if (rec_NeutralHadEnergy > cut3) {
0383         resNeutralHadEnergy = (rec_NeutralHadEnergy - true_NeutralHadEnergy) / rec_NeutralHadEnergy;
0384         plot3 = true;
0385       }
0386       if (true_NeutralEmEnergy > cut4) {
0387         resNeutralEmEnergy = (rec_NeutralEmEnergy - true_NeutralEmEnergy) / true_NeutralEmEnergy;
0388         plot4 = true;
0389       }
0390       if (true_NeutralEnergy > cut5) {
0391         resNeutralEnergy = (rec_NeutralEnergy - true_NeutralEnergy) / true_NeutralEnergy;
0392         plot5 = true;
0393       }
0394 
0395       //double deltaEta = algo_->deltaEta(&pfj, truth);
0396       //double deltaPhi = algo_->deltaPhi(&pfj, truth);
0397 
0398       // Print outliers for further debugging
0399       if ((resPt > 0.2 && true_pt > 100.) || (resPt < -0.5 && true_pt > 100.)) {
0400         //if ( ( true_pt > 50. &&
0401         //     ( ( truth->eta()>3.0 && rec_eta-truth->eta() < -0.1 ) ||
0402         //       ( truth->eta()<-3.0 && rec_eta-truth->eta() > 0.1 ) ))) {
0403         std::cout << "Entry " << entry_ << " resPt = " << resPt << " resCharged  " << resChargedHadEnergy
0404                   << " resNeutralHad  " << resNeutralHadEnergy << " resNeutralEm  " << resNeutralEmEnergy
0405                   << " pT (T/R) " << true_pt << "/" << rec_pt << " Eta (T/R) " << truth->eta() << "/" << rec_eta
0406                   << " Phi (T/R) " << truth->phi() << "/" << rec_phi << std::endl;
0407 
0408         // check overlapping PF jets
0409         const reco::PFJet* pfoj = nullptr;
0410         double dRo = 1E9;
0411         for (unsigned j = 0; j < pfJets.size(); j++) {
0412           const reco::PFJet& pfo = pfJets[j];
0413           if (j != i && algo_->deltaR(&pfj, &pfo) < dRo && pfo.pt() > 0.25 * pfj.pt()) {
0414             dRo = algo_->deltaR(&pfj, &pfo);
0415             pfoj = &pfo;
0416           }
0417         }
0418 
0419         // Check overlapping Gen Jet
0420         math::XYZTLorentzVector overlappinGenJet(0., 0., 0., 0.);
0421         const reco::GenJet* genoj = nullptr;
0422         double dRgo = 1E9;
0423         for (unsigned j = 0; j < genJets.size(); j++) {
0424           const reco::GenJet* gjo = &(genJets[j]);
0425           if (gjo != truth && algo_->deltaR(truth, gjo) < dRgo && gjo->pt() > 0.25 * truth->pt()) {
0426             dRgo = algo_->deltaR(truth, gjo);
0427             genoj = gjo;
0428           }
0429         }
0430 
0431         if (dRo < 0.8 && dRgo < 0.8 && algo_->deltaR(genoj, pfoj) < 2. * deltaRMax_)
0432           std::cout << "Excess probably due to overlapping jets (DR = " << algo_->deltaR(genoj, pfoj) << "),"
0433                     << " at DeltaR(T/R) = " << dRgo << "/" << dRo << " with pT(T/R) " << genoj->pt() << "/"
0434                     << pfoj->pt() << " and Eta (T/R) " << genoj->eta() << "/" << pfoj->eta() << " and Phi (T/R) "
0435                     << genoj->phi() << "/" << pfoj->phi() << std::endl;
0436       }
0437 
0438       if (std::abs(resPt) > std::abs(resPtMax_))
0439         resPtMax_ = resPt;
0440       if (std::abs(resChargedHadEnergy) > std::abs(resChargedHadEnergyMax_))
0441         resChargedHadEnergyMax_ = resChargedHadEnergy;
0442       if (std::abs(resNeutralHadEnergy) > std::abs(resNeutralHadEnergyMax_))
0443         resNeutralHadEnergyMax_ = resNeutralHadEnergy;
0444       if (std::abs(resNeutralEmEnergy) > std::abs(resNeutralEmEnergyMax_))
0445         resNeutralEmEnergyMax_ = resNeutralEmEnergy;
0446       if (debug_) {
0447         cout << i << "  =========PFJet Pt " << rec_pt << " eta " << rec_eta << " phi " << rec_phi
0448              << " Charged Had Energy " << rec_ChargedHadEnergy << " Neutral Had Energy " << rec_NeutralHadEnergy
0449              << " Neutral elm Energy " << rec_NeutralEmEnergy << endl;
0450         cout << " matching Gen Jet Pt " << true_pt << " eta " << truth->eta() << " phi " << truth->phi()
0451              << " Charged Had Energy " << true_ChargedHadEnergy << " Neutral Had Energy " << true_NeutralHadEnergy
0452              << " Neutral elm Energy " << true_NeutralEmEnergy << endl;
0453         printPFJet(&pfj);
0454         //      cout<<pfj.print()<<endl;
0455         printGenJet(truth);
0456         //cout <<truth->print()<<endl;
0457 
0458         cout << "==============deltaR " << deltaR << "  resPt " << resPt << " resChargedHadEnergy "
0459              << resChargedHadEnergy << " resNeutralHadEnergy " << resNeutralHadEnergy << " resNeutralEmEnergy "
0460              << resNeutralEmEnergy << endl;
0461       }
0462 
0463       if (plot1) {
0464         if (rec_eta > 0.)
0465           hDEtavsEta->Fill(true_eta, rec_eta - true_eta);
0466         else
0467           hDEtavsEta->Fill(true_eta, -rec_eta + true_eta);
0468         hDPhivsEta->Fill(true_eta, rec_phi - true_phi);
0469 
0470         hRPtvsEta->Fill(true_eta, resPt);
0471         hNCHvsEta->Fill(true_eta, rec_ChargedMultiplicity);
0472         hNCH0vsEta->Fill(true_eta, chMult[0]);
0473         hNCH1vsEta->Fill(true_eta, chMult[1]);
0474         hNCH2vsEta->Fill(true_eta, chMult[2]);
0475         hNCH3vsEta->Fill(true_eta, chMult[3]);
0476         hNCH4vsEta->Fill(true_eta, chMult[4]);
0477         hNCH5vsEta->Fill(true_eta, chMult[5]);
0478         hNCH6vsEta->Fill(true_eta, chMult[6]);
0479         hNCH7vsEta->Fill(true_eta, chMult[7]);
0480       }
0481       if (plot2)
0482         hRCHEvsEta->Fill(true_eta, resChargedHadEnergy);
0483       if (plot5)
0484         hRNeutvsEta->Fill(true_eta, resNeutralEnergy);
0485       if (plot6) {
0486         hRHCALvsEta->Fill(true_eta, resHCALEnergy);
0487         hRNEUTvsEta->Fill(true_eta, resNEUTEnergy);
0488       }
0489       if (plot7) {
0490         hRHONLvsEta->Fill(true_eta, resHCALEnergy);
0491         hRNONLvsEta->Fill(true_eta, resNEUTEnergy);
0492       }
0493 
0494       // fill histograms for relative delta quantitites of matched jets
0495       // delta Pt or E quantities for Barrel
0496       if (Barrel) {
0497         if (plot1) {
0498           hBRPt->Fill(resPt);
0499           if (pt_denom > 20. && pt_denom < 40.)
0500             hBRPt20_40->Fill(resPt);
0501           if (pt_denom > 40. && pt_denom < 60.)
0502             hBRPt40_60->Fill(resPt);
0503           if (pt_denom > 60. && pt_denom < 80.)
0504             hBRPt60_80->Fill(resPt);
0505           if (pt_denom > 80. && pt_denom < 100.)
0506             hBRPt80_100->Fill(resPt);
0507           if (pt_denom > 100. && pt_denom < 150.)
0508             hBRPt100_150->Fill(resPt);
0509           if (pt_denom > 150. && pt_denom < 200.)
0510             hBRPt150_200->Fill(resPt);
0511           if (pt_denom > 200. && pt_denom < 250.)
0512             hBRPt200_250->Fill(resPt);
0513           if (pt_denom > 250. && pt_denom < 300.)
0514             hBRPt250_300->Fill(resPt);
0515           if (pt_denom > 300. && pt_denom < 400.)
0516             hBRPt300_400->Fill(resPt);
0517           if (pt_denom > 400. && pt_denom < 500.)
0518             hBRPt400_500->Fill(resPt);
0519           if (pt_denom > 500. && pt_denom < 750.)
0520             hBRPt500_750->Fill(resPt);
0521           if (pt_denom > 750. && pt_denom < 1250.)
0522             hBRPt750_1250->Fill(resPt);
0523           if (pt_denom > 1250. && pt_denom < 2000.)
0524             hBRPt1250_2000->Fill(resPt);
0525           if (pt_denom > 2000. && pt_denom < 5000.)
0526             hBRPt2000_5000->Fill(resPt);
0527           hBNCH->Fill(rec_ChargedMultiplicity);
0528           hBNCH0vsPt->Fill(pt_denom, chMult[0]);
0529           hBNCH1vsPt->Fill(pt_denom, chMult[1]);
0530           hBNCH2vsPt->Fill(pt_denom, chMult[2]);
0531           hBNCH3vsPt->Fill(pt_denom, chMult[3]);
0532           hBNCH4vsPt->Fill(pt_denom, chMult[4]);
0533           hBNCH5vsPt->Fill(pt_denom, chMult[5]);
0534           hBNCH6vsPt->Fill(pt_denom, chMult[6]);
0535           hBNCH7vsPt->Fill(pt_denom, chMult[7]);
0536           hBNCHvsPt->Fill(pt_denom, rec_ChargedMultiplicity);
0537           if (rec_eta > 0.)
0538             hBDEtavsPt->Fill(pt_denom, rec_eta - true_eta);
0539           else
0540             hBDEtavsPt->Fill(pt_denom, -rec_eta + true_eta);
0541           hBDPhivsPt->Fill(pt_denom, rec_phi - true_phi);
0542         }
0543         if (plot2)
0544           hBRCHE->Fill(resChargedHadEnergy);
0545         if (plot3)
0546           hBRNHE->Fill(resNeutralHadEnergy);
0547         if (plot4)
0548           hBRNEE->Fill(resNeutralEmEnergy);
0549         if (plot5)
0550           hBRneut->Fill(resNeutralEnergy);
0551         if (plot1)
0552           hBRPtvsPt->Fill(pt_denom, resPt);
0553         if (plot2)
0554           hBRCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
0555         if (plot3)
0556           hBRNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
0557         if (plot4)
0558           hBRNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
0559         if (plot5)
0560           hBRneutvsPt->Fill(pt_denom, resNeutralEnergy);
0561         if (plot6) {
0562           hBRHCALvsP->Fill(true_E, resHCALEnergy);
0563           hBRNEUTvsP->Fill(true_E, resNEUTEnergy);
0564         }
0565         if (plot7) {
0566           hBRHONLvsP->Fill(true_E, resHCALEnergy);
0567           hBRNONLvsP->Fill(true_E, resNEUTEnergy);
0568         }
0569       }
0570       // delta Pt or E quantities for Endcap
0571       if (Endcap) {
0572         if (plot1) {
0573           hERPt->Fill(resPt);
0574           if (pt_denom > 20. && pt_denom < 40.)
0575             hERPt20_40->Fill(resPt);
0576           if (pt_denom > 40. && pt_denom < 60.)
0577             hERPt40_60->Fill(resPt);
0578           if (pt_denom > 60. && pt_denom < 80.)
0579             hERPt60_80->Fill(resPt);
0580           if (pt_denom > 80. && pt_denom < 100.)
0581             hERPt80_100->Fill(resPt);
0582           if (pt_denom > 100. && pt_denom < 150.)
0583             hERPt100_150->Fill(resPt);
0584           if (pt_denom > 150. && pt_denom < 200.)
0585             hERPt150_200->Fill(resPt);
0586           if (pt_denom > 200. && pt_denom < 250.)
0587             hERPt200_250->Fill(resPt);
0588           if (pt_denom > 250. && pt_denom < 300.)
0589             hERPt250_300->Fill(resPt);
0590           if (pt_denom > 300. && pt_denom < 400.)
0591             hERPt300_400->Fill(resPt);
0592           if (pt_denom > 400. && pt_denom < 500.)
0593             hERPt400_500->Fill(resPt);
0594           if (pt_denom > 500. && pt_denom < 750.)
0595             hERPt500_750->Fill(resPt);
0596           if (pt_denom > 750. && pt_denom < 1250.)
0597             hERPt750_1250->Fill(resPt);
0598           if (pt_denom > 1250. && pt_denom < 2000.)
0599             hERPt1250_2000->Fill(resPt);
0600           if (pt_denom > 2000. && pt_denom < 5000.)
0601             hERPt2000_5000->Fill(resPt);
0602           hENCH->Fill(rec_ChargedMultiplicity);
0603           hENCHvsPt->Fill(pt_denom, rec_ChargedMultiplicity);
0604           hENCH0vsPt->Fill(pt_denom, chMult[0]);
0605           hENCH1vsPt->Fill(pt_denom, chMult[1]);
0606           hENCH2vsPt->Fill(pt_denom, chMult[2]);
0607           hENCH3vsPt->Fill(pt_denom, chMult[3]);
0608           hENCH4vsPt->Fill(pt_denom, chMult[4]);
0609           hENCH5vsPt->Fill(pt_denom, chMult[5]);
0610           hENCH6vsPt->Fill(pt_denom, chMult[6]);
0611           hENCH7vsPt->Fill(pt_denom, chMult[7]);
0612           if (rec_eta > 0.)
0613             hEDEtavsPt->Fill(pt_denom, rec_eta - true_eta);
0614           else
0615             hEDEtavsPt->Fill(pt_denom, -rec_eta + true_eta);
0616           hEDPhivsPt->Fill(pt_denom, rec_phi - true_phi);
0617         }
0618         if (plot2)
0619           hERCHE->Fill(resChargedHadEnergy);
0620         if (plot3)
0621           hERNHE->Fill(resNeutralHadEnergy);
0622         if (plot4)
0623           hERNEE->Fill(resNeutralEmEnergy);
0624         if (plot5)
0625           hERneut->Fill(resNeutralEnergy);
0626         if (plot1)
0627           hERPtvsPt->Fill(pt_denom, resPt);
0628         if (plot2)
0629           hERCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
0630         if (plot3)
0631           hERNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
0632         if (plot4)
0633           hERNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
0634         if (plot5)
0635           hERneutvsPt->Fill(pt_denom, resNeutralEnergy);
0636         if (plot6) {
0637           hERHCALvsP->Fill(true_E, resHCALEnergy);
0638           hERNEUTvsP->Fill(true_E, resNEUTEnergy);
0639         }
0640         if (plot7) {
0641           hERHONLvsP->Fill(true_E, resHCALEnergy);
0642           hERNONLvsP->Fill(true_E, resNEUTEnergy);
0643         }
0644       }
0645       // delta Pt or E quantities for Forward
0646       if (Forward) {
0647         if (plot1) {
0648           hFRPt->Fill(resPt);
0649           if (pt_denom > 20. && pt_denom < 40.)
0650             hFRPt20_40->Fill(resPt);
0651           if (pt_denom > 40. && pt_denom < 60.)
0652             hFRPt40_60->Fill(resPt);
0653           if (pt_denom > 60. && pt_denom < 80.)
0654             hFRPt60_80->Fill(resPt);
0655           if (pt_denom > 80. && pt_denom < 100.)
0656             hFRPt80_100->Fill(resPt);
0657           if (pt_denom > 100. && pt_denom < 150.)
0658             hFRPt100_150->Fill(resPt);
0659           if (pt_denom > 150. && pt_denom < 200.)
0660             hFRPt150_200->Fill(resPt);
0661           if (pt_denom > 200. && pt_denom < 250.)
0662             hFRPt200_250->Fill(resPt);
0663           if (pt_denom > 250. && pt_denom < 300.)
0664             hFRPt250_300->Fill(resPt);
0665           if (pt_denom > 300. && pt_denom < 400.)
0666             hFRPt300_400->Fill(resPt);
0667           if (pt_denom > 400. && pt_denom < 500.)
0668             hFRPt400_500->Fill(resPt);
0669           if (pt_denom > 500. && pt_denom < 750.)
0670             hFRPt500_750->Fill(resPt);
0671           if (pt_denom > 750. && pt_denom < 1250.)
0672             hFRPt750_1250->Fill(resPt);
0673           if (pt_denom > 1250. && pt_denom < 2000.)
0674             hFRPt1250_2000->Fill(resPt);
0675           if (pt_denom > 2000. && pt_denom < 5000.)
0676             hFRPt2000_5000->Fill(resPt);
0677           if (rec_eta > 0.)
0678             hFDEtavsPt->Fill(pt_denom, rec_eta - true_eta);
0679           else
0680             hFDEtavsPt->Fill(pt_denom, -rec_eta + true_eta);
0681           hFDPhivsPt->Fill(pt_denom, rec_phi - true_phi);
0682         }
0683         if (plot2)
0684           hFRCHE->Fill(resChargedHadEnergy);
0685         if (plot3)
0686           hFRNHE->Fill(resNeutralHadEnergy);
0687         if (plot4)
0688           hFRNEE->Fill(resNeutralEmEnergy);
0689         if (plot5)
0690           hFRneut->Fill(resNeutralEnergy);
0691         if (plot1)
0692           hFRPtvsPt->Fill(pt_denom, resPt);
0693         if (plot2)
0694           hFRCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
0695         if (plot3)
0696           hFRNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
0697         if (plot4)
0698           hFRNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
0699         if (plot5)
0700           hFRneutvsPt->Fill(pt_denom, resNeutralEnergy);
0701         if (plot6) {
0702           hFRHCALvsP->Fill(true_E, resHCALEnergy);
0703           hFRNEUTvsP->Fill(true_E, resNEUTEnergy);
0704         }
0705         if (plot7) {
0706           hFRHONLvsP->Fill(true_E, resHCALEnergy);
0707           hFRNONLvsP->Fill(true_E, resNEUTEnergy);
0708         }
0709       }
0710     }  // end case deltaR < deltaRMax
0711 
0712   }  // i loop on pf Jets
0713 
0714   // Increment counter
0715   entry_++;
0716 }
0717 
0718 void PFJetBenchmark::gettrue(const reco::GenJet* truth,
0719                              double& true_ChargedHadEnergy,
0720                              double& true_NeutralHadEnergy,
0721                              double& true_NeutralEmEnergy) {
0722   std::vector<const GenParticle*> mcparts = truth->getGenConstituents();
0723   true_NeutralEmEnergy = 0.;
0724   true_ChargedHadEnergy = 0.;
0725   true_NeutralHadEnergy = 0.;
0726   // for each MC particle in turn
0727   for (unsigned i = 0; i < mcparts.size(); i++) {
0728     const GenParticle* mcpart = mcparts[i];
0729     int PDG = std::abs(mcpart->pdgId());
0730     double e = mcpart->energy();
0731     switch (PDG) {  // start PDG switch
0732       case 22:      // photon
0733         true_NeutralEmEnergy += e;
0734         break;
0735       case 211:   // pi
0736       case 321:   // K
0737       case 2212:  // p
0738       case 11:    //electrons (until recognised)
0739         true_ChargedHadEnergy += e;
0740         break;
0741       case 310:   // K_S0
0742       case 130:   // K_L0
0743       case 3122:  // Lambda0
0744       case 2112:  // n0
0745         true_NeutralHadEnergy += e;
0746       default:
0747         break;
0748     }  // end PDG switch
0749   }  // end loop on constituents.
0750 }
0751 
0752 void PFJetBenchmark::printPFJet(const reco::PFJet* pfj) {
0753   cout << setiosflags(ios::right);
0754   cout << setiosflags(ios::fixed);
0755   cout << setprecision(3);
0756 
0757   cout << "PFJet  p/px/py/pz/pt: " << pfj->p() << "/" << pfj->px() << "/" << pfj->py() << "/" << pfj->pz() << "/"
0758        << pfj->pt() << endl
0759        << "    eta/phi: " << pfj->eta() << "/" << pfj->phi() << endl
0760        << "    PFJet specific:" << std::endl
0761        << "      charged/neutral hadrons energy: " << pfj->chargedHadronEnergy() << '/' << pfj->neutralHadronEnergy()
0762        << endl
0763        << "      charged/neutral em energy: " << pfj->chargedEmEnergy() << '/' << pfj->neutralEmEnergy() << endl
0764        << "      charged muon energy: " << pfj->chargedMuEnergy() << '/' << endl
0765        << "      charged/neutral multiplicity: " << pfj->chargedMultiplicity() << '/' << pfj->neutralMultiplicity()
0766        << endl;
0767 
0768   // And print the constituents
0769   std::cout << pfj->print() << std::endl;
0770 
0771   cout << resetiosflags(ios::right | ios::fixed);
0772 }
0773 
0774 void PFJetBenchmark::printGenJet(const reco::GenJet* truth) {
0775   std::vector<const GenParticle*> mcparts = truth->getGenConstituents();
0776   cout << "GenJet p/px/py/pz/pt: " << truth->p() << '/' << truth->px() << '/' << truth->py() << '/' << truth->pz()
0777        << '/' << truth->pt() << endl
0778        << "    eta/phi: " << truth->eta() << '/' << truth->phi() << endl
0779        << "    # of constituents: " << mcparts.size() << endl;
0780   cout << "    constituents:" << endl;
0781   for (unsigned i = 0; i < mcparts.size(); i++) {
0782     const GenParticle* mcpart = mcparts[i];
0783     cout << "      #" << i << "  PDG code:" << mcpart->pdgId() << ", p/pt/eta/phi: " << mcpart->p() << '/'
0784          << mcpart->pt() << '/' << mcpart->eta() << '/' << mcpart->phi() << endl;
0785   }
0786 }
[ 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!