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