calorimeter/calo_cdb/genStatus.cc

0001 #include "genStatus.h"
0002
0003 #include "geometry_constants.h"
0004
0005 // -- sPHENIX includes --
0006 #include <calobase/TowerInfoDefs.h>
0007
0008 #include <cdbobjects/CDBTTree.h>
0009
0010 #include <emcnoisytowerfinder/emcNoisyTowerFinder.h>
0011
0012 // -- root includes --
0013 #include <TFile.h>
0014
0015 // c++ includes --
0016 #include <filesystem>
0017 #include <fstream>
0018 #include <iomanip>
0019 #include <iostream>
0020 #include <map>
0021
0022 void GenStatus::setRunDataset(const std::string &input)
0023 {
0024   std::string basename = std::filesystem::path(input).filename().stem().string();
0025   m_run = basename.substr(0, basename.find('_'));
0026   m_dataset = basename.substr(basename.find('_') + 1, basename.size() - basename.find('_'));
0027 }
0028
0029 int GenStatus::readHists(const std::string &input)
0030 {
0031   // Create an input stream
0032   std::ifstream file(input);
0033
0034   // Check if the file was successfully opened
0035   if (!file.is_open())
0036   {
0037     std::cout << "Failed to open file list: " << input << std::endl;
0038     return 1;
0039   }
0040
0041   std::cout << "Reading Hists" << std::endl;
0042   std::cout << "======================================" << std::endl;
0043
0044   h_CaloValid_cemc_etaphi_badChi2 = std::make_unique<TProfile2D>("cemc_etaphi_badChi2", "", cemc_bins_eta, 0, cemc_bins_eta, cemc_bins_phi, 0, cemc_bins_phi);
0045   h_CaloValid_ihcal_etaphi_badChi2 = std::make_unique<TProfile2D>("ihcal_etaphi_badChi2", "", hcal_bins_eta, 0, hcal_bins_eta, hcal_bins_phi, 0, hcal_bins_phi);
0046   h_CaloValid_ohcal_etaphi_badChi2 = std::make_unique<TProfile2D>("ohcal_etaphi_badChi2", "", hcal_bins_eta, 0, hcal_bins_eta, hcal_bins_phi, 0, hcal_bins_phi);
0047
0048   h_CaloValid_cemc_etaphi_time_raw = std::make_unique<TProfile2D>("cemc_etaphi_time_raw", "", cemc_bins_eta, 0, cemc_bins_eta, cemc_bins_phi, 0, cemc_bins_phi);
0049   h_CaloValid_ihcal_etaphi_time_raw = std::make_unique<TProfile2D>("ihcal_etaphi_time_raw", "", hcal_bins_eta, 0, hcal_bins_eta, hcal_bins_phi, 0, hcal_bins_phi);
0050   h_CaloValid_ohcal_etaphi_time_raw = std::make_unique<TProfile2D>("ohcal_etaphi_time_raw", "", hcal_bins_eta, 0, hcal_bins_eta, hcal_bins_phi, 0, hcal_bins_phi);
0051
0052   std::string line;
0053   std::map<std::string, int> ctr;
0054   while (std::getline(file, line))
0055   {
0056     ++ctr["total_files"];
0057
0058     std::cout << "Reading File: " << line << std::endl;
0059     std::unique_ptr<TFile> tf(TFile::Open(line.c_str()));
0060     if (!tf || tf->IsZombie())
0061     {
0062       std::cout << "Error: Could not open ROOT file: " << line << std::endl;
0063       continue;  // Indicate an error
0064     }
0065
0066     ++ctr["successfully_opened_files"];
0067
0068     auto *h = dynamic_cast<TProfile2D *>(tf->Get("h_CaloValid_cemc_etaphi_badChi2"));
0069
0070     if (h)
0071     {
0072       h_CaloValid_cemc_etaphi_badChi2->Add(h);
0073       ++ctr["h_CaloValid_cemc_etaphi_badChi2"];
0074     }
0075
0076     h = dynamic_cast<TProfile2D *>(tf->Get("h_CaloValid_ihcal_etaphi_badChi2"));
0077
0078     if (h)
0079     {
0080       h_CaloValid_ihcal_etaphi_badChi2->Add(h);
0081       ++ctr["h_CaloValid_ihcal_etaphi_badChi2"];
0082     }
0083
0084     h = dynamic_cast<TProfile2D *>(tf->Get("h_CaloValid_ohcal_etaphi_badChi2"));
0085
0086     if (h)
0087     {
0088       h_CaloValid_ohcal_etaphi_badChi2->Add(h);
0089       ++ctr["h_CaloValid_ohcal_etaphi_badChi2"];
0090     }
0091
0092     h = dynamic_cast<TProfile2D *>(tf->Get("h_CaloValid_cemc_etaphi_time_raw"));
0093
0094     if (h)
0095     {
0096       h_CaloValid_cemc_etaphi_time_raw->Add(h);
0097       ++ctr["h_CaloValid_cemc_etaphi_time_raw"];
0098     }
0099
0100     h = dynamic_cast<TProfile2D *>(tf->Get("h_CaloValid_ihcal_etaphi_time_raw"));
0101
0102     if (h)
0103     {
0104       h_CaloValid_ihcal_etaphi_time_raw->Add(h);
0105       ++ctr["h_CaloValid_ihcal_etaphi_time_raw"];
0106     }
0107
0108     h = dynamic_cast<TProfile2D *>(tf->Get("h_CaloValid_ohcal_etaphi_time_raw"));
0109
0110     if (h)
0111     {
0112       h_CaloValid_ohcal_etaphi_time_raw->Add(h);
0113       ++ctr["h_CaloValid_ohcal_etaphi_time_raw"];
0114     }
0115
0116     tf->Close();
0117   }
0118
0119   std::cout << "===============================" << std::endl;
0120   std::cout << "Stats" << std::endl;
0121   std::cout << "Successfully opened files: " << ctr["successfully_opened_files"] << ", " << ctr["successfully_opened_files"] * 100. / ctr["total_files"] << " %" << std::endl;
0122   for (const auto &[name, value] : ctr)
0123   {
0124     if (name.starts_with("h_CaloValid"))
0125     {
0126       std::cout << "Hist: " << name << ", Found: " << value << ", " << value * 100. / ctr["successfully_opened_files"] << " %" << std::endl;
0127     }
0128   }
0129   std::cout << "===============================" << std::endl;
0130
0131   // Close the file
0132   file.close();
0133
0134   return 0;
0135 }
0136
0137 void GenStatus::histToCaloCDBTree(const std::string &outputfile, const std::string &fieldName, int icalo, TProfile2D *hist)
0138 {
0139   unsigned int neta;
0140   unsigned int nphi;
0141
0142   if (icalo != 0 && icalo != 1)
0143   {
0144     return;
0145   }
0146
0147   if (icalo == 0)
0148   {
0149     neta = CaloGeometry::CEMC_ETA_BINS;
0150     nphi = CaloGeometry::CEMC_PHI_BINS;
0151   }
0152   if (icalo == 1)
0153   {
0154     neta = CaloGeometry::HCAL_ETA_BINS;
0155     nphi = CaloGeometry::HCAL_PHI_BINS;
0156   }
0157
0158   std::unique_ptr<CDBTTree> cdbttree = std::make_unique<CDBTTree>(outputfile);
0159
0160   double mean = 0;
0161   int count = 0;
0162
0163   for (unsigned int ie = 0; ie < neta; ie++)
0164   {
0165     for (unsigned int ip = 0; ip < nphi; ip++)
0166     {
0167       unsigned int key;
0168       if (icalo == 0)
0169       {
0170         key = TowerInfoDefs::encode_emcal(ie, ip);
0171       }
0172       if (icalo == 1)
0173       {
0174         key = TowerInfoDefs::encode_hcal(ie, ip);
0175       }
0176       float val = hist->GetBinContent(ie + 1, ip + 1);
0177       cdbttree->SetFloatValue(key, fieldName, val);
0178       mean += val;
0179       count++;
0180     }
0181   }
0182
0183   std::cout << "Writing " << outputfile << "   with mean=" << mean / count << std::endl;
0184   cdbttree->Commit();
0185   cdbttree->WriteCDBTTree();
0186 }
0187
0188 void GenStatus::analyze(const std::string &outputDir)
0189 {
0190   std::string detector = "CEMC";
0191   // fracBadChi2
0192   std::string payloadName = outputDir + "/" + detector + "_hotTowers_fracBadChi2" + "_" + m_dataset + "_" + m_run + ".root";
0193   if (h_CaloValid_cemc_etaphi_badChi2)
0194   {
0195     histToCaloCDBTree(payloadName, "fraction", 0, h_CaloValid_cemc_etaphi_badChi2.get());
0196   }
0197   // time
0198   payloadName = outputDir + "/" + detector + "_meanTime" + "_" + m_dataset + "_" + m_run + ".root";
0199   if (h_CaloValid_cemc_etaphi_time_raw)
0200   {
0201     histToCaloCDBTree(payloadName, "time", 0, h_CaloValid_cemc_etaphi_time_raw.get());
0202   }
0203
0204   detector = "HCALIN";
0205   // fracBadChi2
0206   payloadName = outputDir + "/" + detector + "_hotTowers_fracBadChi2" + "_" + m_dataset + "_" + m_run + ".root";
0207   if (h_CaloValid_ihcal_etaphi_badChi2)
0208   {
0209     histToCaloCDBTree(payloadName, "fraction", 1, h_CaloValid_ihcal_etaphi_badChi2.get());
0210   }
0211   // time
0212   payloadName = outputDir + "/" + detector + "_meanTime" + "_" + m_dataset + "_" + m_run + ".root";
0213   if (h_CaloValid_ihcal_etaphi_time_raw)
0214   {
0215     histToCaloCDBTree(payloadName, "time", 1, h_CaloValid_ihcal_etaphi_time_raw.get());
0216   }
0217
0218   detector = "HCALOUT";
0219   // fracBadChi2
0220   payloadName = outputDir + "/" + detector + "_hotTowers_fracBadChi2" + "_" + m_dataset + "_" + m_run + ".root";
0221   if (h_CaloValid_ohcal_etaphi_badChi2)
0222   {
0223     histToCaloCDBTree(payloadName, "fraction", 1, h_CaloValid_ohcal_etaphi_badChi2.get());
0224   }
0225   // time
0226   payloadName = outputDir + "/" + detector + "_meanTime" + "_" + m_dataset + "_" + m_run + ".root";
0227   if (h_CaloValid_ohcal_etaphi_time_raw)
0228   {
0229     histToCaloCDBTree(payloadName, "time", 1, h_CaloValid_ohcal_etaphi_time_raw.get());
0230   }
0231 }
0232
0233 void GenStatus::process(const std::string &input, const std::string &output)
0234 {
0235   std::cout << "#############################" << std::endl;
0236   std::cout << "Run Parameters" << std::endl;
0237   std::cout << "input: " << input << std::endl;
0238   std::cout << "output: " << output << std::endl;
0239   std::cout << "#############################" << std::endl;
0240
0241   setRunDataset(input);
0242
0243   std::cout << "Processing: Run: " << m_run << ", Dataset: " << m_dataset << std::endl;
0244
0245   std::stringstream datasetDir;
0246   datasetDir.str("");
0247   datasetDir << output << "/" << m_run << "_" << m_dataset;
0248
0249   std::string outputDir = datasetDir.str();
0250
0251   std::string hotMapFile = "EMCalHotMap_" + m_dataset + "_" + m_run + ".root";
0252   std::string hotMapOutput = datasetDir.str() + "/" + hotMapFile;
0253
0254   datasetDir << "/QA";
0255
0256   // create output & QA directory
0257   std::filesystem::create_directories(datasetDir.str());
0258
0259   std::string s_input = input;
0260   std::string hotMapOutputQA = datasetDir.str() + "/" + hotMapFile;
0261
0262   // merges individal qa into one per run
0263   readHists(input);
0264   analyze(outputDir);
0265
0266   std::unique_ptr<emcNoisyTowerFinder> calo = std::make_unique<emcNoisyTowerFinder>();
0267   calo->FindHot(s_input, hotMapOutput, "h_CaloValid_cemc_etaphi");
0268
0269   if (std::filesystem::exists(hotMapOutput))
0270   {
0271     std::filesystem::rename(hotMapOutput, hotMapOutputQA);
0272   }
0273   else
0274   {
0275     std::cout << "ERROR: EMCal Hot Map FAILED to Create." << std::endl;
0276   }
0277 }