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 <iostream>
0019 #include <map>
0020
0021 void GenStatus::setRunDataset(const std::string &input)
0022 {
0023   std::string basename = std::filesystem::path(input).filename().stem().string();
0024   m_run = basename.substr(0, basename.find('_'));
0025   m_dataset = basename.substr(basename.find('_') + 1, basename.size() - basename.find('_'));
0026 }
0027
0028 int GenStatus::readHists(const std::string &input)
0029 {
0030   // Create an input stream
0031   std::ifstream file(input);
0032
0033   // Check if the file was successfully opened
0034   if (!file.is_open())
0035   {
0036     std::cout << "Failed to open file list: " << input << std::endl;
0037     return 1;
0038   }
0039
0040   std::cout << "Reading Hists" << std::endl;
0041   std::cout << "======================================" << std::endl;
0042
0043   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);
0044   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);
0045   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);
0046
0047   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);
0048   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);
0049   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);
0050
0051   h_CaloFittingQA_cemc_etaphi_ZScrosscalib = std::make_unique<TProfile2D>("cemc_etaphi_ZScrosscalib", "", cemc_bins_eta, 0, cemc_bins_eta, cemc_bins_phi, 0, cemc_bins_phi);
0052   h_CaloFittingQA_ihcal_etaphi_ZScrosscalib = std::make_unique<TProfile2D>("ihcal_etaphi_ZScrosscalib", "", hcal_bins_eta, 0, hcal_bins_eta, hcal_bins_phi, 0, hcal_bins_phi);
0053   h_CaloFittingQA_ohcal_etaphi_ZScrosscalib = std::make_unique<TProfile2D>("ohcal_etaphi_ZScrosscalib", "", hcal_bins_eta, 0, hcal_bins_eta, hcal_bins_phi, 0, hcal_bins_phi);
0054
0055   std::ofstream outfile(m_CaloValid_list);
0056
0057   if (!outfile.is_open())
0058   {
0059     std::cout << "Error: Unable to open file " << m_CaloValid_list << " for writing." << std::endl;
0060     return 1;
0061   }
0062
0063   std::string line;
0064   std::map<std::string, int> ctr;
0065   while (std::getline(file, line))
0066   {
0067     ++ctr["total_files"];
0068
0069     std::cout << "Reading File: " << line << std::endl;
0070     std::unique_ptr<TFile> tf(TFile::Open(line.c_str()));
0071     if (!tf || tf->IsZombie())
0072     {
0073       std::cout << "Error: Could not open ROOT file: " << line << std::endl;
0074       continue;  // Indicate an error
0075     }
0076
0077     ++ctr["successfully_opened_files"];
0078
0079     if (line.find("HIST_CALOQA") != std::string::npos)
0080     {
0081       outfile << line << std::endl;
0082     }
0083
0084     auto *h = dynamic_cast<TProfile2D *>(tf->Get("h_CaloValid_cemc_etaphi_badChi2"));
0085
0086     if (h)
0087     {
0088       h_CaloValid_cemc_etaphi_badChi2->Add(h);
0089       ++ctr["h_CaloValid_cemc_etaphi_badChi2"];
0090     }
0091
0092     h = dynamic_cast<TProfile2D *>(tf->Get("h_CaloValid_ihcal_etaphi_badChi2"));
0093
0094     if (h)
0095     {
0096       h_CaloValid_ihcal_etaphi_badChi2->Add(h);
0097       ++ctr["h_CaloValid_ihcal_etaphi_badChi2"];
0098     }
0099
0100     h = dynamic_cast<TProfile2D *>(tf->Get("h_CaloValid_ohcal_etaphi_badChi2"));
0101
0102     if (h)
0103     {
0104       h_CaloValid_ohcal_etaphi_badChi2->Add(h);
0105       ++ctr["h_CaloValid_ohcal_etaphi_badChi2"];
0106     }
0107
0108     h = dynamic_cast<TProfile2D *>(tf->Get("h_CaloValid_cemc_etaphi_time_raw"));
0109
0110     if (h)
0111     {
0112       h_CaloValid_cemc_etaphi_time_raw->Add(h);
0113       ++ctr["h_CaloValid_cemc_etaphi_time_raw"];
0114     }
0115
0116     h = dynamic_cast<TProfile2D *>(tf->Get("h_CaloValid_ihcal_etaphi_time_raw"));
0117
0118     if (h)
0119     {
0120       h_CaloValid_ihcal_etaphi_time_raw->Add(h);
0121       ++ctr["h_CaloValid_ihcal_etaphi_time_raw"];
0122     }
0123
0124     h = dynamic_cast<TProfile2D *>(tf->Get("h_CaloValid_ohcal_etaphi_time_raw"));
0125
0126     if (h)
0127     {
0128       h_CaloValid_ohcal_etaphi_time_raw->Add(h);
0129       ++ctr["h_CaloValid_ohcal_etaphi_time_raw"];
0130     }
0131
0132     h = dynamic_cast<TProfile2D *>(tf->Get("h_CaloFittingQA_cemc_etaphi_ZScrosscalib"));
0133
0134     if (h)
0135     {
0136       h_CaloFittingQA_cemc_etaphi_ZScrosscalib->Add(h);
0137       ++ctr["h_CaloFittingQA_cemc_etaphi_ZScrosscalib"];
0138     }
0139
0140     h = dynamic_cast<TProfile2D *>(tf->Get("h_CaloFittingQA_ihcal_etaphi_ZScrosscalib"));
0141
0142     if (h)
0143     {
0144       h_CaloFittingQA_ihcal_etaphi_ZScrosscalib->Add(h);
0145       ++ctr["h_CaloFittingQA_ihcal_etaphi_ZScrosscalib"];
0146     }
0147
0148     h = dynamic_cast<TProfile2D *>(tf->Get("h_CaloFittingQA_ohcal_etaphi_ZScrosscalib"));
0149
0150     if (h)
0151     {
0152       h_CaloFittingQA_ohcal_etaphi_ZScrosscalib->Add(h);
0153       ++ctr["h_CaloFittingQA_ohcal_etaphi_ZScrosscalib"];
0154     }
0155
0156     tf->Close();
0157   }
0158
0159   std::cout << "===============================" << std::endl;
0160   std::cout << "Stats" << std::endl;
0161   std::cout << "Successfully opened files: " << ctr["successfully_opened_files"] << ", " << ctr["successfully_opened_files"] * 100. / ctr["total_files"] << " %" << std::endl;
0162   for (const auto &[name, value] : ctr)
0163   {
0164     if (name.starts_with("h_Calo"))
0165     {
0166       std::cout << "Hist: " << name << ", Found: " << value << ", " << value * 100. / ctr["successfully_opened_files"] << " %" << std::endl;
0167     }
0168   }
0169   std::cout << "===============================" << std::endl;
0170
0171   // Close the file
0172   file.close();
0173   outfile.close();
0174
0175   // Check if the close/flush was successful
0176   if (outfile.fail())
0177   {
0178     std::cout << "Error: Failed to flush/close the file properly: " << m_CaloValid_list << std::endl;
0179     return 1;
0180   }
0181
0182   return 0;
0183 }
0184
0185 void GenStatus::histToCaloCDBTree(const std::string &outputfile, const std::string &fieldName, int icalo, TProfile2D *hist)
0186 {
0187   unsigned int neta;
0188   unsigned int nphi;
0189
0190   if (icalo != 0 && icalo != 1)
0191   {
0192     return;
0193   }
0194
0195   if (icalo == 0)
0196   {
0197     neta = CaloGeometry::CEMC_ETA_BINS;
0198     nphi = CaloGeometry::CEMC_PHI_BINS;
0199   }
0200   if (icalo == 1)
0201   {
0202     neta = CaloGeometry::HCAL_ETA_BINS;
0203     nphi = CaloGeometry::HCAL_PHI_BINS;
0204   }
0205
0206   std::unique_ptr<CDBTTree> cdbttree = std::make_unique<CDBTTree>(outputfile);
0207
0208   double mean = 0;
0209   int count = 0;
0210
0211   for (unsigned int ie = 0; ie < neta; ie++)
0212   {
0213     for (unsigned int ip = 0; ip < nphi; ip++)
0214     {
0215       unsigned int key;
0216       if (icalo == 0)
0217       {
0218         key = TowerInfoDefs::encode_emcal(ie, ip);
0219       }
0220       if (icalo == 1)
0221       {
0222         key = TowerInfoDefs::encode_hcal(ie, ip);
0223       }
0224       float val = hist->GetBinContent(ie + 1, ip + 1);
0225       cdbttree->SetFloatValue(key, fieldName, val);
0226       mean += val;
0227       count++;
0228     }
0229   }
0230
0231   std::cout << "Writing " << outputfile << "   with mean=" << mean / count << std::endl;
0232   cdbttree->Commit();
0233   cdbttree->WriteCDBTTree();
0234 }
0235
0236 void GenStatus::analyze(const std::string &outputDir)
0237 {
0238   std::string detector = "CEMC";
0239   // fracBadChi2
0240   std::string payloadName = outputDir + "/" + detector + "_hotTowers_fracBadChi2" + "_" + m_dataset + "_" + m_run + ".root";
0241   if (h_CaloValid_cemc_etaphi_badChi2->GetEntries())
0242   {
0243     histToCaloCDBTree(payloadName, "fraction", 0, h_CaloValid_cemc_etaphi_badChi2.get());
0244   }
0245   // time
0246   payloadName = outputDir + "/" + detector + "_meanTime" + "_" + m_dataset + "_" + m_run + ".root";
0247   if (h_CaloValid_cemc_etaphi_time_raw->GetEntries())
0248   {
0249     histToCaloCDBTree(payloadName, "time", 0, h_CaloValid_cemc_etaphi_time_raw.get());
0250   }
0251   // ZSCrossCalib
0252   payloadName = outputDir + "/" + detector + "_ZSCrossCalib" + "_" + m_dataset + "_" + m_run + ".root";
0253   if (h_CaloFittingQA_cemc_etaphi_ZScrosscalib->GetEntries())
0254   {
0255     histToCaloCDBTree(payloadName, "ratio", 0, h_CaloFittingQA_cemc_etaphi_ZScrosscalib.get());
0256   }
0257
0258   detector = "HCALIN";
0259   // fracBadChi2
0260   payloadName = outputDir + "/" + detector + "_hotTowers_fracBadChi2" + "_" + m_dataset + "_" + m_run + ".root";
0261   if (h_CaloValid_ihcal_etaphi_badChi2->GetEntries())
0262   {
0263     histToCaloCDBTree(payloadName, "fraction", 1, h_CaloValid_ihcal_etaphi_badChi2.get());
0264   }
0265   // time
0266   payloadName = outputDir + "/" + detector + "_meanTime" + "_" + m_dataset + "_" + m_run + ".root";
0267   if (h_CaloValid_ihcal_etaphi_time_raw->GetEntries())
0268   {
0269     histToCaloCDBTree(payloadName, "time", 1, h_CaloValid_ihcal_etaphi_time_raw.get());
0270   }
0271   // ZSCrossCalib
0272   payloadName = outputDir + "/" + detector + "_ZSCrossCalib" + "_" + m_dataset + "_" + m_run + ".root";
0273   if (h_CaloFittingQA_ihcal_etaphi_ZScrosscalib->GetEntries())
0274   {
0275     histToCaloCDBTree(payloadName, "ratio", 1, h_CaloFittingQA_ihcal_etaphi_ZScrosscalib.get());
0276   }
0277
0278   detector = "HCALOUT";
0279   // fracBadChi2
0280   payloadName = outputDir + "/" + detector + "_hotTowers_fracBadChi2" + "_" + m_dataset + "_" + m_run + ".root";
0281   if (h_CaloValid_ohcal_etaphi_badChi2->GetEntries())
0282   {
0283     histToCaloCDBTree(payloadName, "fraction", 1, h_CaloValid_ohcal_etaphi_badChi2.get());
0284   }
0285   // time
0286   payloadName = outputDir + "/" + detector + "_meanTime" + "_" + m_dataset + "_" + m_run + ".root";
0287   if (h_CaloValid_ohcal_etaphi_time_raw->GetEntries())
0288   {
0289     histToCaloCDBTree(payloadName, "time", 1, h_CaloValid_ohcal_etaphi_time_raw.get());
0290   }
0291   // ZSCrossCalib
0292   payloadName = outputDir + "/" + detector + "_ZSCrossCalib" + "_" + m_dataset + "_" + m_run + ".root";
0293   if (h_CaloFittingQA_ohcal_etaphi_ZScrosscalib->GetEntries())
0294   {
0295     histToCaloCDBTree(payloadName, "ratio", 1, h_CaloFittingQA_ohcal_etaphi_ZScrosscalib.get());
0296   }
0297 }
0298
0299 void GenStatus::process(const std::string &input, const std::string &output)
0300 {
0301   std::cout << "#############################" << std::endl;
0302   std::cout << "Run Parameters" << std::endl;
0303   std::cout << "input: " << input << std::endl;
0304   std::cout << "output: " << output << std::endl;
0305   std::cout << "#############################" << std::endl;
0306
0307   setRunDataset(input);
0308
0309   std::cout << "Processing: Run: " << m_run << ", Dataset: " << m_dataset << std::endl;
0310
0311   std::stringstream datasetDir;
0312   datasetDir.str("");
0313   datasetDir << output << "/" << m_run << "_" << m_dataset;
0314
0315   std::string outputDir = datasetDir.str();
0316
0317   std::string hotMapFile = "EMCalHotMap_" + m_dataset + "_" + m_run + ".root";
0318   std::string hotMapOutput = datasetDir.str() + "/" + hotMapFile;
0319
0320   datasetDir << "/QA";
0321
0322   // create output & QA directory
0323   std::filesystem::create_directories(datasetDir.str());
0324
0325   std::string hotMapOutputQA = datasetDir.str() + "/" + hotMapFile;
0326
0327   // merges individal qa into one per run
0328   int ret = readHists(input);
0329   if (ret)
0330   {
0331     return;
0332   }
0333
0334   analyze(outputDir);
0335
0336   std::unique_ptr<emcNoisyTowerFinder> calo = std::make_unique<emcNoisyTowerFinder>();
0337   calo->FindHot(m_CaloValid_list, hotMapOutput, "h_CaloValid_cemc_etaphi");
0338
0339   if (std::filesystem::exists(hotMapOutput))
0340   {
0341     std::filesystem::rename(hotMapOutput, hotMapOutputQA);
0342   }
0343   else
0344   {
0345     std::cout << "ERROR: EMCal Hot Map FAILED to Create." << std::endl;
0346   }
0347 }