- Timestamp:
- Aug 23, 2016, 2:37:13 PM (8 years ago)
- Branches:
- ImprovedOutputFile, Timing, dual_readout, llp, master
- Children:
- 20ca0cd, df315dc
- Parents:
- 9f8d4e7
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
examples/Validation.C
r9f8d4e7 r6563d4d 49 49 double ptmin; 50 50 double ptmax; 51 double xmin; 52 double xmax; 51 53 TString obj; 52 54 53 55 resolPlot(); 54 56 resolPlot(double ptdown, double ptup, TString object); 55 void set(double ptdown, double ptup, TString object );57 void set(double ptdown, double ptup, TString object, double xmin = 0, double xmax = 2); 56 58 void print(){std::cout << ptmin << std::endl;} 57 59 }; … … 67 69 } 68 70 69 void resolPlot::set(double ptdown, double ptup, TString object ){71 void resolPlot::set(double ptdown, double ptup, TString object, double xmin, double xmax){ 70 72 ptmin = ptdown; 71 73 ptmax = ptup; 72 74 obj = object; 73 75 74 cenResolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_cen", obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_cen", 500, 0, 2);75 fwdResolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_fwd", obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_fwd", 500, 0.4, 0.4);76 77 } 78 79 void HistogramsCollection(std::vector<resolPlot> *histos, double ptmin, double ptmax, TString obj )76 cenResolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_cen", obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_cen", 500, xmin, xmax); 77 fwdResolHist = new TH1D(obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_fwd", obj+"_delta_pt_"+Form("%4.2f",ptmin)+"_"+Form("%4.2f",ptmax)+"_fwd", 500, xmin, xmax); 78 79 } 80 81 void HistogramsCollection(std::vector<resolPlot> *histos, double ptmin, double ptmax, TString obj, double xmin = 0, double xmax = 2) 80 82 { 81 83 double width; … … 89 91 ptdown = TMath::Power(10,ptmin + i * width ); 90 92 ptup = TMath::Power(10,ptmin + (i+1) * width ); 91 ptemp.set(ptdown, ptup, obj );93 ptemp.set(ptdown, ptup, obj, xmin, xmax); 92 94 histos->push_back(ptemp); 93 95 } … … 362 364 } 363 365 } 366 367 void GetMetres(std::vector<resolPlot> *histos, TClonesArray *branchScalarHT, TClonesArray *branchMet, TClonesArray *branchJet, ExRootTreeReader *treeReader) 368 { 369 370 Long64_t allEntries = treeReader->GetEntries(); 371 372 cout << "** Computing resolution of " << branchMet->GetName() << " vs " << branchScalarHT->GetName() << endl; 373 374 MissingET *met; 375 ScalarHT *scalarHT; 376 377 Long64_t entry; 378 379 Int_t bin; 380 Double_t ht; 381 382 // Loop over all events 383 for(entry = 0; entry < allEntries; ++entry) 384 { 385 // Load selected branches with data from specified event 386 treeReader->ReadEntry(entry); 387 388 if(entry%10000 == 0) cout << "Event number: "<< entry <<endl; 389 390 if (branchJet->GetEntriesFast() > 1) 391 { 392 met = (MissingET*) branchMet->At(0); 393 scalarHT = (ScalarHT*) branchScalarHT->At(0); 394 ht = scalarHT->HT; 395 396 for (bin = 0; bin < Nbins; bin++) 397 { 398 if(ht > histos->at(bin).ptmin && ht < histos->at(bin).ptmax ) 399 { 400 histos->at(bin).cenResolHist->Fill(met->P4().Px()); 401 histos->at(bin).fwdResolHist->Fill(met->P4().Py()); 402 } 403 } 404 } 405 } 406 } 407 408 364 409 std::pair<Double_t, Double_t> GausFit(TH1* hist) 365 410 { … … 374 419 375 420 376 TGraphErrors EresGraph(std::vector<resolPlot> *histos, bool central )421 TGraphErrors EresGraph(std::vector<resolPlot> *histos, bool central, bool rms = false) 377 422 { 378 423 Int_t bin; … … 389 434 std::cout << " entries : " << histos->at(bin).cenResolHist->GetEntries() << std::endl; 390 435 std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).cenResolHist); 391 gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first); 392 gr.SetPointError(count,0, sigvalues.second); 436 if (rms == true) 437 { 438 gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.second); 439 gr.SetPointError(count,0, sigvalues.second); // to correct 440 } 441 else 442 { 443 gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first); 444 gr.SetPointError(count,0, sigvalues.second); 445 } 393 446 count++; 394 447 } … … 551 604 TClonesArray *branchPFJet = treeReader->UseBranch("Jet"); 552 605 TClonesArray *branchCaloJet = treeReader->UseBranch("CaloJet"); 606 TClonesArray *branchScalarHT = treeReader->UseBranch("ScalarHT"); 607 TClonesArray *branchMet = treeReader->UseBranch("MissingET"); 553 608 554 609 #ifdef ELECTRON … … 811 866 gr_pfjets.SetName("pfJet"); 812 867 813 814 // PFJets Energy Resolution 868 // CaloJets Energy Resolution 815 869 std::vector<resolPlot> plots_calojets; 816 870 HistogramsCollection(&plots_calojets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "CaloJet"); … … 819 873 gr_calojets.SetName("caloJet"); 820 874 875 // MET Resolution vs HT 876 std::vector<resolPlot> plots_met; 877 HistogramsCollection(&plots_met, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "MET", -200, 200); 878 GetMetres( &plots_met, branchScalarHT, branchMet, branchPFJet, treeReader); 879 TGraphErrors gr_met = EresGraph(&plots_met, true); 880 gr_calojets.SetName("MET"); 881 821 882 // Canvas 822 883 … … 882 943 C_btag1->SaveAs(btagEff+".eps"); 883 944 C_jet->SaveAs(jetRes+".eps"); 945 946 TString metRes = "MetRes"; 947 TCanvas *C_met = new TCanvas(metRes,metRes, 1000, 500); 948 949 TMultiGraph *mg_met = new TMultiGraph(metRes,metRes); 950 TLegend *leg_met = new TLegend(0.52,0.7,0.9,0.9); 951 952 addGraph(mg_met, &gr_met, leg_met, 3); 953 954 mg_met->Draw("ACX"); 955 leg_met->Draw(); 956 957 DrawAxis(mg_met, leg_met, 100); 958 959 C_met->SaveAs(metRes+".eps"); 960 884 961 885 962 #endif … … 934 1011 plots_pfjets.at(bin).cenResolHist->Write(); 935 1012 plots_calojets.at(bin).cenResolHist->Write(); 1013 plots_met.at(bin).cenResolHist->Write(); 936 1014 } 937 1015 histos_btag.first->Write(); … … 942 1020 C_tautag1->Write(); 943 1021 C_jet->Write(); 1022 C_met->Write(); 944 1023 #endif 945 1024
Note:
See TracChangeset
for help on using the changeset viewer.