- Timestamp:
- Aug 26, 2016, 5:27:30 PM (8 years ago)
- Branches:
- ImprovedOutputFile, Timing, dual_readout, llp, master
- Children:
- d091310
- Parents:
- eab6d46 (diff), caa953a (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the(diff)
links above to see all the changes relative to each parent. - File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
examples/Validation.cpp
reab6d46 r8f0b34c 343 343 } 344 344 } 345 346 347 template<typename T> 348 void GetPtres(std::vector<resolPlot> *histos, TClonesArray *branchReco, TClonesArray *branchParticle, int pdgID, ExRootTreeReader *treeReader) 349 { 350 Long64_t allEntries = treeReader->GetEntries(); 351 352 cout << "** Computing pt resolution of " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl; 353 354 GenParticle *particle; 355 T* recoObj; 356 357 TLorentzVector recoMomentum, genMomentum, bestGenMomentum; 358 359 Float_t deltaR; 360 Float_t pt, eta; 361 Long64_t entry; 362 363 Int_t i, j, bin; 364 365 // Loop over all events 366 for(entry = 0; entry < allEntries; ++entry) 367 { 368 // Load selected branches with data from specified event 369 treeReader->ReadEntry(entry); 370 371 // Loop over all reconstructed jets in event 372 for(i = 0; i < branchReco->GetEntriesFast(); ++i) 373 { 374 recoObj = (T*) branchReco->At(i); 375 recoMomentum = recoObj->P4(); 376 377 deltaR = 999; 378 379 // Loop over all hard partons in event 380 for(j = 0; j < branchParticle->GetEntriesFast(); ++j) 381 { 382 particle = (GenParticle*) branchParticle->At(j); 383 if (particle->PID == pdgID && particle->Status == 1) 384 { 385 genMomentum = particle->P4(); 386 387 // this is simply to avoid warnings from initial state particle 388 // having infite rapidity ... 389 if(genMomentum.Px() == 0 && genMomentum.Py() == 0) continue; 390 391 // take the closest parton candidate 392 if(genMomentum.DeltaR(recoMomentum) < deltaR) 393 { 394 deltaR = genMomentum.DeltaR(recoMomentum); 395 bestGenMomentum = genMomentum; 396 } 397 } 398 } 399 400 if(deltaR < 0.3) 401 { 402 pt = bestGenMomentum.Pt(); 403 eta = TMath::Abs(bestGenMomentum.Eta()); 404 405 for (bin = 0; bin < Nbins; bin++) 406 { 407 if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 2.5) 408 { 409 if (eta < 1.5) {histos->at(bin).cenResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());} 410 else if (eta < 2.5) {histos->at(bin).fwdResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());} 411 } 412 } 413 } 414 } 415 } 416 } 417 418 345 419 void GetJetsEres(std::vector<resolPlot> *histos, TClonesArray *branchJet, TClonesArray *branchGenJet, ExRootTreeReader *treeReader) 346 420 { … … 398 472 if(deltaR < 0.25) 399 473 { 400 pt = genJetMomentum. Pt();474 pt = genJetMomentum.E(); 401 475 eta = TMath::Abs(genJetMomentum.Eta()); 402 476 … … 405 479 if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 1.5) 406 480 { 407 histos->at(bin).cenResolHist->Fill(jetMomentum. Pt()/bestGenJetMomentum.Pt());481 histos->at(bin).cenResolHist->Fill(jetMomentum.E()/bestGenJetMomentum.E()); 408 482 } 409 483 else if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 2.5) 410 484 { 411 histos->at(bin).fwdResolHist->Fill(jetMomentum. Pt()/bestGenJetMomentum.Pt());485 histos->at(bin).fwdResolHist->Fill(jetMomentum.E()/bestGenJetMomentum.E()); 412 486 } 413 414 487 } 415 488 } … … 476 549 TF1 *f1 = new TF1("f1", "gaus", hist->GetMean()-2*hist->GetRMS(), hist->GetMean()+2*hist->GetRMS()); 477 550 hist->Fit("f1","RQ"); 478 Double_t sig = f1->GetParameter(2); 479 Double_t sigErr = f1->GetParError(2); 551 552 TF1 *f2 = new TF1("f2", "gaus", f1->GetParameter(1) - 2*f1->GetParameter(2), f1->GetParameter(1) + 2*f1->GetParameter(2)); 553 hist->Fit("f2","RQ"); 554 555 Double_t sig = f2->GetParameter(2); 556 Double_t sigErr = f2->GetParError(2); 557 480 558 delete f1; 559 delete f2; 481 560 return make_pair (sig, sigErr); 482 561 } … … 610 689 } 611 690 612 void DrawAxis(TMultiGraph *mg, TLegend *leg, double max )691 void DrawAxis(TMultiGraph *mg, TLegend *leg, double max, int type = 0) 613 692 { 614 693 mg->SetMinimum(0.); … … 616 695 mg->GetXaxis()->SetLimits(ptrangemin,ptrangemax); 617 696 mg->GetYaxis()->SetTitle("resolution"); 618 mg->GetXaxis()->SetTitle("p_{T} [GeV]"); 697 if (type == 0) mg->GetXaxis()->SetTitle("E [GeV]"); 698 else mg->GetXaxis()->SetTitle("p_{T} [GeV]"); 619 699 mg->GetYaxis()->SetTitleSize(0.07); 620 700 mg->GetXaxis()->SetTitleSize(0.07); … … 641 721 } 642 722 643 void DrawAxis(TH1D *h, TLegend *leg, int type )723 void DrawAxis(TH1D *h, TLegend *leg, int type = 0) 644 724 { 645 725 646 726 h->GetYaxis()->SetRangeUser(0,1.0); 647 if (type == 0) h->GetXaxis()->SetTitle(" p_{T}[GeV]");648 else h->GetXaxis()->SetTitle(" #eta");727 if (type == 0) h->GetXaxis()->SetTitle("E [GeV]"); 728 else h->GetXaxis()->SetTitle("p_{T} [GeV]"); 649 729 h->GetYaxis()->SetTitle("efficiency"); 650 730 h->GetYaxis()->SetTitleSize(0.07); … … 782 862 histos_el.first->Draw("same ]["); 783 863 addHist(histos_el.first, leg_el1, 1); 784 DrawAxis(histos_eltrack.first, leg_el1, 0);864 DrawAxis(histos_eltrack.first, leg_el1,1); 785 865 786 866 leg_el1->Draw(); … … 797 877 addHist(histos_el.second, leg_el2, 1); 798 878 799 DrawAxis(histos_eltrack.second, leg_el2, 0);879 DrawAxis(histos_eltrack.second, leg_el2, 1); 800 880 leg_el2->Draw(); 801 881 … … 850 930 std::pair <TH1D*,TH1D*> histos_mutrack = GetEff<Track>(branchTrackMuon, branchParticleMuon, "muonTrack", muID, treeReaderMuon); 851 931 852 // Muon EnergyResolution932 // Muon Pt Resolution 853 933 std::vector<resolPlot> plots_mu; 854 934 HistogramsCollection(&plots_mu, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muons"); 855 Get Eres<Muon>(&plots_mu, branchMuon, branchParticleMuon, muID, treeReaderMuon);935 GetPtres<Muon>(&plots_mu, branchMuon, branchParticleMuon, muID, treeReaderMuon); 856 936 TGraphErrors gr_mu = EresGraph(&plots_mu, true); 857 937 TGraphErrors gr_muFwd = EresGraph(&plots_mu, false); … … 862 942 std::vector<resolPlot> plots_mutrack; 863 943 HistogramsCollection(&plots_mutrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muonsTracks"); 864 Get Eres<Track>(&plots_mutrack, branchTrackMuon, branchParticleMuon, muID, treeReaderMuon);944 GetPtres<Track>(&plots_mutrack, branchTrackMuon, branchParticleMuon, muID, treeReaderMuon); 865 945 TGraphErrors gr_mutrack = EresGraph(&plots_mutrack, true); 866 946 TGraphErrors gr_mutrackFwd = EresGraph(&plots_mutrack, false); … … 884 964 addHist(histos_mu.first, leg_mu1, 1); 885 965 886 DrawAxis(histos_mutrack.first, leg_mu1, 0);966 DrawAxis(histos_mutrack.first, leg_mu1, 1); 887 967 888 968 leg_mu1->Draw(); … … 919 999 leg_mu->Draw(); 920 1000 921 DrawAxis(mg_mu, leg_mu, 0.3 );1001 DrawAxis(mg_mu, leg_mu, 0.3, 1); 922 1002 923 1003 C_mu->cd(2); … … 933 1013 leg_muFwd->Draw(); 934 1014 935 DrawAxis(mg_muFwd, leg_muFwd, 0.3 );1015 DrawAxis(mg_muFwd, leg_muFwd, 0.3, 1); 936 1016 937 1017 C_mu1->Print("delphes_validation.pdf","pdf"); … … 985 1065 addHist(histos_ph.first, leg_ph1, 1); 986 1066 987 DrawAxis(histos_phtower.first, leg_ph1, 0);1067 DrawAxis(histos_phtower.first, leg_ph1, 1); 988 1068 leg_ph1->Draw(); 989 1069 … … 1020 1100 leg_ph->Draw(); 1021 1101 1022 DrawAxis(mg_ph, leg_ph, 0. 3);1102 DrawAxis(mg_ph, leg_ph, 0.1); 1023 1103 1024 1104 C_ph->cd(2); … … 1034 1114 leg_phFwd->Draw(); 1035 1115 1036 DrawAxis(mg_phFwd, leg_phFwd, 0. 3);1116 DrawAxis(mg_phFwd, leg_phFwd, 0.1); 1037 1117 1038 1118 C_ph1->Print("delphes_validation.pdf","pdf"); … … 1091 1171 addHist(histos_btag.first, leg_btag1, 0); 1092 1172 1093 DrawAxis(histos_btag.first, leg_btag1, 0);1173 DrawAxis(histos_btag.first, leg_btag1, 1); 1094 1174 leg_btag1->Draw(); 1095 1175 … … 1103 1183 addHist(histos_btag.second, leg_btag2, 0); 1104 1184 1105 DrawAxis(histos_btag.second, leg_btag2, 0);1185 DrawAxis(histos_btag.second, leg_btag2, 1); 1106 1186 leg_btag2->Draw(); 1107 1187 C_btag1->cd(0); … … 1119 1199 addHist(histos_tautag.first, leg_tautag1, 0); 1120 1200 1121 DrawAxis(histos_tautag.first, leg_tautag1, 0);1201 DrawAxis(histos_tautag.first, leg_tautag1, 1); 1122 1202 leg_tautag1->Draw(); 1123 1203 … … 1131 1211 addHist(histos_tautag.second, leg_tautag2, 0); 1132 1212 1133 DrawAxis(histos_tautag.second, leg_tautag2, 0);1213 DrawAxis(histos_tautag.second, leg_tautag2, 1); 1134 1214 leg_tautag2->Draw(); 1135 1215 C_tautag1->cd(0); … … 1149 1229 addGraph(mg_jet, &gr_pfjets, leg_jet, 1); 1150 1230 1151 mg_jet->Draw("A CX");1231 mg_jet->Draw("ALX"); 1152 1232 leg_jet->Draw(); 1153 1233 1154 DrawAxis(mg_jet, leg_jet, 0. 25);1234 DrawAxis(mg_jet, leg_jet, 0.5); 1155 1235 1156 1236 C_jet->cd(2); 1157 1237 TMultiGraph *mg_jetFwd = new TMultiGraph(jetResFwd,jetResFwd); 1158 1238 TLegend *leg_jetFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax); 1159 leg_jetFwd->SetHeader("#splitline{jets}{ |#eta| < 1.5}");1239 leg_jetFwd->SetHeader("#splitline{jets}{1.5 < |#eta| < 2.5}"); 1160 1240 leg_jetFwd->AddEntry("","",""); 1161 1241 … … 1163 1243 addGraph(mg_jetFwd, &gr_pfjetsFwd, leg_jetFwd, 1); 1164 1244 1165 mg_jetFwd->Draw("A CX");1245 mg_jetFwd->Draw("ALX"); 1166 1246 leg_jetFwd->Draw(); 1167 1247 1168 DrawAxis(mg_jetFwd, leg_jetFwd, 0. 25);1248 DrawAxis(mg_jetFwd, leg_jetFwd, 0.5); 1169 1249 1170 1250 TString metRes = "MetRes"; … … 1212 1292 C_el1->Write(); 1213 1293 C_el2->Write(); 1294 1295 for (int bin = 0; bin < Nbins; bin++) 1296 { 1297 plots_mu.at(bin).cenResolHist->Write(); 1298 plots_mutrack.at(bin).cenResolHist->Write(); 1299 plots_mu.at(bin).fwdResolHist->Write(); 1300 plots_mutrack.at(bin).fwdResolHist->Write(); 1301 } 1214 1302 1215 1303 histos_mu.first->Write();
Note:
See TracChangeset
for help on using the changeset viewer.