Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LivePlot-Energy-EMC.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file LivePlot-Energy-EMC.C
1 #include <iostream>
2 #include <cstdlib>
3 #include <memory>
4 #include <sstream>
5 #include <string>
6 #include <cstring>
7 
8 #include "TFile.h"
9 #include "TTree.h"
10 #include "TCanvas.h"
11 #include "TLegend.h"
12 
13 #define NELEMS(arr) (sizeof(arr)/sizeof(arr[0]))
14 
15 /* Directory where data is stored for plots */
16 const char *const data_directory =
17  "/sphenix/user/gregtom3/data/Summer2018/ECAL_energy_studies";
18 /* The energy levels we have in GeV */
19 const static int energy_levels[] = { 1, 2, 5, 10, 20 };
20 
21 /* The particle types we have */
22 enum particle_type { electron, pion };
23 /* The detectors we have */
24 enum detector { cemc, eemc, femc };
25 
26 TTree *load_tree(const char *const file_name, const char *const tree_name);
27 void fill_histogram(TH1F * const h, TTree * const t, const Float_t min_value,
28  const bool normalize);
29 void display_histogram(TH1F * const h_pion, TH1F * const h_electron,
30  const char *const title, const char *const label);
31 char *generate_name(const particle_type p, const int particle_energy_gev,
32  const detector d);
33 char *generate_file_path(const particle_type p, const int particle_energy_gev,
34  const detector d);
35 char *generate_title(const int particle_energy_gev, const detector d);
36 char *generate_label(const int particle_energy_gev, const detector d);
37 
38 void LivePlot_Energy_EMC()
39 {
40 
41  /*
42  * sPHENIX Style
43  */
44 
45  gROOT->LoadMacro
46  ("/sphenix/user/gregtom3/SBU/research/macros/macros/sPHENIXStyle/sPhenixStyle.C");
47  SetsPhenixStyle();
48 
49  /*
50  * Base Histogram (Recreated from Matching Plots)
51  */
52 
53  TH1F *h_base = new TH1F("h_base", "", 25, 0.0, 2.5);
54  TH1F *h_base_e = (TH1F *) h_base->Clone();
55  TH1F *h_base_p = (TH1F *) h_base->Clone();
56  h_base_e->SetLineColor(kRed);
57  h_base_p->SetLineColor(kBlue);
58 
59  /* iterate through energy levels and create plots for Pions and Electrons
60  * in the CEMC and EEMC detectors */
61  for (size_t i = 0; i < NELEMS(energy_levels); ++i) {
62  /* CEMC */
63  TH1F *const h_pion_cemc = h_base_p->Clone();
64  TH1F *const h_electron_cemc = h_base_e->Clone();
65  h_pion_cemc->SetName(generate_name
66  (pion, energy_levels[i], cemc));
67  h_electron_cemc->SetName(generate_name
68  (electron, energy_levels[i], cemc));
69 
70  TTree *const t_pion_cemc =
71  load_tree(generate_file_path(pion, energy_levels[i], cemc),
72  "ntp_cluster");
73  TTree *const t_electron_cemc =
74  load_tree(generate_file_path
75  (electron, energy_levels[i], cemc),
76  "ntp_cluster");
77 
78  fill_histogram(h_pion_cemc, t_pion_cemc, 0.3, true);
79  fill_histogram(h_electron_cemc, t_electron_cemc, 0.3, true);
80  display_histogram(h_pion_cemc, h_electron_cemc,
81  generate_title(energy_levels[i], cemc),
82  generate_label(energy_levels[i], cemc));
83 
84  /* EEMC */
85  TH1F *const h_pion_eemc = h_base_p->Clone();
86  TH1F *const h_electron_eemc = h_base_e->Clone();
87  h_pion_eemc->SetName(generate_name
88  (pion, energy_levels[i], eemc));
89  h_electron_eemc->SetName(generate_name
90  (electron, energy_levels[i], eemc));
91 
92  TTree *const t_pion_eemc =
93  load_tree(generate_file_path(pion, energy_levels[i], eemc),
94  "ntp_cluster");
95  TTree *const t_electron_eemc =
96  load_tree(generate_file_path
97  (electron, energy_levels[i], eemc),
98  "ntp_cluster");
99 
100  fill_histogram(h_pion_eemc, t_pion_eemc, 0.3, true);
101  fill_histogram(h_electron_eemc, t_electron_eemc, 0.3, true);
102  display_histogram(h_pion_eemc, h_electron_eemc,
103  generate_title(energy_levels[i], eemc),
104  generate_label(energy_levels[i], eemc));
105 
106  /* FEMC */
107  TH1F *const h_pion_femc = h_base_p->Clone();
108  TH1F *const h_electron_femc = h_base_e->Clone();
109  h_pion_femc->SetName(generate_name
110  (pion, energy_levels[i], femc));
111  h_electron_femc->SetName(generate_name
112  (electron, energy_levels[i], femc));
113 
114  TTree *const t_pion_femc =
115  load_tree(generate_file_path(pion, energy_levels[i], femc),
116  "ntp_cluster");
117  TTree *const t_electron_femc =
118  load_tree(generate_file_path
119  (electron, energy_levels[i], femc),
120  "ntp_cluster");
121 
122  fill_histogram(h_pion_femc, t_pion_femc, 0.3, true);
123  fill_histogram(h_electron_femc, t_electron_femc, 0.3, true);
124  display_histogram(h_pion_femc, h_electron_femc,
125  generate_title(energy_levels[i], femc),
126  generate_label(energy_levels[i], femc));
127 
128  }
129 
130 }
131 
132 TTree *load_tree(const char *const file_name, const char *const tree_name)
133 {
134  return (TTree *) (new TFile(file_name, "READ"))->Get(tree_name);
135 }
136 
137 /*
138  * This function, when fed a histogram, tree,
139  * a floating point min_value variable, and a boolean, will fill each entry
140  * inside the specificed branch into the histogram, so long as each entry
141  * is a floating point number GREATER than the min_value. If normalize is
142  * set to 'true', the histogram will be normalize based on its number of
143  * entries. The axes titles are furthermore assumed to be generic and have
144  * been already set.
145  */
146 void fill_histogram(TH1F * const h, TTree * const t, const Float_t min_value,
147  const bool normalize)
148 {
149  Float_t measured_energy;
150  Float_t true_energy;
151  t->SetBranchAddress("e", &measured_energy);
152  t->SetBranchAddress("ge", &true_energy);
153  Int_t nentries = Int_t(t->GetEntries());
154 
155  for (Int_t i = 0; i < nentries; ++i) {
156  if (t->LoadTree(i) < 0)
157  break;
158 
159  t->GetEntry(i);
160  if (measured_energy > min_value && true_energy > 0.1)
161  h->Fill(measured_energy / true_energy);
162  }
163  if (normalize)
164  h->Scale(1 / h->GetEntries());
165 
166  h->SetXTitle("em_cluster_e / true_e");
167  h->SetYTitle("entries / #scale[0.5]{#sum} entries ");
168 }
169 
170 void display_histogram(TH1F * const h_pion, TH1F * const h_electron,
171  const char *const title, const char *const label)
172 {
173  TCanvas *cPNG = new TCanvas(label, title, 600, 400);
174 
175  h_pion->GetYaxis()->SetRangeUser(0.0001, 1);
176  h_electron->GetYaxis()->SetRangeUser(0.0001, 1);
177  gPad->SetLogy();
178  h_pion->Draw();
179  h_electron->Draw("SAME");
180  gPad->RedrawAxis();
181 
182  auto legend = new TLegend(0.70, 0.9, 0.95, 0.65, label);
183  legend->AddEntry(h_pion, "Pions", "l");
184  legend->AddEntry(h_electron, "Electrons", "l");
185  legend->Draw();
186 }
187 
188 char *strdup(const char *s)
189 {
190  char *const t = new char[strlen(s) + 1];
191  return strcpy(t, s);
192 }
193 
194 char *generate_name(const particle_type p, const int particle_energy_gev,
195  const detector d)
196 {
197  std::stringstream name;
198  name << "Histogram";
199  switch (p) {
200  case electron:
201  name << "_Electron";
202  break;
203  case pion:
204  name << "_Pion";
205  break;
206  }
207 
208  name << "_" << particle_energy_gev << "GeV";
209  switch (d) {
210  case cemc:
211  name << "_CEMC";
212  break;
213  case eemc:
214  name << "_EEMC";
215  break;
216  }
217 
218  return strdup(name.str().c_str());
219 }
220 
221 char *generate_file_path(const particle_type p,
222  const int particle_energy_gev, const detector d)
223 {
224  std::stringstream path;
225  path << data_directory;
226 
227  switch (p) {
228  case electron:
229  path << "/Electrons/Electrons";
230  break;
231  case pion:
232  path << "/Pions/Pions";
233  break;
234  }
235 
236  path << particle_energy_gev;
237  switch (d) {
238  case cemc:
239  path << "C";
240  break;
241  case eemc:
242  path << "E";
243  break;
244  case femc:
245  path << "F";
246  break;
247  }
248 
249  path << ".root";
250 
251  return strdup(path.str().c_str());
252 }
253 
254 char *generate_title(const int particle_energy_gev, const detector d)
255 {
256  return generate_label(particle_energy_gev, d);
257 }
258 
259 char *generate_label(const int particle_energy_gev, const detector d)
260 {
261  std::stringstream label;
262 
263  switch (d) {
264  case cemc:
265  label << "CEMC ";
266  break;
267  case eemc:
268  label << "EEMC ";
269  break;
270  case femc:
271  label << "FEMC ";
272  break;
273  }
274  label << particle_energy_gev << "GeV";
275 
276  return strdup(label.str().c_str());
277 }