File indexing completed on 2025-08-05 08:12:18
0001
0002 #include <string>
0003 #include <iostream>
0004 #include <iomanip>
0005 #include <fstream>
0006 #include <unordered_set>
0007
0008
0009 #include <TSystem.h>
0010 #include <TROOT.h>
0011 #include <TH2F.h>
0012 #include <TFile.h>
0013 #include <TMath.h>
0014 #include <TDatime.h>
0015 #include <TGraph.h>
0016
0017 #include <calobase/TowerInfoDefs.h>
0018 #include <cdbobjects/CDBTTree.h>
0019 #include <ffamodules/CDBInterface.h>
0020 #include <phool/recoConsts.h>
0021
0022 using std::cout;
0023 using std::cerr;
0024 using std::endl;
0025 using std::string;
0026 using std::stringstream;
0027 using std::ofstream;
0028 using std::vector;
0029 using std::pair;
0030 using std::min;
0031 using std::max;
0032 using std::unordered_set;
0033
0034 R__LOAD_LIBRARY(libcalo_io.so)
0035 R__LOAD_LIBRARY(libcdbobjects.so)
0036 R__LOAD_LIBRARY(libffamodules.so)
0037 R__LOAD_LIBRARY(libphool.so)
0038
0039 namespace myAnalysis {
0040
0041 Int_t init(const string &input);
0042 void init_hists();
0043 Int_t readFile(const string &input);
0044
0045 Int_t process_event();
0046 void finalize(const string &i_output);
0047
0048 TH1F* hBadTowersHot;
0049
0050 TH1F* hSigmaNeverHot;
0051 TH1F* hSigmaHot;
0052 TH1F* hSigmaFreqHot;
0053 TH1F* hSigmaTypeA;
0054 TH1F* hSigmaTypeB;
0055 TH1F* hSigmaTypeC;
0056
0057 vector<TH1F*> hHotTowerSigma_vec;
0058
0059 vector<pair<UInt_t,string>> runs;
0060 vector<UInt_t> runsMissingHotMap;
0061 recoConsts* rc;
0062
0063 UInt_t m_ntowers = 24576;
0064 UInt_t m_bins_phi = 256;
0065 UInt_t m_bins_eta = 96;
0066
0067 UInt_t m_bins_sigma = 1000;
0068 Float_t m_sigma_low = -50;
0069 Float_t m_sigma_high = 50;
0070
0071 Int_t m_maxRuns = 0;
0072
0073
0074 UInt_t m_threshold = 400;
0075
0076 string m_detector = "CEMC";
0077
0078
0079
0080 string m_calibName_hotMap = m_detector + "_BadTowerMap";
0081 string m_fieldname_hotMap = "status";
0082 string m_sigma_hotMap = m_detector + "_sigma";
0083
0084 unordered_set<UInt_t> hotTypeA;
0085 unordered_set<UInt_t> hotTypeB;
0086 unordered_set<UInt_t> hotTypeC;
0087 }
0088
0089 void myAnalysis::init_hists() {
0090
0091 hBadTowersHot = new TH1F("hBadTowersHot", "Bad Towers: Hot; Tower Index; Runs", m_ntowers, -0.5, m_ntowers-0.5);
0092
0093 hSigmaHot = new TH1F("hSigmaHot", "Towers (Flagged Hot #leq 50% of Runs); sigma; Counts", m_bins_sigma, m_sigma_low, m_sigma_high);
0094 hSigmaNeverHot = new TH1F("hSigmaNeverHot", "Towers (Flagged Hot for 0% of Runs); sigma; Counts", m_bins_sigma, m_sigma_low, m_sigma_high);
0095 hSigmaFreqHot = new TH1F("hSigmaFreqHot", "Towers (Flagged Hot > 50% of Runs); sigma; Counts", m_bins_sigma, m_sigma_low, m_sigma_high);
0096 hSigmaTypeA = new TH1F("hSigmaTypeA", "Towers (Type A); sigma; Counts", m_bins_sigma, m_sigma_low, m_sigma_high);
0097 hSigmaTypeB = new TH1F("hSigmaTypeB", "Towers (Type B); sigma; Counts", m_bins_sigma, m_sigma_low, m_sigma_high);
0098 hSigmaTypeC = new TH1F("hSigmaTypeC", "Towers (Type C); sigma; Counts", m_bins_sigma, m_sigma_low, m_sigma_high);
0099
0100 stringstream name, title;
0101 for (UInt_t i = 0; i < m_ntowers; ++i) {
0102 UInt_t key = (m_detector == "CEMC") ? TowerInfoDefs::encode_emcal(i) : TowerInfoDefs::encode_hcal(i);
0103 UInt_t etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
0104 UInt_t phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
0105 name.str("");
0106 title.str("");
0107
0108 name << "hHotTowerSigma_" << phibin << "_" << etabin;
0109 title << "Hot Tower Sigma: (" << phibin << "," << etabin << "); sigma; Runs";
0110
0111 auto h = new TH1F(name.str().c_str(), title.str().c_str(), m_bins_sigma, m_sigma_low, m_sigma_high);
0112
0113 hHotTowerSigma_vec.push_back(h);
0114 }
0115 }
0116
0117 Int_t myAnalysis::init(const string &input) {
0118 rc = recoConsts::instance();
0119 rc->set_StringFlag("CDB_GLOBALTAG", "ProdA_2024");
0120 CDBInterface::instance()->Verbosity(0);
0121
0122 std::ifstream file(input);
0123
0124
0125 if (!file.is_open()) {
0126 cerr << "Failed to open file: " << input << endl;
0127 return 1;
0128 }
0129
0130 string line;
0131
0132
0133 while (std::getline(file, line)) {
0134 std::istringstream lineStream(line);
0135
0136 UInt_t runnumber;
0137 string timestamp;
0138 string date;
0139 string time;
0140 char comma;
0141
0142 if (lineStream >> runnumber >> comma >> date >> time) {
0143 timestamp = date+" "+time;
0144 runs.push_back(make_pair(runnumber,timestamp));
0145 }
0146 else {
0147 cerr << "Failed to parse line: " << line << endl;
0148 return 1;
0149 }
0150 }
0151
0152
0153 file.close();
0154
0155 init_hists();
0156
0157 return 0;
0158 }
0159
0160 Int_t myAnalysis::readFile(const string &input) {
0161 std::ifstream file(input);
0162
0163
0164 if (!file.is_open()) {
0165 cerr << "Failed to open file: " << input << endl;
0166 return 1;
0167 }
0168
0169 string line;
0170
0171
0172 std::getline(file, line);
0173
0174
0175 while (std::getline(file, line)) {
0176 std::istringstream lineStream(line);
0177
0178 char hotType;
0179 UInt_t phibin;
0180 UInt_t etabin;
0181 UInt_t towerIndex;
0182 char comma;
0183
0184 if (lineStream >> hotType >> comma >> phibin >> comma >> etabin) {
0185 towerIndex = (m_detector == "CEMC") ? TowerInfoDefs::decode_emcal(TowerInfoDefs::encode_emcal(etabin, phibin)) :
0186 TowerInfoDefs::decode_hcal(TowerInfoDefs::encode_hcal(etabin, phibin));
0187
0188 if(hotType == 'A') hotTypeA.insert(towerIndex);
0189 else if(hotType == 'B') hotTypeB.insert(towerIndex);
0190 else if(hotType == 'C') hotTypeC.insert(towerIndex);
0191 else {
0192 cout << "ERROR: Unknown Hot Tower Type: " << hotType << endl;
0193 return 2;
0194 }
0195 cout << "Hot Type: " << hotType << ", iphi: " << phibin << ", ieta: " << etabin << ", towerIndex: " << towerIndex << endl;
0196 }
0197 else {
0198 cerr << "Failed to parse line: " << line << endl;
0199 return 1;
0200 }
0201 }
0202
0203
0204 file.close();
0205
0206 return 0;
0207 }
0208
0209 Int_t myAnalysis::process_event() {
0210
0211 UInt_t bin_run = 1;
0212 Float_t sigma_min = 9999;
0213 Float_t sigma_max = 0;
0214 UInt_t ctr_overwrite = 0;
0215 unordered_set<UInt_t> badSigma;
0216
0217 for(auto p : runs) {
0218 if (m_maxRuns && bin_run > m_maxRuns) break;
0219
0220 UInt_t run = p.first;
0221 string timestamp = p.second;
0222 TDatime d;
0223
0224 cout << "Run: " << run << ", Timestamp: " << timestamp << endl;
0225
0226 rc->set_uint64Flag("TIMESTAMP", run);
0227
0228 Bool_t hasHotMap = true;
0229 string calibdir = CDBInterface::instance()->getUrl(m_calibName_hotMap);
0230
0231 std::unique_ptr<CDBTTree> m_cdbttree_hotMap = nullptr;
0232
0233 if (!calibdir.empty()) {
0234 m_cdbttree_hotMap = std::make_unique<CDBTTree>(calibdir);
0235 }
0236 else {
0237 cout << "Run: " << run << ", Missing: " << m_calibName_hotMap << endl;
0238 runsMissingHotMap.push_back(run);
0239 hasHotMap = false;
0240 }
0241
0242 Bool_t hasBadSigma = false;
0243 UInt_t acceptance = 0;
0244 UInt_t bad_ctr[4] = {0};
0245
0246 for (UInt_t channel = 0; channel < m_ntowers; ++channel) {
0247 UInt_t key = (m_detector == "CEMC") ? TowerInfoDefs::encode_emcal(channel) : TowerInfoDefs::encode_hcal(channel);
0248 UInt_t etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
0249 UInt_t phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
0250
0251 Int_t hotMap_val = 0;
0252 if (hasHotMap) {
0253 hotMap_val = m_cdbttree_hotMap->GetIntValue(key, m_fieldname_hotMap);
0254 }
0255
0256
0257
0258
0259
0260
0261 if(hotMap_val > 3) {
0262 cout << "ERROR: Unknown hotMap value: " << hotMap_val << endl;
0263 return 1;
0264 }
0265
0266 Float_t sigma_val = 0;
0267 if(!hasBadSigma && hasHotMap) {
0268 sigma_val = m_cdbttree_hotMap->GetFloatValue(key, m_sigma_hotMap);
0269 if(std::isnan(sigma_val) || std::isinf(sigma_val)) {
0270 cout << "WARNING: sigma does not exist for channel: " << channel << ", key: " << key << endl;
0271 hasBadSigma = true;
0272 badSigma.insert(run);
0273 }
0274 else {
0275 hHotTowerSigma_vec[channel]->Fill(sigma_val);
0276 }
0277 }
0278
0279 if (hasHotMap && hotMap_val != 0) {
0280
0281 if(hotMap_val == 2) {
0282 hBadTowersHot->Fill(channel);
0283 if(!hasBadSigma) {
0284 sigma_min = min(sigma_min, sigma_val);
0285 if(sigma_val < 100) sigma_max = max(sigma_max, sigma_val);
0286 }
0287 }
0288 }
0289 }
0290 ++bin_run;
0291 }
0292
0293 cout << "###################" << endl;
0294 cout << "Runs with Bad Sigma" << endl;
0295 for(auto run : badSigma) {
0296 cout << "Run: " << run << endl;
0297 }
0298
0299 cout << "###################" << endl;
0300
0301
0302 CDBInterface::instance()->Print();
0303
0304 UInt_t runs_processed = bin_run-runsMissingHotMap.size()-1;
0305 cout << "Sigma Min: " << sigma_min << ", Sigma Max: " << sigma_max << endl;
0306 cout << "Runs with bad sigma (nan or inf): " << badSigma.size() << endl;
0307 cout << "Runs Missing HotMap: " << runsMissingHotMap.size() << endl;
0308 cout << "Runs Processed: " << runs_processed << ", " << runs_processed*100./(bin_run-1) << " %" << endl;
0309
0310 return 0;
0311 }
0312
0313 void myAnalysis::finalize(const string &i_output) {
0314
0315 m_threshold = min(m_threshold, (UInt_t)(runs.size() - runsMissingHotMap.size()) / 2);
0316 if(m_maxRuns) m_threshold = min(m_threshold, (UInt_t)(m_maxRuns - runsMissingHotMap.size()) / 2);
0317 cout << "threshold: " << m_threshold << endl;
0318
0319 TFile output(i_output.c_str(), "recreate");
0320 output.cd();
0321
0322 for (UInt_t i = 0; i < m_ntowers; ++i) {
0323 if (hBadTowersHot->GetBinContent(i + 1) >= m_threshold) {
0324 hSigmaFreqHot->Add(hHotTowerSigma_vec[i]);
0325
0326 if(hotTypeA.find(i) != hotTypeA.end()) hSigmaTypeA->Add(hHotTowerSigma_vec[i]);
0327 if(hotTypeB.find(i) != hotTypeB.end()) hSigmaTypeB->Add(hHotTowerSigma_vec[i]);
0328 if(hotTypeC.find(i) != hotTypeC.end()) hSigmaTypeC->Add(hHotTowerSigma_vec[i]);
0329 }
0330 else hSigmaHot->Add(hHotTowerSigma_vec[i]);
0331 }
0332
0333 output.cd();
0334 hSigmaNeverHot->Write();
0335 hSigmaHot->Write();
0336 hSigmaFreqHot->Write();
0337 hSigmaTypeA->Write();
0338 hSigmaTypeB->Write();
0339 hSigmaTypeC->Write();
0340
0341 output.Close();
0342 }
0343
0344 void HotTowerTypeAnalysis(const string &input,
0345 const string &hotTypeList,
0346 const string &output = "hotType.root",
0347 const Int_t maxRuns = 0) {
0348
0349 cout << "#############################" << endl;
0350 cout << "Input Parameters" << endl;
0351 cout << "input: " << input << endl;
0352 cout << "Hot Type List: " << hotTypeList << endl;
0353 cout << "output: " << output << endl;
0354
0355 myAnalysis::m_maxRuns = maxRuns;
0356
0357 if (myAnalysis::m_detector == "HCALIN" || myAnalysis::m_detector == "HCALOUT") {
0358 myAnalysis::m_bins_phi = 64;
0359 myAnalysis::m_bins_eta = 24;
0360 myAnalysis::m_ntowers = myAnalysis::m_bins_phi * myAnalysis::m_bins_eta;
0361 }
0362
0363 cout << "Detector: " << myAnalysis::m_detector << endl;
0364
0365 Int_t ret = myAnalysis::init(input);
0366 if(ret) return;
0367
0368 ret = myAnalysis::readFile(hotTypeList);
0369 if(ret) return;
0370
0371 if(maxRuns) cout << "Processing: " << min((Int_t) myAnalysis::runs.size(), maxRuns) << " Runs" << endl;
0372 else cout << "Processing: " << myAnalysis::runs.size() << " Runs" << endl;
0373 cout << "#############################" << endl;
0374
0375 ret = myAnalysis::process_event();
0376 if(ret) return;
0377
0378 myAnalysis::finalize(output);
0379
0380 cout << "======================================" << endl;
0381 cout << "done" << endl;
0382 std::quick_exit(0);
0383 }
0384
0385 # ifndef __CINT__
0386 Int_t main(Int_t argc, char* argv[]) {
0387 if(argc < 3 || argc > 5) {
0388 cout << "usage: ./hotAna inputFile hotTypeList [outputFile] [maxRuns]" << endl;
0389 cout << "inputFile: containing list of run numbers" << endl;
0390 cout << "Hot Type List: hot tower list by types." << endl;
0391 cout << "outputFile: location of output file. Default: hotType.root." << endl;
0392 cout << "maxRuns: Max runs to process (useful when testing). Default: 0 (means run all)." << endl;
0393 return 1;
0394 }
0395
0396 string outputFile = "hotType.root";
0397 Int_t maxRuns = 0;
0398
0399 if (argc >= 4) {
0400 outputFile = argv[3];
0401 }
0402 if (argc >= 5) {
0403 maxRuns = atoi(argv[4]);
0404 }
0405
0406 HotTowerTypeAnalysis(argv[1], argv[2], outputFile, maxRuns);
0407
0408 return 0;
0409 }
0410 # endif
