Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
eic_sphenix_test_smearing.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file eic_sphenix_test_smearing.C
1 /* STL includes */
2 #include <iostream>
3 
4 /* ROOT includes */
5 #include <TROOT.h>
6 //#include <gSystem.h>
7 #include <TFile.h>
8 #include <TTree.h>
9 #include <TH2F.h>
10 #include <THnSparse.h>
11 
12 /* eicsmear includes */
13 //#include <eicsmear/erhic/EventPythia.h>
14 //#include <eicsmear/erhic/ParticleMC.h>
15 
16 using namespace std;
17 
18 int
19 eic_sphenix_test_smearing( TString filename_output,
20  TString filename_mc,
21  TString filename_mc_smeared = "",
22  bool debug = false )
23 {
24  /* Uncomment this line when running without compilation */
25  gSystem->Load("libeicsmear");
26 
27  /* Open input files. */
28  TFile *file_mc = new TFile(filename_mc, "OPEN");
29  TFile *file_mc_smeared = new TFile(filename_mc_smeared, "OPEN");
30 
31  /* Get trees from files. */
32  TTree *tree = (TTree*)file_mc->Get("EICTree");
33  TTree *tree_smeared = (TTree*)file_mc_smeared->Get("Smeared");
34 
35  /* Output file. */
36  //TFile *file_out = new TFile(filename_output, "RECREATE");
37 
38  /* Add friend to match branches in trees. */
39  tree->AddFriend(tree_smeared);
40 
41  erhic::EventPythia *event = NULL;
42  Smear::Event *eventS = NULL;
43 
44  tree->SetBranchAddress("event", &event);
45  tree->SetBranchAddress("eventS", &eventS);
46 
47  /* Create histogram */
48  TH2F* h_e_smeared_vs_eta = new TH2F("h_e_smeared_vs_eta","Energy Smeared vs True Pseudorapidity",100,-5,5,70,0,35);
49  h_e_smeared_vs_eta->GetXaxis()->SetTitle("#eta_{true}");
50  h_e_smeared_vs_eta->GetYaxis()->SetTitle("E_{smeared} (GeV)");
51 
52  TH2F* h_e_smeared_vs_true = new TH2F("h_e_smeared_vs_true","Energy Smeared vs True",60,0,30,70,0,35);
53  h_e_smeared_vs_true->GetXaxis()->SetTitle("E_{true} (GeV)");
54  h_e_smeared_vs_true->GetYaxis()->SetTitle("E_{smeared} (GeV)");
55 
56  TH1F* h_eta = new TH1F("h_eta", ";#eta;dN/d#eta", 100, -5, 5);
57  TH1F* h_eta_accept = (TH1F*)h_eta->Clone("h_eta_accept");
58  h_eta_accept->SetLineColor(kGreen+4);
59 
60  TH1F* h_e_eref_true = new TH1F("h_e_eref_true","True reference energy",300,0,30);
61  h_e_eref_true->GetXaxis()->SetTitle("E_{true} (GeV)");
62  h_e_eref_true->GetYaxis()->SetTitle("# entries");
63  h_e_eref_true->SetLineColor(kRed);
64 
65  TH1F* h_e_eref_smeared = new TH1F("h_e_eref_smeared","Smeared reference energy",300,0,30);
66  h_e_eref_smeared->GetXaxis()->SetTitle("E_{smeared} (GeV)");
67  h_e_eref_smeared->GetYaxis()->SetTitle("# entries");
68  h_e_eref_smeared->SetLineColor(kBlue);
69 
70  /* Loop over all events in tree. */
71  unsigned max_event = tree->GetEntries();
72 
73  for ( unsigned ievent = 0; ievent < max_event; ievent++ )
74  {
75  if ( ievent%1000 == 0 )
76  cout << "Processing event " << ievent << endl;
77 
78  /* load event */
79  tree->GetEntry(ievent);
80 
81  /* Cut on EVENT kinematics */
82  float y = event->GetTrueY();
83  if ( y > 0.99 || y < 0.01 )
84  continue;
85 
86  float energy = event->ScatteredLepton()->GetE();
87  float energy_smeared = eventS->ScatteredLepton()->GetE();
88 
89  float eta = event->ScatteredLepton()->GetEta();
90 
91  /* Fill histograms */
92  h_e_smeared_vs_eta->Fill(eta,energy_smeared);
93  h_e_smeared_vs_true->Fill(energy,energy_smeared);
94 
95  h_eta->Fill(eta);
96  if ( energy_smeared > 0 )
97  h_eta_accept->Fill( eta );
98 
99  /* Fill histograms if truth energy within range around reference energy */
100  float eref = 19.05;
101  float erange = 0.1;
102  if ( energy > (eref-erange/2.) && energy < (eref+erange/2.) )
103  {
104  h_e_eref_true->Fill(energy);
105  h_e_eref_smeared->Fill(energy_smeared);
106  }
107 
108  } // end loop over events
109 
110  float underflow = tree->GetEntry(0);
111  float overflow = tree->GetEntry(max_event + 1);
112  cout << "underflow: " << underflow << endl;
113  cout << "overflow: " << overflow << endl;
114 
115  /* Write histograms. */
116  //h_e_smeared_vs_true->Write();
117  //h_e_smeared_vs_eta->Write();
118 
119  /*draw histograms */
120  TCanvas *c1 = new TCanvas();
121  h_e_smeared_vs_true->DrawClone("COLZ");
122  gPad->RedrawAxis();
123 
124  TCanvas *c2 = new TCanvas();
125  h_e_smeared_vs_eta->DrawClone("COLZ");
126  gPad->RedrawAxis();
127 
128  TCanvas *c3 = new TCanvas();
129  h_e_eref_smeared->Draw();
130  h_e_eref_smeared->Fit("gaus");
131  h_e_eref_true->Draw("sames");
132  gPad->RedrawAxis();
133 
134  TCanvas *c4 = new TCanvas();
135  h_eta_accept->Divide(h_eta);
136  h_eta_accept->Draw();
137  gPad->RedrawAxis();
138 
139  /* Close output file. */
140  //file_out->Close();
141 
142  return 0;
143 }