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