Changeset 8f0b34c in git for examples/Validation.cpp

Timestamp:

Aug 26, 2016, 5:27:30 PM (8 years ago)

Author:

Michele Selvaggi <michele.selvaggi@…>

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.

Message:

Merge branch 'dev_01' of github.com:delphes/delphes into dev_01

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)
 {
@@ -398 +472 @@
       if(deltaR < 0.25)
       {
-        pt  = genJetMomentum.Pt();
+        pt  = genJetMomentum.E();
         eta = TMath::Abs(genJetMomentum.Eta());
 
@@ -405 +479 @@
             if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 1.5)
             {
-                histos->at(bin).cenResolHist->Fill(jetMomentum.Pt()/bestGenJetMomentum.Pt());
+                histos->at(bin).cenResolHist->Fill(jetMomentum.E()/bestGenJetMomentum.E());
             }
             else if(pt > histos->at(bin).ptmin && pt < histos->at(bin).ptmax && eta < 2.5)
             {
-                histos->at(bin).fwdResolHist->Fill(jetMomentum.Pt()/bestGenJetMomentum.Pt());
+                histos->at(bin).fwdResolHist->Fill(jetMomentum.E()/bestGenJetMomentum.E());
             }
-
         }
       }
@@ -476 +549 @@
   TF1 *f1 = new TF1("f1", "gaus", hist->GetMean()-2*hist->GetRMS(), hist->GetMean()+2*hist->GetRMS());
   hist->Fit("f1","RQ");
-  Double_t sig = f1->GetParameter(2);
-  Double_t sigErr = f1->GetParError(2);
+
+  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);
+
   delete f1;
+  delete f2;
   return make_pair (sig, sigErr);
 }
@@ -610 +689 @@
 }
 
-void DrawAxis(TMultiGraph *mg, TLegend *leg, double max)
+void DrawAxis(TMultiGraph *mg, TLegend *leg, double max, int type = 0)
 {
   mg->SetMinimum(0.);
@@ -616 +695 @@
   mg->GetXaxis()->SetLimits(ptrangemin,ptrangemax);
   mg->GetYaxis()->SetTitle("resolution");
-  mg->GetXaxis()->SetTitle("p_{T} [GeV]");
+  if (type == 0) mg->GetXaxis()->SetTitle("E [GeV]");
+  else mg->GetXaxis()->SetTitle("p_{T} [GeV]");
   mg->GetYaxis()->SetTitleSize(0.07);
   mg->GetXaxis()->SetTitleSize(0.07);
@@ -641 +721 @@
 }
 
