Changes in / [8f0b34c:eab6d46] in git

File:

• examples/Validation.cpp (modified) (24 diffs)

examples/Validation.cpp

@@ -343 +343 @@
   }
 }
-
-
-template<typename T>
-void GetPtres(std::vector<resolPlot> *histos, TClonesArray *branchReco, TClonesArray *branchParticle, int pdgID, ExRootTreeReader *treeReader)
-{
-  Long64_t allEntries = treeReader->GetEntries();
-
-  cout << "** Computing pt resolution of " << branchReco->GetName() << " induced by " << branchParticle->GetName() << " with PID " << pdgID << endl;
-
-  GenParticle *particle;
-  T* recoObj;
-
-  TLorentzVector recoMomentum, genMomentum, bestGenMomentum;
-
-  Float_t deltaR;
-  Float_t pt, eta;
-  Long64_t entry;
-
-  Int_t i, j, bin;
-
-  // Loop over all events
-  for(entry = 0; entry < allEntries; ++entry)
-  {
-    // Load selected branches with data from specified event
-    treeReader->ReadEntry(entry);
-
-    // Loop over all reconstructed jets in event
-    for(i = 0; i < branchReco->GetEntriesFast(); ++i)
-    {
-      recoObj = (T*) branchReco->At(i);
-      recoMomentum = recoObj->P4();
-
-      deltaR = 999;
-
-     // Loop over all hard partons in event
-     for(j = 0; j < branchParticle->GetEntriesFast(); ++j)
-     {
-        particle = (GenParticle*) branchParticle->At(j);
-        if (particle->PID == pdgID && particle->Status == 1)
-        {
-          genMomentum = particle->P4();
-
-          // this is simply to avoid warnings from initial state particle
-          // having infite rapidity ...
-          if(genMomentum.Px() == 0 && genMomentum.Py() == 0) continue;
-
-          // take the closest parton candidate
-          if(genMomentum.DeltaR(recoMomentum) < deltaR)
-          {
-             deltaR = genMomentum.DeltaR(recoMomentum);
-             bestGenMomentum = genMomentum;
-          }
-        }
-      }
-
-      if(deltaR < 0.3)
-      {
-        pt  = bestGenMomentum.Pt();
-        eta = TMath::Abs(bestGenMomentum.Eta());
-
-        for (bin = 0; bin < Nbins; bin++)
-        {
-          if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 2.5)
-          {
-            if (eta < 1.5) {histos->at(bin).cenResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());}
-            else if (eta < 2.5) {histos->at(bin).fwdResolHist->Fill(recoMomentum.Pt()/bestGenMomentum.Pt());}
-          }
-        }
-      }
-    }
-  }
-}
-
-
 void GetJetsEres(std::vector<resolPlot> *histos, TClonesArray *branchJet, TClonesArray *branchGenJet, ExRootTreeReader *treeReader)
 {
@@ -472 +398 @@
       if(deltaR < 0.25)
       {
-        pt  = genJetMomentum.E();
+        pt  = genJetMomentum.Pt();
         eta = TMath::Abs(genJetMomentum.Eta());
 
@@ -479 +405 @@
             if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 1.5)
             {
-                histos->at(bin).cenResolHist->Fill(jetMomentum.E()/bestGenJetMomentum.E());
+                histos->at(bin).cenResolHist->Fill(jetMomentum.Pt()/bestGenJetMomentum.Pt());
             }
             else if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 2.5)
             {
-                histos->at(bin).fwdResolHist->Fill(jetMomentum.E()/bestGenJetMomentum.E());
+                histos->at(bin).fwdResolHist->Fill(jetMomentum.Pt()/bestGenJetMomentum.Pt());
             }
+
         }
       }
@@ -549 +476 @@
   TF1 *f1 = new TF1("f1", "gaus", hist->GetMean()-2*hist->GetRMS(), hist->GetMean()+2*hist->GetRMS());
   hist->Fit("f1","RQ");
-
-  TF1 *f2 = new TF1("f2", "gaus", f1->GetParameter(1) - 2*f1->GetParameter(2), f1->GetParameter(1) + 2*f1->GetParameter(2));
-  hist->Fit("f2","RQ");  
-
-  Double_t sig = f2->GetParameter(2);
-  Double_t sigErr = f2->GetParError(2);
-
+  Double_t sig = f1->GetParameter(2);
+  Double_t sigErr = f1->GetParError(2);
   delete f1;
