CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PFMETBenchmark.cc
Go to the documentation of this file.
1 #include "RecoParticleFlow/Benchmark/interface/PFMETBenchmark.h"
2 //#include "TCanvas.h"
3 #include "TProfile.h"
4 #include "TF1.h"
5 #include "TH1F.h"
6 
7 // preprocessor macro for booking 1d histos with DQMStore -or- bare Root
8 #define BOOK1D(name,title,nbinsx,lowx,highx) \
9  h##name = dbe_ ? dbe_->book1D(#name,title,nbinsx,lowx,highx)->getTH1F() \
10  : new TH1F(#name,title,nbinsx,lowx,highx)
11 
12 // preprocessor macro for booking 2d histos with DQMStore -or- bare Root
13 #define BOOK2D(name,title,nbinsx,lowx,highx,nbinsy,lowy,highy) \
14  h##name = dbe_ ? dbe_->book2D(#name,title,nbinsx,lowx,highx,nbinsy,lowy,highy)->getTH2F() \
15  : new TH2F(#name,title,nbinsx,lowx,highx,nbinsy,lowy,highy)
16 
17 // all versions OK
18 // preprocesor macro for setting axis titles
19 #define SETAXES(name,xtitle,ytitle) \
20  h##name->GetXaxis()->SetTitle(xtitle); h##name->GetYaxis()->SetTitle(ytitle)
21 
22 
23 /*#define SET2AXES(name,xtitle,ytitle) \
24  hE##name->GetXaxis()->SetTitle(xtitle); hE##name->GetYaxis()->SetTitle(ytitle); hB##name->GetXaxis()->SetTitle(xtitle); hB##name->GetYaxis()->SetTitle(ytitle)
25 */
26 
27 #define PT (plotAgainstReco_)?"reconstructed P_{T}" :"generated P_{T}"
28 
29 using namespace reco;
30 using namespace std;
31 
32 class MonitorElement;
33 
34 PFMETBenchmark::PFMETBenchmark() : file_(0) {}
35 
36 PFMETBenchmark::~PFMETBenchmark() {
37  if(file_) file_->Close();
38 }
39 
40 void PFMETBenchmark::write() {
41  // Store the DAQ Histograms
42  if (outputFile_.size() != 0) {
43  if (dbe_)
44  dbe_->save(outputFile_.c_str());
45  // use bare Root if no DQM (FWLite applications)
46  else if (file_) {
47  file_->Write(outputFile_.c_str());
48  cout << "Benchmark output written to file " << outputFile_.c_str() << endl;
49  file_->Close();
50  }
51  }
52  else
53  cout << "No output file specified ("<<outputFile_<<"). Results will not be saved!" << endl;
54 
55 }
56 
57 void PFMETBenchmark::setup(
58  string Filename,
59  bool debug,
60  bool plotAgainstReco,
61  string benchmarkLabel_,
62  DQMStore * dbe_store) {
63  debug_ = debug;
64  plotAgainstReco_ = plotAgainstReco;
65  outputFile_=Filename;
66  file_ = NULL;
67  dbe_ = dbe_store;
68  // print parameters
69  //cout<< "PFMETBenchmark Setup parameters =============================================="<<endl;
70  cout << "Filename to write histograms " << Filename<<endl;
71  cout << "PFMETBenchmark debug " << debug_<< endl;
72  cout << "plotAgainstReco " << plotAgainstReco_ << endl;
73  cout << "benchmarkLabel " << benchmarkLabel_ << endl;
74 
75  // Book histogram
76 
77  // Establish DQM Store
78  string path = "PFTask/Benchmarks/"+ benchmarkLabel_ + "/";
79  if (plotAgainstReco) path += "Reco"; else path += "Gen";
80  if (dbe_) {
81  dbe_->setCurrentFolder(path.c_str());
82  }
83  else {
84  file_ = new TFile(outputFile_.c_str(), "recreate");
85 // TTree * tr = new TTree("PFTast");
86 // tr->Branch("Benchmarks/ParticleFlow")
87  cout << "Info: DQM is not available to provide data storage service. Using TFile to save histograms. "<<endl;
88  }
89 
90  // delta Pt or E quantities for Barrel
91  BOOK1D(MEX,"Particle Flow",400,-200,200);
92  BOOK1D(DeltaMEX,"Particle Flow",400,-200,200);
93  BOOK1D(DeltaMET,"Particle Flow",400,-200,200);
94  BOOK1D(DeltaPhi,"Particle Flow", 1000, -3.2, 3.2);
95  BOOK1D(DeltaSET,"Particle Flow",400,-200,200);
96  BOOK2D(SETvsDeltaMET,"Particle Flow",200, 0.0, 1000.0, 400, -200.0, 200.0);
97  BOOK2D(SETvsDeltaSET,"Particle Flow",200, 0.0, 1000.0, 400, -200.0, 200.0);
98  profileSETvsSETresp = new TProfile("#DeltaPFSET / trueSET vs trueSET", "", 200, 0.0, 1000.0, -1.0, 1.0);
99  profileMETvsMETresp = new TProfile("#DeltaPFMET / trueMET vs trueMET", "", 50, 0.0, 200.0, -1.0, 1.0);
100 
101  BOOK1D(CaloMEX,"Calorimeter",400,-200,200);
102  BOOK1D(DeltaCaloMEX,"Calorimeter",400,-200,200);
103  BOOK1D(DeltaCaloMET,"Calorimeter",400,-200,200);
104  BOOK1D(DeltaCaloPhi,"Calorimeter", 1000, -3.2, 3.2);
105  BOOK1D(DeltaCaloSET,"Calorimeter",400,-200,200);
106  BOOK2D(CaloSETvsDeltaCaloMET,"Calorimeter",200, 0.0, 1000.0, 400, -200.0, 200.0);
107  BOOK2D(CaloSETvsDeltaCaloSET,"Calorimeter",200, 0.0, 1000.0, 400, -200.0, 200.0);
108  profileCaloSETvsCaloSETresp = new TProfile("#DeltaCaloSET / trueSET vs trueSET", "", 200, 0.0, 1000.0, -1.0, 1.0);
109  profileCaloMETvsCaloMETresp = new TProfile("#DeltaCaloMET / trueMET vs trueMET", "", 200, 0.0, 200.0, -1.0, 1.0);
110 
111  BOOK2D(DeltaPFMETvstrueMET,"Particle Flow", 500, 0.0, 1000.0, 400, -200.0, 200.0);
112  BOOK2D(DeltaCaloMETvstrueMET,"Calorimeter", 500, 0.0, 1000.0, 400, -200.0, 200.0);
113  BOOK2D(DeltaPFPhivstrueMET,"Particle Flow", 500, 0.0, 1000.0, 400, -3.2, 3.2);
114  BOOK2D(DeltaCaloPhivstrueMET,"Calorimeter", 500, 0.0, 1000.0, 400, -3.2, 3.2);
115  BOOK2D(CaloMETvstrueMET,"Calorimeter", 500, 0.0, 1000.0, 500, 0.0, 1000.0);
116  BOOK2D(PFMETvstrueMET,"Particle Flow", 500, 0.0, 1000.0, 500, 0.0, 1000.0);
117 
118  BOOK2D(DeltaCaloMEXvstrueSET,"Calorimeter",200, 0.0, 1000.0, 400, -200.0, 200.0);
119  BOOK2D(DeltaPFMEXvstrueSET,"Particle Flow",200, 0.0, 1000.0, 400, -200.0, 200.0);
120 
121  BOOK1D(TrueMET,"MC truth", 500, 0.0, 1000.0);
122  BOOK1D(CaloMET,"Calorimeter", 500, 0.0, 1000.0);
123  BOOK1D(PFMET,"Particle Flow", 500, 0.0, 1000.0);
124 
125  BOOK1D(TCMEX,"Track Corrected",400,-200,200);
126  BOOK1D(DeltaTCMEX,"Track Corrected",400,-200,200);
127  BOOK1D(DeltaTCMET,"Track Corrected",400,-200,200);
128  BOOK1D(DeltaTCPhi,"Track Corrected", 1000, -3.2, 3.2);
129  BOOK1D(DeltaTCSET,"Track Corrected",400,-200,200);
130  BOOK2D(TCSETvsDeltaTCMET,"Track Corrected",200, 0.0, 1000.0, 400, -200.0, 200.0);
131  BOOK2D(TCSETvsDeltaTCSET,"Track Corrected",200, 0.0, 1000.0, 400, -200.0, 200.0);
132  profileTCSETvsTCSETresp = new TProfile("#DeltaTCSET / trueSET vs trueSET", "", 200, 0.0, 1000.0, -1.0, 1.0);
133  profileTCMETvsTCMETresp = new TProfile("#DeltaTCMET / trueMET vs trueMET", "", 200, 0.0, 200.0, -1.0, 1.0);
134 
135  BOOK1D(TCMET,"Track Corrected",500,0,1000);
136  BOOK2D(DeltaTCMETvstrueMET,"Track Corrected", 500, 0.0, 1000.0, 400, -200.0, 200.0);
137  BOOK2D(DeltaTCPhivstrueMET,"Track Corrected", 500, 0.0, 1000.0, 400, -3.2, 3.2);
138 
139  BOOK2D(DeltaTCMEXvstrueSET,"Track Corrected", 200, 0.0, 1000.0, 400, -200.0, 200.0);
140  BOOK2D(TCMETvstrueMET,"Track Corrected", 200, 0.0, 1000.0, 500, 0.0, 1000.0);
141 
142  //BOOK1D(meanPF, "Mean PFMEX", 100, 0.0, 1600.0);
143  //BOOK1D(meanCalo, "Mean CaloMEX", 100, 0.0, 1600.0);
144  //BOOK1D(sigmaPF, "#sigma(PFMEX)", 100, 0.0, 1600.0);
145  //BOOK1D(sigmaCalo, "#sigma(CaloMEX)", 100, 0.0, 1600.0);
146  //BOOK1D(rmsPF, "RMS(PFMEX)", 100, 0.0, 1600.0);
147  //BOOK1D(rmsCalo, "RMS(CaloMEX)", 100, 0.0, 1600.0);
148 
149  // Set Axis Titles
150  // delta Pt or E quantities for Barrel and Endcap
151  SETAXES(MEX, "MEX", "Events");
152  SETAXES(DeltaMEX, "#DeltaMEX", "Events");
153  SETAXES(DeltaMET, "#DeltaMET", "Events");
154  SETAXES(DeltaPhi, "#Delta#phi", "Events");
155  SETAXES(DeltaSET, "#DeltaSET", "Events");
156  SETAXES(SETvsDeltaMET, "SET", "#DeltaMET");
157  SETAXES(SETvsDeltaSET, "SET", "#DeltaSET");
158 
159  SETAXES(DeltaPFMETvstrueMET, "trueMET", "#DeltaMET");
160  SETAXES(DeltaCaloMETvstrueMET, "trueMET", "#DeltaCaloMET");
161  SETAXES(DeltaPFPhivstrueMET, "trueMET", "#DeltaPFPhi");
162  SETAXES(DeltaCaloPhivstrueMET, "trueMET", "#DeltaCaloPhi");
163  SETAXES(CaloMETvstrueMET, "trueMET", "CaloMET");
164  SETAXES(PFMETvstrueMET, "trueMET", "PFMET");
165  SETAXES(DeltaCaloMEXvstrueSET, "trueSET", "#DeltaCaloMEX");
166  SETAXES(DeltaPFMEXvstrueSET, "trueSET", "#DeltaPFMEX");
167  SETAXES(TrueMET, "trueMET", "Events");
168  SETAXES(CaloMET, "CaloMET", "Events");
169  SETAXES(PFMET, "PFMET", "Events");
170  SETAXES(TCMET, "TCMET", "Events");
171 
172  SETAXES(CaloMEX, "MEX", "Events");
173  SETAXES(DeltaCaloMEX, "#DeltaMEX", "Events");
174  SETAXES(DeltaCaloMET, "#DeltaMET", "Events");
175  SETAXES(DeltaCaloPhi, "#Delta#phi", "Events");
176  SETAXES(DeltaCaloSET, "#DeltaSET", "Events");
177  SETAXES(CaloSETvsDeltaCaloMET, "SET", "#DeltaMET");
178  SETAXES(CaloSETvsDeltaCaloSET, "SET", "#DeltaSET");
179 
180  SETAXES(TCMEX, "MEX", "Events");
181  SETAXES(DeltaTCMEX, "#DeltaMEX", "Events");
182  SETAXES(DeltaTCMET, "#DeltaMET", "Events");
183  SETAXES(DeltaTCPhi, "#Delta#phi", "Events");
184  SETAXES(DeltaTCSET, "#DeltaSET", "Events");
185  SETAXES(TCSETvsDeltaTCMET, "SET", "#DeltaMET");
186  SETAXES(TCSETvsDeltaTCSET, "SET", "#DeltaSET");
187 
188  SETAXES(DeltaTCMETvstrueMET, "trueMET", "#DeltaTCMET");
189  SETAXES(DeltaTCPhivstrueMET, "trueMET", "#DeltaTCPhi");
190 
191  SETAXES(DeltaTCMEXvstrueSET, "trueSET", "#DeltaTCMEX");
192  SETAXES(TCMETvstrueMET, "trueMET", "TCMET");
193 }
194 
195 
196 //void PFMETBenchmark::process(const reco::PFMETCollection& pfMets, const reco::GenMETCollection& genMets) {
197 void PFMETBenchmark::process( const reco::PFMETCollection& pfMets,
198  const reco::GenParticleCollection& genParticleList,
199  const reco::CaloMETCollection& caloMets,
200  const reco::METCollection& tcMets ) {
201  calculateQuantities(pfMets, genParticleList, caloMets, tcMets);
202  if (debug_) {
203  cout << " =========PFMET " << rec_met << ", " << rec_phi << endl;
204  cout << " =========GenMET " << true_met << ", " << true_phi << endl;
205  cout << " =========CaloMET " << calo_met << ", " << calo_phi << endl;
206  cout << " =========TCMET " << tc_met << ", " << tc_phi << endl;
207  }
208  // fill histograms
209  hTrueMET->Fill(true_met);
210  // delta Pt or E quantities
211  // PF
212  hDeltaMET->Fill( rec_met - true_met );
213  hMEX->Fill( rec_mex );
214  hMEX->Fill( rec_mey );
215  hDeltaMEX->Fill( rec_mex - true_mex );
216  hDeltaMEX->Fill( rec_mey - true_mey );
217  hDeltaPhi->Fill( rec_phi - true_phi );
218  hDeltaSET->Fill( rec_set - true_set );
219  if( true_met > 5.0 ) hSETvsDeltaMET->Fill( rec_set, rec_met - true_met );
220  else hSETvsDeltaMET->Fill( rec_set, rec_mex - true_mex );
221  hSETvsDeltaSET->Fill( rec_set, rec_set - true_set );
222  if( true_met > 5.0 ) profileMETvsMETresp->Fill(true_met, (rec_met-true_met)/true_met);
223  profileSETvsSETresp->Fill(true_set, (rec_set-true_set)/true_set);
224 
225  hDeltaPFMETvstrueMET->Fill(true_met, rec_met - true_met);
226  hDeltaPFPhivstrueMET->Fill(true_met, rec_phi - true_phi);
227  hPFMETvstrueMET->Fill(true_met, rec_met);
228  hPFMET->Fill(rec_met);
229 
230  hDeltaPFMEXvstrueSET->Fill(true_set, rec_mex - true_mex);
231  hDeltaPFMEXvstrueSET->Fill(true_set, rec_mey - true_mey);
232 
233  // Calo
234  hDeltaCaloMET->Fill( calo_met - true_met );
235  hCaloMEX->Fill( calo_mex );
236  hCaloMEX->Fill( calo_mey );
237  hDeltaCaloMEX->Fill( calo_mex - true_mex);
238  hDeltaCaloMEX->Fill( calo_mey - true_mey);
239  hDeltaCaloPhi->Fill( calo_phi - true_phi );
240  hDeltaCaloSET->Fill( calo_set - true_set );
241  if( true_met > 5.0 ) hCaloSETvsDeltaCaloMET->Fill( calo_set, calo_met - true_met );
242  else hCaloSETvsDeltaCaloMET->Fill( calo_set, calo_mex - true_mex);
243  hCaloSETvsDeltaCaloSET->Fill( calo_set, calo_set - true_set );
244  if( true_met > 5.0 ) profileCaloMETvsCaloMETresp->Fill(true_met, (calo_met-true_met)/true_met);
245  profileCaloSETvsCaloSETresp->Fill(true_set, (calo_set-true_set)/true_set);
246 
247  hDeltaCaloMETvstrueMET->Fill(true_met, calo_met - true_met);
248  hDeltaCaloPhivstrueMET->Fill(true_met, calo_phi - true_phi);
249  hCaloMETvstrueMET->Fill(true_met, calo_met);
250  hCaloMET->Fill(calo_met);
251 
252  hDeltaCaloMEXvstrueSET->Fill(true_set, calo_mex - true_mex);
253  hDeltaCaloMEXvstrueSET->Fill(true_set, calo_mey - true_mey);
254 
255  // TC
256  hDeltaTCMET->Fill( tc_met - true_met );
257  hTCMET->Fill( tc_met );
258  hTCMEX->Fill( tc_mex );
259  hTCMEX->Fill( tc_mey );
260  hDeltaTCMEX->Fill( tc_mex - true_mex);
261  hDeltaTCMEX->Fill( tc_mey - true_mey);
262  hDeltaTCPhi->Fill( tc_phi - true_phi );
263  hDeltaTCSET->Fill( tc_set - true_set );
264  if( true_met > 5.0 ) hTCSETvsDeltaTCMET->Fill( tc_set, tc_met - true_met );
265  else hTCSETvsDeltaTCMET->Fill( tc_set, tc_mex - true_mex);
266  hTCSETvsDeltaTCSET->Fill( tc_set, tc_set - true_set );
267  if( true_met > 5.0 ) profileTCMETvsTCMETresp->Fill(true_met, (tc_met-true_met)/true_met);
268  profileTCSETvsTCSETresp->Fill(true_set, (tc_set-true_set)/true_set);
269 
270  hDeltaTCMETvstrueMET->Fill( true_met, tc_met - true_met);
271  hDeltaTCPhivstrueMET->Fill( true_met, tc_phi - true_phi);
272  hDeltaTCMEXvstrueSET->Fill( true_set, tc_mex - true_mex);
273  hDeltaTCMEXvstrueSET->Fill( true_set, tc_mey - true_mey);
274  hTCMETvstrueMET->Fill( true_met, tc_met );
275 
276 }
277 
278 void PFMETBenchmark::calculateQuantities( const reco::PFMETCollection& pfMets,
279  const reco::GenParticleCollection& genParticleList,
280  const reco::CaloMETCollection& caloMets,
281  const reco::METCollection& tcMets)
282 {
283  const reco::PFMET& pfm = pfMets[0];
284  const reco::CaloMET& cm = caloMets[0];
285  const reco::MET& tcm = tcMets[0];
286 
287  double trueMEY = 0.0;
288  double trueMEX = 0.0;;
289  true_set = 0.0;;
290 
291  // for( genParticle = genParticleList.begin(); genParticle != genParticleList.end(); genParticle++ )
292  for( unsigned i = 0; i < genParticleList.size(); i++ ) {
293  if( genParticleList[i].status() == 1 && fabs(genParticleList[i].eta()) < 5.0 ) {
294  if( std::abs(genParticleList[i].pdgId()) == 12 ||
295  std::abs(genParticleList[i].pdgId()) == 14 ||
296  std::abs(genParticleList[i].pdgId()) == 16 ||
297  std::abs(genParticleList[i].pdgId()) < 7 ||
298  std::abs(genParticleList[i].pdgId()) == 21 ) {
299  trueMEX += genParticleList[i].px();
300  trueMEY += genParticleList[i].py();
301  } else {
302  true_set += genParticleList[i].pt();
303  }
304  }
305  }
306  true_mex = -trueMEX;
307  true_mey = -trueMEY;
308  true_met = sqrt( trueMEX*trueMEX + trueMEY*trueMEY );
309  true_phi = atan2(trueMEY,trueMEX);
310  rec_met = pfm.pt();
311  rec_mex = pfm.px();
312  rec_mex = pfm.py();
313  rec_phi = pfm.phi();
314  rec_set = pfm.sumEt();
315  calo_met = cm.pt();
316  calo_mex = cm.px();
317  calo_mey = cm.py();
318  calo_phi = cm.phi();
319  calo_set = cm.sumEt();
320  tc_met = tcm.pt();
321  tc_mex = tcm.px();
322  tc_mey = tcm.py();
323  tc_phi = tcm.phi();
324  tc_set = tcm.sumEt();
325 
326  if (debug_) {
327  cout << " =========PFMET " << rec_met << ", " << rec_phi << endl;
328  cout << " =========trueMET " << true_met << ", " << true_phi << endl;
329  cout << " =========CaloMET " << calo_met << ", " << calo_phi << endl;
330  cout << " =========TCMET " << tc_met << ", " << tc_phi << endl;
331  }
332 }
333 
334 void PFMETBenchmark::calculateQuantities( const reco::PFMETCollection& pfMets,
335  const reco::GenParticleCollection& genParticleList,
336  const reco::CaloMETCollection& caloMets,
337  const reco::METCollection& tcMets,
338  const std::vector<reco::CaloJet> caloJets,
339  const std::vector<reco::CaloJet> corr_caloJets)
340 {
341  const reco::PFMET& pfm = pfMets[0];
342  const reco::CaloMET& cm = caloMets[0];
343  const reco::MET& tcm = tcMets[0];
344 
345  double trueMEY = 0.0;
346  double trueMEX = 0.0;;
347  true_set = 0.0;;
348 
349  // for( genParticle = genParticleList.begin(); genParticle != genParticleList.end(); genParticle++ )
350  for( unsigned i = 0; i < genParticleList.size(); i++ ) {
351  if( genParticleList[i].status() == 1 && fabs(genParticleList[i].eta()) < 5.0 ) {
352  if( std::abs(genParticleList[i].pdgId()) == 12 ||
353  std::abs(genParticleList[i].pdgId()) == 14 ||
354  std::abs(genParticleList[i].pdgId()) == 16 ||
355  std::abs(genParticleList[i].pdgId()) < 7 ||
356  std::abs(genParticleList[i].pdgId()) == 21 ) {
357  trueMEX += genParticleList[i].px();
358  trueMEY += genParticleList[i].py();
359  } else {
360  true_set += genParticleList[i].pt();
361  }
362  }
363  }
364  true_mex = -trueMEX;
365  true_mey = -trueMEY;
366  true_met = sqrt( trueMEX*trueMEX + trueMEY*trueMEY );
367  true_phi = atan2(trueMEY,trueMEX);
368  rec_met = pfm.pt();
369  rec_mex = pfm.px();
370  rec_mex = pfm.py();
371  rec_phi = pfm.phi();
372  rec_set = pfm.sumEt();
373 
374  // propagation of the JEC to the caloMET:
375  double caloJetCorPX = 0.0;
376  double caloJetCorPY = 0.0;
377 
378  for(unsigned int caloJetc=0;caloJetc<caloJets.size();++caloJetc)
379  {
380  //std::cout << "caloJets[" << caloJetc << "].pt() = " << caloJets[caloJetc].pt() << std::endl;
381  //std::cout << "caloJets[" << caloJetc << "].phi() = " << caloJets[caloJetc].phi() << std::endl;
382  //std::cout << "caloJets[" << caloJetc << "].eta() = " << caloJets[caloJetc].eta() << std::endl;
383  //}
384  for(unsigned int corr_caloJetc=0;corr_caloJetc<corr_caloJets.size();++corr_caloJetc)
385  {
386  //std::cout << "corr_caloJets[" << corr_caloJetc << "].pt() = " << corr_caloJets[corr_caloJetc].pt() << std::endl;
387  //std::cout << "corr_caloJets[" << corr_caloJetc << "].phi() = " << corr_caloJets[corr_caloJetc].phi() << std::endl;
388  //std::cout << "corr_caloJets[" << corr_caloJetc << "].eta() = " << corr_caloJets[corr_caloJetc].eta() << std::endl;
389  //}
390  Float_t DeltaPhi = corr_caloJets[corr_caloJetc].phi() - caloJets[caloJetc].phi();
391  Float_t DeltaEta = corr_caloJets[corr_caloJetc].eta() - caloJets[caloJetc].eta();
392  Float_t DeltaR2 = DeltaPhi*DeltaPhi + DeltaEta*DeltaEta;
393  if( DeltaR2 < 0.0001 && caloJets[caloJetc].pt() > 20.0 )
394  {
395  caloJetCorPX += (corr_caloJets[corr_caloJetc].px() - caloJets[caloJetc].px());
396  caloJetCorPY += (corr_caloJets[corr_caloJetc].py() - caloJets[caloJetc].py());
397  }
398  }
399  }
400  double corr_calomet=sqrt((cm.px()-caloJetCorPX)*(cm.px()-caloJetCorPX)+(cm.py()-caloJetCorPY)*(cm.py()-caloJetCorPY));
401  calo_met = corr_calomet;
402  calo_mex = cm.px()-caloJetCorPX;
403  calo_mey = cm.py()-caloJetCorPY;
404  calo_phi = atan2((cm.py()-caloJetCorPY),(cm.px()-caloJetCorPX));
405  //calo_met = cm.pt();
406  //calo_mex = cm.px();
407  //calo_mey = cm.py();
408  //calo_phi = cm.phi();
409 
410  calo_set = cm.sumEt();
411  tc_met = tcm.pt();
412  tc_mex = tcm.px();
413  tc_mey = tcm.py();
414  tc_phi = tcm.phi();
415  tc_set = tcm.sumEt();
416 
417  if (debug_) {
418  cout << " =========PFMET " << rec_met << ", " << rec_phi << endl;
419  cout << " =========trueMET " << true_met << ", " << true_phi << endl;
420  cout << " =========CaloMET " << calo_met << ", " << calo_phi << endl;
421  cout << " =========TCMET " << tc_met << ", " << tc_phi << endl;
422  }
423 }
424 
425 //double fitf(double *x, double *par)
426 //{
427 // double fitval = sqrt( par[0]*par[0] +
428 // par[1]*par[1]*(x[0]-par[3]) +
429 // par[2]*par[2]*(x[0]-par[3])*(x[0]-par[3]) );
430 // return fitval;
431 //}
432 
433 void PFMETBenchmark::analyse()
434 {
435  //Define fit functions and histograms
436  //TF1* func1 = new TF1("fit1", fitf, 0, 40, 4);
437  //TF1* func2 = new TF1("fit2", fitf, 0, 40, 4);
438  //TF1* func3 = new TF1("fit3", fitf, 0, 40, 4);
439  //TF1* func4 = new TF1("fit4", fitf, 0, 40, 4);
440 
441  //fit gaussian to Delta MET corresponding to different slices in MET, store fit values (mean,sigma) in histos
446 
453 
454  // Make the MET resolution versus SET plot
455  /*
456  TCanvas* canvas_MetResVsRecoSet = new TCanvas("MetResVsRecoSet", "MET Sigma vs Reco SET", 500,500);
457  hsigmaPF->SetStats(0);
458  func1->SetLineColor(1);
459  func1->SetParNames("Noise", "Stochastic", "Constant", "Offset");
460  func1->SetParameters(10.0, 0.8, 0.1, 100.0);
461  hsigmaPF->Fit("fit1", "", "", 100.0, 900.0);
462  func2->SetLineColor(2);
463  func2->SetParNames("Noise", "Stochastic", "Constant", "Offset");
464  func2->SetParameters(10.0, 0.8, 0.1, 100.0);
465  hsigmaCalo->Fit("fit2", "", "", 100.0, 900.0);
466  func3->SetLineColor(4);
467  func3->SetParNames("Noise", "Stochastic", "Constant", "Offset");
468  func3->SetParameters(10.0, 0.8, 0.1, 100.0);
469  hrmsPF->Fit("fit3", "", "", 100.0, 900.0);
470  func4->SetLineColor(6);
471  func4->SetParNames("Noise", "Stochastic", "Constant", "Offset");
472  func4->SetParameters(10.0, 0.8, 0.1, 100.0);
473  hrmsCalo->Fit("fit4", "", "", 100.0, 900.0);
474  (hsigmaPF->GetYaxis())->SetRangeUser( 0.0, 50.0);
475  hsigmaPF->SetLineWidth(2);
476  hsigmaPF->SetLineColor(1);
477  hsigmaPF->Draw();
478  hsigmaCalo->SetLineWidth(2);
479  hsigmaCalo->SetLineColor(2);
480  hsigmaCalo->Draw("SAME");
481  hrmsPF->SetLineWidth(2);
482  hrmsPF->SetLineColor(4);
483  hrmsPF->Draw("SAME");
484  hrmsCalo->SetLineWidth(2);
485  hrmsCalo->SetLineColor(6);
486  hrmsCalo->Draw("SAME");
487  */
488 
489  // Make the SET response versus SET plot
490  /*
491  TCanvas* canvas_SetRespVsTrueSet = new TCanvas("SetRespVsTrueSet", "SET Response vs True SET", 500,500);
492  profileSETvsSETresp->SetStats(0);
493  profileSETvsSETresp->SetStats(0);
494  (profileSETvsSETresp->GetYaxis())->SetRangeUser(-1.0, 1.0);
495  profileSETvsSETresp->SetLineWidth(2);
496  profileSETvsSETresp->SetLineColor(4);
497  profileSETvsSETresp->Draw();
498  profileCaloSETvsCaloSETresp->SetLineWidth(2);
499  profileCaloSETvsCaloSETresp->SetLineColor(2);
500  profileCaloSETvsCaloSETresp->Draw("SAME");
501  */
502 
503  // Make the MET response versus MET plot
504  /*
505  TCanvas* canvas_MetRespVsTrueMet = new TCanvas("MetRespVsTrueMet", "MET Response vs True MET", 500,500);
506  profileMETvsMETresp->SetStats(0);
507  profileMETvsMETresp->SetStats(0);
508  (profileMETvsMETresp->GetYaxis())->SetRangeUser(-1.0, 1.0);
509  profileMETvsMETresp->SetLineWidth(2);
510  profileMETvsMETresp->SetLineColor(4);
511  profileMETvsMETresp->Draw();
512  profileCaloMETvsCaloMETresp->SetLineWidth(2);
513  profileCaloMETvsCaloMETresp->SetLineColor(2);
514  profileCaloMETvsCaloMETresp->Draw("SAME");
515  */
516 
517  //print the resulting plots to file
518  /*
519  canvas_MetResVsRecoSet->Print("MetResVsRecoSet.ps");
520  canvas_SetRespVsTrueSet->Print("SetRespVsTrueSet.ps");
521  canvas_MetRespVsTrueMet->Print("MetRespVsTrueMet.ps");
522  */
523 }
524 
525 //void PFMETBenchmark::FitSlicesInY(TH2F* h, TH1F* mean, TH1F* sigma, bool doGausFit, int type )
526 //{
527 // TAxis *fXaxis = h->GetXaxis();
528 // TAxis *fYaxis = h->GetYaxis();
529 // Int_t nbins = fXaxis->GetNbins();
530 //
531 // //std::cout << "nbins = " << nbins << std::endl;
532 //
533 // Int_t binmin = 1;
534 // Int_t binmax = nbins;
535 // TString option = "QNR";
536 // TString opt = option;
537 // opt.ToLower();
538 // Float_t ngroup = 1;
539 // ngroup = 1;
540 //
541 // //std::cout << "fYaxis->GetXmin() = " << fYaxis->GetXmin() << std::endl;
542 // //std::cout << "fYaxis->GetXmax() = " << fYaxis->GetXmax() << std::endl;
543 //
544 // //default is to fit with a gaussian
545 // TF1 *f1 = 0;
546 // if (f1 == 0)
547 // {
548 // //f1 = (TF1*)gROOT->GetFunction("gaus");
549 // if (f1 == 0) f1 = new TF1("gaus","gaus", fYaxis->GetXmin(), fYaxis->GetXmax());
550 // else f1->SetRange( fYaxis->GetXmin(), fYaxis->GetXmax());
551 // }
552 // Int_t npar = f1->GetNpar();
553 //
554 // //std::cout << "npar = " << npar << std::endl;
555 //
556 // if (npar <= 0) return;
557 //
558 // Double_t *parsave = new Double_t[npar];
559 // f1->GetParameters(parsave);
560 //
562 // Int_t ipar;
563 // TH1F **hlist = new TH1F*[npar];
564 // char *name = new char[2000];
565 // char *title = new char[2000];
566 // const TArrayD *bins = fXaxis->GetXbins();
567 // for( ipar=0; ipar < npar; ipar++ )
568 // {
569 // if( ipar == 1 )
570 // if( type == 1 ) sprintf(name,"meanPF");
571 // else sprintf(name,"meanCalo");
572 // else
573 // if( doGausFit )
574 // if( type == 1 ) sprintf(name,"sigmaPF");
575 // else sprintf(name,"sigmaCalo");
576 // else
577 // if( type == 1 ) sprintf(name,"rmsPF");
578 // else sprintf(name,"rmsCalo");
579 // if( type == 1 ) sprintf(title,"Particle Flow");
580 // else sprintf(title,"Calorimeter");
581 // delete gDirectory->FindObject(name);
582 // if (bins->fN == 0)
583 // hlist[ipar] = new TH1F(name,title, nbins, fXaxis->GetXmin(), fXaxis->GetXmax());
584 // else
585 // hlist[ipar] = new TH1F(name,title, nbins,bins->fArray);
586 // hlist[ipar]->GetXaxis()->SetTitle(fXaxis->GetTitle());
587 // }
588 // sprintf(name,"test_chi2");
589 // delete gDirectory->FindObject(name);
590 // TH1F *hchi2 = new TH1F(name,"chisquare", nbins, fXaxis->GetXmin(), fXaxis->GetXmax());
591 // hchi2->GetXaxis()->SetTitle(fXaxis->GetTitle());
592 //
593 // //Loop on all bins in X, generate a projection along Y
594 // Int_t bin;
595 // Int_t nentries;
596 // for( bin = (Int_t) binmin; bin <= (Int_t) binmax; bin += (Int_t)ngroup )
597 // {
598 // TH1F *hpy = (TH1F*) h->ProjectionY("_temp", (Int_t) bin, (Int_t) (bin + ngroup - 1), "e");
599 // if(hpy == 0) continue;
600 // nentries = Int_t( hpy->GetEntries() );
601 // if(nentries == 0 ) {delete hpy; continue;}
602 // f1->SetParameters(parsave);
603 // hpy->Fit( f1, opt.Data() );
604 // Int_t npfits = f1->GetNumberFitPoints();
605 // //cout << "bin = " << bin << "; Npfits = " << npfits << "; npar = " << npar << endl;
606 // if( npfits > npar )
607 // {
608 // Int_t biny = bin + (Int_t)ngroup/2;
609 // for( ipar=0; ipar < npar; ipar++ )
610 // {
611 // if( doGausFit ) hlist[ipar]->Fill( fXaxis->GetBinCenter(biny), f1->GetParameter(ipar) );
612 // else hlist[ipar]->Fill( fXaxis->GetBinCenter(biny), hpy->GetRMS() );
613 // //cout << "bin[" << bin << "]: RMS = " << hpy->GetRMS() << "; sigma = " << f1->GetParameter(ipar) << endl;
614 // hlist[ipar]->SetBinError( biny, f1->GetParError(ipar) );
615 // }
616 // hchi2->Fill( fXaxis->GetBinCenter(biny), f1->GetChisquare()/(npfits-npar) );
617 // }
618 // delete hpy;
619 // ngroup += ngroup*0.2;//0.1 //used for non-uniform binning
620 // }
621 // *mean = *hlist[1];
622 // *sigma = *hlist[2];
623 // cout << "Entries = " << hlist[0]->GetEntries() << endl;
624 //}
625 
626 double
627 PFMETBenchmark::mpi_pi(double angle) {
628 
629  const double pi = 3.14159265358979323;
630  const double pi2 = pi*2.;
631  while(angle>pi) angle -= pi2;
632  while(angle<-pi) angle += pi2;
633  return angle;
634 
635 }
PFMETBenchmark::profileTCSETvsTCSETresp
TProfile * profileTCSETvsTCSETresp
Definition: PFMETBenchmark.h:141
reco::GenParticleCollection
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
Definition: GenParticleFwd.h:10
PFMETBenchmark::profileTCMETvsTCMETresp
TProfile * profileTCMETvsTCMETresp
Definition: PFMETBenchmark.h:142
i
int i
Definition: DBlmapReader.cc:9
PFMETBenchmark::hCaloSETvsDeltaCaloSET
TH2F * hCaloSETvsDeltaCaloSET
Definition: PFMETBenchmark.h:133
PFMETBenchmark::hDeltaPFPhivstrueMET
TH2F * hDeltaPFPhivstrueMET
Definition: PFMETBenchmark.h:118
PFMETBenchmark::hDeltaTCPhivstrueMET
TH2F * hDeltaTCPhivstrueMET
Definition: PFMETBenchmark.h:129
PFMETBenchmark::hDeltaTCMETvstrueMET
TH2F * hDeltaTCMETvstrueMET
Definition: PFMETBenchmark.h:128
PFMETBenchmark::profileMETvsMETresp
TProfile * profileMETvsMETresp
Definition: PFMETBenchmark.h:107
PFMETBenchmark::hDeltaCaloMETvstrueMET
TH2F * hDeltaCaloMETvstrueMET
Definition: PFMETBenchmark.h:117
PFMETBenchmark::~PFMETBenchmark
virtual ~PFMETBenchmark()
Definition: PFMETBenchmark.cc:36
PFMETBenchmark::hCaloMETvstrueMET
TH2F * hCaloMETvstrueMET
Definition: PFMETBenchmark.h:120
BOOK1D
#define BOOK1D(name, title, nbinsx, lowx, highx)
Definition: PFMETBenchmark.cc:8
PFMETBenchmark::calculateQuantities
void calculateQuantities(const reco::PFMETCollection &, const reco::GenParticleCollection &, const reco::CaloMETCollection &, const reco::METCollection &)
Definition: PFMETBenchmark.cc:278
PFMETBenchmark::hDeltaCaloPhivstrueMET
TH2F * hDeltaCaloPhivstrueMET
Definition: PFMETBenchmark.h:119
DQMStore::save
void save(const std::string &filename, const std::string &path="", const std::string &pattern="", const std::string &rewrite="", SaveReferenceTag ref=SaveWithReference, int minStatus=dqm::qstatus::STATUS_OK, const std::string &fileupdate="RECREATE")
Definition: DQMStore.cc:2113
pi2
const double pi2
Definition: Thrust.cc:5
abs
#define abs(x)
Definition: mlp_lapack.h:159
PFMETBenchmark::profileSETvsSETresp
TProfile * profileSETvsSETresp
Definition: PFMETBenchmark.h:106
NULL
#define NULL
Definition: scimark2.h:8
SETAXES
#define SETAXES(name, xtitle, ytitle)
Definition: PFMETBenchmark.cc:19
PFMETBenchmark::hCaloSETvsDeltaCaloMET
TH2F * hCaloSETvsDeltaCaloMET
Definition: PFMETBenchmark.h:134
eta
T eta() const
Definition: Basic3DVectorLD.h:172
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:105
PFMETBenchmark::outputFile_
std::string outputFile_
Definition: PFMETBenchmark.h:161
reco::METCollection
std::vector< reco::MET > METCollection
collection of MET objects
Definition: METCollection.h:23
BOOK2D
#define BOOK2D(name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)
Definition: PFMETBenchmark.cc:13
PFMETBenchmark::mpi_pi
double mpi_pi(double angle)
Definition: PFMETBenchmark.cc:627
scaleCards.path
list path
Definition: scaleCards.py:51
PFMETBenchmark::process
void process(const reco::PFMETCollection &, const reco::GenParticleCollection &, const reco::CaloMETCollection &, const reco::METCollection &)
Definition: PFMETBenchmark.cc:197
reco::MET::sumEt
double sumEt() const
Definition: MET.h:48
reco::MET
Definition: MET.h:32
PFMETBenchmark::hDeltaPFMETvstrueMET
TH2F * hDeltaPFMETvstrueMET
Definition: PFMETBenchmark.h:116
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
PFMETBenchmark::hDeltaTCMEXvstrueSET
TH2F * hDeltaTCMEXvstrueSET
Definition: PFMETBenchmark.h:151
benchmarkLabel_
string benchmarkLabel_
Definition: PFJetBenchmarkAnalyzer.cc:75
PFMETBenchmark::hDeltaCaloMEXvstrueSET
TH2F * hDeltaCaloMEXvstrueSET
Definition: PFMETBenchmark.h:122
PFMETBenchmark::hDeltaPFMEXvstrueSET
TH2F * hDeltaPFMEXvstrueSET
Definition: PFMETBenchmark.h:123
PFMETBenchmark::profileCaloSETvsCaloSETresp
TProfile * profileCaloSETvsCaloSETresp
Definition: PFMETBenchmark.h:131
plotAgainstReco
bool plotAgainstReco
Definition: PFJetBenchmarkAnalyzer.cc:72
reco::CaloMETCollection
std::vector< reco::CaloMET > CaloMETCollection
collection of CaloMET objects
Definition: CaloMETCollection.h:21
reco::LeafCandidate::px
virtual double px() const
x coordinate of momentum vector
Definition: LeafCandidate.h:131
reco::LeafCandidate::pt
virtual double pt() const
transverse momentum
Definition: LeafCandidate.cc:127
reco::PFMETCollection
std::vector< reco::PFMET > PFMETCollection
collection of PFMET objects
Definition: PFMETCollection.h:21
PFMETBenchmark::setup
void setup(std::string Filename, bool debug, bool plotAgainstReco=0, std::string benchmarkLabel_="ParticleFlow", DQMStore *dbe_store=NULL)
Definition: PFMETBenchmark.cc:57
PFMETBenchmark::hTCSETvsDeltaTCSET
TH2F * hTCSETvsDeltaTCSET
Definition: PFMETBenchmark.h:143
PFMETBenchmark::hTCMETvstrueMET
TH2F * hTCMETvstrueMET
Definition: PFMETBenchmark.h:152
PFMETBenchmark::dbe_
DQMStore * dbe_
Definition: PFMETBenchmark.h:189
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
pi
double pi
Definition: Vx3DHLTAnalyzer.h:63
PFMETBenchmark::profileCaloMETvsCaloMETresp
TProfile * profileCaloMETvsCaloMETresp
Definition: PFMETBenchmark.h:132
benchmark_cfg.pdgId
tuple pdgId
Definition: benchmark_cfg.py:18
ntuplemaker.status
tuple status
Definition: ntuplemaker.py:245
debug
#define debug
Definition: MEtoEDMFormat.h:34
PFMETBenchmark::hTCSETvsDeltaTCMET
TH2F * hTCSETvsDeltaTCMET
Definition: PFMETBenchmark.h:144
reco::LeafCandidate::phi
virtual double phi() const
momentum azimuthal angle
Definition: LeafCandidate.h:139
reco::LeafCandidate::py
virtual double py() const
y coordinate of momentum vector
Definition: LeafCandidate.h:133
DQMStore::setCurrentFolder
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:429
PFMETBenchmark::hPFMETvstrueMET
TH2F * hPFMETvstrueMET
Definition: PFMETBenchmark.h:121
angle
T angle(T x1, T y1, T z1, T x2, T y2, T z2)
Definition: angle.h:11
Generated for CMSSW Reference Manual by doxygen 1.8.5