Changes in / [9f8d4e7:df315dc] in git

File:

• examples/Validation.C (modified) (12 diffs)

examples/Validation.C

@@ -49 +49 @@
     double ptmin;
     double ptmax;
+    double xmin;
+    double xmax;
     TString obj;
 
     resolPlot();
     resolPlot(double ptdown, double ptup, TString object);
-    void set(double ptdown, double ptup, TString object);
+    void set(double ptdown, double ptup, TString object, double xmin = 0, double xmax = 2);
     void print(){std::cout << ptmin << std::endl;}
 };
@@ -67 +69 @@
 }
 
-void resolPlot::set(double ptdown, double ptup, TString object){
+void resolPlot::set(double ptdown, double ptup, TString object, double xmin, double xmax){
     ptmin = ptdown;
     ptmax = ptup;
     obj = object;
 
-    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);
-    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);
-
-}
-
-void HistogramsCollection(std::vector<resolPlot> *histos, double ptmin, double ptmax, TString obj)
+    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);
+    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);
+
+}
+
+void HistogramsCollection(std::vector<resolPlot> *histos, double ptmin, double ptmax, TString obj, double xmin = 0, double xmax = 2)
 {
     double width;
@@ -89 +91 @@
         ptdown = TMath::Power(10,ptmin + i * width );
         ptup = TMath::Power(10,ptmin + (i+1) * width );
-        ptemp.set(ptdown, ptup, obj);
+        ptemp.set(ptdown, ptup, obj, xmin, xmax);
         histos->push_back(ptemp);
     }
@@ -362 +364 @@
   }
 }
+
+void GetMetres(std::vector<resolPlot> *histos, TClonesArray *branchScalarHT, TClonesArray *branchMet, TClonesArray *branchJet, ExRootTreeReader *treeReader)
+{
+
+  Long64_t allEntries = treeReader->GetEntries();
+
+  cout << "** Computing resolution of " << branchMet->GetName() << " vs " << branchScalarHT->GetName() << endl;
+
+  MissingET *met;
+  ScalarHT *scalarHT;
+
+  Long64_t entry;
+
+  Int_t bin;
+  Double_t ht;
+
+  // Loop over all events
+  for(entry = 0; entry < allEntries; ++entry)
+  {
+    // Load selected branches with data from specified event
+    treeReader->ReadEntry(entry);
+
+    if(entry%10000 == 0) cout << "Event number: "<< entry <<endl;
+
+    if (branchJet->GetEntriesFast() > 1)
+    {
+      met = (MissingET*) branchMet->At(0);
+      scalarHT = (ScalarHT*) branchScalarHT->At(0);
+      ht = scalarHT->HT;
+
+      for (bin = 0; bin < Nbins; bin++)
+      {
+        if(ht > histos->at(bin).ptmin && ht < histos->at(bin).ptmax ) 
+        {
+          histos->at(bin).cenResolHist->Fill(met->P4().Px());
+          histos->at(bin).fwdResolHist->Fill(met->P4().Py());
+        }
+      }
+    }
+  }
+}
+
+
 std::pair<Double_t, Double_t> GausFit(TH1* hist)
 {
@@ -374 +419 @@
 
 
-TGraphErrors EresGraph(std::vector<resolPlot> *histos, bool central)
+TGraphErrors EresGraph(std::vector<resolPlot> *histos, bool central, bool rms = false)
 {
     Int_t bin;
@@ -389 +434 @@
             std::cout << " entries : " << histos->at(bin).cenResolHist->GetEntries() << std::endl;
             std::pair<Double_t, Double_t> sigvalues = GausFit(histos->at(bin).cenResolHist);
-            gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first);
-            gr.SetPointError(count,0, sigvalues.second);
+            if (rms == true) 
+            {
+              gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.second);
+              gr.SetPointError(count,0, sigvalues.second); // to correct
+            }
+            else 
+            {
+              gr.SetPoint(count,(histos->at(bin).ptmin+histos->at(bin).ptmax)/2.0, sigvalues.first);
+              gr.SetPointError(count,0, sigvalues.second);
+            }
             count++;
         }
@@ -551 +604 @@
   TClonesArray *branchPFJet = treeReader->UseBranch("Jet");
   TClonesArray *branchCaloJet = treeReader->UseBranch("CaloJet");
+  TClonesArray *branchScalarHT = treeReader->UseBranch("ScalarHT");
+  TClonesArray *branchMet = treeReader->UseBranch("MissingET");
 
 #ifdef ELECTRON
@@ -811 +866 @@
   gr_pfjets.SetName("pfJet");
 
-
-  // PFJets Energy Resolution
+  // CaloJets Energy Resolution
   std::vector<resolPlot> plots_calojets;
   HistogramsCollection(&plots_calojets, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "CaloJet");
@@ -819 +873 @@
   gr_calojets.SetName("caloJet");
 
+  // MET Resolution vs HT
+  std::vector<resolPlot> plots_met;
+  HistogramsCollection(&plots_met, TMath::Log10(ptrangemin), TMath::Log10(ptrangemax), "MET", -200, 200);
+  GetMetres( &plots_met, branchScalarHT, branchMet, branchPFJet, treeReader);
+  TGraphErrors gr_met = EresGraph(&plots_met, true);
+  gr_calojets.SetName("MET");
+
   // Canvas
 
@@ -882 +943 @@
   C_btag1->SaveAs(btagEff+".eps");
   C_jet->SaveAs(jetRes+".eps");
+
+  TString metRes = "MetRes";
+  TCanvas *C_met = new TCanvas(metRes,metRes, 1000, 500);
+
+  TMultiGraph *mg_met = new TMultiGraph(metRes,metRes);
+  TLegend *leg_met = new TLegend(0.52,0.7,0.9,0.9);
+
+  addGraph(mg_met, &gr_met, leg_met, 3);
+
+  mg_met->Draw("ACX");
+  leg_met->Draw();
+
+  DrawAxis(mg_met, leg_met, 100);
+
+  C_met->SaveAs(metRes+".eps");
+
 
 #endif
@@ -934 +1011 @@
       plots_pfjets.at(bin).cenResolHist->Write();
       plots_calojets.at(bin).cenResolHist->Write();
+      plots_met.at(bin).cenResolHist->Write();
   }
   histos_btag.first->Write();
@@ -942 +1020 @@
   C_tautag1->Write();
   C_jet->Write();
+  C_met->Write();
 #endif