macros/scattered-electron_analysis/LivePlot-Energy-EMC.C

0001 #include <iostream>
0002 #include <cstdlib>
0003 #include <memory>
0004 #include <sstream>
0005 #include <string>
0006 #include <cstring>
0007
0008 #include "TFile.h"
0009 #include "TTree.h"
0010 #include "TCanvas.h"
0011 #include "TLegend.h"
0012
0013 #define NELEMS(arr) (sizeof(arr)/sizeof(arr[0]))
0014
0015 /* Directory where data is stored for plots */
0016 const char *const data_directory =
0017     "/sphenix/user/gregtom3/data/Summer2018/ECAL_energy_studies";
0018 /* The energy levels we have in GeV */
0019 const static int energy_levels[] = { 1, 2, 5, 10, 20 };
0020
0021 /* The particle types we have */
0022 enum particle_type { electron, pion };
0023 /* The detectors we have */
0024 enum detector { cemc, eemc, femc };
0025
0026 TTree *load_tree(const char *const file_name, const char *const tree_name);
0027 void fill_histogram(TH1F * const h, TTree * const t, const Float_t min_value,
0028             const bool normalize);
0029 void display_histogram(TH1F * const h_pion, TH1F * const h_electron,
0030                const char *const title, const char *const label);
0031 char *generate_name(const particle_type p, const int particle_energy_gev,
0032             const detector d);
0033 char *generate_file_path(const particle_type p, const int particle_energy_gev,
0034              const detector d);
0035 char *generate_title(const int particle_energy_gev, const detector d);
0036 char *generate_label(const int particle_energy_gev, const detector d);
0037
0038 void LivePlot_Energy_EMC()
0039 {
0040
0041     /*
0042      * sPHENIX Style
0043      */
0044
0045     gROOT->LoadMacro
0046         ("/sphenix/user/gregtom3/SBU/research/macros/macros/sPHENIXStyle/sPhenixStyle.C");
0047     SetsPhenixStyle();
0048
0049     /*
0050      * Base Histogram (Recreated from Matching Plots)
0051      */
0052
0053     TH1F *h_base = new TH1F("h_base", "", 25, 0.0, 2.5);
0054     TH1F *h_base_e = (TH1F *) h_base->Clone();
0055     TH1F *h_base_p = (TH1F *) h_base->Clone();
0056     h_base_e->SetLineColor(kRed);
0057     h_base_p->SetLineColor(kBlue);
0058
0059     /* iterate through energy levels and create plots for Pions and Electrons
0060      * in the CEMC and EEMC detectors */
0061     for (size_t i = 0; i < NELEMS(energy_levels); ++i) {
0062         /* CEMC */
0063         TH1F *const h_pion_cemc = h_base_p->Clone();
0064         TH1F *const h_electron_cemc = h_base_e->Clone();
0065         h_pion_cemc->SetName(generate_name
0066                      (pion, energy_levels[i], cemc));
0067         h_electron_cemc->SetName(generate_name
0068                      (electron, energy_levels[i], cemc));
0069
0070         TTree *const t_pion_cemc =
0071             load_tree(generate_file_path(pion, energy_levels[i], cemc),
0072                   "ntp_cluster");
0073         TTree *const t_electron_cemc =
0074             load_tree(generate_file_path
0075                   (electron, energy_levels[i], cemc),
0076                   "ntp_cluster");
0077
0078         fill_histogram(h_pion_cemc, t_pion_cemc, 0.3, true);
0079         fill_histogram(h_electron_cemc, t_electron_cemc, 0.3, true);
0080         display_histogram(h_pion_cemc, h_electron_cemc,
0081                   generate_title(energy_levels[i], cemc),
0082                   generate_label(energy_levels[i], cemc));
0083
0084         /* EEMC */
0085         TH1F *const h_pion_eemc = h_base_p->Clone();
0086         TH1F *const h_electron_eemc = h_base_e->Clone();
0087         h_pion_eemc->SetName(generate_name
0088                      (pion, energy_levels[i], eemc));
0089         h_electron_eemc->SetName(generate_name
0090                      (electron, energy_levels[i], eemc));
0091
0092         TTree *const t_pion_eemc =
0093             load_tree(generate_file_path(pion, energy_levels[i], eemc),
0094                   "ntp_cluster");
0095         TTree *const t_electron_eemc =
0096             load_tree(generate_file_path
0097                   (electron, energy_levels[i], eemc),
0098                   "ntp_cluster");
0099
0100         fill_histogram(h_pion_eemc, t_pion_eemc, 0.3, true);
0101         fill_histogram(h_electron_eemc, t_electron_eemc, 0.3, true);
0102         display_histogram(h_pion_eemc, h_electron_eemc,
0103                   generate_title(energy_levels[i], eemc),
0104                   generate_label(energy_levels[i], eemc));
0105
0106         /* FEMC */
0107         TH1F *const h_pion_femc = h_base_p->Clone();
0108         TH1F *const h_electron_femc = h_base_e->Clone();
0109         h_pion_femc->SetName(generate_name
0110                      (pion, energy_levels[i], femc));
0111         h_electron_femc->SetName(generate_name
0112                      (electron, energy_levels[i], femc));
0113
0114         TTree *const t_pion_femc =
0115             load_tree(generate_file_path(pion, energy_levels[i], femc),
0116                   "ntp_cluster");
0117         TTree *const t_electron_femc =
0118             load_tree(generate_file_path
0119                   (electron, energy_levels[i], femc),
0120                   "ntp_cluster");
0121
0122         fill_histogram(h_pion_femc, t_pion_femc, 0.3, true);
0123         fill_histogram(h_electron_femc, t_electron_femc, 0.3, true);
0124         display_histogram(h_pion_femc, h_electron_femc,
0125                   generate_title(energy_levels[i], femc),
0126                   generate_label(energy_levels[i], femc));
0127
0128     }
0129
0130 }
0131
0132 TTree *load_tree(const char *const file_name, const char *const tree_name)
0133 {
0134     return (TTree *) (new TFile(file_name, "READ"))->Get(tree_name);
0135 }
0136
0137 /*
0138  * This function, when fed a histogram, tree,
0139  * a floating point min_value variable, and a boolean, will fill each entry
0140  * inside the specificed branch into the histogram, so long as each entry
0141  * is a floating point number GREATER than the min_value. If normalize is
0142  * set to 'true', the histogram will be normalize based on its number of
0143  * entries. The axes titles are furthermore assumed to be generic and have
0144  * been already set.
0145  */
0146 void fill_histogram(TH1F * const h, TTree * const t, const Float_t min_value,
0147             const bool normalize)
0148 {
0149     Float_t measured_energy;
0150     Float_t true_energy;
0151     t->SetBranchAddress("e", &measured_energy);
0152     t->SetBranchAddress("ge", &true_energy);
0153     Int_t nentries = Int_t(t->GetEntries());
0154
0155     for (Int_t i = 0; i < nentries; ++i) {
0156         if (t->LoadTree(i) < 0)
0157             break;
0158
0159         t->GetEntry(i);
0160         if (measured_energy > min_value && true_energy > 0.1)
0161             h->Fill(measured_energy / true_energy);
0162     }
0163     if (normalize)
0164         h->Scale(1 / h->GetEntries());
0165
0166     h->SetXTitle("em_cluster_e / true_e");
0167     h->SetYTitle("entries / #scale[0.5]{#sum} entries      ");
0168 }
0169
0170 void display_histogram(TH1F * const h_pion, TH1F * const h_electron,
0171                const char *const title, const char *const label)
0172 {
0173     TCanvas *cPNG = new TCanvas(label, title, 600, 400);
0174
0175     h_pion->GetYaxis()->SetRangeUser(0.0001, 1);
0176     h_electron->GetYaxis()->SetRangeUser(0.0001, 1);
0177     gPad->SetLogy();
0178     h_pion->Draw();
0179     h_electron->Draw("SAME");
0180     gPad->RedrawAxis();
0181
0182     auto legend = new TLegend(0.70, 0.9, 0.95, 0.65, label);
0183     legend->AddEntry(h_pion, "Pions", "l");
0184     legend->AddEntry(h_electron, "Electrons", "l");
0185     legend->Draw();
0186 }
0187
0188 char *strdup(const char *s)
0189 {
0190     char *const t = new char[strlen(s) + 1];
0191     return strcpy(t, s);
0192 }
0193
0194 char *generate_name(const particle_type p, const int particle_energy_gev,
0195             const detector d)
0196 {
0197     std::stringstream name;
0198     name << "Histogram";
0199     switch (p) {
0200     case electron:
0201         name << "_Electron";
0202         break;
0203     case pion:
0204         name << "_Pion";
0205         break;
0206     }
0207
0208     name << "_" << particle_energy_gev << "GeV";
0209     switch (d) {
0210     case cemc:
0211         name << "_CEMC";
0212         break;
0213     case eemc:
0214         name << "_EEMC";
0215         break;
0216     }
0217
0218     return strdup(name.str().c_str());
0219 }
0220
0221 char *generate_file_path(const particle_type p,
0222              const int particle_energy_gev, const detector d)
0223 {
0224     std::stringstream path;
0225     path << data_directory;
0226
0227     switch (p) {
0228     case electron:
0229         path << "/Electrons/Electrons";
0230         break;
0231     case pion:
0232         path << "/Pions/Pions";
0233         break;
0234     }
0235
0236     path << particle_energy_gev;
0237     switch (d) {
0238     case cemc:
0239         path << "C";
0240         break;
0241     case eemc:
0242         path << "E";
0243         break;
0244     case femc:
0245         path << "F";
0246         break;
0247     }
0248
0249     path << ".root";
0250
0251     return strdup(path.str().c_str());
0252 }
0253
0254 char *generate_title(const int particle_energy_gev, const detector d)
0255 {
0256     return generate_label(particle_energy_gev, d);
0257 }
0258
0259 char *generate_label(const int particle_energy_gev, const detector d)
0260 {
0261     std::stringstream label;
0262
0263     switch (d) {
0264     case cemc:
0265         label << "CEMC ";
0266         break;
0267     case eemc:
0268         label << "EEMC ";
0269         break;
0270     case femc:
0271         label << "FEMC ";
0272         break;
0273     }
0274     label << particle_energy_gev << "GeV";
0275
0276     return strdup(label.str().c_str());
0277 }