Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:20 

0001 #include "RecoParticleFlow/Benchmark/interface/PFMETBenchmark.h"
0002 
0003 // preprocessor macro for booking 1d histos with DQMStore -or- bare Root
0004 #define BOOK1D(name, title, nbinsx, lowx, highx) \
0005   h##name =                                      \
0006       dbe_ ? dbe_->book1D(#name, title, nbinsx, lowx, highx)->getTH1F() : new TH1F(#name, title, nbinsx, lowx, highx)
0007 
0008 // preprocessor macro for booking 2d histos with DQMStore -or- bare Root
0009 #define BOOK2D(name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)                              \
0010   h##name = dbe_ ? dbe_->book2D(#name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)->getTH2F() \
0011                  : new TH2F(#name, title, nbinsx, lowx, highx, nbinsy, lowy, highy)
0012 
0013 // all versions OK
0014 // preprocesor macro for setting axis titles
0015 #define SETAXES(name, xtitle, ytitle)    \
0016   h##name->GetXaxis()->SetTitle(xtitle); \
0017   h##name->GetYaxis()->SetTitle(ytitle)
0018 
0019 /*#define SET2AXES(name,xtitle,ytitle) \
0020   hE##name->GetXaxis()->SetTitle(xtitle); hE##name->GetYaxis()->SetTitle(ytitle);  hB##name->GetXaxis()->SetTitle(xtitle); hB##name->GetYaxis()->SetTitle(ytitle)
0021 */
0022 
0023 #define PT (plotAgainstReco_) ? "reconstructed P_{T}" : "generated P_{T}"
0024 
0025 using namespace reco;
0026 using namespace std;
0027 
0028 PFMETBenchmark::PFMETBenchmark() : file_(nullptr) {}
0029 
0030 PFMETBenchmark::~PFMETBenchmark() {
0031   if (file_)
0032     file_->Close();
0033 }
0034 
0035 void PFMETBenchmark::write() {
0036   // Store the DAQ Histograms
0037   if (!outputFile_.empty()) {
0038     if (dbe_)
0039       dbe_->save(outputFile_);
0040     // use bare Root if no DQM (FWLite applications)
0041     else if (file_) {
0042       file_->Write(outputFile_.c_str());
0043       cout << "Benchmark output written to file " << outputFile_.c_str() << endl;
0044       file_->Close();
0045     }
0046   } else
0047     cout << "No output file specified (" << outputFile_ << "). Results will not be saved!" << endl;
0048 }
0049 
0050 void PFMETBenchmark::setup(
0051     string Filename, bool debug, bool plotAgainstReco, string benchmarkLabel_, DQMStore* dbe_store) {
0052   debug_ = debug;
0053   plotAgainstReco_ = plotAgainstReco;
0054   outputFile_ = Filename;
0055   file_ = nullptr;
0056   dbe_ = dbe_store;
0057   // print parameters
0058   //cout<< "PFMETBenchmark Setup parameters =============================================="<<endl;
0059   cout << "Filename to write histograms " << Filename << endl;
0060   cout << "PFMETBenchmark debug " << debug_ << endl;
0061   cout << "plotAgainstReco " << plotAgainstReco_ << endl;
0062   cout << "benchmarkLabel " << benchmarkLabel_ << endl;
0063 
0064   // Book histogram
0065 
0066   // Establish DQM Store
0067   string path = "PFTask/Benchmarks/" + benchmarkLabel_ + "/";
0068   if (plotAgainstReco)
0069     path += "Reco";
0070   else
0071     path += "Gen";
0072   if (dbe_) {
0073     dbe_->setCurrentFolder(path);
0074   } else {
0075     file_ = new TFile(outputFile_.c_str(), "recreate");
0076     //    TTree * tr = new TTree("PFTast");
0077     //    tr->Branch("Benchmarks/ParticleFlow")
0078     cout << "Info: DQM is not available to provide data storage service. Using TFile to save histograms. " << endl;
0079   }
0080 
0081   // delta Pt or E quantities for Barrel
0082   BOOK1D(MEX, "Particle Flow", 400, -200, 200);
0083   BOOK1D(DeltaMEX, "Particle Flow", 400, -200, 200);
0084   BOOK1D(DeltaMET, "Particle Flow", 400, -200, 200);
0085   BOOK1D(DeltaPhi, "Particle Flow", 1000, -3.2, 3.2);
0086   BOOK1D(DeltaSET, "Particle Flow", 400, -200, 200);
0087   BOOK2D(SETvsDeltaMET, "Particle Flow", 200, 0.0, 1000.0, 400, -200.0, 200.0);
0088   BOOK2D(SETvsDeltaSET, "Particle Flow", 200, 0.0, 1000.0, 400, -200.0, 200.0);
0089   profileSETvsSETresp = new TProfile("#DeltaPFSET / trueSET vs trueSET", "", 200, 0.0, 1000.0, -1.0, 1.0);
0090   profileMETvsMETresp = new TProfile("#DeltaPFMET / trueMET vs trueMET", "", 50, 0.0, 200.0, -1.0, 1.0);
0091 
0092   BOOK1D(CaloMEX, "Calorimeter", 400, -200, 200);
0093   BOOK1D(DeltaCaloMEX, "Calorimeter", 400, -200, 200);
0094   BOOK1D(DeltaCaloMET, "Calorimeter", 400, -200, 200);
0095   BOOK1D(DeltaCaloPhi, "Calorimeter", 1000, -3.2, 3.2);
0096   BOOK1D(DeltaCaloSET, "Calorimeter", 400, -200, 200);
0097   BOOK2D(CaloSETvsDeltaCaloMET, "Calorimeter", 200, 0.0, 1000.0, 400, -200.0, 200.0);
0098   BOOK2D(CaloSETvsDeltaCaloSET, "Calorimeter", 200, 0.0, 1000.0, 400, -200.0, 200.0);
0099   profileCaloSETvsCaloSETresp = new TProfile("#DeltaCaloSET / trueSET vs trueSET", "", 200, 0.0, 1000.0, -1.0, 1.0);
0100   profileCaloMETvsCaloMETresp = new TProfile("#DeltaCaloMET / trueMET vs trueMET", "", 200, 0.0, 200.0, -1.0, 1.0);
0101 
0102   BOOK2D(DeltaPFMETvstrueMET, "Particle Flow", 500, 0.0, 1000.0, 400, -200.0, 200.0);
0103   BOOK2D(DeltaCaloMETvstrueMET, "Calorimeter", 500, 0.0, 1000.0, 400, -200.0, 200.0);
0104   BOOK2D(DeltaPFPhivstrueMET, "Particle Flow", 500, 0.0, 1000.0, 400, -3.2, 3.2);
0105   BOOK2D(DeltaCaloPhivstrueMET, "Calorimeter", 500, 0.0, 1000.0, 400, -3.2, 3.2);
0106   BOOK2D(CaloMETvstrueMET, "Calorimeter", 500, 0.0, 1000.0, 500, 0.0, 1000.0);
0107   BOOK2D(PFMETvstrueMET, "Particle Flow", 500, 0.0, 1000.0, 500, 0.0, 1000.0);
0108 
0109   BOOK2D(DeltaCaloMEXvstrueSET, "Calorimeter", 200, 0.0, 1000.0, 400, -200.0, 200.0);
0110   BOOK2D(DeltaPFMEXvstrueSET, "Particle Flow", 200, 0.0, 1000.0, 400, -200.0, 200.0);
0111 
0112   BOOK1D(TrueMET, "MC truth", 500, 0.0, 1000.0);
0113   BOOK1D(CaloMET, "Calorimeter", 500, 0.0, 1000.0);
0114   BOOK1D(PFMET, "Particle Flow", 500, 0.0, 1000.0);
0115 
0116   BOOK1D(TCMEX, "Track Corrected", 400, -200, 200);
0117   BOOK1D(DeltaTCMEX, "Track Corrected", 400, -200, 200);
0118   BOOK1D(DeltaTCMET, "Track Corrected", 400, -200, 200);
0119   BOOK1D(DeltaTCPhi, "Track Corrected", 1000, -3.2, 3.2);
0120   BOOK1D(DeltaTCSET, "Track Corrected", 400, -200, 200);
0121   BOOK2D(TCSETvsDeltaTCMET, "Track Corrected", 200, 0.0, 1000.0, 400, -200.0, 200.0);
0122   BOOK2D(TCSETvsDeltaTCSET, "Track Corrected", 200, 0.0, 1000.0, 400, -200.0, 200.0);
0123   profileTCSETvsTCSETresp = new TProfile("#DeltaTCSET / trueSET vs trueSET", "", 200, 0.0, 1000.0, -1.0, 1.0);
0124   profileTCMETvsTCMETresp = new TProfile("#DeltaTCMET / trueMET vs trueMET", "", 200, 0.0, 200.0, -1.0, 1.0);
0125 
0126   BOOK1D(TCMET, "Track Corrected", 500, 0, 1000);
0127   BOOK2D(DeltaTCMETvstrueMET, "Track Corrected", 500, 0.0, 1000.0, 400, -200.0, 200.0);
0128   BOOK2D(DeltaTCPhivstrueMET, "Track Corrected", 500, 0.0, 1000.0, 400, -3.2, 3.2);
0129 
0130   BOOK2D(DeltaTCMEXvstrueSET, "Track Corrected", 200, 0.0, 1000.0, 400, -200.0, 200.0);
0131   BOOK2D(TCMETvstrueMET, "Track Corrected", 200, 0.0, 1000.0, 500, 0.0, 1000.0);
0132 
0133   //BOOK1D(meanPF,    "Mean PFMEX", 100, 0.0, 1600.0);
0134   //BOOK1D(meanCalo,  "Mean CaloMEX", 100, 0.0, 1600.0);
0135   //BOOK1D(sigmaPF,   "#sigma(PFMEX)", 100, 0.0, 1600.0);
0136   //BOOK1D(sigmaCalo, "#sigma(CaloMEX)", 100, 0.0, 1600.0);
0137   //BOOK1D(rmsPF,     "RMS(PFMEX)", 100, 0.0, 1600.0);
0138   //BOOK1D(rmsCalo,   "RMS(CaloMEX)", 100, 0.0, 1600.0);
0139 
0140   // Set Axis Titles
0141   // delta Pt or E quantities for Barrel and Endcap
0142   SETAXES(MEX, "MEX", "Events");
0143   SETAXES(DeltaMEX, "#DeltaMEX", "Events");
0144   SETAXES(DeltaMET, "#DeltaMET", "Events");
0145   SETAXES(DeltaPhi, "#Delta#phi", "Events");
0146   SETAXES(DeltaSET, "#DeltaSET", "Events");
0147   SETAXES(SETvsDeltaMET, "SET", "#DeltaMET");
0148   SETAXES(SETvsDeltaSET, "SET", "#DeltaSET");
0149 
0150   SETAXES(DeltaPFMETvstrueMET, "trueMET", "#DeltaMET");
0151   SETAXES(DeltaCaloMETvstrueMET, "trueMET", "#DeltaCaloMET");
0152   SETAXES(DeltaPFPhivstrueMET, "trueMET", "#DeltaPFPhi");
0153   SETAXES(DeltaCaloPhivstrueMET, "trueMET", "#DeltaCaloPhi");
0154   SETAXES(CaloMETvstrueMET, "trueMET", "CaloMET");
0155   SETAXES(PFMETvstrueMET, "trueMET", "PFMET");
0156   SETAXES(DeltaCaloMEXvstrueSET, "trueSET", "#DeltaCaloMEX");
0157   SETAXES(DeltaPFMEXvstrueSET, "trueSET", "#DeltaPFMEX");
0158   SETAXES(TrueMET, "trueMET", "Events");
0159   SETAXES(CaloMET, "CaloMET", "Events");
0160   SETAXES(PFMET, "PFMET", "Events");
0161   SETAXES(TCMET, "TCMET", "Events");
0162 
0163   SETAXES(CaloMEX, "MEX", "Events");
0164   SETAXES(DeltaCaloMEX, "#DeltaMEX", "Events");
0165   SETAXES(DeltaCaloMET, "#DeltaMET", "Events");
0166   SETAXES(DeltaCaloPhi, "#Delta#phi", "Events");
0167   SETAXES(DeltaCaloSET, "#DeltaSET", "Events");
0168   SETAXES(CaloSETvsDeltaCaloMET, "SET", "#DeltaMET");
0169   SETAXES(CaloSETvsDeltaCaloSET, "SET", "#DeltaSET");
0170 
0171   SETAXES(TCMEX, "MEX", "Events");
0172   SETAXES(DeltaTCMEX, "#DeltaMEX", "Events");
0173   SETAXES(DeltaTCMET, "#DeltaMET", "Events");
0174   SETAXES(DeltaTCPhi, "#Delta#phi", "Events");
0175   SETAXES(DeltaTCSET, "#DeltaSET", "Events");
0176   SETAXES(TCSETvsDeltaTCMET, "SET", "#DeltaMET");
0177   SETAXES(TCSETvsDeltaTCSET, "SET", "#DeltaSET");
0178 
0179   SETAXES(DeltaTCMETvstrueMET, "trueMET", "#DeltaTCMET");
0180   SETAXES(DeltaTCPhivstrueMET, "trueMET", "#DeltaTCPhi");
0181 
0182   SETAXES(DeltaTCMEXvstrueSET, "trueSET", "#DeltaTCMEX");
0183   SETAXES(TCMETvstrueMET, "trueMET", "TCMET");
0184 }
0185 
0186 //void PFMETBenchmark::process(const reco::PFMETCollection& pfMets, const reco::GenMETCollection& genMets) {
0187 void PFMETBenchmark::process(const reco::PFMETCollection& pfMets,
0188                              const reco::GenParticleCollection& genParticleList,
0189                              const reco::CaloMETCollection& caloMets,
0190                              const reco::METCollection& tcMets) {
0191   calculateQuantities(pfMets, genParticleList, caloMets, tcMets);
0192   if (debug_) {
0193     cout << "  =========PFMET  " << rec_met << ", " << rec_phi << endl;
0194     cout << "  =========GenMET " << true_met << ", " << true_phi << endl;
0195     cout << "  =========CaloMET " << calo_met << ", " << calo_phi << endl;
0196     cout << "  =========TCMET " << tc_met << ", " << tc_phi << endl;
0197   }
0198   // fill histograms
0199   hTrueMET->Fill(true_met);
0200   // delta Pt or E quantities
0201   // PF
0202   hDeltaMET->Fill(rec_met - true_met);
0203   hMEX->Fill(rec_mex);
0204   hMEX->Fill(rec_mey);
0205   hDeltaMEX->Fill(rec_mex - true_mex);
0206   hDeltaMEX->Fill(rec_mey - true_mey);
0207   hDeltaPhi->Fill(rec_phi - true_phi);
0208   hDeltaSET->Fill(rec_set - true_set);
0209   if (true_met > 5.0)
0210     hSETvsDeltaMET->Fill(rec_set, rec_met - true_met);
0211   else
0212     hSETvsDeltaMET->Fill(rec_set, rec_mex - true_mex);
0213   hSETvsDeltaSET->Fill(rec_set, rec_set - true_set);
0214   if (true_met > 5.0)
0215     profileMETvsMETresp->Fill(true_met, (rec_met - true_met) / true_met);
0216   profileSETvsSETresp->Fill(true_set, (rec_set - true_set) / true_set);
0217 
0218   hDeltaPFMETvstrueMET->Fill(true_met, rec_met - true_met);
0219   hDeltaPFPhivstrueMET->Fill(true_met, rec_phi - true_phi);
0220   hPFMETvstrueMET->Fill(true_met, rec_met);
0221   hPFMET->Fill(rec_met);
0222 
0223   hDeltaPFMEXvstrueSET->Fill(true_set, rec_mex - true_mex);
0224   hDeltaPFMEXvstrueSET->Fill(true_set, rec_mey - true_mey);
0225 
0226   // Calo
0227   hDeltaCaloMET->Fill(calo_met - true_met);
0228   hCaloMEX->Fill(calo_mex);
0229   hCaloMEX->Fill(calo_mey);
0230   hDeltaCaloMEX->Fill(calo_mex - true_mex);
0231   hDeltaCaloMEX->Fill(calo_mey - true_mey);
0232   hDeltaCaloPhi->Fill(calo_phi - true_phi);
0233   hDeltaCaloSET->Fill(calo_set - true_set);
0234   if (true_met > 5.0)
0235     hCaloSETvsDeltaCaloMET->Fill(calo_set, calo_met - true_met);
0236   else
0237     hCaloSETvsDeltaCaloMET->Fill(calo_set, calo_mex - true_mex);
0238   hCaloSETvsDeltaCaloSET->Fill(calo_set, calo_set - true_set);
0239   if (true_met > 5.0)
0240     profileCaloMETvsCaloMETresp->Fill(true_met, (calo_met - true_met) / true_met);
0241   profileCaloSETvsCaloSETresp->Fill(true_set, (calo_set - true_set) / true_set);
0242 
0243   hDeltaCaloMETvstrueMET->Fill(true_met, calo_met - true_met);
0244   hDeltaCaloPhivstrueMET->Fill(true_met, calo_phi - true_phi);
0245   hCaloMETvstrueMET->Fill(true_met, calo_met);
0246   hCaloMET->Fill(calo_met);
0247 
0248   hDeltaCaloMEXvstrueSET->Fill(true_set, calo_mex - true_mex);
0249   hDeltaCaloMEXvstrueSET->Fill(true_set, calo_mey - true_mey);
0250 
0251   // TC
0252   hDeltaTCMET->Fill(tc_met - true_met);
0253   hTCMET->Fill(tc_met);
0254   hTCMEX->Fill(tc_mex);
0255   hTCMEX->Fill(tc_mey);
0256   hDeltaTCMEX->Fill(tc_mex - true_mex);
0257   hDeltaTCMEX->Fill(tc_mey - true_mey);
0258   hDeltaTCPhi->Fill(tc_phi - true_phi);
0259   hDeltaTCSET->Fill(tc_set - true_set);
0260   if (true_met > 5.0)
0261     hTCSETvsDeltaTCMET->Fill(tc_set, tc_met - true_met);
0262   else
0263     hTCSETvsDeltaTCMET->Fill(tc_set, tc_mex - true_mex);
0264   hTCSETvsDeltaTCSET->Fill(tc_set, tc_set - true_set);
0265   if (true_met > 5.0)
0266     profileTCMETvsTCMETresp->Fill(true_met, (tc_met - true_met) / true_met);
0267   profileTCSETvsTCSETresp->Fill(true_set, (tc_set - true_set) / true_set);
0268 
0269   hDeltaTCMETvstrueMET->Fill(true_met, tc_met - true_met);
0270   hDeltaTCPhivstrueMET->Fill(true_met, tc_phi - true_phi);
0271   hDeltaTCMEXvstrueSET->Fill(true_set, tc_mex - true_mex);
0272   hDeltaTCMEXvstrueSET->Fill(true_set, tc_mey - true_mey);
0273   hTCMETvstrueMET->Fill(true_met, tc_met);
0274 }
0275 
0276 void PFMETBenchmark::calculateQuantities(const reco::PFMETCollection& pfMets,
0277                                          const reco::GenParticleCollection& genParticleList,
0278                                          const reco::CaloMETCollection& caloMets,
0279                                          const reco::METCollection& tcMets) {
0280   const reco::PFMET& pfm = pfMets[0];
0281   const reco::CaloMET& cm = caloMets[0];
0282   const reco::MET& tcm = tcMets[0];
0283 
0284   double trueMEY = 0.0;
0285   double trueMEX = 0.0;
0286   ;
0287   true_set = 0.0;
0288   ;
0289 
0290   //  for( genParticle = genParticleList.begin(); genParticle != genParticleList.end(); genParticle++ )
0291   for (unsigned i = 0; i < genParticleList.size(); i++) {
0292     if (genParticleList[i].status() == 1 && fabs(genParticleList[i].eta()) < 5.0) {
0293       if (std::abs(genParticleList[i].pdgId()) == 12 || std::abs(genParticleList[i].pdgId()) == 14 ||
0294           std::abs(genParticleList[i].pdgId()) == 16 || std::abs(genParticleList[i].pdgId()) < 7 ||
0295           std::abs(genParticleList[i].pdgId()) == 21) {
0296         trueMEX += genParticleList[i].px();
0297         trueMEY += genParticleList[i].py();
0298       } else {
0299         true_set += genParticleList[i].pt();
0300       }
0301     }
0302   }
0303   true_mex = -trueMEX;
0304   true_mey = -trueMEY;
0305   true_met = sqrt(trueMEX * trueMEX + trueMEY * trueMEY);
0306   true_phi = atan2(trueMEY, trueMEX);
0307   rec_met = pfm.pt();
0308   rec_mex = pfm.px();
0309   rec_mex = pfm.py();
0310   rec_phi = pfm.phi();
0311   rec_set = pfm.sumEt();
0312   calo_met = cm.pt();
0313   calo_mex = cm.px();
0314   calo_mey = cm.py();
0315   calo_phi = cm.phi();
0316   calo_set = cm.sumEt();
0317   tc_met = tcm.pt();
0318   tc_mex = tcm.px();
0319   tc_mey = tcm.py();
0320   tc_phi = tcm.phi();
0321   tc_set = tcm.sumEt();
0322 
0323   if (debug_) {
0324     cout << "  =========PFMET  " << rec_met << ", " << rec_phi << endl;
0325     cout << "  =========trueMET " << true_met << ", " << true_phi << endl;
0326     cout << "  =========CaloMET " << calo_met << ", " << calo_phi << endl;
0327     cout << "  =========TCMET " << tc_met << ", " << tc_phi << endl;
0328   }
0329 }
0330 
0331 void PFMETBenchmark::calculateQuantities(const reco::PFMETCollection& pfMets,
0332                                          const reco::GenParticleCollection& genParticleList,
0333                                          const reco::CaloMETCollection& caloMets,
0334                                          const reco::METCollection& tcMets,
0335                                          const std::vector<reco::CaloJet>& caloJets,
0336                                          const std::vector<reco::CaloJet>& corr_caloJets) {
0337   const reco::PFMET& pfm = pfMets[0];
0338   const reco::CaloMET& cm = caloMets[0];
0339   const reco::MET& tcm = tcMets[0];
0340 
0341   double trueMEY = 0.0;
0342   double trueMEX = 0.0;
0343   ;
0344   true_set = 0.0;
0345   ;
0346 
0347   //  for( genParticle = genParticleList.begin(); genParticle != genParticleList.end(); genParticle++ )
0348   for (unsigned i = 0; i < genParticleList.size(); i++) {
0349     if (genParticleList[i].status() == 1 && fabs(genParticleList[i].eta()) < 5.0) {
0350       if (std::abs(genParticleList[i].pdgId()) == 12 || std::abs(genParticleList[i].pdgId()) == 14 ||
0351           std::abs(genParticleList[i].pdgId()) == 16 || std::abs(genParticleList[i].pdgId()) < 7 ||
0352           std::abs(genParticleList[i].pdgId()) == 21) {
0353         trueMEX += genParticleList[i].px();
0354         trueMEY += genParticleList[i].py();
0355       } else {
0356         true_set += genParticleList[i].pt();
0357       }
0358     }
0359   }
0360   true_mex = -trueMEX;
0361   true_mey = -trueMEY;
0362   true_met = sqrt(trueMEX * trueMEX + trueMEY * trueMEY);
0363   true_phi = atan2(trueMEY, trueMEX);
0364   rec_met = pfm.pt();
0365   rec_mex = pfm.px();
0366   rec_mex = pfm.py();
0367   rec_phi = pfm.phi();
0368   rec_set = pfm.sumEt();
0369 
0370   // propagation of the JEC to the caloMET:
0371   double caloJetCorPX = 0.0;
0372   double caloJetCorPY = 0.0;
0373 
0374   for (unsigned int caloJetc = 0; caloJetc < caloJets.size(); ++caloJetc) {
0375     //std::cout << "caloJets[" << caloJetc << "].pt() = " << caloJets[caloJetc].pt() << std::endl;
0376     //std::cout << "caloJets[" << caloJetc << "].phi() = " << caloJets[caloJetc].phi() << std::endl;
0377     //std::cout << "caloJets[" << caloJetc << "].eta() = " << caloJets[caloJetc].eta() << std::endl;
0378     //}
0379     for (unsigned int corr_caloJetc = 0; corr_caloJetc < corr_caloJets.size(); ++corr_caloJetc) {
0380       //std::cout << "corr_caloJets[" << corr_caloJetc << "].pt() = " << corr_caloJets[corr_caloJetc].pt() << std::endl;
0381       //std::cout << "corr_caloJets[" << corr_caloJetc << "].phi() = " << corr_caloJets[corr_caloJetc].phi() << std::endl;
0382       //std::cout << "corr_caloJets[" << corr_caloJetc << "].eta() = " << corr_caloJets[corr_caloJetc].eta() << std::endl;
0383       //}
0384       Float_t DeltaPhi = corr_caloJets[corr_caloJetc].phi() - caloJets[caloJetc].phi();
0385       Float_t DeltaEta = corr_caloJets[corr_caloJetc].eta() - caloJets[caloJetc].eta();
0386       Float_t DeltaR2 = DeltaPhi * DeltaPhi + DeltaEta * DeltaEta;
0387       if (DeltaR2 < 0.0001 && caloJets[caloJetc].pt() > 20.0) {
0388         caloJetCorPX += (corr_caloJets[corr_caloJetc].px() - caloJets[caloJetc].px());
0389         caloJetCorPY += (corr_caloJets[corr_caloJetc].py() - caloJets[caloJetc].py());
0390       }
0391     }
0392   }
0393   double corr_calomet =
0394       sqrt((cm.px() - caloJetCorPX) * (cm.px() - caloJetCorPX) + (cm.py() - caloJetCorPY) * (cm.py() - caloJetCorPY));
0395   calo_met = corr_calomet;
0396   calo_mex = cm.px() - caloJetCorPX;
0397   calo_mey = cm.py() - caloJetCorPY;
0398   calo_phi = atan2((cm.py() - caloJetCorPY), (cm.px() - caloJetCorPX));
0399   //calo_met = cm.pt();
0400   //calo_mex = cm.px();
0401   //calo_mey = cm.py();
0402   //calo_phi = cm.phi();
0403 
0404   calo_set = cm.sumEt();
0405   tc_met = tcm.pt();
0406   tc_mex = tcm.px();
0407   tc_mey = tcm.py();
0408   tc_phi = tcm.phi();
0409   tc_set = tcm.sumEt();
0410 
0411   if (debug_) {
0412     cout << "  =========PFMET  " << rec_met << ", " << rec_phi << endl;
0413     cout << "  =========trueMET " << true_met << ", " << true_phi << endl;
0414     cout << "  =========CaloMET " << calo_met << ", " << calo_phi << endl;
0415     cout << "  =========TCMET " << tc_met << ", " << tc_phi << endl;
0416   }
0417 }
0418 
0419 //double fitf(double *x, double *par)
0420 //{
0421 //  double fitval = sqrt( par[0]*par[0] +
0422 //          par[1]*par[1]*(x[0]-par[3]) +
0423 //          par[2]*par[2]*(x[0]-par[3])*(x[0]-par[3]) );
0424 //  return fitval;
0425 //}
0426 
0427 void PFMETBenchmark::analyse() {
0428   //Define fit functions and histograms
0429   //TF1* func1 = new TF1("fit1", fitf, 0, 40, 4);
0430   //TF1* func2 = new TF1("fit2", fitf, 0, 40, 4);
0431   //TF1* func3 = new TF1("fit3", fitf, 0, 40, 4);
0432   //TF1* func4 = new TF1("fit4", fitf, 0, 40, 4);
0433 
0434   //fit gaussian to Delta MET corresponding to different slices in MET, store fit values (mean,sigma) in histos
0435   ////FitSlicesInY(hSETvsDeltaMET, hmeanPF, hrmsPF, false, 1); //set option flag for RMS or gaussian
0436   ////FitSlicesInY(hSETvsDeltaMET, hmeanPF, hsigmaPF, true, 1); //set option flag for RMS or gaussian
0437   ////FitSlicesInY(hCaloSETvsDeltaCaloMET, hmeanCalo, hrmsCalo, false, 2);
0438   ////FitSlicesInY(hCaloSETvsDeltaCaloMET, hmeanCalo, hsigmaCalo, true, 2);
0439 
0440   ////SETAXES(meanPF,    "SET", "Mean(MEX)");
0441   ////SETAXES(meanCalo,  "SET", "Mean(MEX)");
0442   ////SETAXES(sigmaPF,   "SET", "#sigma(MEX)");
0443   ////SETAXES(sigmaCalo, "SET", "#sigma(MEX)");
0444   ////SETAXES(rmsPF,     "SET", "RMS(MEX)");
0445   ////SETAXES(rmsCalo,   "SET", "RMS(MEX)");
0446 
0447   // Make the MET resolution versus SET plot
0448   /*
0449   TCanvas* canvas_MetResVsRecoSet = new TCanvas("MetResVsRecoSet", "MET Sigma vs Reco SET", 500,500);
0450   hsigmaPF->SetStats(0); 
0451   func1->SetLineColor(1); 
0452   func1->SetParNames("Noise", "Stochastic", "Constant", "Offset");
0453   func1->SetParameters(10.0, 0.8, 0.1, 100.0);
0454   hsigmaPF->Fit("fit1", "", "", 100.0, 900.0);
0455   func2->SetLineColor(2); 
0456   func2->SetParNames("Noise", "Stochastic", "Constant", "Offset");
0457   func2->SetParameters(10.0, 0.8, 0.1, 100.0);
0458   hsigmaCalo->Fit("fit2", "", "", 100.0, 900.0);
0459   func3->SetLineColor(4); 
0460   func3->SetParNames("Noise", "Stochastic", "Constant", "Offset");
0461   func3->SetParameters(10.0, 0.8, 0.1, 100.0);
0462   hrmsPF->Fit("fit3", "", "", 100.0, 900.0);
0463   func4->SetLineColor(6); 
0464   func4->SetParNames("Noise", "Stochastic", "Constant", "Offset");
0465   func4->SetParameters(10.0, 0.8, 0.1, 100.0);
0466   hrmsCalo->Fit("fit4", "", "", 100.0, 900.0);
0467   (hsigmaPF->GetYaxis())->SetRangeUser( 0.0, 50.0);
0468   hsigmaPF->SetLineWidth(2); 
0469   hsigmaPF->SetLineColor(1); 
0470   hsigmaPF->Draw();
0471   hsigmaCalo->SetLineWidth(2);
0472   hsigmaCalo->SetLineColor(2);
0473   hsigmaCalo->Draw("SAME");
0474   hrmsPF->SetLineWidth(2);
0475   hrmsPF->SetLineColor(4);
0476   hrmsPF->Draw("SAME");  
0477   hrmsCalo->SetLineWidth(2);
0478   hrmsCalo->SetLineColor(6);
0479   hrmsCalo->Draw("SAME");
0480   */
0481 
0482   // Make the SET response versus SET plot
0483   /*
0484   TCanvas* canvas_SetRespVsTrueSet = new TCanvas("SetRespVsTrueSet", "SET Response vs True SET", 500,500);
0485   profileSETvsSETresp->SetStats(0); 
0486   profileSETvsSETresp->SetStats(0); 
0487   (profileSETvsSETresp->GetYaxis())->SetRangeUser(-1.0, 1.0);
0488   profileSETvsSETresp->SetLineWidth(2); 
0489   profileSETvsSETresp->SetLineColor(4); 
0490   profileSETvsSETresp->Draw();
0491   profileCaloSETvsCaloSETresp->SetLineWidth(2); 
0492   profileCaloSETvsCaloSETresp->SetLineColor(2); 
0493   profileCaloSETvsCaloSETresp->Draw("SAME");
0494   */
0495 
0496   // Make the MET response versus MET plot
0497   /*
0498   TCanvas* canvas_MetRespVsTrueMet = new TCanvas("MetRespVsTrueMet", "MET Response vs True MET", 500,500);
0499   profileMETvsMETresp->SetStats(0); 
0500   profileMETvsMETresp->SetStats(0); 
0501   (profileMETvsMETresp->GetYaxis())->SetRangeUser(-1.0, 1.0);
0502   profileMETvsMETresp->SetLineWidth(2); 
0503   profileMETvsMETresp->SetLineColor(4); 
0504   profileMETvsMETresp->Draw();
0505   profileCaloMETvsCaloMETresp->SetLineWidth(2); 
0506   profileCaloMETvsCaloMETresp->SetLineColor(2); 
0507   profileCaloMETvsCaloMETresp->Draw("SAME");
0508   */
0509 
0510   //print the resulting plots to file
0511   /*
0512   canvas_MetResVsRecoSet->Print("MetResVsRecoSet.ps");
0513   canvas_SetRespVsTrueSet->Print("SetRespVsTrueSet.ps");
0514   canvas_MetRespVsTrueMet->Print("MetRespVsTrueMet.ps");  
0515   */
0516 }
0517 
0518 //void PFMETBenchmark::FitSlicesInY(TH2F* h, TH1F* mean, TH1F* sigma, bool doGausFit, int type )
0519 //{
0520 //  TAxis *fXaxis = h->GetXaxis();
0521 //  TAxis *fYaxis = h->GetYaxis();
0522 //  Int_t nbins  = fXaxis->GetNbins();
0523 //
0524 //  //std::cout << "nbins = " << nbins << std::endl;
0525 //
0526 //  Int_t binmin = 1;
0527 //  Int_t binmax = nbins;
0528 //  TString option = "QNR";
0529 //  TString opt = option;
0530 //  opt.ToLower();
0531 //  Float_t ngroup = 1;
0532 //  ngroup = 1;
0533 //
0534 //  //std::cout << "fYaxis->GetXmin() = " << fYaxis->GetXmin() << std::endl;
0535 //  //std::cout << "fYaxis->GetXmax() = " << fYaxis->GetXmax() << std::endl;
0536 //
0537 //  //default is to fit with a gaussian
0538 //  TF1 *f1 = 0;
0539 //  if (f1 == 0)
0540 //    {
0541 //      //f1 = (TF1*)gROOT->GetFunction("gaus");
0542 //      if (f1 == 0) f1 = new TF1("gaus","gaus", fYaxis->GetXmin(), fYaxis->GetXmax());
0543 //      else         f1->SetRange( fYaxis->GetXmin(), fYaxis->GetXmax());
0544 //    }
0545 //  Int_t npar = f1->GetNpar();
0546 //
0547 //  //std::cout << "npar = " << npar << std::endl;
0548 //
0549 //  if (npar <= 0) return;
0550 //
0551 //  Double_t *parsave = new Double_t[npar];
0552 //  f1->GetParameters(parsave);
0553 //
0554 ////  //Create one histogram for each function parameter
0555 //  Int_t ipar;
0556 //  TH1F **hlist = new TH1F*[npar];
0557 //  char *name   = new char[2000];
0558 //  char *title  = new char[2000];
0559 //  const TArrayD *bins = fXaxis->GetXbins();
0560 //  for( ipar=0; ipar < npar; ipar++ )
0561 //    {
0562 //      if( ipar == 1 )
0563 //  if( type == 1 )   sprintf(name,"meanPF");
0564 //  else              sprintf(name,"meanCalo");
0565 //      else
0566 //  if( doGausFit )
0567 //    if( type == 1 ) sprintf(name,"sigmaPF");
0568 //    else            sprintf(name,"sigmaCalo");
0569 //  else
0570 //    if( type == 1 ) sprintf(name,"rmsPF");
0571 //    else            sprintf(name,"rmsCalo");
0572 //      if( type == 1 )     sprintf(title,"Particle Flow");
0573 //      else                sprintf(title,"Calorimeter");
0574 //      delete gDirectory->FindObject(name);
0575 //      if (bins->fN == 0)
0576 //  hlist[ipar] = new TH1F(name,title, nbins, fXaxis->GetXmin(), fXaxis->GetXmax());
0577 //      else
0578 //  hlist[ipar] = new TH1F(name,title, nbins,bins->fArray);
0579 //      hlist[ipar]->GetXaxis()->SetTitle(fXaxis->GetTitle());
0580 //    }
0581 //  sprintf(name,"test_chi2");
0582 //  delete gDirectory->FindObject(name);
0583 //  TH1F *hchi2 = new TH1F(name,"chisquare", nbins, fXaxis->GetXmin(), fXaxis->GetXmax());
0584 //  hchi2->GetXaxis()->SetTitle(fXaxis->GetTitle());
0585 //
0586 //  //Loop on all bins in X, generate a projection along Y
0587 //  Int_t bin;
0588 //  Int_t nentries;
0589 //  for( bin = (Int_t) binmin; bin <= (Int_t) binmax; bin += (Int_t)ngroup )
0590 //    {
0591 //      TH1F *hpy = (TH1F*) h->ProjectionY("_temp", (Int_t) bin, (Int_t) (bin + ngroup - 1), "e");
0592 //      if(hpy == 0) continue;
0593 //      nentries = Int_t( hpy->GetEntries() );
0594 //      if(nentries == 0 ) {delete hpy; continue;}
0595 //      f1->SetParameters(parsave);
0596 //      hpy->Fit( f1, opt.Data() );
0597 //      Int_t npfits = f1->GetNumberFitPoints();
0598 //      //cout << "bin = " << bin << "; Npfits = " << npfits << "; npar = " << npar << endl;
0599 //      if( npfits > npar )
0600 //  {
0601 //    Int_t biny = bin + (Int_t)ngroup/2;
0602 //    for( ipar=0; ipar < npar; ipar++ )
0603 //      {
0604 //        if( doGausFit ) hlist[ipar]->Fill( fXaxis->GetBinCenter(biny), f1->GetParameter(ipar) );
0605 //        else            hlist[ipar]->Fill( fXaxis->GetBinCenter(biny), hpy->GetRMS() );
0606 //        //cout << "bin[" << bin << "]: RMS = " << hpy->GetRMS() << "; sigma = " << f1->GetParameter(ipar) << endl;
0607 //        hlist[ipar]->SetBinError( biny, f1->GetParError(ipar) );
0608 //      }
0609 //    hchi2->Fill( fXaxis->GetBinCenter(biny), f1->GetChisquare()/(npfits-npar) );
0610 //  }
0611 //      delete hpy;
0612 //      ngroup += ngroup*0.2;//0.1  //used for non-uniform binning
0613 //    }
0614 //  *mean = *hlist[1];
0615 //  *sigma = *hlist[2];
0616 //  cout << "Entries = " << hlist[0]->GetEntries() << endl;
0617 //}
0618 
0619 double PFMETBenchmark::mpi_pi(double angle) {
0620   const double pi = 3.14159265358979323;
0621   const double pi2 = pi * 2.;
0622   while (angle > pi)
0623     angle -= pi2;
0624   while (angle < -pi)
0625     angle += pi2;
0626   return angle;
0627 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!