-void DrawAxis(TH1D *h, TLegend *leg, int type)
+void DrawAxis(TH1D *h, TLegend *leg, int type = 0)
 {
 
   h->GetYaxis()->SetRangeUser(0,1.0);
-  if (type == 0) h->GetXaxis()->SetTitle("p_{T} [GeV]");
-  else h->GetXaxis()->SetTitle("#eta");
+  if (type == 0) h->GetXaxis()->SetTitle("E [GeV]");
+  else h->GetXaxis()->SetTitle("p_{T} [GeV]");
   h->GetYaxis()->SetTitle("efficiency");
   h->GetYaxis()->SetTitleSize(0.07);
@@ -782 +862 @@
   histos_el.first->Draw("same ][");
   addHist(histos_el.first, leg_el1, 1);
-  DrawAxis(histos_eltrack.first, leg_el1, 0);
+  DrawAxis(histos_eltrack.first, leg_el1,1);
 
   leg_el1->Draw();
@@ -797 +877 @@
   addHist(histos_el.second, leg_el2, 1);
 
-  DrawAxis(histos_eltrack.second, leg_el2, 0);
+  DrawAxis(histos_eltrack.second, leg_el2, 1);
   leg_el2->Draw();
 
@@ -850 +930 @@
   std::pair <TH1D*,TH1D*> histos_mutrack = GetEff<Track>(branchTrackMuon, branchParticleMuon, "muonTrack", muID, treeReaderMuon);
 
-  // Muon Energy Resolution
+  // Muon Pt Resolution
   std::vector<resolPlot> plots_mu;
   HistogramsCollection(&plots_mu, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muons");
-  GetEres<Muon>(&plots_mu, branchMuon, branchParticleMuon, muID, treeReaderMuon);
+  GetPtres<Muon>(&plots_mu, branchMuon, branchParticleMuon, muID, treeReaderMuon);
   TGraphErrors gr_mu = EresGraph(&plots_mu, true);
   TGraphErrors gr_muFwd = EresGraph(&plots_mu, false);
@@ -862 +942 @@
   std::vector<resolPlot> plots_mutrack;
   HistogramsCollection(&plots_mutrack, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "muonsTracks");
-  GetEres<Track>(&plots_mutrack, branchTrackMuon, branchParticleMuon, muID, treeReaderMuon);
+  GetPtres<Track>(&plots_mutrack, branchTrackMuon, branchParticleMuon, muID, treeReaderMuon);
   TGraphErrors gr_mutrack = EresGraph(&plots_mutrack, true);
   TGraphErrors gr_mutrackFwd = EresGraph(&plots_mutrack, false);
@@ -884 +964 @@
   addHist(histos_mu.first, leg_mu1, 1);
 
-  DrawAxis(histos_mutrack.first, leg_mu1, 0);
+  DrawAxis(histos_mutrack.first, leg_mu1, 1);
 
   leg_mu1->Draw();
@@ -919 +999 @@
   leg_mu->Draw();
 
-  DrawAxis(mg_mu, leg_mu, 0.3);
+  DrawAxis(mg_mu, leg_mu, 0.3, 1);
 
   C_mu->cd(2);
@@ -933 +1013 @@
   leg_muFwd->Draw();
 
-  DrawAxis(mg_muFwd, leg_muFwd, 0.3);
+  DrawAxis(mg_muFwd, leg_muFwd, 0.3, 1);
 
   C_mu1->Print("delphes_validation.pdf","pdf");
@@ -985 +1065 @@
   addHist(histos_ph.first, leg_ph1, 1);
 
-  DrawAxis(histos_phtower.first, leg_ph1, 0);
+  DrawAxis(histos_phtower.first, leg_ph1, 1);
   leg_ph1->Draw();
 
@@ -1020 +1100 @@
   leg_ph->Draw();
 
-  DrawAxis(mg_ph, leg_ph, 0.3);
+  DrawAxis(mg_ph, leg_ph, 0.1);
 
   C_ph->cd(2);
@@ -1034 +1114 @@
   leg_phFwd->Draw();
 
-  DrawAxis(mg_phFwd, leg_phFwd, 0.3);
+  DrawAxis(mg_phFwd, leg_phFwd, 0.1);
 
   C_ph1->Print("delphes_validation.pdf","pdf");
@@ -1091 +1171 @@
   addHist(histos_btag.first, leg_btag1, 0);
 
-  DrawAxis(histos_btag.first, leg_btag1, 0);
+  DrawAxis(histos_btag.first, leg_btag1, 1);
   leg_btag1->Draw();
 
@@ -1103 +1183 @@
   addHist(histos_btag.second, leg_btag2, 0);
 
-  DrawAxis(histos_btag.second, leg_btag2, 0);
+  DrawAxis(histos_btag.second, leg_btag2, 1);
   leg_btag2->Draw();
   C_btag1->cd(0);
@@ -1119 +1199 @@
   addHist(histos_tautag.first, leg_tautag1, 0);
 
-  DrawAxis(histos_tautag.first, leg_tautag1, 0);
+  DrawAxis(histos_tautag.first, leg_tautag1, 1);
   leg_tautag1->Draw();
 
@@ -1131 +1211 @@
   addHist(histos_tautag.second, leg_tautag2, 0);
 
-  DrawAxis(histos_tautag.second, leg_tautag2, 0);
+  DrawAxis(histos_tautag.second, leg_tautag2, 1);
   leg_tautag2->Draw();
   C_tautag1->cd(0);
@@ -1149 +1229 @@
   addGraph(mg_jet, &gr_pfjets, leg_jet, 1);
 
-  mg_jet->Draw("ACX");
+  mg_jet->Draw("ALX");
   leg_jet->Draw();
 
-  DrawAxis(mg_jet, leg_jet, 0.25);
+  DrawAxis(mg_jet, leg_jet, 0.5);
 
   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}{|#eta| < 1.5}");
+  leg_jetFwd->SetHeader("#splitline{jets}{1.5 < |#eta| < 2.5}");
   leg_jetFwd->AddEntry("","","");
 
@@ -1163 +1243 @@
   addGraph(mg_jetFwd, &gr_pfjetsFwd, leg_jetFwd, 1);
 
-  mg_jetFwd->Draw("ACX");
+  mg_jetFwd->Draw("ALX");
   leg_jetFwd->Draw();
 
-  DrawAxis(mg_jetFwd, leg_jetFwd, 0.25);
+  DrawAxis(mg_jetFwd, leg_jetFwd, 0.5);
 
   TString metRes = "MetRes";
@@ -1212 +1292 @@
   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();