Eres fit

ecal/scripts/rec_emcal_resolution_analysis.cxx

+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot energy resolution
+// 2021/01/20
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "dd4pod/CalorimeterHitCollection.h"
+#include "dd4pod/TrackerHitCollection.h"
+#include "eicd/CalorimeterHitCollection.h"
+#include "eicd/CalorimeterHitData.h"
+#include "eicd/ClusterCollection.h"
+#include "eicd/ClusterData.h"
+
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TMath.h"
+#include "TH1.h"
+#include "TH1D.h"
+
+// If root cannot find a dd4pod header make sure to add the include dir to ROOT_INCLUDE_PATH:
+//  export ROOT_INCLUDE_PATH=$HOME/include:$ROOT_INCLUDE_PATH
+//  export ROOT_INCLUDE_PATH=/home/jihee/stow/development/include:$ROOT_INCLUDE_PATH
+
+using ROOT::RDataFrame;
+using namespace ROOT::VecOps;
+
+void rec_emcal_resolution_analysis(const char* input_fname = "rec_emcal_uniform_electrons.root")
+{
+  // Setting for graphs
+  gROOT->SetStyle("Plain");
+  gStyle->SetOptFit(1);
+  gStyle->SetLineWidth(2);
+  gStyle->SetPadTickX(1);
+  gStyle->SetPadTickY(1);
+  gStyle->SetPadGridX(1);
+  gStyle->SetPadGridY(1);
+  gStyle->SetPadLeftMargin(0.14);
+  gStyle->SetPadRightMargin(0.14);
+
+  ROOT::EnableImplicitMT();
+  ROOT::RDataFrame d0("events", input_fname);
+
+  // Number of Clusters
+  auto ncluster = [] (const std::vector<eic::ClusterData>& evt) {return (int) evt.size(); };
+
+  // Cluster Energy [GeV]
+  auto clusterE = [] (const std::vector<eic::ClusterData>& evt) {
+  	std::vector<double> result;
+  	for (const auto& i: evt)
+    		result.push_back(i.energy);
+  return result;
+  };
+
+  // Thrown Energy [GeV]
+  auto E_thr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+  	std::vector<double> result;
+    	result.push_back(TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass));
+  return result;
+  };
+
+  // Energy Resolution = delta(E)/E 
+  // delta(E) = (Erec - Ethr)
+  // E        = Ethr
+  auto E_res = [] (const std::vector<double>& Erec, const std::vector<double>& Ethr) {
+        std::vector<double> result;
+        for (const auto& E1 : Ethr) {
+                for (const auto& E2 : Erec) {
+                        result.push_back((E2-E1)/E1);
+                }
+        }
+  return result;
+  };
+
+  // Define variables
+  auto d1 = d0.Define("ncluster", ncluster, {"EcalClusters"})
+	      .Define("clusterE", clusterE, {"EcalClusters"})
+	      .Define("E_thr",    E_thr,    {"mcparticles2"})
+	      .Define("E_res",    E_res,    {"clusterE","E_thr"})
+	      ;
+
+  // Define Cuts
+  // d2: Select Events with One Cluster
+  // d3: Select Events which determined their center of cluster not being on the edge of detector in addition to d2 cut
+  auto d2 = d1.Filter("ncluster==1");
+  auto d3 = d2.Filter([=] (const std::vector<eic::ClusterData>& evt) {
+                                        for(const auto& i: evt) {
+                                                auto pos_x = i.position.x;
+                                                auto pos_y = i.position.y;
+                                                auto radial_dist = TMath::Sqrt(pos_x*pos_x + pos_y*pos_y);
+                                                if (radial_dist > 8.0 && radial_dist < 30.0)
+                                                        return true;
+                                        }
+                                   return false;}, {"EcalClusters"});
+  
+  // Define histograms
+  auto hEresNC1      = d2.Histo1D({"hEresNC1",      "Energy Resolution CE vs Ethr w/ NC1; #DeltaE/E; Events",       100, -0.5, 0.5}, "E_res"); 
+  auto hEresNC1EGCUT = d3.Histo1D({"hEresNC1EGCUT", "Energy Resolution CE vs Ethr w/ NC1 EGCUT; #DeltaE/E; Events", 100, -0.5, 0.5}, "E_res");
+ 
+  // Energy Resolution Graphs
+  // Divided into serveral energy bin group
+  // Do gaussian git to extract mean and sigma
+  int    nbin      = 12;                       // number of energy bin group
+  double maxEnergy = 30.0;                     // maximum of energy in dataset
+  double multiple  = maxEnergy / (double)nbin; // energy interval between energy bin group
+  double maximumE;                             // high end in each energy bin group
+  double minimumE;                             // low end in each energy bin group
+  double mean[nbin];                           // mean from gaussian fit
+  double sigma[nbin];                          // sigma from gaussian fit
+  double dmean[nbin];                          // error on mean
+  double dsigma[nbin];                         // error on sigma
+  double gEbin[nbin];                          // center of energy in each group
+  double invSqrtgEbin[nbin];                   // inverse sqrt center of energy in each group
+  double resolution[nbin];                     // resolution = sigma/mean
+  double resolutionError[nbin];                // error on resolution
+  double resolutionErrorSigma;                 // error on resolution sigma term
+  double resolutionErrorMean;                  // error on resolution mean term
+  double xerror[nbin];                         // set to 0.0 for all
+  // Loop over each energy bin group
+  // Define hitogram and do gaussian fit
+  // Save fit results; mean, sigma, error on mean & sigam
+  // Calculate resolution and resolution error
+  for(int n=1; n <= nbin; n++) {
+	maximumE          = n * multiple;
+	minimumE          = maximumE - multiple;
+  	gEbin[n-1]        = (maximumE + minimumE) / 2.0;
+	invSqrtgEbin[n-1] = 1.0 / TMath::Sqrt(gEbin[n-1]); 
+	// Define histogram
+	auto h = d3.Filter([=] (const std::vector<eic::ClusterData>& evt) {
+					for (const auto& i: evt) {
+                                        	auto eng     = i.energy;
+						if (eng >= minimumE && eng < maximumE)
+                                                	return true;
+                                        }
+                                        return false;}, {"EcalClusters"})
+		   .Histo1D({"h", "Energy Resolution; #DeltaE/E; Events", 100, -0.5, 0.5}, "E_res");
+	// Gaussian fit
+  	h->Fit("gaus","W,E");
+  	// Save fit results
+	// how to access how many events in histogram: h->GetEntries();
+	if(h->GetFunction("gaus")) {
+  		mean[n-1]            = TMath::Abs(h->GetFunction("gaus")->GetParameter(1));
+  		dmean[n-1]           = h->GetFunction("gaus")->GetParError(1);
+  		sigma[n-1]           = h->GetFunction("gaus")->GetParameter(2);
+  		dsigma[n-1]          = h->GetFunction("gaus")->GetParError(2);
+  		resolution[n-1]      = sigma[n-1] / mean[n-1];
+		resolutionErrorSigma = dsigma[n-1] / mean[n-1];
+		resolutionErrorMean  = dmean[n-1] * sigma[n-1];
+		resolutionError[n-1] = TMath::Sqrt(resolutionErrorSigma*resolutionErrorSigma +  resolutionErrorMean*resolutionErrorMean);
+     		xerror[n-1]          = 0.0;
+  	}
+	// Draw histogram, but
+	// Just save fit results
+	//TCanvas *c = new TCanvas("c", "c", 500, 500);
+        //h->GetYaxis()->SetTitleOffset(1.4);
+        //h->SetLineWidth(2);
+        //h->SetLineColor(kBlue);
+        //h->Fit("gaus","W,E");
+        //h->DrawClone();
+        //c->SaveAs("results/h.png");
+  }
+
+  // Generate graphs using gaussian fit results
+  // Cluster energy VS Simga
+  TGraphErrors* gCEVsSigmaNC1EGCUT = new TGraphErrors(nbin,gEbin,sigma,xerror,dsigma);
+  TCanvas *c1 = new TCanvas("c1","c1",500,500);
+  gCEVsSigmaNC1EGCUT->SetTitle("Cluster Energy vs #sigma");
+  gCEVsSigmaNC1EGCUT->GetYaxis()->SetTitleOffset(1.4);
+  gCEVsSigmaNC1EGCUT->GetXaxis()->SetTitle("Cluster Energy [GeV]");
+  gCEVsSigmaNC1EGCUT->GetYaxis()->SetTitle("#sigma");
+  gCEVsSigmaNC1EGCUT->GetYaxis()->SetRangeUser(0.0,0.05);
+  gCEVsSigmaNC1EGCUT->SetMarkerStyle(21);
+  gCEVsSigmaNC1EGCUT->DrawClone("AEP");
+  c1->SaveAs("results/emcal_electrons_gCEVsSigmaNC1EGCUT.png");
+  c1->SaveAs("results/emcal_electrons_gCEVsSigmaNC1EGCUT.pdf");
+  // Cluster energy VS Mean
+  TGraphErrors* gCEVsMeanNC1EGCUT = new TGraphErrors(nbin,gEbin,mean,xerror,dmean);
+  TCanvas *c2 = new TCanvas("c2","c2",500,500);
+  gCEVsMeanNC1EGCUT->SetTitle("Cluster Energy vs #mu");
+  gCEVsMeanNC1EGCUT->GetYaxis()->SetTitleOffset(1.4);
+  gCEVsMeanNC1EGCUT->GetXaxis()->SetTitle("Cluster Energy [GeV]");
+  gCEVsMeanNC1EGCUT->GetYaxis()->SetTitle("#mu");
+  gCEVsMeanNC1EGCUT->GetYaxis()->SetRangeUser(0.0,0.1);
+  gCEVsMeanNC1EGCUT->SetMarkerStyle(21);
+  gCEVsMeanNC1EGCUT->DrawClone("AEP");
+  c2->SaveAs("results/emcal_electrons_gCEVsMeanNC1EGCUT.png");
+  c2->SaveAs("results/emcal_electrons_gCEVsMeanNC1EGCUT.pdf");
+  // Inverse of sqrt(cluster energy) VS resolution
+  TGraphErrors* gInvSqCEVsSigmaENC1EGCUT = new TGraphErrors(nbin,invSqrtgEbin,resolution,xerror,resolutionError);
+  TCanvas *c3 = new TCanvas("c3","c3",500,500);
+  gInvSqCEVsSigmaENC1EGCUT->SetTitle("#sigma/E [%] vs 1/sqrt(Cluster Energy)");
+  gInvSqCEVsSigmaENC1EGCUT->GetYaxis()->SetTitleOffset(1.4);
+  gInvSqCEVsSigmaENC1EGCUT->GetXaxis()->SetTitle("1/sqrt(Cluster Energy)");
+  gInvSqCEVsSigmaENC1EGCUT->GetYaxis()->SetTitle("#sigma/E [%]");
+  gInvSqCEVsSigmaENC1EGCUT->GetYaxis()->SetRangeUser(0.0,1.5);
+  gInvSqCEVsSigmaENC1EGCUT->GetXaxis()->SetRangeUser(0.0,1.5);
+  gInvSqCEVsSigmaENC1EGCUT->SetMarkerStyle(21);
+  gInvSqCEVsSigmaENC1EGCUT->DrawClone("AEP");
+  c3->SaveAs("results/emcal_electrons_gInvSqCEVsSimgaENC1EGCUT.png");
+  c3->SaveAs("results/emcal_electrons_gInvSqCEVsSigmaENC1EGCUT.pdf");
+  
+  // Fit Function
+  // sigma/E = S/sqrt(E) + N/E + C; added in quadrature
+  // S: Stochastic term
+  // N: Niose term
+  // C: Constant term
+  TF1 *f1 = new TF1("f1","sqrt([0]*[0] + pow([1]/sqrt(x),2) + pow([2]/x,2))", 0.5, 30.0);
+  f1->SetParNames("C", "S", "N");    // param names
+  f1->SetParameters(0.30,2.80,0.12); // initial param values
+  f1->SetLineWidth(2);
+  f1->SetLineColor(kRed);
+
+  // Cluster energy VS resolution
+  TGraphErrors* gCEVsSigmaENC1EGCUT = new TGraphErrors(nbin,gEbin,resolution,xerror,resolutionError);
+  TCanvas *c4 = new TCanvas("c4","c4",500,500);
+  gCEVsSigmaENC1EGCUT->SetTitle("Cluster Energy vs #sigma/E [%]");
+  gCEVsSigmaENC1EGCUT->GetYaxis()->SetTitleOffset(1.4);
+  gCEVsSigmaENC1EGCUT->GetXaxis()->SetTitle("Cluster Energy [GeV]");
+  gCEVsSigmaENC1EGCUT->GetYaxis()->SetTitle("#sigma/E [%]");
+  gCEVsSigmaENC1EGCUT->GetYaxis()->SetRangeUser(0.0,1.5);
+  gCEVsSigmaENC1EGCUT->SetMarkerStyle(21);
+  gCEVsSigmaENC1EGCUT->Fit("f1","W,E");
+  gCEVsSigmaENC1EGCUT->DrawClone("AEP");
+  c4->SaveAs("results/emcal_electrons_gCEVsSimgaENC1EGCUT.png");
+  c4->SaveAs("results/emcal_electrons_gCEVsSigmaENC1EGCUT.pdf");
+
+  // Event Counts
+  auto nevents_thrown      = d1.Count();
+  auto nevents_cluster1    = d2.Count();
+  auto nevents_cluster1cut = d3.Count();
+
+  std::cout << "Number of Events: " << (*nevents_cluster1) << " / " << (*nevents_thrown) << " = single cluster events/all \n";
+  std::cout << "Number of Events: " << (*nevents_cluster1cut) << " / " << (*nevents_thrown) << " = single cluster events that passed edge cut/all \n";
+}