-  delete f2;
   return make_pair (sig, sigErr);
 }
@@ -689 +610 @@
 }
 
-void DrawAxis(TMultiGraph *mg, TLegend *leg, double max, int type = 0)
+void DrawAxis(TMultiGraph *mg, TLegend *leg, double max)
 {
   mg->SetMinimum(0.);
@@ -695 +616 @@
   mg->GetXaxis()->SetLimits(ptrangemin,ptrangemax);
   mg->GetYaxis()->SetTitle("resolution");
-  if (type == 0) mg->GetXaxis()->SetTitle("E [GeV]");
-  else mg->GetXaxis()->SetTitle("p_{T} [GeV]");
+  mg->GetXaxis()->SetTitle("p_{T} [GeV]");
   mg->GetYaxis()->SetTitleSize(0.07);
   mg->GetXaxis()->SetTitleSize(0.07);
@@ -721 +641 @@
 }
 
-void DrawAxis(TH1D *h, TLegend *leg, int type = 0)
+void DrawAxis(TH1D *h, TLegend *leg, int type)
 {
 
   h->GetYaxis()->SetRangeUser(0,1.0);
-  if (type == 0) h->GetXaxis()->SetTitle("E [GeV]");
-  else h->GetXaxis()->SetTitle("p_{T} [GeV]");
+  if (type == 0) h->GetXaxis()->SetTitle("p_{T} [GeV]");
+  else h->GetXaxis()->SetTitle("#eta");
   h->GetYaxis()->SetTitle("efficiency");
   h->GetYaxis()->SetTitleSize(0.07);
@@ -862 +782 @@
   histos_el.first->Draw("same ][");
   addHist(histos_el.first, leg_el1, 1);
-  DrawAxis(histos_eltrack.first, leg_el1,1);
+  DrawAxis(histos_eltrack.first, leg_el1, 0);
 
   leg_el1->Draw();
@@ -877 +797 @@
   addHist(histos_el.second, leg_el2, 1);
 
-  DrawAxis(histos_eltrack.second, leg_el2, 1);
+  DrawAxis(histos_eltrack.second, leg_el2, 0);
   leg_el2->Draw();
 
@@ -930 +850 @@
   std::pair <TH1D*,TH1D*> histos_mutrack = GetEff<Track>(branchTrackMuon, branchParticleMuon, "muonTrack", muID, treeReaderMuon);
 
-  // Muon Pt Resolution
+  // Muon Energy Resolution
   std::vector<resolPlot> plots_mu;
   HistogramsCollection(&plots_mu, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muons");
-  GetPtres<Muon>(&plots_mu, branchMuon, branchParticleMuon, muID, treeReaderMuon);
+  GetEres<Muon>(&plots_mu, branchMuon, branchParticleMuon, muID, treeReaderMuon);
   TGraphErrors gr_mu = EresGraph(&plots_mu, true);
   TGraphErrors gr_muFwd = EresGraph(&plots_mu, false);
@@ -942 +862 @@
   std::vector<resolPlot> plots_mutrack;
   HistogramsCollection(&plots_mutrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muonsTracks");
-  GetPtres<Track>(&plots_mutrack, branchTrackMuon, branchParticleMuon, muID, treeReaderMuon);
+  GetEres<Track>(&plots_mutrack, branchTrackMuon, branchParticleMuon, muID, treeReaderMuon);
   TGraphErrors gr_mutrack = EresGraph(&plots_mutrack, true);
   TGraphErrors gr_mutrackFwd = EresGraph(&plots_mutrack, false);
@@ -964 +884 @@
   addHist(histos_mu.first, leg_mu1, 1);
 
-  DrawAxis(histos_mutrack.first, leg_mu1, 1);
+  DrawAxis(histos_mutrack.first, leg_mu1, 0);
 
   leg_mu1->Draw();
@@ -999 +919 @@
   leg_mu->Draw();
 
-  DrawAxis(mg_mu, leg_mu, 0.3, 1);
+  DrawAxis(mg_mu, leg_mu, 0.3);
 
   C_mu->cd(2);
@@ -1013 +933 @@
   leg_muFwd->Draw();
 
-  DrawAxis(mg_muFwd, leg_muFwd, 0.3, 1);
+  DrawAxis(mg_muFwd, leg_muFwd, 0.3);
 
   C_mu1->Print("delphes_validation.pdf","pdf");
@@ -1065 +985 @@
   addHist(histos_ph.first, leg_ph1, 1);
 
-  DrawAxis(histos_phtower.first, leg_ph1, 1);
+  DrawAxis(histos_phtower.first, leg_ph1, 0);
   leg_ph1->Draw();
 
@@ -1100 +1020 @@
   leg_ph->Draw();
 
-  DrawAxis(mg_ph, leg_ph, 0.1);
+  DrawAxis(mg_ph, leg_ph, 0.3);
 
   C_ph->cd(2);
@@ -1114 +1034 @@
   leg_phFwd->Draw();
 
-  DrawAxis(mg_phFwd, leg_phFwd, 0.1);
+  DrawAxis(mg_phFwd, leg_phFwd, 0.3);
 
   C_ph1->Print("delphes_validation.pdf","pdf");
@@ -1171 +1091 @@
   addHist(histos_btag.first, leg_btag1, 0);
 
-  DrawAxis(histos_btag.first, leg_btag1, 1);
+  DrawAxis(histos_btag.first, leg_btag1, 0);
   leg_btag1->Draw();
 
@@ -1183 +1103 @@
   addHist(histos_btag.second, leg_btag2, 0);
 
-  DrawAxis(histos_btag.second, leg_btag2, 1);
+  DrawAxis(histos_btag.second, leg_btag2, 0);
   leg_btag2->Draw();
   C_btag1->cd(0);
@@ -1199 +1119 @@
   addHist(histos_tautag.first, leg_tautag1, 0);
 
-  DrawAxis(histos_tautag.first, leg_tautag1, 1);
+  DrawAxis(histos_tautag.first, leg_tautag1, 0);
   leg_tautag1->Draw();
 
@@ -1211 +1131 @@
   addHist(histos_tautag.second, leg_tautag2, 0);
 
-  DrawAxis(histos_tautag.second, leg_tautag2, 1);
+  DrawAxis(histos_tautag.second, leg_tautag2, 0);
   leg_tautag2->Draw();
   C_tautag1->cd(0);
@@ -1229 +1149 @@
   addGraph(mg_jet, &gr_pfjets, leg_jet, 1);
 
-  mg_jet->Draw("ALX");
+  mg_jet->Draw("ACX");
   leg_jet->Draw();
 
-  DrawAxis(mg_jet, leg_jet, 0.5);
+  DrawAxis(mg_jet, leg_jet, 0.25);
 
   C_jet->cd(2);
   TMultiGraph *mg_jetFwd = new TMultiGraph(jetResFwd,jetResFwd);
   TLegend *leg_jetFwd = new TLegend(resLegXmin,resLegYmin,resLegXmax,resLegYmax);
-  leg_jetFwd->SetHeader("#splitline{jets}{1.5 < |#eta| < 2.5}");
+  leg_jetFwd->SetHeader("#splitline{jets}{|#eta| < 1.5}");
   leg_jetFwd->AddEntry("","","");
 
@@ -1243 +1163 @@
   addGraph(mg_jetFwd, &gr_pfjetsFwd, leg_jetFwd, 1);
 
-  mg_jetFwd->Draw("ALX");
+  mg_jetFwd->Draw("ACX");
   leg_jetFwd->Draw();
 
-  DrawAxis(mg_jetFwd, leg_jetFwd, 0.5);
+  DrawAxis(mg_jetFwd, leg_jetFwd, 0.25);
 
   TString metRes = "MetRes";
@@ -1292 +1212 @@
   C_el1->Write();
   C_el2->Write();
-
-  for (int bin = 0; bin < Nbins; bin++)
-  {
-    plots_mu.at(bin).cenResolHist->Write();
-    plots_mutrack.at(bin).cenResolHist->Write();
-    plots_mu.at(bin).fwdResolHist->Write();
-    plots_mutrack.at(bin).fwdResolHist->Write();
-  }
 
   histos_mu.first->Write();