CMS 3D CMS Logo

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