EMCal-Hot-Towers/macro/display-hot.C

0001 // -- c++ includes --
0002 #include <string>
0003 #include <iostream>
0004 #include <iomanip>
0005 #include <fstream>
0006 #include <filesystem>
0007
0008 // -- root includes --
0009 #include <TH2F.h>
0010 #include <TF1.h>
0011 #include <TFile.h>
0012 #include <TLine.h>
0013 #include <TLegend.h>
0014 #include <TLatex.h>
0015 #include <TCanvas.h>
0016
0017 // -- Tower Info
0018 #include <calobase/TowerInfoDefs.h>
0019
0020 // -- sPHENIX Style
0021 #include "sPhenixStyle.C"
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::pair;
0034
0035 R__LOAD_LIBRARY(libcalo_io.so)
0036
0037 namespace myAnalysis {
0038     void plots(const string& i_input, const string &hotMap, const string &outputDir);
0039     Int_t readHotTowerIndexFile(const string &hotList);
0040
0041     vector<pair<UInt_t,UInt_t>> hotTowerIndex;
0042 }
0043
0044 Int_t myAnalysis::readHotTowerIndexFile(const string &hotList) {
0045
0046   cout << "Reading Hot and Ref tower indices" << endl;
0047
0048   std::ifstream file(hotList);
0049
0050   // Check if the file was successfully opened
0051   if (!file.is_open()) {
0052       cerr << "Failed to open file: " << hotList << endl;
0053       return 1;
0054   }
0055
0056   string line;
0057
0058   // skip header
0059   std::getline(file, line);
0060
0061   // loop over each run
0062   while (std::getline(file, line)) {
0063       std::istringstream lineStream(line);
0064
0065       UInt_t towerHotIndex;
0066       UInt_t towerRefIndex;
0067       Char_t comma;
0068
0069       if (lineStream >> towerHotIndex >> comma >> towerRefIndex) {
0070           hotTowerIndex.push_back(std::make_pair(towerHotIndex, towerRefIndex));
0071           cout << "Hot: " << towerHotIndex << ", Ref: " << towerRefIndex << endl;
0072       }
0073       else {
0074           cerr << "Failed to parse line: " << line << endl;
0075           return 1;
0076       }
0077   }
0078
0079   // Close the file
0080   file.close();
0081
0082   return 0;
0083 }
0084
0085 void myAnalysis::plots(const string& i_input, const string &hotMap, const string &outputDir) {
0086     TFile input(i_input.c_str());
0087     TFile inputHot(hotMap.c_str());
0088
0089     TCanvas* c1 = new TCanvas();
0090     c1->SetTickx();
0091     c1->SetTicky();
0092
0093     c1->SetCanvasSize(1500, 1000);
0094     c1->SetLeftMargin(.13);
0095     c1->SetRightMargin(.05);
0096
0097     string outputEnergyDir = outputDir + "/energy/";
0098     string output = outputEnergyDir + "energy.pdf";
0099
0100     c1->Print((output + "[").c_str(), "pdf portrait");
0101
0102     c1->cd();
0103
0104     TLatex l1;
0105     TLatex l2;
0106     // l1.SetTextSize(0.05);
0107
0108     auto hHotMap = (TH1F*)inputHot.Get("hBadTowersHot");
0109     auto hRunStatus = (TH1F*)inputHot.Get("hRunStatus");
0110
0111     UInt_t nRuns = hRunStatus->GetBinContent(1);
0112     cout << "Total Runs: " << nRuns << endl;
0113
0114     stringstream s;
0115
0116     UInt_t i = 0;
0117     // initialize histograms
0118     for(auto idx : hotTowerIndex) {
0119
0120         UInt_t key    = TowerInfoDefs::encode_emcal(idx.first);
0121         UInt_t etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
0122         UInt_t phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
0123
0124         string name = "Hot/HotTower_"+to_string(phibin)+"_"+to_string(etabin);
0125         string lHot = "Hot (Status = 2): (" + to_string(phibin) + "," + to_string(etabin) + ")";
0126         auto hHot   = (TH1F*)input.Get(name.c_str());
0127         Int_t overFlowHot = hHot->GetBinContent(hHot->GetNbinsX()+1);
0128         if(hHot->Integral()) hHot->Scale(1./hHot->Integral());
0129
0130         name            = "HotComplement/HotTowerComplement_"+to_string(phibin)+"_"+to_string(etabin);
0131         string lHotComp = "Hot (Status #neq 2): (" + to_string(phibin) + "," + to_string(etabin) + ")";
0132         auto hHotComp   = (TH1F*)input.Get(name.c_str());
0133         Int_t overFlowHotComp = hHotComp->GetBinContent(hHotComp->GetNbinsX()+1);
0134         if(hHotComp->Integral()) hHotComp->Scale(1./hHotComp->Integral());
0135
0136         UInt_t nRunsHot = hHotMap->GetBinContent(idx.first+1);
0137
0138         key    = TowerInfoDefs::encode_emcal(idx.second);
0139         etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
0140         phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
0141
0142         name        = "Ref/RefTower_"+to_string(phibin)+"_"+to_string(etabin)+"_"+to_string(i);
0143         string lRef = "Ref (Status = 0): (" + to_string(phibin) + "," + to_string(etabin) + ")";
0144         auto hRef   = (TH1F*)input.Get(name.c_str());
0145         Int_t overFlowRef = hRef->GetBinContent(hRef->GetNbinsX()+1);
0146         if(hRef->Integral()) hRef->Scale(1./hRef->Integral());
0147
0148         gPad->SetLogy();
0149         hHot->SetLineColor(kRed);
0150         hHot->GetYaxis()->SetTitle("Normalized Counts");
0151         hHot->GetYaxis()->SetTitleOffset(1.2);
0152         hHot->Draw("HIST");
0153         hHotComp->SetLineColor(kGreen+2);
0154         hHotComp->Draw("HIST same");
0155         hRef->SetLineColor(kBlue);
0156         hRef->Draw("HIST same");
0157
0158         name = "#splitline{#color[2]{Overflow: " + to_string(overFlowHot) + "}}{#splitline{" +
0159           "#color[3]{Overflow: " + to_string(overFlowHotComp) + "}}{" +
0160           "#color[4]{Overflow: " + to_string(overFlowRef) + "}}}";
0161
0162         l1.DrawLatexNDC(0.57,0.7, name.c_str());
0163
0164         s.str("");
0165         s << "#color[2]{Runs: " << nRunsHot << ", (" << (Int_t)(nRunsHot*10000./nRuns)/100. << " %)}";
0166
0167         l2.DrawLatexNDC(0.2,0.85, s.str().c_str());
0168
0169         TLegend *leg = new TLegend(0.52,.75,0.72,.92);
0170         leg->SetFillStyle(0);
0171         leg->AddEntry(hHot,lHot.c_str(),"f");
0172         leg->AddEntry(hHotComp,lHotComp.c_str(),"f");
0173         leg->AddEntry(hRef,lRef.c_str(),"f");
0174         leg->Draw("same");
0175
0176         c1->Print((outputEnergyDir + string(hHot->GetName()) + ".png").c_str());
0177         c1->Print(output.c_str(), "pdf portrait");
0178
0179         gPad->SetLogy(0);
0180         hHot->Rebin(10);
0181         hHotComp->Rebin(10);
0182         hRef->Rebin(10);
0183         auto h = (TH1F*)hHot->Clone("h");
0184         hHot->Divide(hRef);
0185         hHot->GetYaxis()->SetTitle("Ratio");
0186         hHot->SetLineColor(kBlue);
0187         hHot->Draw("HIST");
0188
0189         h->Divide(hHotComp);
0190         h->SetLineColor(kGreen+2);
0191         h->Draw("HIST same");
0192
0193         Float_t high = max(hHot->GetMaximum(), h->GetMaximum());
0194         hHot->GetYaxis()->SetRangeUser(0, high);
0195
0196         TLine* line = new TLine(0, 1, 1.6e4, 1);
0197         line->SetLineColor(kBlack);
0198         line->SetLineStyle(9);
0199         // line->SetLineWidth(1);
0200         line->Draw("same");
0201
0202         leg = new TLegend(0.6,.75,0.8,.92);
0203         leg->SetFillStyle(0);
0204         leg->AddEntry(hHot,"Hot (Status = 2) / Ref (Status = 0)","l");
0205         leg->AddEntry(h,"Hot (Status = 2) / Hot (Status #neq 2)","l");
0206         leg->SetTextSize(0.03);
0207         leg->Draw("same");
0208
0209         c1->Print((outputEnergyDir + string(hHot->GetName()) + "-ratio.png").c_str());
0210         c1->Print(output.c_str(), "pdf portrait");
0211
0212         ++i;
0213     }
0214
0215     c1->Print((output + "]").c_str(), "pdf portrait");
0216
0217     inputHot.Close();
0218     input.Close();
0219 }
0220
0221 void display_hot(const string &input, const string &hotList, const string &hotMap, const string &outputDir=".") {
0222     cout << "#############################" << endl;
0223     cout << "Run Parameters" << endl;
0224     cout << "input: "  << input << endl;
0225     cout << "hotList: "  << hotList << endl;
0226     cout << "hotMap: "  << hotMap << endl;
0227     cout << "outputDir: " << outputDir << endl;
0228     cout << "#############################" << endl;
0229
0230     // set sPHENIX plotting style
0231     SetsPhenixStyle();
0232
0233     Int_t ret = myAnalysis::readHotTowerIndexFile(hotList);
0234     if(ret) return;
0235
0236     std::filesystem::create_directories(outputDir + "/energy");
0237
0238     myAnalysis::plots(input, hotMap, outputDir);
0239 }
0240
0241 # ifndef __CINT__
0242 Int_t main(Int_t argc, char* argv[]) {
0243 if(argc < 4 || argc > 5){
0244         cout << "usage: ./plots-hot input hotList hotMap [outputDir]" << endl;
0245         cout << "input: input root file" << endl;
0246         cout << "hotList: list of hot and reference towers" << endl;
0247         cout << "outputDir: output directory" << endl;
0248         return 1;
0249     }
0250
0251     string outputDir  = ".";
0252
0253     if(argc >= 5) {
0254         outputDir = argv[4];
0255     }
0256
0257     display_hot(argv[1], argv[2], argv[3], outputDir);
0258
0259     cout << "======================================" << endl;
0260     cout << "done" << endl;
0261     return 0;
0262 }
0263 # endif