CMS 3D CMS Logo

PFJetBenchmark.cc
Go to the documentation of this file.
3 
7 
8 // preprocessor macro for booking 1d histos with DQMStore -or- bare Root
9 #define BOOK1D(name,title,nbinsx,lowx,highx) \
10  h##name = dbe_ ? dbe_->book1D(#name,title,nbinsx,lowx,highx)->getTH1F() \
11  : new TH1F(#name,title,nbinsx,lowx,highx)
12 
13 // preprocessor macro for booking 2d histos with DQMStore -or- bare Root
14 #define BOOK2D(name,title,nbinsx,lowx,highx,nbinsy,lowy,highy) \
15  h##name = dbe_ ? dbe_->book2D(#name,title,nbinsx,lowx,highx,nbinsy,lowy,highy)->getTH2F() \
16  : new TH2F(#name,title,nbinsx,lowx,highx,nbinsy,lowy,highy)
17 
18 //macros for building barrel and endcap histos with one call
19 #define DBOOK1D(name,title,nbinsx,lowx,highx) \
20  BOOK1D(B##name,"Barrel "#title,nbinsx,lowx,highx); BOOK1D(E##name,"Endcap "#title,nbinsx,lowx,highx); BOOK1D(F##name,"Forward "#title,nbinsx,lowx,highx);
21 #define DBOOK2D(name,title,nbinsx,lowx,highx,nbinsy,lowy,highy) \
22  BOOK2D(B##name,"Barrel "#title,nbinsx,lowx,highx,nbinsy,lowy,highy); BOOK2D(E##name,"Endcap "#title,nbinsx,lowx,highx,nbinsy,lowy,highy); BOOK2D(F##name,"Forward "#title,nbinsx,lowx,highx,nbinsy,lowy,highy);
23 
24 // all versions OK
25 // preprocesor macro for setting axis titles
26 #define SETAXES(name,xtitle,ytitle) \
27  h##name->GetXaxis()->SetTitle(xtitle); h##name->GetYaxis()->SetTitle(ytitle)
28 
29 //macro for setting the titles for barrel and endcap together
30 #define DSETAXES(name,xtitle,ytitle) \
31  SETAXES(B##name,xtitle,ytitle);SETAXES(E##name,xtitle,ytitle);SETAXES(F##name,xtitle,ytitle);
32 /*#define SET2AXES(name,xtitle,ytitle) \
33  hE##name->GetXaxis()->SetTitle(xtitle); hE##name->GetYaxis()->SetTitle(ytitle); hB##name->GetXaxis()->SetTitle(xtitle); hB##name->GetYaxis()->SetTitle(ytitle)
34 */
35 
36 #define PT (plotAgainstReco_)?"reconstructed P_{T}" :"generated P_{T}"
37 #define P (plotAgainstReco_)?"generated P" :"generated P"
38 
39 using namespace reco;
40 using namespace std;
41 
42 PFJetBenchmark::PFJetBenchmark() : file_(nullptr), entry_(0) {}
43 
45  if(file_) file_->Close();
46 }
47 
49  // Store the DAQ Histograms
50  if (!outputFile_.empty()) {
51  if (dbe_)
52  dbe_->save(outputFile_);
53  // use bare Root if no DQM (FWLite applications)
54  else if (file_) {
55  file_->Write(outputFile_.c_str());
56  cout << "Benchmark output written to file " << outputFile_.c_str() << endl;
57  file_->Close();
58  }
59  }
60  else
61  cout << "No output file specified ("<<outputFile_<<"). Results will not be saved!" << endl;
62 
63 }
64 
66  string Filename,
67  bool debug,
68  bool plotAgainstReco,
69  bool onlyTwoJets,
70  double deltaRMax,
71  string benchmarkLabel_,
72  double recPt,
73  double maxEta,
74  DQMStore * dbe_store) {
75  debug_ = debug;
79  outputFile_=Filename;
80  recPt_cut = recPt;
82  file_ = nullptr;
83  dbe_ = dbe_store;
84  // print parameters
85  cout<< "PFJetBenchmark Setup parameters =============================================="<<endl;
86  cout << "Filename to write histograms " << Filename<<endl;
87  cout << "PFJetBenchmark debug " << debug_<< endl;
88  cout << "plotAgainstReco " << plotAgainstReco_ << endl;
89  cout << "onlyTwoJets " << onlyTwoJets_ << endl;
90  cout << "deltaRMax " << deltaRMax << endl;
91  cout << "benchmarkLabel " << benchmarkLabel_ << endl;
92  cout << "recPt_cut " << recPt_cut << endl;
93  cout << "maxEta_cut " << maxEta_cut << endl;
94 
95  // Book histogram
96 
97  // Establish DQM Store
98  string path = "PFTask/Benchmarks/"+ benchmarkLabel_ + "/";
99  if (plotAgainstReco) path += "Reco"; else path += "Gen";
100  if (dbe_) {
101  dbe_->setCurrentFolder(path);
102  }
103  else {
104  file_ = new TFile(outputFile_.c_str(), "recreate");
105 // TTree * tr = new TTree("PFTast");
106 // tr->Branch("Benchmarks/ParticleFlow")
107  cout << "Info: DQM is not available to provide data storage service. Using TFile to save histograms. "<<endl;
108  }
109  // Jets inclusive distributions (Pt > 20 or specified recPt GeV)
110  char cutString[35];
111  sprintf(cutString,"Jet multiplicity P_{T}>%4.1f GeV", recPt_cut);
112  BOOK1D(Njets,cutString,50, 0, 50);
113 
114  BOOK1D(jetsPt,"Jets P_{T} Distribution",100, 0, 500);
115 
116  sprintf(cutString,"Jets #eta Distribution |#eta|<%4.1f", maxEta_cut);
117  BOOK1D(jetsEta,cutString,100, -5, 5);
118 
119  BOOK2D(RPtvsEta,"#DeltaP_{T}/P_{T} vs #eta",200, -5., 5., 200,-1,1);
120  BOOK2D(RNeutvsEta,"R_{Neutral} vs #eta",200, -5., 5., 200,-1,1);
121  BOOK2D(RNEUTvsEta,"R_{HCAL+ECAL} vs #eta",200, -5., 5., 200,-1,1);
122  BOOK2D(RNONLvsEta,"R_{HCAL+ECAL - Hcal Only} vs #eta",200, -5., 5., 200,-1,1);
123  BOOK2D(RHCALvsEta,"R_{HCAL} vs #eta",200, -5., 5., 200,-1,1);
124  BOOK2D(RHONLvsEta,"R_{HCAL only} vs #eta",200, -5., 5., 200,-1,1);
125  BOOK2D(RCHEvsEta,"R_{Charged} vs #eta",200, -5., 5., 200,-1,1);
126  BOOK2D(NCHvsEta,"N_{Charged} vs #eta",200, -5., 5., 200,0.,200.);
127  BOOK2D(NCH0vsEta,"N_{Charged} vs #eta, iter 0",200, -5., 5., 200,0.,200.);
128  BOOK2D(NCH1vsEta,"N_{Charged} vs #eta, iter 1",200, -5., 5., 200,0.,200.);
129  BOOK2D(NCH2vsEta,"N_{Charged} vs #eta, iter 2",200, -5., 5., 200,0.,200.);
130  BOOK2D(NCH3vsEta,"N_{Charged} vs #eta, iter 3",200, -5., 5., 200,0.,200.);
131  BOOK2D(NCH4vsEta,"N_{Charged} vs #eta, iter 4",200, -5., 5., 200,0.,200.);
132  BOOK2D(NCH5vsEta,"N_{Charged} vs #eta, iter 5",200, -5., 5., 200,0.,200.);
133  BOOK2D(NCH6vsEta,"N_{Charged} vs #eta, iter 6",200, -5., 5., 200,0.,200.);
134  // delta Pt or E quantities for Barrel
135  DBOOK1D(RPt,#DeltaP_{T}/P_{T},80,-1,1);
136  DBOOK1D(RCHE,#DeltaE/E (charged had),80,-2,2);
137  DBOOK1D(RNHE,#DeltaE/E (neutral had),80,-2,2);
138  DBOOK1D(RNEE,#DeltaE/E (neutral em),80,-2,2);
139  DBOOK1D(Rneut,#DeltaE/E (neutral),80,-2,2);
140  DBOOK1D(NCH, #N_{charged},200,0.,200.);
141  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
142  DBOOK2D(RCHEvsPt,#DeltaE/E (charged had) vs P_{T},250, 0, 500, 120,-1,2);
143  DBOOK2D(RNHEvsPt,#DeltaE/E (neutral had) vs P_{T},250, 0, 500, 120,-1,2);
144  DBOOK2D(RNEEvsPt,#DeltaE/E (neutral em) vs P_{T},250, 0, 500, 120,-1,2);
145  DBOOK2D(RneutvsPt,#DeltaE/E (neutral) vs P_{T},250, 0, 500, 120,-1,2);
146  DBOOK2D(NCHvsPt,N_{charged} vs P_{T},250,0,500,200,0.,200.);
147  DBOOK2D(NCH0vsPt, N_{charged} vs P_{T} iter 0,250,0,500,200,0.,200.);
148  DBOOK2D(NCH1vsPt, N_{charged} vs P_{T} iter 1,250,0,500,200,0.,200.);
149  DBOOK2D(NCH2vsPt, N_{charged} vs P_{T} iter 2,250,0,500,200,0.,200.);
150  DBOOK2D(NCH3vsPt, N_{charged} vs P_{T} iter 3,250,0,500,200,0.,200.);
151  DBOOK2D(NCH4vsPt, N_{charged} vs P_{T} iter 4,250,0,500,200,0.,200.);
152  DBOOK2D(NCH5vsPt, N_{charged} vs P_{T} iter 5,250,0,500,200,0.,200.);
153  DBOOK2D(NCH6vsPt, N_{charged} vs P_{T} iter 6,250,0,500,200,0.,200.);
154 
155 
156 
157  DBOOK2D(RNEUTvsP,#DeltaE/E (ECAL+HCAL) vs P,250, 0, 1000, 150,-1.5,1.5);
158  DBOOK2D(RNONLvsP,#DeltaE/E (ECAL+HCAL-only) vs P,250, 0, 1000, 150,-1.5,1.5);
159  DBOOK2D(RHCALvsP,#DeltaE/E (HCAL) vs P,250, 0, 1000, 150,-1.5,1.5);
160  DBOOK2D(RHONLvsP,#DeltaE/E (HCAL only) vs P,250, 0, 1000, 150,-1.5,1.5);
161  DBOOK1D(RPt20_40,#DeltaP_{T}/P_{T},80,-1,1);
162  DBOOK1D(RPt40_60,#DeltaP_{T}/P_{T},80,-1,1);
163  DBOOK1D(RPt60_80,#DeltaP_{T}/P_{T},80,-1,1);
164  DBOOK1D(RPt80_100,#DeltaP_{T}/P_{T},80,-1,1);
165  DBOOK1D(RPt100_150,#DeltaP_{T}/P_{T},80,-1,1);
166  DBOOK1D(RPt150_200,#DeltaP_{T}/P_{T},80,-1,1);
167  DBOOK1D(RPt200_250,#DeltaP_{T}/P_{T},80,-1,1);
168  DBOOK1D(RPt250_300,#DeltaP_{T}/P_{T},80,-1,1);
169  DBOOK1D(RPt300_400,#DeltaP_{T}/P_{T},160,-1,1);
170  DBOOK1D(RPt400_500,#DeltaP_{T}/P_{T},160,-1,1);
171  DBOOK1D(RPt500_750,#DeltaP_{T}/P_{T},160,-1,1);
172  DBOOK1D(RPt750_1250,#DeltaP_{T}/P_{T},160,-1,1);
173  DBOOK1D(RPt1250_2000,#DeltaP_{T}/P_{T},160,-1,1);
174  DBOOK1D(RPt2000_5000,#DeltaP_{T}/P_{T},160,-1,1);
175 
176  DBOOK2D(DEtavsPt,#Delta#eta vs P_{T},1000,0,2000,500,-0.5,0.5);
177  DBOOK2D(DPhivsPt,#Delta#phi vs P_{T},1000,0,2000,500,-0.5,0.5);
178  BOOK2D(DEtavsEta,"#Delta#eta vs P_{T}",1000,-5.,+5.,500,-0.5,0.5);
179  BOOK2D(DPhivsEta,"#Delta#phi vs P_{T}",1000,-5.,+5.,500,-0.5,0.5);
180 
181  // Set Axis Titles
182 
183  // Jets inclusive distributions (Pt > 20 GeV)
184  SETAXES(Njets,"","Multiplicity");
185  SETAXES(jetsPt, PT, "Number of Events");
186  SETAXES(jetsEta, "#eta", "Number of Events");
187  SETAXES(RNeutvsEta, "#eta", "#DeltaE/E (Neutral)");
188  SETAXES(RNEUTvsEta, "#eta", "#DeltaE/E (ECAL+HCAL)");
189  SETAXES(RNONLvsEta, "#eta", "#DeltaE/E (ECAL+HCAL-only)");
190  SETAXES(RHCALvsEta, "#eta", "#DeltaE/E (HCAL)");
191  SETAXES(RHONLvsEta, "#eta", "#DeltaE/E (HCAL Only)");
192  SETAXES(RCHEvsEta, "#eta", "#DeltaE/E (Charged)");
193  SETAXES(RPtvsEta, "#eta", "#DeltaP_{T}/P_{T}");
194  SETAXES(DEtavsEta, "#eta", "#Delta#eta");
195  SETAXES(DPhivsEta,"#eta", "#Delta#phi");
196  // delta Pt or E quantities for Barrel and Endcap
197  DSETAXES(RPt, "#DeltaP_{T}/P_{T}", "Events");
198  DSETAXES(RPt20_40, "#DeltaP_{T}/P_{T}", "Events");
199  DSETAXES(RPt40_60, "#DeltaP_{T}/P_{T}", "Events");
200  DSETAXES(RPt60_80, "#DeltaP_{T}/P_{T}", "Events");
201  DSETAXES(RPt80_100, "#DeltaP_{T}/P_{T}", "Events");
202  DSETAXES(RPt100_150, "#DeltaP_{T}/P_{T}", "Events");
203  DSETAXES(RPt150_200, "#DeltaP_{T}/P_{T}", "Events");
204  DSETAXES(RPt200_250, "#DeltaP_{T}/P_{T}", "Events");
205  DSETAXES(RPt250_300, "#DeltaP_{T}/P_{T}", "Events");
206  DSETAXES(RPt300_400, "#DeltaP_{T}/P_{T}", "Events");
207  DSETAXES(RPt400_500, "#DeltaP_{T}/P_{T}", "Events");
208  DSETAXES(RPt500_750, "#DeltaP_{T}/P_{T}", "Events");
209  DSETAXES(RPt750_1250, "#DeltaP_{T}/P_{T}", "Events");
210  DSETAXES(RPt1250_2000, "#DeltaP_{T}/P_{T}", "Events");
211  DSETAXES(RPt2000_5000, "#DeltaP_{T}/P_{T}", "Events");
212  DSETAXES(RCHE, "#DeltaE/E(charged had)", "Events");
213  DSETAXES(RNHE, "#DeltaE/E(neutral had)", "Events");
214  DSETAXES(RNEE, "#DeltaE/E(neutral em)", "Events");
215  DSETAXES(Rneut, "#DeltaE/E(neutral)", "Events");
216  DSETAXES(RPtvsPt, PT, "#DeltaP_{T}/P_{T}");
217  DSETAXES(RCHEvsPt, PT, "#DeltaE/E(charged had)");
218  DSETAXES(RNHEvsPt, PT, "#DeltaE/E(neutral had)");
219  DSETAXES(RNEEvsPt, PT, "#DeltaE/E(neutral em)");
220  DSETAXES(RneutvsPt, PT, "#DeltaE/E(neutral)");
221  DSETAXES(RHCALvsP, P, "#DeltaE/E(HCAL)");
222  DSETAXES(RHONLvsP, P, "#DeltaE/E(HCAL-only)");
223  DSETAXES(RNEUTvsP, P, "#DeltaE/E(ECAL+HCAL)");
224  DSETAXES(RNONLvsP, P, "#DeltaE/E(ECAL+HCAL-only)");
225  DSETAXES(DEtavsPt, PT, "#Delta#eta");
226  DSETAXES(DPhivsPt, PT, "#Delta#phi");
227 
228 }
229 
230 
232  // loop over reco pf jets
233  resPtMax_ = 0.;
237  int NPFJets = 0;
238 
239  for(unsigned i=0; i<pfJets.size(); i++) {
240 
241  // Count the number of jets with a larger energy
242  unsigned highJets = 0;
243  for(unsigned j=0; j<pfJets.size(); j++) {
244  if ( j != i && pfJets[j].pt() > pfJets[i].pt() ) highJets++;
245  }
246  if ( onlyTwoJets_ && highJets > 1 ) continue;
247 
248 
249  const reco::PFJet& pfj = pfJets[i];
250  double rec_pt = pfj.pt();
251  double rec_eta = pfj.eta();
252  double rec_phi = pfj.phi();
253 
254  // skip PFjets with pt < recPt_cut GeV
255  if (rec_pt<recPt_cut and recPt_cut != -1.) continue;
256  // skip PFjets with eta > maxEta_cut
257  if (fabs(rec_eta)>maxEta_cut and maxEta_cut != -1.) continue;
258 
259  NPFJets++;
260 
261  // fill inclusive PFjet distribution pt > 20 GeV
262  hNjets->Fill(NPFJets);
263  hjetsPt->Fill(rec_pt);
264  hjetsEta->Fill(rec_eta);
265 
266  // separate Barrel PFJets from Endcap PFJets
267  bool Barrel = false;
268  bool Endcap = false;
269  bool Forward = false;
270  if (std::abs(rec_eta) < 1.4 ) Barrel = true;
271  if (std::abs (rec_eta) > 1.6 && std::abs (rec_eta) < 2.4 ) Endcap = true;
272  if (std::abs (rec_eta) > 2.5 && std::abs (rec_eta) < 2.9 ) Forward = true;
273  if (std::abs (rec_eta) > 3.1 && std::abs (rec_eta) < 4.7 ) Forward = true;
274 
275  // do only barrel for now
276  // if(!Barrel) continue;
277 
278  // look for the closets gen Jet : truth
279  const GenJet *truth = algo_->matchByDeltaR(&pfj,&genJets);
280  if(!truth) continue;
281  double deltaR = algo_->deltaR(&pfj, truth);
282  // check deltaR is small enough
283  if(deltaR < deltaRMax_ || (abs(rec_eta)>2.5 && deltaR < 0.2) || deltaRMax_ == -1.0 ) {//start case deltaR < deltaRMax
284 
285  // generate histograms comparing the reco and truth candidate (truth = closest in delta-R)
286  // get the quantities to place on the denominator and/or divide by
287  double pt_denom;
288  double true_E = truth->p();
289  double true_pt = truth->pt();
290  double true_eta = truth->eta();
291  double true_phi = truth->phi();
292 
293  if (plotAgainstReco_) {pt_denom = rec_pt;}
294  else {pt_denom = true_pt;}
295  // get true specific quantities
296  double true_ChargedHadEnergy;
297  double true_NeutralHadEnergy;
298  double true_NeutralEmEnergy;
299  gettrue (truth, true_ChargedHadEnergy, true_NeutralHadEnergy, true_NeutralEmEnergy);
300  double true_NeutralEnergy = true_NeutralHadEnergy + true_NeutralEmEnergy;
301  double rec_ChargedHadEnergy = pfj.chargedHadronEnergy();
302  double rec_NeutralHadEnergy = pfj.neutralHadronEnergy();
303  double rec_NeutralEmEnergy = pfj.neutralEmEnergy();
304  double rec_NeutralEnergy = rec_NeutralHadEnergy + rec_NeutralEmEnergy;
305  double rec_ChargedMultiplicity = pfj.chargedMultiplicity();
306  std::vector <PFCandidatePtr> constituents = pfj.getPFConstituents ();
307  std::vector <unsigned int> chMult(9, static_cast<unsigned int>(0));
308  for (unsigned ic = 0; ic < constituents.size (); ++ic) {
309  if ( constituents[ic]->particleId() > 3 ) continue;
310  reco::TrackRef trackRef = constituents[ic]->trackRef();
311  if ( trackRef.isNull() ) {
312  //std::cout << "Warning in entry " << entry_
313  // << " : a track with Id " << constituents[ic]->particleId()
314  // << " has no track ref.." << std::endl;
315  continue;
316  }
317  unsigned int iter = trackRef->algo()>6 ? 6: trackRef->algo();
318  ++(chMult[iter]);
319  }
320 
321  bool plot1 = false;
322  bool plot2 = false;
323  bool plot3 = false;
324  bool plot4 = false;
325  bool plot5 = false;
326  bool plot6 = false;
327  bool plot7 = false;
328  double cut1 = 0.0001;
329  double cut2 = 0.0001;
330  double cut3 = 0.0001;
331  double cut4 = 0.0001;
332  double cut5 = 0.0001;
333  double cut6 = 0.0001;
334  double cut7 = 0.0001;
335  double resPt =0.;
336  double resChargedHadEnergy= 0.;
337  double resNeutralHadEnergy= 0.;
338  double resNeutralEmEnergy= 0.;
339  double resNeutralEnergy= 0.;
340 
341  double resHCALEnergy = 0.;
342  double resNEUTEnergy = 0.;
343  if ( rec_NeutralHadEnergy > cut6 && rec_ChargedHadEnergy < cut1 ) {
344  double true_NEUTEnergy = true_NeutralHadEnergy + true_NeutralEmEnergy;
345  double true_HCALEnergy = true_NEUTEnergy - rec_NeutralEmEnergy;
346  double rec_NEUTEnergy = rec_NeutralHadEnergy+rec_NeutralEmEnergy;
347  double rec_HCALEnergy = rec_NeutralHadEnergy;
348  resHCALEnergy = (rec_HCALEnergy-true_HCALEnergy)/rec_HCALEnergy;
349  resNEUTEnergy = (rec_NEUTEnergy-true_NEUTEnergy)/rec_NEUTEnergy;
350  if ( rec_NeutralEmEnergy > cut7 ) {
351  plot6 = true;
352  } else {
353  plot7 = true;
354  }
355  }
356 
357  // get relative delta quantities (protect against division by zero!)
358  if (true_pt > 0.0001){
359  resPt = (rec_pt -true_pt)/true_pt ;
360  plot1 = true;}
361  if (true_ChargedHadEnergy > cut1){
362  resChargedHadEnergy = (rec_ChargedHadEnergy- true_ChargedHadEnergy)/true_ChargedHadEnergy;
363  plot2 = true;}
364  if (true_NeutralHadEnergy > cut2){
365  resNeutralHadEnergy = (rec_NeutralHadEnergy- true_NeutralHadEnergy)/true_NeutralHadEnergy;
366  plot3 = true;}
367  else
368  if (rec_NeutralHadEnergy > cut3){
369  resNeutralHadEnergy = (rec_NeutralHadEnergy- true_NeutralHadEnergy)/rec_NeutralHadEnergy;
370  plot3 = true;}
371  if (true_NeutralEmEnergy > cut4){
372  resNeutralEmEnergy = (rec_NeutralEmEnergy- true_NeutralEmEnergy)/true_NeutralEmEnergy;
373  plot4 = true;}
374  if (true_NeutralEnergy > cut5){
375  resNeutralEnergy = (rec_NeutralEnergy- true_NeutralEnergy)/true_NeutralEnergy;
376  plot5 = true;}
377 
378  //double deltaEta = algo_->deltaEta(&pfj, truth);
379  //double deltaPhi = algo_->deltaPhi(&pfj, truth);
380 
381  // Print outliers for further debugging
382  if ( ( resPt > 0.2 && true_pt > 100. ) ||
383  ( resPt < -0.5 && true_pt > 100. ) ) {
384  //if ( ( true_pt > 50. &&
385  // ( ( truth->eta()>3.0 && rec_eta-truth->eta() < -0.1 ) ||
386  // ( truth->eta()<-3.0 && rec_eta-truth->eta() > 0.1 ) ))) {
387  std::cout << "Entry " << entry_
388  << " resPt = " << resPt
389  <<" resCharged " << resChargedHadEnergy
390  <<" resNeutralHad " << resNeutralHadEnergy
391  << " resNeutralEm " << resNeutralEmEnergy
392  << " pT (T/R) " << true_pt << "/" << rec_pt
393  << " Eta (T/R) " << truth->eta() << "/" << rec_eta
394  << " Phi (T/R) " << truth->phi() << "/" << rec_phi
395  << std::endl;
396 
397  // check overlapping PF jets
398  const reco::PFJet* pfoj = nullptr;
399  double dRo = 1E9;
400  for(unsigned j=0; j<pfJets.size(); j++) {
401  const reco::PFJet& pfo = pfJets[j];
402  if ( j != i && algo_->deltaR(&pfj,&pfo) < dRo && pfo.pt() > 0.25*pfj.pt()) {
403  dRo = algo_->deltaR(&pfj,&pfo);
404  pfoj = &pfo;
405  }
406  }
407 
408  // Check overlapping Gen Jet
409  math::XYZTLorentzVector overlappinGenJet(0.,0.,0.,0.);
410  const reco::GenJet* genoj = nullptr;
411  double dRgo = 1E9;
412  for(unsigned j=0; j<genJets.size(); j++) {
413  const reco::GenJet* gjo = &(genJets[j]);
414  if ( gjo != truth && algo_->deltaR(truth,gjo) < dRgo && gjo->pt() > 0.25*truth->pt() ) {
415  dRgo = algo_->deltaR(truth,gjo);
416  genoj = gjo;
417  }
418  }
419 
420  if ( dRo < 0.8 && dRgo < 0.8 && algo_->deltaR(genoj,pfoj) < 2.*deltaRMax_ )
421  std::cout << "Excess probably due to overlapping jets (DR = " << algo_->deltaR(genoj,pfoj) << "),"
422  << " at DeltaR(T/R) = " << dRgo << "/" << dRo
423  << " with pT(T/R) " << genoj->pt() << "/" << pfoj->pt()
424  << " and Eta (T/R) " << genoj->eta() << "/" << pfoj->eta()
425  << " and Phi (T/R) " << genoj->phi() << "/" << pfoj->phi()
426  << std::endl;
427  }
428 
429  if(std::abs(resPt) > std::abs(resPtMax_)) resPtMax_ = resPt;
430  if(std::abs(resChargedHadEnergy) > std::abs(resChargedHadEnergyMax_) ) resChargedHadEnergyMax_ = resChargedHadEnergy;
431  if(std::abs(resNeutralHadEnergy) > std::abs(resNeutralHadEnergyMax_) ) resNeutralHadEnergyMax_ = resNeutralHadEnergy;
432  if(std::abs(resNeutralEmEnergy) > std::abs(resNeutralEmEnergyMax_) ) resNeutralEmEnergyMax_ = resNeutralEmEnergy;
433  if (debug_) {
434  cout << i <<" =========PFJet Pt "<< rec_pt
435  << " eta " << rec_eta
436  << " phi " << rec_phi
437  << " Charged Had Energy " << rec_ChargedHadEnergy
438  << " Neutral Had Energy " << rec_NeutralHadEnergy
439  << " Neutral elm Energy " << rec_NeutralEmEnergy << endl;
440  cout << " matching Gen Jet Pt " << true_pt
441  << " eta " << truth->eta()
442  << " phi " << truth->phi()
443  << " Charged Had Energy " << true_ChargedHadEnergy
444  << " Neutral Had Energy " << true_NeutralHadEnergy
445  << " Neutral elm Energy " << true_NeutralEmEnergy << endl;
446  printPFJet(&pfj);
447  // cout<<pfj.print()<<endl;
448  printGenJet(truth);
449  //cout <<truth->print()<<endl;
450 
451  cout << "==============deltaR " << deltaR << " resPt " << resPt
452  << " resChargedHadEnergy " << resChargedHadEnergy
453  << " resNeutralHadEnergy " << resNeutralHadEnergy
454  << " resNeutralEmEnergy " << resNeutralEmEnergy
455  << endl;
456  }
457 
458 
459  if(plot1) {
460  if ( rec_eta > 0. )
461  hDEtavsEta->Fill(true_eta,rec_eta-true_eta);
462  else
463  hDEtavsEta->Fill(true_eta,-rec_eta+true_eta);
464  hDPhivsEta->Fill(true_eta,rec_phi-true_phi);
465 
466  hRPtvsEta->Fill(true_eta, resPt);
467  hNCHvsEta->Fill(true_eta, rec_ChargedMultiplicity);
468  hNCH0vsEta->Fill(true_eta,chMult[0]);
469  hNCH1vsEta->Fill(true_eta,chMult[1]);
470  hNCH2vsEta->Fill(true_eta,chMult[2]);
471  hNCH3vsEta->Fill(true_eta,chMult[3]);
472  hNCH4vsEta->Fill(true_eta,chMult[4]);
473  hNCH5vsEta->Fill(true_eta,chMult[5]);
474  hNCH6vsEta->Fill(true_eta,chMult[6]);
475  hNCH7vsEta->Fill(true_eta,chMult[7]);
476  }
477  if(plot2)hRCHEvsEta->Fill(true_eta, resChargedHadEnergy);
478  if(plot5)hRNeutvsEta->Fill(true_eta, resNeutralEnergy);
479  if(plot6) {
480  hRHCALvsEta->Fill(true_eta, resHCALEnergy);
481  hRNEUTvsEta->Fill(true_eta, resNEUTEnergy);
482  }
483  if(plot7) {
484  hRHONLvsEta->Fill(true_eta, resHCALEnergy);
485  hRNONLvsEta->Fill(true_eta, resNEUTEnergy);
486  }
487 
488  // fill histograms for relative delta quantitites of matched jets
489  // delta Pt or E quantities for Barrel
490  if (Barrel){
491  if(plot1) {
492  hBRPt->Fill (resPt);
493  if ( pt_denom > 20. && pt_denom < 40. ) hBRPt20_40->Fill (resPt);
494  if ( pt_denom > 40. && pt_denom < 60. ) hBRPt40_60->Fill (resPt);
495  if ( pt_denom > 60. && pt_denom < 80. ) hBRPt60_80->Fill (resPt);
496  if ( pt_denom > 80. && pt_denom < 100. ) hBRPt80_100->Fill (resPt);
497  if ( pt_denom > 100. && pt_denom < 150. ) hBRPt100_150->Fill (resPt);
498  if ( pt_denom > 150. && pt_denom < 200. ) hBRPt150_200->Fill (resPt);
499  if ( pt_denom > 200. && pt_denom < 250. ) hBRPt200_250->Fill (resPt);
500  if ( pt_denom > 250. && pt_denom < 300. ) hBRPt250_300->Fill (resPt);
501  if ( pt_denom > 300. && pt_denom < 400. ) hBRPt300_400->Fill (resPt);
502  if ( pt_denom > 400. && pt_denom < 500. ) hBRPt400_500->Fill (resPt);
503  if ( pt_denom > 500. && pt_denom < 750. ) hBRPt500_750->Fill (resPt);
504  if ( pt_denom > 750. && pt_denom < 1250. ) hBRPt750_1250->Fill (resPt);
505  if ( pt_denom > 1250. && pt_denom < 2000. ) hBRPt1250_2000->Fill (resPt);
506  if ( pt_denom > 2000. && pt_denom < 5000. ) hBRPt2000_5000->Fill (resPt);
507  hBNCH->Fill(rec_ChargedMultiplicity);
508  hBNCH0vsPt->Fill(pt_denom,chMult[0]);
509  hBNCH1vsPt->Fill(pt_denom,chMult[1]);
510  hBNCH2vsPt->Fill(pt_denom,chMult[2]);
511  hBNCH3vsPt->Fill(pt_denom,chMult[3]);
512  hBNCH4vsPt->Fill(pt_denom,chMult[4]);
513  hBNCH5vsPt->Fill(pt_denom,chMult[5]);
514  hBNCH6vsPt->Fill(pt_denom,chMult[6]);
515  hBNCH7vsPt->Fill(pt_denom,chMult[7]);
516  hBNCHvsPt->Fill(pt_denom,rec_ChargedMultiplicity);
517  if ( rec_eta > 0. )
518  hBDEtavsPt->Fill(pt_denom,rec_eta-true_eta);
519  else
520  hBDEtavsPt->Fill(pt_denom,-rec_eta+true_eta);
521  hBDPhivsPt->Fill(pt_denom,rec_phi-true_phi);
522  }
523  if(plot2)hBRCHE->Fill(resChargedHadEnergy);
524  if(plot3)hBRNHE->Fill(resNeutralHadEnergy);
525  if(plot4)hBRNEE->Fill(resNeutralEmEnergy);
526  if(plot5)hBRneut->Fill(resNeutralEnergy);
527  if(plot1)hBRPtvsPt->Fill(pt_denom, resPt);
528  if(plot2)hBRCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
529  if(plot3)hBRNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
530  if(plot4)hBRNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
531  if(plot5)hBRneutvsPt->Fill(pt_denom, resNeutralEnergy);
532  if(plot6) {
533  hBRHCALvsP->Fill(true_E, resHCALEnergy);
534  hBRNEUTvsP->Fill(true_E, resNEUTEnergy);
535  }
536  if(plot7) {
537  hBRHONLvsP->Fill(true_E, resHCALEnergy);
538  hBRNONLvsP->Fill(true_E, resNEUTEnergy);
539  }
540 
541  }
542  // delta Pt or E quantities for Endcap
543  if (Endcap){
544  if(plot1) {
545  hERPt->Fill (resPt);
546  if ( pt_denom > 20. && pt_denom < 40. ) hERPt20_40->Fill (resPt);
547  if ( pt_denom > 40. && pt_denom < 60. ) hERPt40_60->Fill (resPt);
548  if ( pt_denom > 60. && pt_denom < 80. ) hERPt60_80->Fill (resPt);
549  if ( pt_denom > 80. && pt_denom < 100. ) hERPt80_100->Fill (resPt);
550  if ( pt_denom > 100. && pt_denom < 150. ) hERPt100_150->Fill (resPt);
551  if ( pt_denom > 150. && pt_denom < 200. ) hERPt150_200->Fill (resPt);
552  if ( pt_denom > 200. && pt_denom < 250. ) hERPt200_250->Fill (resPt);
553  if ( pt_denom > 250. && pt_denom < 300. ) hERPt250_300->Fill (resPt);
554  if ( pt_denom > 300. && pt_denom < 400. ) hERPt300_400->Fill (resPt);
555  if ( pt_denom > 400. && pt_denom < 500. ) hERPt400_500->Fill (resPt);
556  if ( pt_denom > 500. && pt_denom < 750. ) hERPt500_750->Fill (resPt);
557  if ( pt_denom > 750. && pt_denom < 1250. ) hERPt750_1250->Fill (resPt);
558  if ( pt_denom > 1250. && pt_denom < 2000. ) hERPt1250_2000->Fill (resPt);
559  if ( pt_denom > 2000. && pt_denom < 5000. ) hERPt2000_5000->Fill (resPt);
560  hENCH->Fill(rec_ChargedMultiplicity);
561  hENCHvsPt->Fill(pt_denom,rec_ChargedMultiplicity);
562  hENCH0vsPt->Fill(pt_denom,chMult[0]);
563  hENCH1vsPt->Fill(pt_denom,chMult[1]);
564  hENCH2vsPt->Fill(pt_denom,chMult[2]);
565  hENCH3vsPt->Fill(pt_denom,chMult[3]);
566  hENCH4vsPt->Fill(pt_denom,chMult[4]);
567  hENCH5vsPt->Fill(pt_denom,chMult[5]);
568  hENCH6vsPt->Fill(pt_denom,chMult[6]);
569  hENCH7vsPt->Fill(pt_denom,chMult[7]);
570  if ( rec_eta > 0. )
571  hEDEtavsPt->Fill(pt_denom,rec_eta-true_eta);
572  else
573  hEDEtavsPt->Fill(pt_denom,-rec_eta+true_eta);
574  hEDPhivsPt->Fill(pt_denom,rec_phi-true_phi);
575  }
576  if(plot2)hERCHE->Fill(resChargedHadEnergy);
577  if(plot3)hERNHE->Fill(resNeutralHadEnergy);
578  if(plot4)hERNEE->Fill(resNeutralEmEnergy);
579  if(plot5)hERneut->Fill(resNeutralEnergy);
580  if(plot1)hERPtvsPt->Fill(pt_denom, resPt);
581  if(plot2)hERCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
582  if(plot3)hERNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
583  if(plot4)hERNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
584  if(plot5)hERneutvsPt->Fill(pt_denom, resNeutralEnergy);
585  if(plot6) {
586  hERHCALvsP->Fill(true_E, resHCALEnergy);
587  hERNEUTvsP->Fill(true_E, resNEUTEnergy);
588  }
589  if(plot7) {
590  hERHONLvsP->Fill(true_E, resHCALEnergy);
591  hERNONLvsP->Fill(true_E, resNEUTEnergy);
592  }
593  }
594  // delta Pt or E quantities for Forward
595  if (Forward){
596  if(plot1) {
597  hFRPt->Fill (resPt);
598  if ( pt_denom > 20. && pt_denom < 40. ) hFRPt20_40->Fill (resPt);
599  if ( pt_denom > 40. && pt_denom < 60. ) hFRPt40_60->Fill (resPt);
600  if ( pt_denom > 60. && pt_denom < 80. ) hFRPt60_80->Fill (resPt);
601  if ( pt_denom > 80. && pt_denom < 100. ) hFRPt80_100->Fill (resPt);
602  if ( pt_denom > 100. && pt_denom < 150. ) hFRPt100_150->Fill (resPt);
603  if ( pt_denom > 150. && pt_denom < 200. ) hFRPt150_200->Fill (resPt);
604  if ( pt_denom > 200. && pt_denom < 250. ) hFRPt200_250->Fill (resPt);
605  if ( pt_denom > 250. && pt_denom < 300. ) hFRPt250_300->Fill (resPt);
606  if ( pt_denom > 300. && pt_denom < 400. ) hFRPt300_400->Fill (resPt);
607  if ( pt_denom > 400. && pt_denom < 500. ) hFRPt400_500->Fill (resPt);
608  if ( pt_denom > 500. && pt_denom < 750. ) hFRPt500_750->Fill (resPt);
609  if ( pt_denom > 750. && pt_denom < 1250. ) hFRPt750_1250->Fill (resPt);
610  if ( pt_denom > 1250. && pt_denom < 2000. ) hFRPt1250_2000->Fill (resPt);
611  if ( pt_denom > 2000. && pt_denom < 5000. ) hFRPt2000_5000->Fill (resPt);
612  if ( rec_eta > 0. )
613  hFDEtavsPt->Fill(pt_denom,rec_eta-true_eta);
614  else
615  hFDEtavsPt->Fill(pt_denom,-rec_eta+true_eta);
616  hFDPhivsPt->Fill(pt_denom,rec_phi-true_phi);
617  }
618  if(plot2)hFRCHE->Fill(resChargedHadEnergy);
619  if(plot3)hFRNHE->Fill(resNeutralHadEnergy);
620  if(plot4)hFRNEE->Fill(resNeutralEmEnergy);
621  if(plot5)hFRneut->Fill(resNeutralEnergy);
622  if(plot1)hFRPtvsPt->Fill(pt_denom, resPt);
623  if(plot2)hFRCHEvsPt->Fill(pt_denom, resChargedHadEnergy);
624  if(plot3)hFRNHEvsPt->Fill(pt_denom, resNeutralHadEnergy);
625  if(plot4)hFRNEEvsPt->Fill(pt_denom, resNeutralEmEnergy);
626  if(plot5)hFRneutvsPt->Fill(pt_denom, resNeutralEnergy);
627  if(plot6) {
628  hFRHCALvsP->Fill(true_E, resHCALEnergy);
629  hFRNEUTvsP->Fill(true_E, resNEUTEnergy);
630  }
631  if(plot7) {
632  hFRHONLvsP->Fill(true_E, resHCALEnergy);
633  hFRNONLvsP->Fill(true_E, resNEUTEnergy);
634  }
635  }
636  } // end case deltaR < deltaRMax
637 
638  } // i loop on pf Jets
639 
640  // Increment counter
641  entry_++;
642 }
643 
644 void PFJetBenchmark::gettrue (const reco::GenJet* truth, double& true_ChargedHadEnergy,
645  double& true_NeutralHadEnergy, double& true_NeutralEmEnergy){
646  std::vector <const GenParticle*> mcparts = truth->getGenConstituents ();
647  true_NeutralEmEnergy = 0.;
648  true_ChargedHadEnergy = 0.;
649  true_NeutralHadEnergy = 0.;
650  // for each MC particle in turn
651  for (unsigned i = 0; i < mcparts.size (); i++) {
652  const GenParticle* mcpart = mcparts[i];
653  int PDG = std::abs( mcpart->pdgId());
654  double e = mcpart->energy();
655  switch(PDG){ // start PDG switch
656  case 22: // photon
657  true_NeutralEmEnergy += e;
658  break;
659  case 211: // pi
660  case 321: // K
661  case 2212: // p
662  case 11: //electrons (until recognised)
663  true_ChargedHadEnergy += e;
664  break;
665  case 310: // K_S0
666  case 130: // K_L0
667  case 3122: // Lambda0
668  case 2112: // n0
669  true_NeutralHadEnergy += e;
670  default:
671  break;
672  } // end PDG switch
673  } // end loop on constituents.
674 }
675 
677  cout<<setiosflags(ios::right);
678  cout<<setiosflags(ios::fixed);
679  cout<<setprecision(3);
680 
681  cout << "PFJet p/px/py/pz/pt: " << pfj->p() << "/" << pfj->px ()
682  << "/" << pfj->py() << "/" << pfj->pz() << "/" << pfj->pt() << endl
683  << " eta/phi: " << pfj->eta () << "/" << pfj->phi () << endl
684  << " PFJet specific:" << std::endl
685  << " charged/neutral hadrons energy: " << pfj->chargedHadronEnergy () << '/' << pfj->neutralHadronEnergy () << endl
686  << " charged/neutral em energy: " << pfj->chargedEmEnergy () << '/' << pfj->neutralEmEnergy () << endl
687  << " charged muon energy: " << pfj->chargedMuEnergy () << '/' << endl
688  << " charged/neutral multiplicity: " << pfj->chargedMultiplicity () << '/' << pfj->neutralMultiplicity () << endl;
689 
690  // And print the constituents
691  std::cout << pfj->print() << std::endl;
692 
693  cout<<resetiosflags(ios::right|ios::fixed);
694 }
695 
696 
698  std::vector <const GenParticle*> mcparts = truth->getGenConstituents ();
699  cout << "GenJet p/px/py/pz/pt: " << truth->p() << '/' << truth->px ()
700  << '/' << truth->py() << '/' << truth->pz() << '/' << truth->pt() << endl
701  << " eta/phi: " << truth->eta () << '/' << truth->phi () << endl
702  << " # of constituents: " << mcparts.size() << endl;
703  cout << " constituents:" << endl;
704  for (unsigned i = 0; i < mcparts.size (); i++) {
705  const GenParticle* mcpart = mcparts[i];
706  cout << " #" << i << " PDG code:" << mcpart->pdgId()
707  << ", p/pt/eta/phi: " << mcpart->p() << '/' << mcpart->pt()
708  << '/' << mcpart->eta() << '/' << mcpart->phi() << endl;
709  }
710 }
711 
int pdgId() const final
PDG identifier.
#define BOOK1D(name, title, nbinsx, lowx, highx)
double eta() const final
momentum pseudorapidity
float chargedEmEnergy() const
chargedEmEnergy
Definition: PFJet.h:142
bool onlyTwoJets
#define PT
double deltaRMax
double px() const final
x coordinate of momentum vector
std::vector< GenJet > GenJetCollection
collection of GenJet objects
#define SETAXES(name, xtitle, ytitle)
double pt() const final
transverse momentum
double maxEta
#define nullptr
int chargedMultiplicity() const
chargedMultiplicity
Definition: PFJet.h:155
Jets made from PFObjects.
Definition: PFJet.h:21
float neutralEmEnergy() const
neutralEmEnergy
Definition: PFJet.h:150
unsigned int entry_
DQMStore * dbe_
virtual ~PFJetBenchmark()
void gettrue(const reco::GenJet *truth, double &true_ChargedHadEnergy, double &true_NeutralHadEnergy, double &true_NeutralEmEnergy)
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
#define BOOK2D(name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)
double pz() const final
z coordinate of momentum vector
static double deltaR(const T *, const U *)
double energy() const final
energy
#define DBOOK2D(name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
Jets made from MC generator particles.
Definition: GenJet.h:24
std::string outputFile_
void setup(std::string Filename, bool debug, bool plotAgainstReco=false, bool onlyTwoJets=true, double deltaRMax=0.1, std::string benchmarkLabel_="ParticleFlow", double recPt=-1, double maxEta=-1, DQMStore *dbe_store=0)
string benchmarkLabel_
PFBenchmarkAlgo * algo_
virtual std::vector< const GenParticle * > getGenConstituents() const
get all constituents
Definition: GenJet.cc:58
bool isNull() const
Checks for null.
Definition: Ref.h:250
int neutralMultiplicity() const
neutralMultiplicity
Definition: PFJet.h:157
void process(const reco::PFJetCollection &, const reco::GenJetCollection &)
std::string print() const override
Print object in details.
Definition: PFJet.cc:83
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
Definition: HCAL.py:1
bool plotAgainstReco
double p() const final
magnitude of momentum vector
double resChargedHadEnergyMax_
static const Collection::value_type * matchByDeltaR(const T *, const Collection *)
#define debug
Definition: HDRShower.cc:19
double py() const final
y coordinate of momentum vector
std::pair< OmniClusterRef, TrackingParticleRef > P
std::vector< PFJet > PFJetCollection
collection of PFJet objects
double resNeutralHadEnergyMax_
fixed size matrix
void printPFJet(const reco::PFJet *)
virtual std::vector< reco::PFCandidatePtr > getPFConstituents() const
get all constituents
Definition: PFJet.cc:52
float neutralHadronEnergy() const
neutralHadronEnergy
Definition: PFJet.h:102
float chargedMuEnergy() const
chargedMuEnergy
Definition: PFJet.h:146
double recPt
double resNeutralEmEnergyMax_
#define DSETAXES(name, xtitle, ytitle)
long double T
double phi() const final
momentum azimuthal angle
void printGenJet(const reco::GenJet *)
#define DBOOK1D(name, title, nbinsx, lowx, highx)
float chargedHadronEnergy() const
chargedHadronEnergy
Definition: PFJet.h:98