Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-06 08:12:55

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 true_energy,
0028             const Float_t min_value, const bool normalize);
0029 void histogram_to_png(TH1F * const h_pion_CEMC, TH1F * const h_electron_CEMC,
0030               TH1F * const h_pion_EEMC, TH1F * const h_electron_EEMC,
0031               TH1F * const h_pion_FEMC, TH1F * const h_electron_FEMC,
0032               const char *const title, const char *const save_file_name,
0033               const char *const cemc_label,
0034               const char *const eemc_label,
0035               const char *const femc_label);
0036 char *generate_name(const particle_type p, const int particle_energy_gev,
0037             const detector d);
0038 char *generate_file_path(const particle_type p, const int particle_energy_gev,
0039              const detector d);
0040 char *generate_save_name(const int particle_energy_gev);
0041 char *generate_title(const int particle_energy_gev);
0042 char *generate_label(const int particle_energy_gev, const detector d);
0043 
0044 void Plot_Measured_Energy_Over_True_Energy_EMC()
0045 {
0046 
0047     /* 
0048      * sPHENIX Style
0049      */
0050 
0051     gROOT->LoadMacro
0052         ("/sphenix/user/gregtom3/SBU/research/macros/macros/sPHENIXStyle/sPhenixStyle.C");
0053     SetsPhenixStyle();
0054 
0055     gROOT->SetBatch(kTRUE);
0056 
0057     /*
0058      * Base Histogram (Recreated from Matching Plots)
0059      */
0060 
0061     TH1F *h_base = new TH1F("h_base", "", 25, 0.0, 2.5);
0062     TH1F *h_base_e = (TH1F *) h_base->Clone();
0063     TH1F *h_base_p = (TH1F *) h_base->Clone();
0064     h_base_e->SetLineColor(kRed);
0065     h_base_p->SetLineColor(kBlue);
0066 
0067     /* iterate through energy levels and create plots for Pions and Electrons
0068      * in the CEMC and EEMC detectors */
0069     for (size_t i = 0; i < NELEMS(energy_levels); ++i) {
0070         /* CEMC */
0071         TH1F *const h_pion_cemc = h_base_p->Clone();
0072         TH1F *const h_electron_cemc = h_base_e->Clone();
0073         h_pion_cemc->SetName(generate_name
0074                      (pion, energy_levels[i], cemc));
0075         h_electron_cemc->SetName(generate_name
0076                      (electron, energy_levels[i], cemc));
0077 
0078         TTree *const t_pion_cemc =
0079             load_tree(generate_file_path(pion, energy_levels[i], cemc),
0080                   "ntp_cluster");
0081         TTree *const t_electron_cemc =
0082             load_tree(generate_file_path
0083                   (electron, energy_levels[i], cemc),
0084                   "ntp_cluster");
0085 
0086         fill_histogram(h_pion_cemc, t_pion_cemc, energy_levels[i], 0.3,
0087                    true);
0088         fill_histogram(h_electron_cemc, t_electron_cemc,
0089                    energy_levels[i], 0.3, true);
0090 
0091         /* EEMC */
0092         TH1F *const h_pion_eemc = h_base_p->Clone();
0093         TH1F *const h_electron_eemc = h_base_e->Clone();
0094         h_pion_eemc->SetName(generate_name
0095                      (pion, energy_levels[i], eemc));
0096         h_electron_eemc->SetName(generate_name
0097                      (electron, energy_levels[i], eemc));
0098 
0099         TTree *const t_pion_eemc =
0100             load_tree(generate_file_path(pion, energy_levels[i], eemc),
0101                   "ntp_cluster");
0102         TTree *const t_electron_eemc =
0103             load_tree(generate_file_path
0104                   (electron, energy_levels[i], eemc),
0105                   "ntp_cluster");
0106 
0107         fill_histogram(h_pion_eemc, t_pion_eemc, energy_levels[i], 0.3,
0108                    true);
0109         fill_histogram(h_electron_eemc, t_electron_eemc,
0110                    energy_levels[i], 0.3, true);
0111 
0112         /* FEMC */
0113         TH1F *const h_pion_femc = h_base_p->Clone();
0114         TH1F *const h_electron_femc = h_base_e->Clone();
0115         h_pion_femc->SetName(generate_name
0116                      (pion, energy_levels[i], femc));
0117         h_electron_femc->SetName(generate_name
0118                      (electron, energy_levels[i], femc));
0119 
0120         TTree *const t_pion_femc =
0121             load_tree(generate_file_path(pion, energy_levels[i], femc),
0122                   "ntp_cluster");
0123         TTree *const t_electron_femc =
0124             load_tree(generate_file_path
0125                   (electron, energy_levels[i], femc),
0126                   "ntp_cluster");
0127 
0128         fill_histogram(h_pion_femc, t_pion_femc, energy_levels[i], 0.3,
0129                    true);
0130         fill_histogram(h_electron_femc, t_electron_femc,
0131                    energy_levels[i], 0.3, true);
0132 
0133         histogram_to_png(h_pion_cemc, h_electron_cemc,
0134                  h_pion_eemc, h_electron_eemc,
0135                  h_pion_femc, h_electron_femc,
0136                  generate_title(energy_levels[i]),
0137                  generate_save_name(energy_levels[i]),
0138                  generate_label(energy_levels[i], cemc),
0139                  generate_label(energy_levels[i], eemc),
0140                  generate_label(energy_levels[i], femc));
0141     }
0142 
0143 }
0144 
0145 TTree *load_tree(const char *const file_name, const char *const tree_name)
0146 {
0147     return (TTree *) (new TFile(file_name, "READ"))->Get(tree_name);
0148 }
0149 
0150 /*
0151  * This function, when fed a histogram, tree,
0152  * a floating point min_value variable, and a boolean, will fill each entry
0153  * inside the specificed branch into the histogram, so long as each entry
0154  * is a floating point number GREATER than the min_value. If normalize is
0155  * set to 'true', the histogram will be normalize based on its number of
0156  * entries. The axes titles are furthermore assumed to be generic and have
0157  * been already set.
0158  */
0159 void fill_histogram(TH1F * const h, TTree * const t, const Float_t true_energy,
0160             const Float_t min_value, const bool normalize)
0161 {
0162     Float_t measured_energy;
0163 //      Float_t true_energy;
0164     t->SetBranchAddress("e", &measured_energy);
0165 //      t->SetBranchAddress("ge", &true_energy);
0166     Float_t true_eta;
0167     t->SetBranchAddress("geta", &true_eta);
0168 
0169     Int_t nentries = Int_t(t->GetEntries());
0170 
0171     for (Int_t i = 0; i < nentries; ++i) {
0172         if (t->LoadTree(i) < 0)
0173             break;
0174 
0175         t->GetEntry(i);
0176         if (((true_eta > -0.5 && true_eta < 0.5)
0177              || (true_eta > -3 && true_eta < -2)
0178              || (true_eta > 2 && true_eta < 3))
0179             && (measured_energy > min_value && true_energy > 0.1))
0180             h->Fill(measured_energy / true_energy);
0181     }
0182     if (normalize)
0183         h->Scale(1 / h->GetEntries());
0184 
0185     h->SetXTitle("E_{cluster} / E_{true}");
0186     h->SetYTitle("entries / #scale[0.5]{#sum} entries      ");
0187 }
0188 
0189 void histogram_to_png(TH1F * const h_pion_CEMC, TH1F * const h_electron_CEMC,
0190               TH1F * const h_pion_EEMC, TH1F * const h_electron_EEMC,
0191               TH1F * const h_pion_FEMC, TH1F * const h_electron_FEMC,
0192               const char *const title, const char *const save_file_name,
0193               const char *const cemc_label,
0194               const char *const eemc_label,
0195               const char *const femc_label)
0196 {
0197     /* 3 for CEMC, EEMC, and FEMC */
0198     const unsigned ndetectors = 3;
0199     TCanvas cPNG("cPNG", title, gStyle->GetCanvasDefW() * 3,
0200              gStyle->GetCanvasDefH());
0201     TImage *img = TImage::Create();
0202 
0203     cPNG.Divide(3, 1);
0204     cPNG.cd(1);
0205     h_pion_CEMC->GetYaxis()->SetRangeUser(0.0001, 1);
0206     h_electron_CEMC->GetYaxis()->SetRangeUser(0.0001, 1);
0207     gPad->SetLogy();
0208     h_pion_CEMC->Draw();
0209     h_electron_CEMC->Draw("SAME");
0210     gPad->RedrawAxis();
0211 
0212     auto cemc_legend = new TLegend(0.70, 0.9, 0.95, 0.65, cemc_label);
0213     cemc_legend->AddEntry(h_pion_CEMC, "Pions", "l");
0214     cemc_legend->AddEntry(h_electron_CEMC, "Electrons", "l");
0215     cemc_legend->Draw();
0216 
0217     cPNG.cd(2);
0218     h_pion_EEMC->GetYaxis()->SetRangeUser(0.0001, 1);
0219     h_electron_EEMC->GetYaxis()->SetRangeUser(0.0001, 1);
0220     gPad->SetLogy();
0221     h_pion_EEMC->Draw();
0222     h_electron_EEMC->Draw("SAME");
0223     gPad->RedrawAxis();
0224 
0225     auto eemc_legend = new TLegend(0.65, 0.9, 0.9, 0.65, eemc_label);
0226     eemc_legend->AddEntry(h_pion_EEMC, "Pions", "l");
0227     eemc_legend->AddEntry(h_electron_EEMC, "Electrons", "l");
0228     eemc_legend->Draw();
0229 
0230     cPNG.cd(3);
0231     h_pion_FEMC->GetYaxis()->SetRangeUser(0.0001, 1);
0232     h_electron_FEMC->GetYaxis()->SetRangeUser(0.0001, 1);
0233     gPad->SetLogy();
0234     h_pion_FEMC->Draw();
0235     h_electron_FEMC->Draw("SAME");
0236     gPad->RedrawAxis();
0237 
0238     auto femc_legend = new TLegend(0.65, 0.9, 0.9, 0.65, femc_label);
0239     femc_legend->AddEntry(h_pion_FEMC, "Pions", "l");
0240     femc_legend->AddEntry(h_electron_FEMC, "Electrons", "l");
0241     femc_legend->Draw();
0242 
0243     img->FromPad(&cPNG);
0244     img->WriteImage(save_file_name);
0245 
0246     delete img;
0247 }
0248 
0249 char *strdup(const char *s)
0250 {
0251     char *const t = new char[strlen(s) + 1];
0252     return strcpy(t, s);
0253 }
0254 
0255 char *generate_name(const particle_type p, const int particle_energy_gev,
0256             const detector d)
0257 {
0258     std::stringstream name;
0259     name << "Histogram";
0260     switch (p) {
0261     case electron:
0262         name << "_Electron";
0263         break;
0264     case pion:
0265         name << "_Pion";
0266         break;
0267     }
0268 
0269     name << "_" << particle_energy_gev << "GeV";
0270     switch (d) {
0271     case cemc:
0272         name << "_CEMC";
0273         break;
0274     case eemc:
0275         name << "_EEMC";
0276         break;
0277     }
0278 
0279     return strdup(name.str().c_str());
0280 }
0281 
0282 char *generate_file_path(const particle_type p,
0283              const int particle_energy_gev, const detector d)
0284 {
0285     std::stringstream path;
0286     path << data_directory;
0287 
0288     switch (p) {
0289     case electron:
0290         path << "/Electrons/Electrons";
0291         break;
0292     case pion:
0293         path << "/Pions/Pions";
0294         break;
0295     }
0296 
0297     path << particle_energy_gev;
0298     switch (d) {
0299     case cemc:
0300         path << "C";
0301         break;
0302     case eemc:
0303         path << "E";
0304         break;
0305     case femc:
0306         path << "F";
0307         break;
0308     }
0309 
0310     path << ".root";
0311 
0312     return strdup(path.str().c_str());
0313 }
0314 
0315 char *generate_save_name(const int particle_energy_gev)
0316 {
0317     std::stringstream name;
0318     name << "Electron-Pion-" << particle_energy_gev <<
0319         "GeV-CEMC-EEMC-FEMC.png";
0320 
0321     return strdup(name.str().c_str());
0322 }
0323 
0324 char *generate_title(const int particle_energy_gev)
0325 {
0326     std::stringstream title;
0327     title << particle_energy_gev << "GeV";
0328 
0329     return strdup(title.str().c_str());
0330 }
0331 
0332 char *generate_label(const int particle_energy_gev, const detector d)
0333 {
0334     std::stringstream label;
0335 
0336     switch (d) {
0337     case cemc:
0338         label << "CEMC ";
0339         break;
0340     case eemc:
0341         label << "EEMC ";
0342         break;
0343     case femc:
0344         label << "FEMC ";
0345         break;
0346     }
0347     label << particle_energy_gev << "GeV";
0348 
0349     return strdup(label.str().c_str());
0350 }