Fork me on GitHub

Changeset 32677fb in git for examples/Validation.cpp


Ignore:
Timestamp:
Aug 26, 2016, 4:53:30 PM (8 years ago)
Author:
Alexandre Mertens <alexandre.mertens@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
d0f661b
Parents:
14fe596
Message:

style changes

File:
1 edited

Legend:

Unmodified
Added
Removed
  • examples/Validation.cpp

    r14fe596 r32677fb  
    343343  }
    344344}
     345
     346
     347template<typename T>
     348void 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
    345419void GetJetsEres(std::vector<resolPlot> *histos, TClonesArray *branchJet, TClonesArray *branchGenJet, ExRootTreeReader *treeReader)
    346420{
     
    398472      if(deltaR < 0.25)
    399473      {
    400         pt  = genJetMomentum.Pt();
     474        pt  = genJetMomentum.E();
    401475        eta = TMath::Abs(genJetMomentum.Eta());
    402476
     
    405479            if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 1.5)
    406480            {
    407                 histos->at(bin).cenResolHist->Fill(jetMomentum.Pt()/bestGenJetMomentum.Pt());
     481                histos->at(bin).cenResolHist->Fill(jetMomentum.E()/bestGenJetMomentum.E());
    408482            }
    409483            else if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 2.5)
    410484            {
    411                 histos->at(bin).fwdResolHist->Fill(jetMomentum.Pt()/bestGenJetMomentum.Pt());
     485                histos->at(bin).fwdResolHist->Fill(jetMomentum.E()/bestGenJetMomentum.E());
    412486            }
    413 
    414487        }
    415488      }
     
    484557
    485558  delete f1;
     559  delete f2;
    486560  return make_pair (sig, sigErr);
    487561}
     
    615689}
    616690
    617 void DrawAxis(TMultiGraph *mg, TLegend *leg, double max)
     691void DrawAxis(TMultiGraph *mg, TLegend *leg, double max, int type = 0)
    618692{
    619693  mg->SetMinimum(0.);
     
    621695  mg->GetXaxis()->SetLimits(ptrangemin,ptrangemax);
    622696  mg->GetYaxis()->SetTitle("resolution");
    623   mg->GetXaxis()->SetTitle("p_{T} [GeV]");
     697  if (type == 0) mg->GetXaxis()->SetTitle("E [GeV]");
     698  else mg->GetXaxis()->SetTitle("p_{T} [GeV]");
    624699  mg->GetYaxis()->SetTitleSize(0.07);
    625700  mg->GetXaxis()->SetTitleSize(0.07);
     
    646721}
    647722
    648 void DrawAxis(TH1D *h, TLegend *leg, int type)
     723void DrawAxis(TH1D *h, TLegend *leg, int type = 0)
    649724{
    650725
    651726  h->GetYaxis()->SetRangeUser(0,1.0);
    652   if (type == 0) h->GetXaxis()->SetTitle("p_{T} [GeV]");
    653   else h->GetXaxis()->SetTitle("#eta");
     727  if (type == 0) h->GetXaxis()->SetTitle("E [GeV]");
     728  else h->GetXaxis()->SetTitle("p_{T} [GeV]");
    654729  h->GetYaxis()->SetTitle("efficiency");
    655730  h->GetYaxis()->SetTitleSize(0.07);
     
    787862  histos_el.first->Draw("same ][");
    788863  addHist(histos_el.first, leg_el1, 1);
    789   DrawAxis(histos_eltrack.first, leg_el1, 0);
     864  DrawAxis(histos_eltrack.first, leg_el1,1);
    790865
    791866  leg_el1->Draw();
     
    802877  addHist(histos_el.second, leg_el2, 1);
    803878
    804   DrawAxis(histos_eltrack.second, leg_el2, 0);
     879  DrawAxis(histos_eltrack.second, leg_el2, 1);
    805880  leg_el2->Draw();
    806881
     
    855930  std::pair <TH1D*,TH1D*> histos_mutrack = GetEff<Track>(branchTrackMuon, branchParticleMuon, "muonTrack", muID, treeReaderMuon);
    856931
    857   // Muon Energy Resolution
     932  // Muon Pt Resolution
    858933  std::vector<resolPlot> plots_mu;
    859934  HistogramsCollection(&plots_mu, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muons");
    860   GetEres<Muon>(&plots_mu, branchMuon, branchParticleMuon, muID, treeReaderMuon);
     935  GetPtres<Muon>(&plots_mu, branchMuon, branchParticleMuon, muID, treeReaderMuon);
    861936  TGraphErrors gr_mu = EresGraph(&plots_mu, true);
    862937  TGraphErrors gr_muFwd = EresGraph(&plots_mu, false);
     
    867942  std::vector<resolPlot> plots_mutrack;
    868943  HistogramsCollection(&plots_mutrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muonsTracks");
    869   GetEres<Track>(&plots_mutrack, branchTrackMuon, branchParticleMuon, muID, treeReaderMuon);
     944  GetPtres<Track>(&plots_mutrack, branchTrackMuon, branchParticleMuon, muID, treeReaderMuon);
    870945  TGraphErrors gr_mutrack = EresGraph(&plots_mutrack, true);
    871946  TGraphErrors gr_mutrackFwd = EresGraph(&plots_mutrack, false);
     
    889964  addHist(histos_mu.first, leg_mu1, 1);
    890965
    891   DrawAxis(histos_mutrack.first, leg_mu1, 0);
     966  DrawAxis(histos_mutrack.first, leg_mu1, 1);
    892967
    893968  leg_mu1->Draw();
     
    924999  leg_mu->Draw();
    9251000
    926   DrawAxis(mg_mu, leg_mu, 0.3);
     1001  DrawAxis(mg_mu, leg_mu, 0.3, 1);
    9271002
    9281003  C_mu->cd(2);
     
    9381013  leg_muFwd->Draw();
    9391014
    940   DrawAxis(mg_muFwd, leg_muFwd, 0.3);
     1015  DrawAxis(mg_muFwd, leg_muFwd, 0.3, 1);
    9411016
    9421017  C_mu1->Print("delphes_validation.pdf","pdf");
     
    9901065  addHist(histos_ph.first, leg_ph1, 1);
    9911066
    992   DrawAxis(histos_phtower.first, leg_ph1, 0);
     1067  DrawAxis(histos_phtower.first, leg_ph1, 1);
    9931068  leg_ph1->Draw();
    9941069
     
    10251100  leg_ph->Draw();
    10261101
    1027   DrawAxis(mg_ph, leg_ph, 0.3);
     1102  DrawAxis(mg_ph, leg_ph, 0.1);
    10281103
    10291104  C_ph->cd(2);
     
    10391114  leg_phFwd->Draw();
    10401115
    1041   DrawAxis(mg_phFwd, leg_phFwd, 0.3);
     1116  DrawAxis(mg_phFwd, leg_phFwd, 0.1);
    10421117
    10431118  C_ph1->Print("delphes_validation.pdf","pdf");
     
    10961171  addHist(histos_btag.first, leg_btag1, 0);
    10971172
    1098   DrawAxis(histos_btag.first, leg_btag1, 0);
     1173  DrawAxis(histos_btag.first, leg_btag1, 1);
    10991174  leg_btag1->Draw();
    11001175
     
    11081183  addHist(histos_btag.second, leg_btag2, 0);
    11091184
    1110   DrawAxis(histos_btag.second, leg_btag2, 0);
     1185  DrawAxis(histos_btag.second, leg_btag2, 1);
    11111186  leg_btag2->Draw();
    11121187  C_btag1->cd(0);
     
    11241199  addHist(histos_tautag.first, leg_tautag1, 0);
    11251200
    1126   DrawAxis(histos_tautag.first, leg_tautag1, 0);
     1201  DrawAxis(histos_tautag.first, leg_tautag1, 1);
    11271202  leg_tautag1->Draw();
    11281203
     
    11361211  addHist(histos_tautag.second, leg_tautag2, 0);
    11371212
    1138   DrawAxis(histos_tautag.second, leg_tautag2, 0);
     1213  DrawAxis(histos_tautag.second, leg_tautag2, 1);
    11391214  leg_tautag2->Draw();
    11401215  C_tautag1->cd(0);
     
    11541229  addGraph(mg_jet, &gr_pfjets, leg_jet, 1);
    11551230
    1156   mg_jet->Draw("ACX");
     1231  mg_jet->Draw("ALX");
    11571232  leg_jet->Draw();
    11581233
    1159   DrawAxis(mg_jet, leg_jet, 0.25);
     1234  DrawAxis(mg_jet, leg_jet, 0.5);
    11601235
    11611236  C_jet->cd(2);
    11621237  TMultiGraph *mg_jetFwd = new TMultiGraph(jetResFwd,jetResFwd);
    11631238  TLegend *leg_jetFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
    1164   leg_jetFwd->SetHeader("#splitline{jets}{|#eta| < 1.5}");
     1239  leg_jetFwd->SetHeader("#splitline{jets}{1.5 < |#eta| < 2.5}");
    11651240  leg_jetFwd->AddEntry("","","");
    11661241
     
    11681243  addGraph(mg_jetFwd, &gr_pfjetsFwd, leg_jetFwd, 1);
    11691244
    1170   mg_jetFwd->Draw("ACX");
     1245  mg_jetFwd->Draw("ALX");
    11711246  leg_jetFwd->Draw();
    11721247
    1173   DrawAxis(mg_jetFwd, leg_jetFwd, 0.25);
     1248  DrawAxis(mg_jetFwd, leg_jetFwd, 0.5);
    11741249
    11751250  TString metRes = "MetRes";
     
    12171292  C_el1->Write();
    12181293  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  }
    12191302
    12201303  histos_mu.first->Write();
Note: See TracChangeset for help on using the changeset viewer.