macros/scattered-electron_analysis/Plot-Measured_Energy-EMC.C

0001 #include <iostream>
0002 #include <cstdlib>
0003 #include <memory>
0004 #include <sstream>
0005 #include <string>
0006 #include <cstring>
0007 #include <vector>
0008 #include <cassert>
0009
0010 #include "TFile.h"
0011 #include "TTree.h"
0012 #include "TCanvas.h"
0013 #include "TLegend.h"
0014 #include "TROOT.h"
0015 #include "TH1F.h"
0016 #include "TColor.h"
0017 #include "TImage.h"
0018
0019 /*
0020  * sPHENIX Style
0021  */
0022 #include "/sphenix/user/gregtom3/SBU/research/macros/macros/sPHENIXStyle/sPhenixStyle.C"
0023
0024 #define NELEMS(arr) (sizeof(arr)/sizeof(arr[0]))
0025
0026 /* The particle types we have */
0027 enum particle_type { electron, pion };
0028 /* The detectors we have */
0029 enum detector { cemc, eemc, femc };
0030
0031 TTree *load_tree(const char *const file_name, const char *const tree_name);
0032 void draw_histogram(TH1F * const h, const particle_type p,
0033             const int energy_level_gev, const detector d);
0034 void fill_histogram(TH1F * const h, TTree * const t, const Float_t min_value);
0035 char *generate_histogram_name(const particle_type p,
0036                   const int particle_energy_gev, const detector d);
0037 char *generate_file_path(const particle_type p, const int particle_energy_gev,
0038              const detector d);
0039 char *generate_save_file_path(const particle_type p);
0040 char *generate_canvas_name(const particle_type p);
0041 char *generate_canvas_title(const particle_type p);
0042 char *generate_legend_header(const particle_type p, const detector d);
0043 char *generate_legend_entry_label(const int particle_energy_gev);
0044
0045 /* Directory where data is stored for plots */
0046 const char *const data_directory =
0047     "/sphenix/user/gregtom3/data/Summer2018/ECAL_energy_studies";
0048 /* The energy levels we have in GeV */
0049 const static int energy_levels[] = { 1, 2, 5, 10, 20 };
0050
0051 /* How to change lines for each energy level */
0052 const static int color_offsets[] = { -10, -9, -7, -4, 0 };
0053
0054 /* Detectors that we want to plot */
0055 const static enum detector detectors[] = { cemc, eemc, femc };
0056
0057 /* Particles we want to plot */
0058 const static enum particle_type particles[] = { pion, electron };
0059
0060 /* The color of the plot corresponding to each particle */
0061 const static Color_t plot_colors[] = { kBlue, kRed, kGreen, kOrange, kViolet };
0062
0063 /* These should not be necessary, but root is buggy */
0064 const int nenergy_levels = NELEMS(energy_levels);
0065 const int ndetectors = NELEMS(detectors);
0066 const int nparticles = NELEMS(particles);
0067
0068 void Plot_Measured_Energy_EMC()
0069 {
0070     /* sanity check */
0071     assert(NELEMS(energy_levels) <= NELEMS(plot_colors));
0072
0073     SetsPhenixStyle();
0074     gROOT->SetBatch(kTRUE);
0075
0076     /*
0077      * Base Histogram (Recreated from Matching Plots)
0078      */
0079     TH1F *const h_base = new TH1F("h_base", "", 100, 0.0, 25);
0080
0081     std::vector < TCanvas * >particle_canvases;
0082     for (int i = 0; i < nparticles; ++i) {
0083         TCanvas *const c =
0084             new TCanvas(generate_canvas_name(particles[i]),
0085                 generate_canvas_title(particles[i]),
0086                 gStyle->GetCanvasDefW() * ndetectors,
0087                 gStyle->GetCanvasDefH());
0088         c->Divide(ndetectors, 1);
0089         particle_canvases.push_back(c);
0090     }
0091
0092     /* iterate through energy levels, detectors, and particles and create plots for each */
0093     for (int e = 0; e < nenergy_levels; ++e)
0094         for (int d = 0; d < ndetectors; ++d)
0095             for (int p = 0; p < nparticles; ++p) {
0096                 /* +1 to account for 1-based indexing */
0097                 particle_canvases[p]->cd(d + 1);
0098                 TH1F *const h = (TH1F *) h_base->Clone();
0099                 h->SetLineColor(plot_colors[e]);
0100                 draw_histogram(h, particles[p],
0101                            energy_levels[e], detectors[d]);
0102             }
0103
0104     /* add legends */
0105     for (int d = 0; d < ndetectors; ++d)
0106         for (int p = 0; p < nparticles; ++p) {
0107             particle_canvases[p]->cd(d + 1);
0108             gPad->RedrawAxis();
0109             if (d == 0) {
0110                 TLegend *const l =
0111                     new TLegend(0.60, .90, .9, 0.625,
0112                         generate_legend_header(particles
0113                                        [p],
0114                                        detectors
0115                                        [d]));
0116                 l->SetTextSize(0.05);
0117                 for (int e = 0; e < nenergy_levels; ++e)
0118                     l->AddEntry(generate_histogram_name
0119                             (particles[p],
0120                              energy_levels[e],
0121                              detectors[d]),
0122                             generate_legend_entry_label
0123                             (energy_levels[e]), "l");
0124
0125                 l->Draw();
0126             } else {
0127                 TLegend *const l =
0128                     new TLegend(0.60, .90, .9, 0.85,
0129                         generate_legend_header(particles
0130                                        [p],
0131                                        detectors
0132                                        [d]));
0133                 l->SetTextSize(0.05);
0134                 l->Draw();
0135             }
0136
0137             gPad->SetLogy();
0138         }
0139
0140     /* save canvases to PNG */
0141     for (int p = 0; p < nparticles; ++p) {
0142         TImage *const img = TImage::Create();
0143         img->FromPad(particle_canvases[p]);
0144         char *const path = generate_save_file_path(particles[p]);
0145         img->WriteImage(path);
0146         delete img;
0147     }
0148 }
0149
0150 void draw_histogram(TH1F * const h, const particle_type p,
0151             const int energy_level_gev, const detector d)
0152 {
0153     h->SetName(generate_histogram_name(p, energy_level_gev, d));
0154
0155     TTree *const t = load_tree(generate_file_path(p, energy_level_gev, d),
0156                    "ntp_cluster");
0157
0158     fill_histogram(h, t, 0.3);
0159     h->Scale(1 / h->GetEntries());
0160     h->GetYaxis()->SetRangeUser(0.0001, 10);
0161     h->SetXTitle("E_{cluster} (GeV)");
0162     h->SetYTitle("entries / #scale[0.5]{#sum} entries      ");
0163
0164     h->Draw("SAME");
0165 }
0166
0167 TTree *load_tree(const char *const file_name, const char *const tree_name)
0168 {
0169     return (TTree *) (new TFile(file_name, "READ"))->Get(tree_name);
0170 }
0171
0172 /*
0173  * This function, when fed a histogram, tree,
0174  * a floating point min_value variable, and a boolean, will fill each entry
0175  * inside the specificed branch into the histogram, so long as each entry
0176  * is a floating point number GREATER than the min_value. If normalize is
0177  * set to 'true', the histogram will be normalize based on its number of
0178  * entries. The axes titles are furthermore assumed to be generic and have
0179  * been already set.
0180  */
0181 void fill_histogram(TH1F * const h, TTree * const t, const Float_t min_value)
0182 {
0183     Float_t measured_energy;
0184     Float_t true_energy;
0185     t->SetBranchAddress("e", &measured_energy);
0186     t->SetBranchAddress("ge", &true_energy);
0187     Int_t nentries = Int_t(t->GetEntries());
0188
0189     for (Int_t i = 0; i < nentries; ++i) {
0190         if (t->LoadTree(i) < 0)
0191             break;
0192
0193         t->GetEntry(i);
0194         if (measured_energy > min_value && true_energy > 0.1)
0195             h->Fill(measured_energy);
0196     }
0197 }
0198
0199 char *strdup(const char *s)
0200 {
0201     char *const t = new char[strlen(s) + 1];
0202     return strcpy(t, s);
0203 }
0204
0205 char *generate_histogram_name(const particle_type p,
0206                   const int particle_energy_gev, const detector d)
0207 {
0208     std::stringstream name;
0209     name << "Histogram";
0210     switch (p) {
0211     case electron:
0212         name << "_Electron";
0213         break;
0214     case pion:
0215         name << "_Pion";
0216         break;
0217     }
0218
0219     name << "_" << particle_energy_gev << "GeV";
0220     switch (d) {
0221     case cemc:
0222         name << "_CEMC";
0223         break;
0224     case eemc:
0225         name << "_EEMC";
0226         break;
0227     }
0228
0229     return strdup(name.str().c_str());
0230 }
0231
0232 char *generate_file_path(const particle_type p,
0233              const int particle_energy_gev, const detector d)
0234 {
0235     std::stringstream path;
0236     path << data_directory;
0237
0238     switch (p) {
0239     case electron:
0240         path << "/Electrons/Electrons";
0241         break;
0242     case pion:
0243         path << "/Pions/Pions";
0244         break;
0245     }
0246
0247     path << particle_energy_gev;
0248     switch (d) {
0249     case cemc:
0250         path << "C";
0251         break;
0252     case eemc:
0253         path << "E";
0254         break;
0255     case femc:
0256         path << "F";
0257         break;
0258     }
0259
0260     path << ".root";
0261
0262     return strdup(path.str().c_str());
0263 }
0264
0265 char *generate_save_file_path(const particle_type p)
0266 {
0267     std::stringstream name;
0268
0269     switch (p) {
0270     case pion:
0271         name << "Pion";
0272         break;
0273     case electron:
0274         name << "Electron";
0275         break;
0276     }
0277
0278     name << "-Measured_Energy";
0279     for (int d = 0; d < ndetectors; ++d)
0280         switch (detectors[d]) {
0281         case cemc:
0282             name << "-CEMC";
0283             break;
0284         case eemc:
0285             name << "-EEMC";
0286             break;
0287         case femc:
0288             name << "-FEMC";
0289             break;
0290         }
0291
0292     name << ".png";
0293
0294     return strdup(name.str().c_str());
0295 }
0296
0297 char *generate_canvas_name(const particle_type p)
0298 {
0299     std::stringstream canvas_name;
0300
0301     switch (p) {
0302     case pion:
0303         canvas_name << "Pion";
0304         break;
0305     case electron:
0306         canvas_name << "Electron";
0307         break;
0308     }
0309
0310     return strdup(canvas_name.str().c_str());
0311 }
0312
0313 char *generate_canvas_title(const particle_type p)
0314 {
0315     std::stringstream canvas_title;
0316
0317     switch (p) {
0318     case pion:
0319         canvas_title << "Pion";
0320         break;
0321     case electron:
0322         canvas_title << "Electron";
0323         break;
0324     }
0325
0326     return strdup(canvas_title.str().c_str());
0327 }
0328
0329 char *generate_legend_header(const particle_type p, const detector d)
0330 {
0331     std::stringstream legend_header;
0332
0333     switch (d) {
0334     case cemc:
0335         legend_header << "CEMC ";
0336         break;
0337     case eemc:
0338         legend_header << "EEMC ";
0339         break;
0340     case femc:
0341         legend_header << "FEMC ";
0342         break;
0343     }
0344
0345     switch (p) {
0346     case pion:
0347         legend_header << "Pions ";
0348         break;
0349     case electron:
0350         legend_header << "Electrons ";
0351         break;
0352     }
0353
0354     return strdup(legend_header.str().c_str());
0355 }
0356
0357 char *generate_legend_entry_label(const int particle_energy_gev)
0358 {
0359     std::stringstream legend_entry;
0360     legend_entry << particle_energy_gev << " GeV";
0361
0362     return strdup(legend_entry.str().c_str());
0363 }