Changeset 2264876 in git
- Timestamp:
- Aug 25, 2016, 4:27:14 PM (8 years ago)
- Branches:
- ImprovedOutputFile, Timing, dual_readout, llp, master
- Children:
- bc58cf5
- Parents:
- 7993cad (diff), 1408174 (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. - Files:
-
- 1 deleted
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
Makefile
r7993cad r2264876 137 137 tmp/examples/Example1.$(ObjSuf): \ 138 138 examples/Example1.cpp \ 139 classes/DelphesClasses.h \ 140 external/ExRootAnalysis/ExRootTreeReader.h \ 141 external/ExRootAnalysis/ExRootTreeWriter.h \ 142 external/ExRootAnalysis/ExRootTreeBranch.h \ 143 external/ExRootAnalysis/ExRootResult.h \ 144 external/ExRootAnalysis/ExRootUtilities.h 145 Validation$(ExeSuf): \ 146 tmp/examples/Validation.$(ObjSuf) 147 148 tmp/examples/Validation.$(ObjSuf): \ 149 examples/Validation.cpp \ 139 150 classes/DelphesClasses.h \ 140 151 external/ExRootAnalysis/ExRootTreeReader.h \ … … 150 161 root2pileup$(ExeSuf) \ 151 162 stdhep2pileup$(ExeSuf) \ 152 Example1$(ExeSuf) 163 Example1$(ExeSuf) \ 164 Validation$(ExeSuf) 153 165 154 166 EXECUTABLE_OBJ += \ … … 159 171 tmp/converters/root2pileup.$(ObjSuf) \ 160 172 tmp/converters/stdhep2pileup.$(ObjSuf) \ 161 tmp/examples/Example1.$(ObjSuf) 173 tmp/examples/Example1.$(ObjSuf) \ 174 tmp/examples/Validation.$(ObjSuf) 162 175 163 176 DelphesHepMC$(ExeSuf): \ … … 183 196 classes/DelphesLHEFReader.h \ 184 197 external/ExRootAnalysis/ExRootTreeWriter.h \ 198 external/ExRootAnalysis/ExRootTreeBranch.h \ 199 external/ExRootAnalysis/ExRootProgressBar.h 200 DelphesROOT$(ExeSuf): \ 201 tmp/readers/DelphesROOT.$(ObjSuf) 202 203 tmp/readers/DelphesROOT.$(ObjSuf): \ 204 readers/DelphesROOT.cpp \ 205 modules/Delphes.h \ 206 classes/DelphesStream.h \ 207 classes/DelphesClasses.h \ 208 classes/DelphesFactory.h \ 209 external/ExRootAnalysis/ExRootTreeWriter.h \ 210 external/ExRootAnalysis/ExRootTreeReader.h \ 185 211 external/ExRootAnalysis/ExRootTreeBranch.h \ 186 212 external/ExRootAnalysis/ExRootProgressBar.h … … 200 226 DelphesHepMC$(ExeSuf) \ 201 227 DelphesLHEF$(ExeSuf) \ 228 DelphesROOT$(ExeSuf) \ 202 229 DelphesSTDHEP$(ExeSuf) 203 230 … … 205 232 tmp/readers/DelphesHepMC.$(ObjSuf) \ 206 233 tmp/readers/DelphesLHEF.$(ObjSuf) \ 234 tmp/readers/DelphesROOT.$(ObjSuf) \ 207 235 tmp/readers/DelphesSTDHEP.$(ObjSuf) 208 236 … … 344 372 modules/PdgCodeFilter.h \ 345 373 modules/BeamSpotFilter.h \ 374 modules/RecoPuFilter.h \ 346 375 modules/Cloner.h \ 347 376 modules/Weighter.h \ … … 351 380 modules/VertexSorter.h \ 352 381 modules/VertexFinder.h \ 382 modules/VertexFinderDA4D.h \ 353 383 modules/ExampleModule.h 354 384 tmp/modules/ModulesDict$(PcmSuf): \ … … 798 828 external/ExRootAnalysis/ExRootFilter.h \ 799 829 external/ExRootAnalysis/ExRootClassifier.h 830 tmp/modules/RecoPuFilter.$(ObjSuf): \ 831 modules/RecoPuFilter.$(SrcSuf) \ 832 modules/RecoPuFilter.h \ 833 classes/DelphesClasses.h \ 834 classes/DelphesFactory.h \ 835 classes/DelphesFormula.h \ 836 external/ExRootAnalysis/ExRootResult.h \ 837 external/ExRootAnalysis/ExRootFilter.h \ 838 external/ExRootAnalysis/ExRootClassifier.h 800 839 tmp/modules/SimpleCalorimeter.$(ObjSuf): \ 801 840 modules/SimpleCalorimeter.$(SrcSuf) \ … … 896 935 modules/VertexFinder.$(SrcSuf) \ 897 936 modules/VertexFinder.h \ 937 classes/DelphesClasses.h \ 938 classes/DelphesFactory.h \ 939 classes/DelphesFormula.h \ 940 classes/DelphesPileUpReader.h \ 941 external/ExRootAnalysis/ExRootResult.h \ 942 external/ExRootAnalysis/ExRootFilter.h \ 943 external/ExRootAnalysis/ExRootClassifier.h 944 tmp/modules/VertexFinderDA4D.$(ObjSuf): \ 945 modules/VertexFinderDA4D.$(SrcSuf) \ 946 modules/VertexFinderDA4D.h \ 898 947 classes/DelphesClasses.h \ 899 948 classes/DelphesFactory.h \ … … 998 1047 tmp/modules/PileUpJetID.$(ObjSuf) \ 999 1048 tmp/modules/PileUpMerger.$(ObjSuf) \ 1049 tmp/modules/RecoPuFilter.$(ObjSuf) \ 1000 1050 tmp/modules/SimpleCalorimeter.$(ObjSuf) \ 1001 1051 tmp/modules/StatusPidFilter.$(ObjSuf) \ … … 1010 1060 tmp/modules/UniqueObjectFinder.$(ObjSuf) \ 1011 1061 tmp/modules/VertexFinder.$(ObjSuf) \ 1062 tmp/modules/VertexFinderDA4D.$(ObjSuf) \ 1012 1063 tmp/modules/VertexSorter.$(ObjSuf) \ 1013 1064 tmp/modules/Weighter.$(ObjSuf) … … 1615 1666 tmp/external/tcl/tclVar.$(ObjSuf) 1616 1667 1668 modules/VertexFinderDA4D.h: \ 1669 classes/DelphesModule.h \ 1670 classes/DelphesClasses.h 1671 @touch $@ 1672 1617 1673 modules/TrackSmearing.h: \ 1618 1674 classes/DelphesModule.h … … 1959 2015 1960 2016 modules/BTagging.h: \ 2017 classes/DelphesModule.h 2018 @touch $@ 2019 2020 modules/RecoPuFilter.h: \ 1961 2021 classes/DelphesModule.h 1962 2022 @touch $@ -
examples/Validation.cpp
r7993cad r2264876 21 21 #include <utility> 22 22 #include <vector> 23 #include <typeinfo> 23 24 24 25 #include "TROOT.h" … … 28 29 #include "TString.h" 29 30 31 #include "TH1.h" 30 32 #include "TH2.h" 33 #include "TMath.h" 34 #include "TStyle.h" 35 #include "TGraph.h" 36 #include "TCanvas.h" 31 37 #include "THStack.h" 32 38 #include "TLegend.h" … … 34 40 #include "TClonesArray.h" 35 41 #include "TLorentzVector.h" 42 #include "TGraphErrors.h" 43 #include "TMultiGraph.h" 36 44 37 45 #include "classes/DelphesClasses.h" … … 47 55 //------------------------------------------------------------------------------ 48 56 49 // Here you can put your analysis macro 50 51 #include "Validation.C" 57 double ptrangemin = 10; 58 double ptrangemax = 10000; 59 static const int Nbins = 20; 60 61 int objStyle = 1; 62 int trackStyle = 7; 63 int towerStyle = 3; 64 65 Color_t objColor = kBlack; 66 Color_t trackColor = kBlack; 67 Color_t towerColor = kBlack; 68 69 double effLegXmin = 0.22; 70 double effLegXmax = 0.7; 71 double effLegYmin = 0.22; 72 double effLegYmax = 0.5; 73 74 double resLegXmin = 0.62; 75 double resLegXmax = 0.9; 76 double resLegYmin = 0.52; 77 double resLegYmax = 0.85; 78 79 double topLeftLegXmin = 0.22; 80 double topLeftLegXmax = 0.7; 81 double topLeftLegYmin = 0.52; 82 double topLeftLegYmax = 0.85; 83 84 85 struct resolPlot 86 { 87 TH1 *cenResolHist; 88 TH1 *fwdResolHist; 89 double ptmin; 90 double ptmax; 91 double xmin; 92 double xmax; 93 TString obj; 94 95 resolPlot(); 96 resolPlot(double ptdown, double ptup, TString object); 97 void set(double ptdown, double ptup, TString object, double xmin = 0, double xmax = 2); 98 void print(){std::cout << ptmin << std::endl;} 99 }; 100 101 102 resolPlot::resolPlot() 103 { 104 } 105 106 resolPlot::resolPlot(double ptdown, double ptup, TString object) 107 { 108 this->set(ptdown,ptup,object); 109 } 110 111 void resolPlot::set(double ptdown, double ptup, TString object, double xmin, double xmax) 112 { 113 ptmin = ptdown; 114 ptmax = ptup; 115 obj = object; 116 117 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", 200, xmin, xmax); 118 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", 200, xmin, xmax); 119 } 120 121 void HistogramsCollection(std::vector<resolPlot> *histos, double ptmin, double ptmax, TString obj, double xmin = 0, double xmax = 2) 122 { 123 double width; 124 double ptdown; 125 double ptup; 126 resolPlot ptemp; 127 128 for (int i = 0; i < Nbins; i++) 129 { 130 width = (ptmax - ptmin) / Nbins; 131 ptdown = TMath::Power(10,ptmin + i * width ); 132 ptup = TMath::Power(10,ptmin + (i+1) * width ); 133 ptemp.set(ptdown, ptup, obj, xmin, xmax); 134 histos->push_back(ptemp); 135 } 136 } 137 138 //------------------------------------------------------------------------------ 139 140 class ExRootResult; 141 class ExRootTreeReader; 142 143 //------------------------------------------------------------------------------ 144 145 void BinLogX(TH1*h) 146 { 147 TAxis *axis = h->GetXaxis(); 148 int bins = axis->GetNbins(); 149 150 Axis_t from = axis->GetXmin(); 151 Axis_t to = axis->GetXmax(); 152 Axis_t width = (to - from) / bins; 153 Axis_t *new_bins = new Axis_t[bins + 1]; 154 155 for (int i = 0; i <= bins; i++) 156 { 157 new_bins[i] = TMath::Power(10, from + i * width); 158 } 159 axis->Set(bins, new_bins); 160 delete new_bins; 161 } 162 163 164 //------------------------------------------------------------------------------ 165 166 template<typename T> 167 std::pair<TH1D*, TH1D*> GetEff(TClonesArray *branchReco, TClonesArray *branchParticle, TString name, int pdgID, ExRootTreeReader *treeReader) 168 { 169 170 cout << "** Computing Efficiency of reconstructing "<< branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl; 171 172 Long64_t allEntries = treeReader->GetEntries(); 173 174 GenParticle *particle; 175 T *recoObj; 176 177 TLorentzVector recoMomentum, genMomentum, bestRecoMomentum; 178 179 Float_t deltaR; 180 Float_t pt, eta; 181 Long64_t entry; 182 183 Int_t i, j; 184 185 TH1D *histGenPtcen = new TH1D(name+" gen spectra Pt",name+" gen spectra cen", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax)); 186 TH1D *histRecoPtcen = new TH1D(name+" reco spectra Pt",name+" reco spectra cen", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax)); 187 TH1D *histGenPtfwd = new TH1D(name+" gen spectra Eta",name+" gen spectra fwd", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax)); 188 TH1D *histRecoPtfwd = new TH1D(name+" reco spectra Eta",name+" reco spectra fwd", Nbins, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax)); 189 190 191 BinLogX(histGenPtcen); 192 BinLogX(histRecoPtcen); 193 BinLogX(histGenPtfwd); 194 BinLogX(histRecoPtfwd); 195 196 // Loop over all events 197 for(entry = 0; entry < allEntries; ++entry) 198 { 199 // Load selected branches with data from specified event 200 treeReader->ReadEntry(entry); 201 202 // Loop over all generated particle in event 203 for(i = 0; i < branchParticle->GetEntriesFast(); ++i) 204 { 205 206 particle = (GenParticle*) branchParticle->At(i); 207 genMomentum = particle->P4(); 208 209 deltaR = 999; 210 211 if (particle->PID == pdgID && genMomentum.Pt() > ptrangemin && genMomentum.Pt() < ptrangemax ) 212 { 213 214 // Loop over all reco object in event 215 for(j = 0; j < branchReco->GetEntriesFast(); ++j) 216 { 217 recoObj = (T*)branchReco->At(j); 218 recoMomentum = recoObj->P4(); 219 // this is simply to avoid warnings from initial state particle 220 // having infite rapidity ... 221 //if(Momentum.Px() == 0 && genMomentum.Py() == 0) continue; 222 223 // take the closest parton candidate 224 if(TMath::Abs(pdgID) == 5) 225 { 226 Jet *jet = (Jet *)recoObj; 227 if(jet->BTag != 1) continue; 228 } 229 if(TMath::Abs(pdgID) == 15) 230 { 231 Jet *jet = (Jet *)recoObj; 232 if(jet->TauTag != 1) continue; 233 } 234 if(genMomentum.DeltaR(recoMomentum) < deltaR) 235 { 236 deltaR = genMomentum.DeltaR(recoMomentum); 237 bestRecoMomentum = recoMomentum; 238 } 239 } 240 241 pt = genMomentum.Pt(); 242 eta = genMomentum.Eta(); 243 244 if (TMath::Abs(eta) < 1.5) 245 { 246 histGenPtcen->Fill(pt); 247 if(deltaR < 0.3) { histRecoPtcen->Fill(pt); } 248 } 249 else if (TMath::Abs(eta) < 2.5) 250 { 251 histGenPtfwd->Fill(pt); 252 if(deltaR < 0.3) { histRecoPtfwd->Fill(pt); } 253 254 } 255 } 256 } 257 } 258 259 260 std::pair<TH1D*,TH1D*> histos; 261 262 histRecoPtcen->Divide(histGenPtcen); 263 histRecoPtfwd->Divide(histGenPtfwd); 264 265 histos.first = histRecoPtcen; 266 histos.second = histRecoPtfwd; 267 268 return histos; 269 } 270 271 template<typename T> 272 void GetEres(std::vector<resolPlot> *histos, TClonesArray *branchReco, TClonesArray *branchParticle, int pdgID, ExRootTreeReader *treeReader) 273 { 274 Long64_t allEntries = treeReader->GetEntries(); 275 276 cout << "** Computing resolution of " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl; 277 278 GenParticle *particle; 279 T* recoObj; 280 281 TLorentzVector recoMomentum, genMomentum, bestGenMomentum; 282 283 Float_t deltaR; 284 Float_t pt, eta; 285 Long64_t entry; 286 287 Int_t i, j, bin; 288 289 // Loop over all events 290 for(entry = 0; entry < allEntries; ++entry) 291 { 292 // Load selected branches with data from specified event 293 treeReader->ReadEntry(entry); 294 295 // Loop over all reconstructed jets in event 296 for(i = 0; i < branchReco->GetEntriesFast(); ++i) 297 { 298 recoObj = (T*) branchReco->At(i); 299 recoMomentum = recoObj->P4(); 300 301 deltaR = 999; 302 303 // Loop over all hard partons in event 304 for(j = 0; j < branchParticle->GetEntriesFast(); ++j) 305 { 306 particle = (GenParticle*) branchParticle->At(j); 307 if (particle->PID == pdgID && particle->Status == 1) 308 { 309 genMomentum = particle->P4(); 310 311 // this is simply to avoid warnings from initial state particle 312 // having infite rapidity ... 313 if(genMomentum.Px() == 0 && genMomentum.Py() == 0) continue; 314 315 // take the closest parton candidate 316 if(genMomentum.DeltaR(recoMomentum) < deltaR) 317 { 318 deltaR = genMomentum.DeltaR(recoMomentum); 319 bestGenMomentum = genMomentum; 320 } 321 } 322 } 323 324 if(deltaR < 0.3) 325 { 326 pt = bestGenMomentum.Pt(); 327 eta = TMath::Abs(bestGenMomentum.Eta()); 328 329 for (bin = 0; bin < Nbins; bin++) 330 { 331 if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta > 0.0 && eta < 2.5) 332 { 333 if (eta < 1.5) {histos->at(bin).cenResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());} 334 else if (eta < 2.5) {histos->at(bin).fwdResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());} 335 } 336 } 337 } 338 } 339 } 340 } 341 void GetJetsEres(std::vector<resolPlot> *histos, TClonesArray *branchJet, TClonesArray *branchGenJet, ExRootTreeReader *treeReader) 342 { 343 344 Long64_t allEntries = treeReader->GetEntries(); 345 346 cout << "** Computing resolution of " << branchJet->GetName() << " induced by " << branchGenJet->GetName() << endl; 347 348 Jet *jet, *genjet; 349 350 TLorentzVector jetMomentum, genJetMomentum, bestGenJetMomentum; 351 352 Float_t deltaR; 353 Float_t pt, eta; 354 Long64_t entry; 355 356 Int_t i, j, bin; 357 358 // Loop over all events 359 for(entry = 0; entry < allEntries; ++entry) 360 { 361 // Load selected branches with data from specified event 362 treeReader->ReadEntry(entry); 363 364 if(entry%10000 == 0) cout << "Event number: "<< entry <<endl; 365 366 // Loop over all reconstructed jets in event 367 for(i = 0; i < TMath::Min(2,branchJet->GetEntriesFast()); ++i) //branchJet->GetEntriesFast(); ++i) 368 { 369 370 jet = (Jet*) branchJet->At(i); 371 jetMomentum = jet->P4(); 372 373 deltaR = 999; 374 375 // Loop over all hard partons in event 376 for(j = 0; j < TMath::Min(2,branchGenJet->GetEntriesFast()); ++j) 377 { 378 genjet = (Jet*) branchGenJet->At(j); 379 380 genJetMomentum = genjet->P4(); 381 382 // this is simply to avoid warnings from initial state particle 383 // having infite rapidity ... 384 if(genJetMomentum.Px() == 0 && genJetMomentum.Py() == 0) continue; 385 386 // take the closest parton candidate 387 if(genJetMomentum.DeltaR(jetMomentum) < deltaR) 388 { 389 deltaR = genJetMomentum.DeltaR(jetMomentum); 390 bestGenJetMomentum = genJetMomentum; 391 } 392 } 393 394 if(deltaR < 0.25) 395 { 396 pt = genJetMomentum.Pt(); 397 eta = TMath::Abs(genJetMomentum.Eta()); 398 399 for (bin = 0; bin < Nbins; bin++) 400 { 401 if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 1.5) 402 { 403 histos->at(bin).cenResolHist->Fill(jetMomentum.Pt()/bestGenJetMomentum.Pt()); 404 } 405 else if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 2.5) 406 { 407 histos->at(bin).fwdResolHist->Fill(jetMomentum.Pt()/bestGenJetMomentum.Pt()); 408 } 409 410 } 411 } 412 } 413 } 414 } 415 416 void GetMetres(std::vector<resolPlot> *histos, TClonesArray *branchScalarHT, TClonesArray *branchMet, TClonesArray *branchJet, ExRootTreeReader *treeReader) 417 { 418 419 Long64_t allEntries = treeReader->GetEntries(); 420 421 cout << "** Computing resolution of " << branchMet->GetName() << " vs " << branchScalarHT->GetName() << endl; 422 423 MissingET *met; 424 ScalarHT *scalarHT; 425 426 Long64_t entry; 427 428 Int_t bin; 429 Double_t ht; 430 431 Jet *jet; 432 TLorentzVector p1, p2; 433 434 // Loop over all events 435 for(entry = 0; entry < allEntries; ++entry) 436 { 437 // Load selected branches with data from specified event 438 treeReader->ReadEntry(entry); 439 440 if(entry%10000 == 0) cout << "Event number: "<< entry <<endl; 441 442 if (branchJet->GetEntriesFast() > 1) 443 { 444 445 jet = (Jet*) branchJet->At(0); 446 p1 = jet->P4(); 447 jet = (Jet*) branchJet->At(1); 448 p2 = jet->P4(); 449 450 met = (MissingET*) branchMet->At(0); 451 scalarHT = (ScalarHT*) branchScalarHT->At(0); 452 ht = scalarHT->HT; 453 454 if(p1.Pt() < 0.75*ht/2) continue; 455 if(p2.Pt() < 0.75*ht/2) continue; 456 457 for (bin = 0; bin < Nbins; bin++) 458 { 459 if(ht > histos->at(bin).ptmin && ht < histos->at(bin).ptmax ) 460 { 461 histos->at(bin).cenResolHist->Fill(met->P4().Px()); 462 histos->at(bin).fwdResolHist->Fill(met->P4().Py()); 463 } 464 } 465 } 466 } 467 } 468 469 470 std::pair<Double_t, Double_t> GausFit(TH1* hist) 471 { 472 TF1 *f1 = new TF1("f1", "gaus", hist->GetMean()-2*hist->GetRMS(), hist->GetMean()+2*hist->GetRMS()); 473 hist->Fit("f1","RQ"); 474 Double_t sig = f1->GetParameter(2); 475 Double_t sigErr = f1->GetParError(2); 476 delete f1; 477 return make_pair (sig, sigErr); 478 } 479 480 481 TGraphErrors EresGraph(std::vector<resolPlot> *histos, bool central, bool rms = false) 482 { 483 Int_t bin; 484 Int_t count = 0; 485 TGraphErrors gr = TGraphErrors(Nbins/2); 486 Double_t sig = 0; 487 Double_t sigErr = 0; 488 for (bin = 0; bin < Nbins; bin++) 489 { 490 if (central == true && histos->at(bin).cenResolHist->GetEntries() > 100) 491 { 492 std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).cenResolHist); 493 if (rms == true) 494 { 495 gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.second); 496 gr.SetPointError(count,0, sigvalues.second); // to correct 497 } 498 else 499 { 500 gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first); 501 gr.SetPointError(count,0, sigvalues.second); 502 } 503 count++; 504 } 505 506 else if (central == false && histos->at(bin).fwdResolHist->GetEntries() > 10) 507 { 508 std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).fwdResolHist); 509 if (rms == true) 510 { 511 gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.second); 512 gr.SetPointError(count,0, sigvalues.second); // to correct 513 } 514 else 515 { 516 gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first); 517 gr.SetPointError(count,0, sigvalues.second); 518 } 519 count++; 520 } 521 522 } 523 return gr; 524 } 525 526 527 //------------------------------------------------------------------------------ 528 529 530 // type 1 : object, 2 : track, 3 : tower 531 532 void addGraph(TMultiGraph *mg, TGraphErrors *gr, TLegend *leg, int type) 533 { 534 535 gr->SetLineWidth(2); 536 537 switch ( type ) 538 { 539 case 1: 540 gr->SetLineColor(objColor); 541 gr->SetLineStyle(objStyle); 542 std::cout << "Adding " << gr->GetName() << std::endl; 543 mg->Add(gr); 544 leg->AddEntry(gr,"Reco","l"); 545 break; 546 547 case 2: 548 gr->SetLineColor(trackColor); 549 gr->SetLineStyle(trackStyle); 550 mg->Add(gr); 551 leg->AddEntry(gr,"Track","l"); 552 break; 553 554 case 3: 555 gr->SetLineColor(towerColor); 556 gr->SetLineStyle(towerStyle); 557 mg->Add(gr); 558 leg->AddEntry(gr,"Tower","l"); 559 break; 560 561 case 0: 562 gr->SetLineColor(objColor); 563 gr->SetLineStyle(objStyle); 564 mg->Add(gr); 565 break; 566 567 default: 568 std::cout << "wrong type, possibles choices are Object, Track and Tower" << std::endl; 569 break; 570 } 571 } 572 573 void addHist(TH1D *h, TLegend *leg, int type) 574 { 575 h->SetLineWidth(2); 576 577 switch ( type ) 578 { 579 case 1: 580 h->SetLineColor(objColor); 581 h->SetLineStyle(objStyle); 582 leg->AddEntry(h,"Reco","l"); 583 break; 584 585 case 2: 586 h->SetLineColor(trackColor); 587 h->SetLineStyle(trackStyle); 588 leg->AddEntry(h,"Track","l"); 589 break; 590 591 case 3: 592 h->SetLineColor(towerColor); 593 h->SetLineStyle(towerStyle); 594 leg->AddEntry(h,"Tower","l"); 595 break; 596 597 case 0: 598 h->SetLineColor(objColor); 599 h->SetLineStyle(objStyle); 600 break; 601 602 default: 603 std::cout << "wrong type, possibles choices are Object, Track and Tower" << std::endl; 604 break; 605 } 606 } 607 608 void DrawAxis(TMultiGraph *mg, TLegend *leg, double max) 609 { 610 mg->SetMinimum(0.); 611 mg->SetMaximum(max); 612 mg->GetXaxis()->SetLimits(ptrangemin,ptrangemax); 613 mg->GetYaxis()->SetTitle("resolution"); 614 mg->GetXaxis()->SetTitle("p_{T} [GeV]"); 615 mg->GetYaxis()->SetTitleSize(0.07); 616 mg->GetXaxis()->SetTitleSize(0.07); 617 mg->GetYaxis()->SetLabelSize(0.06); 618 mg->GetXaxis()->SetLabelSize(0.06); 619 mg->GetYaxis()->SetLabelOffset(0.03); 620 mg->GetYaxis()->SetTitleOffset(1.4); 621 mg->GetXaxis()->SetTitleOffset(1.4); 622 623 mg->GetYaxis()->SetNdivisions(505); 624 625 leg->SetBorderSize(0); 626 leg->SetShadowColor(0); 627 leg->SetFillColor(0); 628 leg->SetFillStyle(0); 629 630 gStyle->SetOptTitle(0); 631 gPad->SetLogx(); 632 gPad->SetBottomMargin(0.2); 633 gPad->SetLeftMargin(0.2); 634 gPad->Modified(); 635 gPad->Update(); 636 637 } 638 639 void DrawAxis(TH1D *h, TLegend *leg, int type) 640 { 641 642 h->GetYaxis()->SetRangeUser(0,1.0); 643 if (type == 0) h->GetXaxis()->SetTitle("p_{T} [GeV]"); 644 else h->GetXaxis()->SetTitle("#eta"); 645 h->GetYaxis()->SetTitle("efficiency"); 646 h->GetYaxis()->SetTitleSize(0.07); 647 h->GetXaxis()->SetTitleSize(0.07); 648 h->GetYaxis()->SetLabelSize(0.06); 649 h->GetXaxis()->SetLabelSize(0.06); 650 h->GetYaxis()->SetLabelOffset(0.03); 651 h->GetYaxis()->SetTitleOffset(1.3); 652 h->GetXaxis()->SetTitleOffset(1.4); 653 654 h->GetYaxis()->SetNdivisions(505); 655 656 leg->SetBorderSize(0); 657 leg->SetShadowColor(0); 658 leg->SetFillColor(0); 659 leg->SetFillStyle(0); 660 661 gStyle->SetOptTitle(0); 662 gStyle->SetOptStat(0); 663 gPad->SetBottomMargin(0.2); 664 gPad->SetLeftMargin(0.2); 665 666 gPad->Modified(); 667 gPad->Update(); 668 669 } 670 671 672 void Validation(const char *inputFileElectron, const char *inputFileMuon, const char *inputFilePhoton, const char *inputFileJet, const char *inputFileBJet, const char *inputFileTauJet, const char *outputFile) 673 { 674 TChain *chainElectron = new TChain("Delphes"); 675 chainElectron->Add(inputFileElectron); 676 ExRootTreeReader *treeReaderElectron = new ExRootTreeReader(chainElectron); 677 678 TChain *chainMuon = new TChain("Delphes"); 679 chainMuon->Add(inputFileMuon); 680 ExRootTreeReader *treeReaderMuon = new ExRootTreeReader(chainMuon); 681 682 TChain *chainPhoton = new TChain("Delphes"); 683 chainPhoton->Add(inputFilePhoton); 684 ExRootTreeReader *treeReaderPhoton = new ExRootTreeReader(chainPhoton); 685 686 TChain *chainJet = new TChain("Delphes"); 687 chainJet->Add(inputFileJet); 688 ExRootTreeReader *treeReaderJet = new ExRootTreeReader(chainJet); 689 690 TChain *chainBJet = new TChain("Delphes"); 691 chainBJet->Add(inputFileBJet); 692 ExRootTreeReader *treeReaderBJet = new ExRootTreeReader(chainBJet); 693 694 TChain *chainTauJet = new TChain("Delphes"); 695 chainTauJet->Add(inputFileTauJet); 696 ExRootTreeReader *treeReaderTauJet = new ExRootTreeReader(chainTauJet); 697 698 TClonesArray *branchParticleElectron = treeReaderElectron->UseBranch("Particle"); 699 TClonesArray *branchTrackElectron = treeReaderElectron->UseBranch("Track"); 700 TClonesArray *branchTowerElectron = treeReaderElectron->UseBranch("Tower"); 701 TClonesArray *branchElectron = treeReaderElectron->UseBranch("Electron"); 702 703 TClonesArray *branchParticleMuon = treeReaderMuon->UseBranch("Particle"); 704 TClonesArray *branchTrackMuon = treeReaderMuon->UseBranch("Track"); 705 TClonesArray *branchMuon = treeReaderMuon->UseBranch("Muon"); 706 707 TClonesArray *branchParticlePhoton = treeReaderPhoton->UseBranch("Particle"); 708 TClonesArray *branchTowerPhoton = treeReaderPhoton->UseBranch("Tower"); 709 TClonesArray *branchPhoton = treeReaderPhoton->UseBranch("Photon"); 710 711 TClonesArray *branchGenJet = treeReaderJet->UseBranch("GenJet"); 712 TClonesArray *branchPFJet = treeReaderJet->UseBranch("Jet"); 713 TClonesArray *branchCaloJet = treeReaderJet->UseBranch("CaloJet"); 714 715 TClonesArray *branchParticleBJet = treeReaderBJet->UseBranch("Particle"); 716 TClonesArray *branchPFBJet = treeReaderBJet->UseBranch("Jet"); 717 718 TClonesArray *branchParticleTauJet = treeReaderTauJet->UseBranch("Particle"); 719 TClonesArray *branchPFTauJet = treeReaderTauJet->UseBranch("Jet"); 720 721 TClonesArray *branchScalarHT = treeReaderJet->UseBranch("ScalarHT"); 722 TClonesArray *branchMet = treeReaderJet->UseBranch("MissingET"); 723 724 /////////////// 725 // Electrons // 726 /////////////// 727 728 // Reconstruction efficiency 729 TString elecs = "Electron"; 730 int elID = 11; 731 std::pair<TH1D*,TH1D*> histos_el = GetEff<Electron>(branchElectron, branchParticleElectron, "Electron", elID, treeReaderElectron); 732 733 // tracking reconstruction efficiency 734 std::pair <TH1D*,TH1D*> histos_eltrack = GetEff<Track>(branchTrackElectron, branchParticleElectron, "electronTrack", elID, treeReaderElectron); 735 736 // Tower reconstruction efficiency 737 std::pair <TH1D*,TH1D*> histos_eltower = GetEff<Tower>(branchTowerElectron, branchParticleElectron, "electronTower", elID, treeReaderElectron); 738 739 // Electron Energy Resolution 740 std::vector<resolPlot> plots_el; 741 HistogramsCollection(&plots_el, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electrons"); 742 GetEres<Electron>(&plots_el, branchElectron, branchParticleElectron, elID, treeReaderElectron); 743 TGraphErrors gr_el = EresGraph(&plots_el, true); 744 TGraphErrors gr_elFwd = EresGraph(&plots_el, false); 745 gr_el.SetName("Electron"); 746 gr_elFwd.SetName("ElectronFwd"); 747 748 // Electron Track Energy Resolution 749 std::vector<resolPlot> plots_eltrack; 750 HistogramsCollection(&plots_eltrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electronsTracks"); 751 GetEres<Track>(&plots_eltrack, branchTrackElectron, branchParticleElectron, elID, treeReaderElectron); 752 TGraphErrors gr_eltrack = EresGraph(&plots_eltrack, true); 753 TGraphErrors gr_eltrackFwd = EresGraph(&plots_eltrack, false); 754 gr_eltrack.SetName("ElectronTracks"); 755 gr_eltrackFwd.SetName("ElectronTracksFwd"); 756 757 // Electron Tower Energy Resolution 758 std::vector<resolPlot> plots_eltower; 759 HistogramsCollection(&plots_eltower, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "electronsTower"); 760 GetEres<Tower>(&plots_eltower, branchTowerElectron, branchParticleElectron, elID, treeReaderElectron); 761 TGraphErrors gr_eltower = EresGraph(&plots_eltower, true); 762 TGraphErrors gr_eltowerFwd = EresGraph(&plots_eltower, false); 763 gr_eltower.SetName("ElectronTower"); 764 gr_eltrackFwd.SetName("ElectronTracksFwd"); 765 766 // Canvases 767 TString elEff = "electronEff"; 768 TCanvas *C_el1 = new TCanvas(elEff,elEff, 1600, 600); 769 C_el1->Divide(2); 770 C_el1->cd(1); 771 TLegend *leg_el1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax); 772 leg_el1->SetHeader("#splitline{electrons}{|#eta| < 1.5}"); 773 leg_el1->AddEntry("","",""); 774 775 gPad->SetLogx(); 776 histos_eltrack.first->Draw("]["); 777 addHist(histos_eltrack.first, leg_el1, 2); 778 histos_el.first->Draw("same ]["); 779 addHist(histos_el.first, leg_el1, 1); 780 DrawAxis(histos_eltrack.first, leg_el1, 0); 781 782 leg_el1->Draw(); 783 784 C_el1->cd(2); 785 TLegend *leg_el2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax); 786 leg_el2->SetHeader("#splitline{electrons}{1.5 < |#eta| < 2.5}"); 787 leg_el2->AddEntry("","",""); 788 789 gPad->SetLogx(); 790 histos_eltrack.second->Draw("]["); 791 addHist(histos_eltrack.second, leg_el2, 2); 792 histos_el.second->Draw("same ]["); 793 addHist(histos_el.second, leg_el2, 1); 794 795 DrawAxis(histos_eltrack.second, leg_el2, 0); 796 leg_el2->Draw(); 797 798 C_el1->cd(0); 799 800 TString elRes = "electronERes"; 801 TString elResFwd = "electronEResForward"; 802 TCanvas *C_el2 = new TCanvas(elRes,elRes, 1600, 600); 803 C_el2->Divide(2); 804 C_el2->cd(1); 805 TMultiGraph *mg_el = new TMultiGraph(elRes,elRes); 806 TLegend *leg_el = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax); 807 leg_el->SetHeader("#splitline{electrons}{|#eta| < 1.5}"); 808 leg_el->AddEntry("","",""); 809 810 addGraph(mg_el, &gr_eltower, leg_el, 3); 811 addGraph(mg_el, &gr_eltrack, leg_el, 2); 812 addGraph(mg_el, &gr_el, leg_el, 1); 813 814 mg_el->Draw("ACX"); 815 leg_el->Draw(); 816 817 DrawAxis(mg_el, leg_el, 0.1); 818 819 C_el2->cd(2); 820 TMultiGraph *mg_elFwd = new TMultiGraph(elResFwd,elResFwd); 821 TLegend *leg_elFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax); 822 leg_elFwd->SetHeader("#splitline{electrons}{1.5 < |#eta| < 2.5}"); 823 leg_elFwd->AddEntry("","",""); 824 825 addGraph(mg_elFwd, &gr_eltowerFwd, leg_elFwd, 3); 826 addGraph(mg_elFwd, &gr_eltrackFwd, leg_elFwd, 2); 827 addGraph(mg_elFwd, &gr_elFwd, leg_elFwd, 1); 828 829 mg_elFwd->Draw("ACX"); 830 leg_elFwd->Draw(); 831 832 DrawAxis(mg_elFwd, leg_elFwd, 0.2); 833 834 C_el1->Print("validation.pdf(","pdf"); 835 C_el2->Print("validation.pdf","pdf"); 836 837 gDirectory->cd(0); 838 839 /////////// 840 // Muons // 841 /////////// 842 843 // Reconstruction efficiency 844 int muID = 13; 845 std::pair<TH1D*,TH1D*> histos_mu = GetEff<Muon>(branchMuon, branchParticleMuon,"Muon", muID, treeReaderMuon); 846 847 // muon tracking reconstruction efficiency 848 std::pair <TH1D*,TH1D*> histos_mutrack = GetEff<Track>(branchTrackMuon, branchParticleMuon, "muonTrack", muID, treeReaderMuon); 849 850 // Muon Energy Resolution 851 std::vector<resolPlot> plots_mu; 852 HistogramsCollection(&plots_mu, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muons"); 853 GetEres<Muon>(&plots_mu, branchMuon, branchParticleMuon, muID, treeReaderMuon); 854 TGraphErrors gr_mu = EresGraph(&plots_mu, true); 855 TGraphErrors gr_muFwd = EresGraph(&plots_mu, false); 856 gr_mu.SetName("Muon"); 857 gr_muFwd.SetName("MuonFwd"); 858 859 // Muon Track Energy Resolution 860 std::vector<resolPlot> plots_mutrack; 861 HistogramsCollection(&plots_mutrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muonsTracks"); 862 GetEres<Track>(&plots_mutrack, branchTrackMuon, branchParticleMuon, muID, treeReaderMuon); 863 TGraphErrors gr_mutrack = EresGraph(&plots_mutrack, true); 864 TGraphErrors gr_mutrackFwd = EresGraph(&plots_mutrack, false); 865 gr_mutrackFwd.SetName("MuonTracksFwd"); 866 867 // Canvas 868 869 TString muEff = "muonEff"; 870 TCanvas *C_mu1 = new TCanvas(muEff,muEff, 1600, 600); 871 C_mu1->Divide(2); 872 C_mu1->cd(1); 873 TLegend *leg_mu1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax); 874 leg_mu1->SetHeader("#splitline{muons}{|#eta| < 1.5}"); 875 leg_mu1->AddEntry("","",""); 876 877 878 gPad->SetLogx(); 879 histos_mutrack.first->Draw("]["); 880 addHist(histos_mutrack.first, leg_mu1, 2); 881 histos_mu.first->Draw("same ]["); 882 addHist(histos_mu.first, leg_mu1, 1); 883 884 DrawAxis(histos_mutrack.first, leg_mu1, 0); 885 886 leg_mu1->Draw(); 887 888 C_mu1->cd(2); 889 TLegend *leg_mu2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax); 890 leg_mu2->SetHeader("#splitline{muons}{1.5 < |#eta| < 2.5}"); 891 leg_mu2->AddEntry("","",""); 892 893 gPad->SetLogx(); 894 histos_mutrack.second->Draw("]["); 895 addHist(histos_mutrack.second, leg_mu2, 2); 896 histos_mu.second->Draw("same ]["); 897 addHist(histos_mu.second, leg_mu2, 1); 898 899 DrawAxis(histos_mutrack.second, leg_mu2, 1); 900 leg_mu2->Draw(); 901 902 TString muRes = "muonERes"; 903 TString muResFwd = "muonEResFwd"; 904 905 TCanvas *C_mu = new TCanvas(muRes,muRes, 1600, 600); 906 C_mu->Divide(2); 907 C_mu->cd(1); 908 TMultiGraph *mg_mu = new TMultiGraph(muRes,muRes); 909 TLegend *leg_mu = new TLegend(topLeftLegXmin,topLeftLegYmin,topLeftLegXmax,topLeftLegYmax); 910 leg_mu->SetHeader("#splitline{muons}{|#eta| < 1.5}"); 911 leg_mu->AddEntry("","",""); 912 913 addGraph(mg_mu, &gr_mutrack, leg_mu, 2); 914 addGraph(mg_mu, &gr_mu, leg_mu, 1); 915 916 mg_mu->Draw("ACX"); 917 leg_mu->Draw(); 918 919 DrawAxis(mg_mu, leg_mu, 0.3); 920 921 C_mu->cd(2); 922 TMultiGraph *mg_muFwd = new TMultiGraph(muResFwd,muResFwd); 923 TLegend *leg_muFwd = new TLegend(topLeftLegXmin,topLeftLegYmin,topLeftLegXmax,topLeftLegYmax); 924 leg_muFwd->SetHeader("#splitline{muons}{1.5 < |#eta| < 2.5}"); 925 leg_muFwd->AddEntry("","",""); 926 927 addGraph(mg_muFwd, &gr_mutrackFwd, leg_muFwd, 2); 928 addGraph(mg_muFwd, &gr_muFwd, leg_muFwd, 1); 929 930 mg_muFwd->Draw("ACX"); 931 leg_muFwd->Draw(); 932 933 DrawAxis(mg_muFwd, leg_muFwd, 0.3); 934 935 //C_mu1->SaveAs(muEff+".eps"); 936 //C_mu->SaveAs(muRes+".eps"); 937 938 C_mu1->Print("validation.pdf","pdf"); 939 C_mu->Print("validation.pdf","pdf"); 940 941 gDirectory->cd(0); 942 943 ///////////// 944 // Photons // 945 ///////////// 946 947 // Reconstruction efficiency 948 int phID = 22; 949 std::pair<TH1D*,TH1D*> histos_ph = GetEff<Electron>(branchPhoton, branchParticlePhoton, "Photon", phID, treeReaderPhoton); 950 std::pair<TH1D*,TH1D*> histos_phtower = GetEff<Electron>(branchTowerPhoton, branchParticlePhoton, "Photon", phID, treeReaderPhoton); 951 952 // Photon Energy Resolution 953 std::vector<resolPlot> plots_ph; 954 HistogramsCollection(&plots_ph, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photons"); 955 GetEres<Photon>(&plots_ph, branchPhoton, branchParticlePhoton, phID, treeReaderPhoton); 956 TGraphErrors gr_ph = EresGraph(&plots_ph, true); 957 TGraphErrors gr_phFwd = EresGraph(&plots_ph, false); 958 gr_ph.SetName("Photon"); 959 gr_phFwd.SetName("PhotonFwd"); 960 961 962 // Photon Tower Energy Resolution 963 std::vector<resolPlot> plots_phtower; 964 HistogramsCollection(&plots_phtower, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "photonsTower"); 965 GetEres<Tower>(&plots_phtower, branchTowerPhoton, branchParticlePhoton, phID, treeReaderPhoton); 966 TGraphErrors gr_phtower = EresGraph(&plots_phtower, true); 967 TGraphErrors gr_phtowerFwd = EresGraph(&plots_phtower, false); 968 gr_phtower.SetName("PhotonTower"); 969 gr_phtowerFwd.SetName("PhotonTowerFwd"); 970 971 // Canvas 972 973 TString phEff = "photonEff"; 974 TCanvas *C_ph1 = new TCanvas(phEff,phEff, 1600, 600); 975 C_ph1->Divide(2); 976 C_ph1->cd(1); 977 TLegend *leg_ph1 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax); 978 leg_ph1->SetHeader("#splitline{photons}{|#eta| < 1.5}"); 979 leg_ph1->AddEntry("","",""); 980 981 982 gPad->SetLogx(); 983 histos_phtower.first->Draw("]["); 984 addHist(histos_phtower.first, leg_ph1, 3); 985 histos_ph.first->Draw("same ]["); 986 addHist(histos_ph.first, leg_ph1, 1); 987 988 DrawAxis(histos_phtower.first, leg_ph1, 0); 989 leg_ph1->Draw(); 990 991 C_ph1->cd(2); 992 TLegend *leg_ph2 = new TLegend(effLegXmin,effLegYmin,effLegXmax,effLegYmax); 993 leg_ph2->SetHeader("#splitline{photons}{1.5 < |#eta| < 2.5}"); 994 leg_ph2->AddEntry("","",""); 995 996 997 gPad->SetLogx(); 998 histos_phtower.second->Draw("]["); 999 addHist(histos_phtower.second, leg_ph2, 3); 1000 histos_ph.second->Draw("same ]["); 1001 addHist(histos_ph.second, leg_ph2, 1); 1002 1003 DrawAxis(histos_phtower.second, leg_ph2, 1); 1004 leg_ph2->Draw(); 1005 1006 C_ph1->SaveAs(phEff+".eps"); 1007 1008 TString phRes = "phERes"; 1009 TString phResFwd = "phEResFwd"; 1010 1011 TCanvas *C_ph = new TCanvas(phRes,phRes, 1600, 600); 1012 C_ph->Divide(2); 1013 C_ph->cd(1); 1014 TMultiGraph *mg_ph = new TMultiGraph(phRes,phRes); 1015 TLegend *leg_ph = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax); 1016 leg_ph->SetHeader("#splitline{photons}{|#eta| < 1.5}"); 1017 leg_ph->AddEntry("","",""); 1018 1019 addGraph(mg_ph, &gr_phtower, leg_ph, 3); 1020 addGraph(mg_ph, &gr_ph, leg_ph, 1); 1021 1022 mg_ph->Draw("ACX"); 1023 leg_ph->Draw(); 1024 1025 DrawAxis(mg_ph, leg_ph, 0.3); 1026 1027 C_ph->cd(2); 1028 TMultiGraph *mg_phFwd = new TMultiGraph(phResFwd,phResFwd); 1029 TLegend *leg_phFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax); 1030 leg_phFwd->SetHeader("#splitline{photons}{1.5 < |#eta| < 2.5}"); 1031 leg_phFwd->AddEntry("","",""); 1032 1033 addGraph(mg_phFwd, &gr_phtowerFwd, leg_phFwd, 3); 1034 addGraph(mg_phFwd, &gr_phFwd, leg_phFwd, 1); 1035 1036 mg_phFwd->Draw("ACX"); 1037 leg_phFwd->Draw(); 1038 1039 DrawAxis(mg_phFwd, leg_phFwd, 0.3); 1040 1041 C_ph->SaveAs(phRes+".eps"); 1042 1043 C_ph1->Print("validation.pdf","pdf"); 1044 C_ph->Print("validation.pdf","pdf"); 1045 1046 gDirectory->cd(0); 1047 1048 ////////// 1049 // Jets // 1050 ////////// 1051 1052 // BJets Reconstruction efficiency 1053 int bID = 5; 1054 std::pair<TH1D*,TH1D*> histos_btag = GetEff<Jet>(branchPFBJet, branchParticleBJet,"BTag", bID, treeReaderBJet); 1055 1056 // TauJets Reconstruction efficiency 1057 int tauID = 15; 1058 std::pair<TH1D*,TH1D*> histos_tautag = GetEff<Jet>(branchPFTauJet, branchParticleTauJet,"TauTag", tauID, treeReaderTauJet); 1059 1060 // PFJets Energy Resolution 1061 std::vector<resolPlot> plots_pfjets; 1062 HistogramsCollection(&plots_pfjets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "PFJet"); 1063 GetJetsEres(&plots_pfjets, branchPFJet, branchGenJet, treeReaderJet); 1064 TGraphErrors gr_pfjets = EresGraph(&plots_pfjets, true); 1065 TGraphErrors gr_pfjetsFwd = EresGraph(&plots_pfjets, false); 1066 gr_pfjets.SetName("pfJet"); 1067 gr_pfjetsFwd.SetName("pfJetFwd"); 1068 1069 // CaloJets Energy Resolution 1070 std::vector<resolPlot> plots_calojets; 1071 HistogramsCollection(&plots_calojets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "CaloJet"); 1072 GetJetsEres(&plots_calojets, branchCaloJet, branchGenJet, treeReaderJet); 1073 TGraphErrors gr_calojets = EresGraph(&plots_calojets, true); 1074 TGraphErrors gr_calojetsFwd = EresGraph(&plots_calojets, false); 1075 gr_calojets.SetName("caloJet"); 1076 gr_calojetsFwd.SetName("caloJetFwd"); 1077 1078 // MET Resolution vs HT 1079 std::vector<resolPlot> plots_met; 1080 HistogramsCollection(&plots_met, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "MET", -500, 500); 1081 GetMetres(&plots_met, branchScalarHT, branchMet, branchPFJet, treeReaderJet); 1082 TGraphErrors gr_met = EresGraph(&plots_met, true); 1083 gr_calojets.SetName("MET"); 1084 1085 // Canvas 1086 TString btagEff = "btagEff"; 1087 TCanvas *C_btag1 = new TCanvas(btagEff,btagEff, 1600, 600); 1088 C_btag1->Divide(2); 1089 C_btag1->cd(1); 1090 TLegend *leg_btag1 = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax); 1091 leg_btag1->SetHeader("#splitline{B-tagging}{|#eta| < 1.5}"); 1092 leg_btag1->AddEntry("","",""); 1093 1094 gPad->SetLogx(); 1095 histos_btag.first->Draw(); 1096 addHist(histos_btag.first, leg_btag1, 0); 1097 1098 DrawAxis(histos_btag.first, leg_btag1, 0); 1099 leg_btag1->Draw(); 1100 1101 C_btag1->cd(2); 1102 TLegend *leg_btag2 = new TLegend(resLegXmin,resLegYmin,resLegXmax+0.05,resLegYmax); 1103 leg_btag2->SetHeader("#splitline{B-tagging}{1.5 < |#eta| < 2.5}"); 1104 leg_btag2->AddEntry("","",""); 1105 1106 gPad->SetLogx(); 1107 histos_btag.second->Draw(); 1108 addHist(histos_btag.second, leg_btag2, 0); 1109 1110 DrawAxis(histos_btag.second, leg_btag2, 0); 1111 leg_btag2->Draw(); 1112 C_btag1->cd(0); 1113 1114 TString tautagEff = "tautagEff"; 1115 TCanvas *C_tautag1 = new TCanvas(tautagEff,tautagEff, 1600, 600); 1116 C_tautag1->Divide(2); 1117 C_tautag1->cd(1); 1118 TLegend *leg_tautag1 = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax); 1119 leg_tautag1->SetHeader("#splitline{#tau-tagging}{|#eta| < 1.5}"); 1120 leg_tautag1->AddEntry("","",""); 1121 1122 gPad->SetLogx(); 1123 histos_tautag.first->Draw(); 1124 addHist(histos_tautag.first, leg_tautag1, 0); 1125 1126 DrawAxis(histos_tautag.first, leg_tautag1, 0); 1127 leg_tautag1->Draw(); 1128 1129 C_tautag1->cd(2); 1130 TLegend *leg_tautag2 = new TLegend(resLegXmin,resLegYmin,resLegXmax+0.05,resLegYmax); 1131 leg_tautag2->SetHeader("#splitline{#tau-tagging}{1.5 < |#eta| < 2.5}"); 1132 leg_tautag2->AddEntry("","",""); 1133 1134 gPad->SetLogx(); 1135 histos_tautag.second->Draw(); 1136 addHist(histos_tautag.second, leg_tautag2, 0); 1137 1138 DrawAxis(histos_tautag.second, leg_tautag2, 0); 1139 leg_tautag2->Draw(); 1140 C_tautag1->cd(0); 1141 1142 TString jetRes = "jetERes"; 1143 TString jetResFwd = "jetEResFwd"; 1144 TCanvas *C_jet = new TCanvas(jetRes,jetRes, 1600, 600); 1145 C_jet->Divide(2); 1146 1147 C_jet->cd(1); 1148 TMultiGraph *mg_jet = new TMultiGraph(jetRes,jetRes); 1149 TLegend *leg_jet = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax); 1150 leg_jet->SetHeader("#splitline{jets}{|#eta| < 1.5}"); 1151 leg_jet->AddEntry("","",""); 1152 1153 addGraph(mg_jet, &gr_calojets, leg_jet, 3); 1154 addGraph(mg_jet, &gr_pfjets, leg_jet, 1); 1155 1156 mg_jet->Draw("ACX"); 1157 leg_jet->Draw(); 1158 1159 DrawAxis(mg_jet, leg_jet, 0.25); 1160 1161 C_jet->cd(2); 1162 TMultiGraph *mg_jetFwd = new TMultiGraph(jetResFwd,jetResFwd); 1163 TLegend *leg_jetFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax); 1164 leg_jetFwd->SetHeader("#splitline{jets}{|#eta| < 1.5}"); 1165 leg_jetFwd->AddEntry("","",""); 1166 1167 addGraph(mg_jetFwd, &gr_calojetsFwd, leg_jetFwd, 3); 1168 addGraph(mg_jetFwd, &gr_pfjetsFwd, leg_jetFwd, 1); 1169 1170 mg_jetFwd->Draw("ACX"); 1171 leg_jetFwd->Draw(); 1172 1173 DrawAxis(mg_jetFwd, leg_jetFwd, 0.25); 1174 1175 C_btag1->SaveAs(btagEff+".eps"); 1176 C_jet->SaveAs(jetRes+".eps"); 1177 1178 TString metRes = "MetRes"; 1179 TCanvas *C_met = new TCanvas(metRes,metRes, 800, 600); 1180 1181 TMultiGraph *mg_met = new TMultiGraph(metRes,metRes); 1182 TLegend *leg_met = new TLegend(topLeftLegXmin,topLeftLegYmin+0.2,topLeftLegXmax-0.2,topLeftLegYmax); 1183 leg_met->SetHeader("E_{T}^{miss}"); 1184 leg_met->AddEntry("","",""); 1185 1186 1187 addGraph(mg_met, &gr_met, leg_met, 0); 1188 1189 mg_met->Draw("ACX"); 1190 leg_met->Draw(); 1191 1192 DrawAxis(mg_met, leg_met, 300); 1193 1194 mg_met->GetXaxis()->SetTitle("H_{T} [GeV]"); 1195 mg_met->GetYaxis()->SetTitle("#sigma(ME_{x}) [GeV]"); 1196 1197 C_met->SaveAs(metRes+".eps"); 1198 1199 C_jet->Print("validation.pdf","pdf"); 1200 C_btag1->Print("validation.pdf","pdf"); 1201 C_tautag1->Print("validation.pdf","pdf"); 1202 C_met->Print("validation.pdf)","pdf"); 1203 1204 TFile *fout = new TFile(outputFile,"recreate"); 1205 1206 for (int bin = 0; bin < Nbins; bin++) 1207 { 1208 plots_el.at(bin).cenResolHist->Write(); 1209 plots_eltrack.at(bin).cenResolHist->Write(); 1210 plots_eltower.at(bin).cenResolHist->Write(); 1211 plots_el.at(bin).fwdResolHist->Write(); 1212 plots_eltrack.at(bin).fwdResolHist->Write(); 1213 plots_eltower.at(bin).fwdResolHist->Write(); 1214 } 1215 1216 histos_el.first->Write(); 1217 histos_el.second->Write(); 1218 histos_eltrack.first->Write(); 1219 histos_eltrack.second->Write(); 1220 histos_eltower.first->Write(); 1221 histos_eltower.second->Write(); 1222 C_el1->Write(); 1223 C_el2->Write(); 1224 1225 histos_mu.first->Write(); 1226 histos_mu.second->Write(); 1227 histos_mutrack.first->Write(); 1228 histos_mutrack.second->Write(); 1229 C_mu1->Write(); 1230 C_mu->Write(); 1231 1232 histos_ph.first->Write(); 1233 histos_ph.second->Write(); 1234 C_ph1->Write(); 1235 C_ph->Write(); 1236 1237 for (int bin = 0; bin < Nbins; bin++) 1238 { 1239 plots_pfjets.at(bin).cenResolHist->Write(); 1240 plots_pfjets.at(bin).fwdResolHist->Write(); 1241 plots_calojets.at(bin).cenResolHist->Write(); 1242 plots_calojets.at(bin).fwdResolHist->Write(); 1243 plots_met.at(bin).cenResolHist->Write(); 1244 } 1245 histos_btag.first->Write(); 1246 histos_btag.second->Write(); 1247 histos_tautag.first->Write(); 1248 histos_tautag.second->Write(); 1249 C_btag1->Write(); 1250 C_tautag1->Write(); 1251 C_jet->Write(); 1252 C_met->Write(); 1253 1254 fout->Write(); 1255 1256 cout << "** Exiting..." << endl; 1257 1258 delete treeReaderElectron; 1259 delete treeReaderMuon; 1260 delete treeReaderPhoton; 1261 delete treeReaderJet; 1262 delete treeReaderBJet; 1263 delete treeReaderTauJet; 1264 delete chainElectron; 1265 delete chainMuon; 1266 delete chainPhoton; 1267 delete chainJet; 1268 delete chainBJet; 1269 delete chainTauJet; 1270 } 52 1271 53 1272 //------------------------------------------------------------------------------ … … 59 1278 if(argc != 3) 60 1279 { 61 cout << " Usage: " << appName << " input_file output_file" << endl; 62 cout << " input_file - input file in ROOT format ('Delphes' tree)," << endl; 1280 cout << " Usage: " << appName << " input_file_electron input_file_muon input_file_photon input_file_jet input_file_bjet input_file_taujet output_file" << endl; 1281 cout << " input_file_electron - input file in ROOT format ('Delphes' tree)," << endl; 1282 cout << " input_file_muon - input file in ROOT format ('Delphes' tree)," << endl; 1283 cout << " input_file_photon - input file in ROOT format ('Delphes' tree)," << endl; 1284 cout << " input_file_jet - input file in ROOT format ('Delphes' tree)," << endl; 1285 cout << " input_file_bjet - input file in ROOT format ('Delphes' tree)," << endl; 1286 cout << " input_file_taujet - input file in ROOT format ('Delphes' tree)," << endl; 63 1287 cout << " output_file - output file in ROOT format" << endl; 64 1288 return 1; … … 71 1295 TApplication app(appName, &appargc, appargv); 72 1296 73 TString inputFile(argv[1]); 74 TString outputFile(argv[2]); 75 76 77 //------------------------------------------------------------------------------ 78 // Here you call your macro's main function 79 80 Validation(inputFile, outputFile); 81 82 //------------------------------------------------------------------------------ 83 84 } 85 86 1297 Validation(argv[1], argv[2], argv[3], argv[4], argv[5], argv[6], argv[7]); 1298 } 1299 1300 -
readers/DelphesROOT.cpp
r7993cad r2264876 167 167 modularDelphes->Clear(); 168 168 treeWriter->Clear(); 169 for(Int_t entry = 0; entry < numberOfEvents ; ++entry)169 for(Int_t entry = 0; entry < numberOfEvents && !interrupted; ++entry) 170 170 { 171 171 172 172 treeReader->ReadEntry(entry); 173 173
Note:
See TracChangeset
for help on using the changeset viewer.