Jet-Study/macros/display-events.C

0001 // -- c++ includes --
0002 #include <string>
0003 #include <iostream>
0004 #include <iomanip>
0005 #include <fstream>
0006
0007 // -- root includes --
0008 #include <TH2F.h>
0009 #include <TF1.h>
0010 #include <TKey.h>
0011 #include <TFile.h>
0012 #include <TLine.h>
0013 #include <TLegend.h>
0014 #include <TLatex.h>
0015 #include <TCanvas.h>
0016
0017 // -- sPHENIX Style
0018 #include "sPhenixStyle.C"
0019
0020 // -- Utils
0021 #include "jetvalidation/JetUtils.h"
0022
0023 using std::cout;
0024 using std::cerr;
0025 using std::endl;
0026 using std::string;
0027 using std::to_string;
0028 using std::vector;
0029 using std::stringstream;
0030 using std::min;
0031 using std::max;
0032 using std::ofstream;
0033 using std::unique_ptr;
0034 using std::unordered_map;
0035
0036 namespace myAnalysis {
0037     void plots(const string &output);
0038     Int_t readFile(const string &input, const string &output = "");
0039
0040     vector<string> files;
0041     string dirname = "CEMC";
0042     Float_t zMin = 1;
0043     Float_t zMax = 100;
0044     UInt_t bins_phi = 64;
0045     Float_t phi_low = -M_PI;
0046     Float_t phi_high = M_PI;
0047     UInt_t bins_eta = 24;
0048     Float_t eta_low = -1.1;
0049     Float_t eta_high = 1.1;
0050     Bool_t savePlots = true;
0051     UInt_t triggerBits = 42;
0052
0053     enum m_event_type {_60_100, _100_200, _200_500, ABOVE_500};
0054     vector<string> labels = {"60 GeV to 100 GeV", "100 GeV to 200 GeV", "200 GeV to 500 GeV", "Above 500 GeV"};
0055 }
0056
0057 Int_t myAnalysis::readFile(const string &input, const string &output) {
0058     std::ifstream file(input);
0059
0060     // Check if the file was successfully opened
0061     if (!file.is_open()) {
0062         cerr << "Failed to open file: " << input << endl;
0063         return 1;
0064     }
0065
0066     string line;
0067     ofstream outputFile;
0068
0069     if(!output.empty()) outputFile.open(output);
0070
0071     UInt_t ctr = 0;
0072     // loop over each run
0073     while (std::getline(file, line)) {
0074         ++ctr;
0075         TFile* tfile = TFile::Open(line.c_str());
0076         if (!tfile || tfile->IsZombie()) {
0077             std::cerr << "Error opening file: " << line << std::endl;
0078             continue;
0079         }
0080
0081         TDirectory* dir = (TDirectory*)tfile->Get(dirname.c_str());
0082         if (!dir) {
0083             std::cerr << "Directory not found: " << dirname << std::endl;
0084             tfile->Close();
0085             continue;
0086         }
0087
0088         if(dir->GetListOfKeys()->GetSize()) {
0089             files.push_back(line);
0090             cout << "Reading file: " << line << endl;
0091             if(!output.empty()) outputFile << line << endl;
0092         }
0093
0094         tfile->Close();
0095     }
0096
0097     if(!output.empty()) outputFile.close();
0098
0099     cout << "Reading: " << files.size() << " files, " << files.size()*100./ctr << " %" << endl;
0100
0101     // Close the file
0102     file.close();
0103
0104     return 0;
0105 }
0106
0107 void myAnalysis::plots(const string &output) {
0108
0109     TCanvas* c1 = new TCanvas();
0110     c1->SetTickx();
0111     c1->SetTicky();
0112
0113     c1->SetCanvasSize(2900, 1000);
0114     c1->SetLeftMargin(.06);
0115     c1->SetRightMargin(.12);
0116     c1->SetTopMargin(.1);
0117     c1->SetBottomMargin(.12);
0118
0119     gStyle->SetOptTitle(1);
0120     gStyle->SetTitleStyle(0);
0121     gStyle->SetTitleFontSize(0.08);
0122     gStyle->SetTitleW(1);
0123     gStyle->SetTitleH(0.08);
0124     gStyle->SetTitleFillColor(0);
0125     gStyle->SetTitleBorderSize(0);
0126     gStyle->SetTitleXOffset(1);
0127     gStyle->SetTitleYOffset(1);
0128
0129     string tag = output.substr(0,output.rfind("."));
0130
0131     ofstream event_log(tag+"-events.csv");
0132     event_log << "Run,Event,Triggers,zvtx,pt1,pt2,frcem,frcoh,maxJetET,dPhi,isDijet" << endl;
0133
0134     vector<string> suffix_vec;
0135
0136     suffix_vec.push_back("-pt-60-100.pdf");
0137     suffix_vec.push_back("-pt-100-200.pdf");
0138     suffix_vec.push_back("-pt-200-500.pdf");
0139     suffix_vec.push_back("-pt-above-500.pdf");
0140
0141     c1->Print((tag + suffix_vec[m_event_type::_60_100]   + "[").c_str(), "pdf portrait");
0142     c1->Print((tag + suffix_vec[m_event_type::_100_200]  + "[").c_str(), "pdf portrait");
0143     c1->Print((tag + suffix_vec[m_event_type::_200_500]  + "[").c_str(), "pdf portrait");
0144     c1->Print((tag + suffix_vec[m_event_type::ABOVE_500] + "[").c_str(), "pdf portrait");
0145
0146     c1->cd();
0147
0148     UInt_t ctr[5] = {0};
0149
0150     TH1::AddDirectory(kFALSE);
0151
0152     // File Loop
0153     for(UInt_t i = 0; i < files.size(); ++i) {
0154         if(i%20 == 0) cout << "Progress: " << i << ", " << (i+1)*100./files.size() << " %" << endl;
0155
0156         string file = files[i];
0157
0158         TFile* tfile = TFile::Open(file.c_str());
0159
0160         TList* keysCEMCBase = ((TDirectory*)tfile->Get("CEMCBase"))->GetListOfKeys();
0161         TList* keysCEMC     = ((TDirectory*)tfile->Get("CEMC"))->GetListOfKeys();
0162         TList* keysIHCal    = ((TDirectory*)tfile->Get("IHCal"))->GetListOfKeys();
0163         TList* keysOHCal    = ((TDirectory*)tfile->Get("OHCal"))->GetListOfKeys();
0164
0165         // Event Loop
0166         for(UInt_t j = 0; j < keysCEMC->GetSize(); ++j) {
0167             string name  = keysCEMC->At(j)->GetName();
0168             vector<string> tokens = JetUtils::split(name,'_');
0169
0170             string run         = tokens[1];
0171             string event       = tokens[2];
0172             vector<string> triggerIndex = JetUtils::split(tokens[3],'-');
0173             Int_t zvtx         = stoi(tokens[4]);
0174             Int_t leadJetPt    = stoi(tokens[5]);
0175             Int_t subLeadJetPt = stoi(tokens[6]);
0176             Float_t frcem      = stof(tokens[7]);
0177             Float_t frcoh      = stof(tokens[8]);
0178             Int_t maxJetET     = stoi(tokens[9]);
0179             Float_t dPhi       = stof(tokens[10]);
0180             Bool_t isDijet     = stoi(tokens[11]);
0181
0182             cout << "Run: " << run << ", Event: " << event
0183                  << ", trigger: " << tokens[3]
0184                  << ", Z: " << zvtx << " cm"
0185                  << ", pt: " << leadJetPt << " GeV" << endl;
0186
0187             event_log << run << "," << event << "," << tokens[3] << "," << zvtx << "," << leadJetPt << "," << subLeadJetPt
0188                       << "," << frcem << "," << frcoh << "," << maxJetET << "," << dPhi << "," << isDijet << endl;
0189
0190             auto hCEMC = (TH2*)tfile->Get(("CEMC/"+name).c_str());
0191
0192             // hCEMC->GetXaxis()->SetLimits(0,64);
0193             // hCEMC->GetXaxis()->SetNdivisions(21, false);
0194             // hCEMC->GetXaxis()->SetLabelSize(0.04);
0195             // hCEMC->GetXaxis()->SetTickSize(0.01);
0196             // hCEMC->GetYaxis()->SetTickSize(0.01);
0197             // hCEMC->GetYaxis()->SetLabelSize(0.04);
0198             // hCEMC->GetYaxis()->SetLimits(0,24);
0199             // hCEMC->GetYaxis()->SetNdivisions(11, false);
0200             hCEMC->GetYaxis()->SetTitleOffset(0.5);
0201
0202             name  = keysCEMCBase->At(j)->GetName();
0203             auto hCEMCBase = (TH2*)tfile->Get(("CEMCBase/"+name).c_str());
0204             hCEMCBase->GetYaxis()->SetTitleOffset(0.5);
0205
0206             name  = keysIHCal->At(j)->GetName();
0207             auto hIHCal = (TH2*)tfile->Get(("IHCal/"+name).c_str());
0208             if(hIHCal->GetMaximum() >= zMin) {
0209                 hIHCal->SetMinimum(zMin);
0210
0211                 hIHCal->GetYaxis()->SetTitleOffset(0.5);
0212             }
0213
0214             name  = keysOHCal->At(j)->GetName();
0215             auto hOHCal = (TH2*)tfile->Get(("OHCal/"+name).c_str());
0216             if(hOHCal->GetMaximum() >= zMin) {
0217                 hOHCal->SetMinimum(zMin);
0218
0219                 hOHCal->GetYaxis()->SetTitleOffset(0.5);
0220             }
0221
0222             m_event_type pt_key = m_event_type::_60_100;
0223
0224             if(leadJetPt >= 60 && leadJetPt < 100) {
0225                 ++ctr[1];
0226             }
0227             else if(leadJetPt >= 100 && leadJetPt < 200) {
0228                 pt_key = m_event_type::_100_200;
0229                 ++ctr[2];
0230             }
0231             else if(leadJetPt >= 200 && leadJetPt < 500) {
0232                 pt_key = m_event_type::_200_500;
0233                 ++ctr[3];
0234             }
0235             else if(leadJetPt >= 500) {
0236                 pt_key = m_event_type::ABOVE_500;
0237                 ++ctr[4];
0238             }
0239
0240             string suffix = suffix_vec[pt_key];
0241
0242             if(savePlots) {
0243                 hCEMCBase->Draw("COLZ1");
0244                 c1->Print((tag+suffix).c_str(), "pdf portrait");
0245
0246                 hCEMCBase->SetMinimum(zMin);
0247                 c1->Print((tag+suffix).c_str(), "pdf portrait");
0248
0249                 hCEMC->Draw("COLZ1");
0250                 c1->Print((tag+suffix).c_str(), "pdf portrait");
0251
0252                 hCEMC->SetMinimum(zMin);
0253                 c1->Print((tag+suffix).c_str(), "pdf portrait");
0254
0255                 if(hIHCal->GetMaximum() >= zMin) {
0256                     hIHCal->Draw("COLZ1");
0257                     c1->Print((tag+suffix).c_str(), "pdf portrait");
0258                 }
0259
0260                 if(hOHCal->GetMaximum() >= zMin) {
0261                     hOHCal->Draw("COLZ1");
0262                     c1->Print((tag+suffix).c_str(), "pdf portrait");
0263                 }
0264             }
0265             ++ctr[0];
0266         }
0267         tfile->Close();
0268     }
0269
0270     event_log.close();
0271
0272     // gPad->SetGrid();
0273     gPad->SetLogy();
0274
0275     c1->Print((tag + suffix_vec[m_event_type::_60_100]   + "]").c_str(), "pdf portrait");
0276     c1->Print((tag + suffix_vec[m_event_type::_100_200]  + "]").c_str(), "pdf portrait");
0277     c1->Print((tag + suffix_vec[m_event_type::_200_500]  + "]").c_str(), "pdf portrait");
0278     c1->Print((tag + suffix_vec[m_event_type::ABOVE_500] + "]").c_str(), "pdf portrait");
0279
0280     // cleaning step, remove extra files
0281     for(UInt_t j = 0; j < labels.size(); ++j) {
0282         string suffix = suffix_vec[j];
0283         if(ctr[j+1] == 0 || !savePlots) {
0284             std::remove((tag+suffix).c_str());
0285         }
0286     }
0287
0288     cout << "----------------------------" << endl;
0289     cout << "Total Events: " << ctr[0] << endl;
0290     for(UInt_t i = 0; i < labels.size(); ++i) {
0291         cout << labels[i] << ": " << ctr[i+1] << " Events" << endl;
0292     }
0293 }
0294
0295 void display_events(const string &input,
0296                     const string &output="plots.pdf",
0297                     Bool_t savePlots=true,
0298                     const string &outputFileList="") {
0299     cout << "#############################" << endl;
0300     cout << "Run Parameters" << endl;
0301     cout << "input: "  << input << endl;
0302     cout << "output: " << output << endl;
0303     cout << "Save Plots: " << savePlots << endl;
0304     cout << "outputFileList: " << outputFileList << endl;
0305     cout << "#############################" << endl;
0306
0307     // set sPHENIX plotting style
0308     SetsPhenixStyle();
0309
0310     myAnalysis::savePlots = savePlots;
0311
0312     myAnalysis::readFile(input, outputFileList);
0313
0314     myAnalysis::plots(output);
0315 }
0316
0317 # ifndef __CINT__
0318 Int_t main(Int_t argc, char* argv[]) {
0319 if(argc < 2 || argc > 5){
0320         cout << "usage: ./display-events input [output] [savePlots] [outputFileList]" << endl;
0321         cout << "input: input list root file" << endl;
0322         cout << "output: output pdf file" << endl;
0323         cout << "savePlots: save event plots" << endl;
0324         cout << "outputFileList: List of files that have data." << endl;
0325         return 1;
0326     }
0327
0328     string output  = "plots.pdf";
0329     Bool_t savePlots = true;
0330     string outputFileList  = "";
0331
0332     if(argc >= 3) {
0333         output = argv[2];
0334     }
0335     if(argc >= 4) {
0336         savePlots = atoi(argv[3]);
0337     }
0338     if(argc >= 5) {
0339         outputFileList = argv[4];
0340     }
0341
0342     display_events(argv[1], output, savePlots, outputFileList);
0343
0344     cout << "======================================" << endl;
0345     cout << "done" << endl;
0346     return 0;
0347 }
0348 # endif