hcal_calib_year2/genFinalCalib/genFinalCalib_pp24.C

0001 #include "tsc_cos_merge.C"  // provides genCdbCorr_HH / genCdbCorr
0002 #include "getCosmicTemps.C" // This macro reads the temp of cosmic runs from DB and it takes few min
0003
0004 #include <TSystem.h>
0005 #include <TFile.h>
0006 #include <TH1F.h>
0007 #include <TString.h>
0008
0009 #include <iostream>
0010 #include <format>
0011 #include <fstream>
0012 #include <sstream>
0013 #include <vector>
0014 #include <string>
0015 #include <algorithm>
0016 #include <cstdlib>
0017 #include <climits>
0018
0019
0020 // Input files needed
0021 // cos_runs_pp_2024.txt, runList_pp_2024_ana509.txt, anaOut.root, cdbFiles/ohcal_cdb_tsc_cos_calib_%d.root
0022
0023 float fetch_ohcal_temp_degC(int run,
0024                 int max_retries = 3,
0025                 int backoff_ms  = 400)
0026 {
0027   // default detector is OHCal
0028   return fetch_hcal_temp_degC(run, "HCALOUT", max_retries, backoff_ms);
0029 }
0030
0031 struct TowerData {
0032   int   tower_ieta;
0033   int   tower_iphi;
0034   float calib_factor;
0035   int   run_range_low;
0036   int   run_range_high;
0037 };
0038
0039 void loadCSV(const std::string& filename, std::vector<TowerData>& data);
0040 std::vector<std::vector<float>> getTowersForRun(const std::vector<TowerData>& data, int run_number);
0041
0042 void genFinalCalib_pp24(const char *runs = "cos_runs_pp_2024.txt")
0043 {
0044   std::ofstream csv_log("beam_cosmic_tempcorr_log_nov22_abscosmic.csv");
0045   csv_log << "beam_run,cosmic_run,t_beam,t_cosmic,correction\n";
0046
0047   gSystem->mkdir("final_ohcal_calib_run2pp_hh_abscosmic_nov23", kTRUE);
0048
0049   // ---- read cosmic runs ----
0050   std::ifstream runlist(runs);
0051   std::vector<int> cosmic_runnumbers;
0052   std::string line;
0053   while (getline(runlist,line)) {
0054     if (line.empty() || line[0]=='#') continue;
0055     cosmic_runnumbers.push_back(stoi(line));
0056   }
0057   sort(cosmic_runnumbers.begin(), cosmic_runnumbers.end());
0058
0059   // ---- read beam runs ----
0060   std::ifstream runlistBeam("runList_pp_2024_ana509.txt");
0061   std::vector<int> runnumbersBeam;
0062   while (getline(runlistBeam, line)) {
0063     if (line.empty() || line[0]=='#') continue;
0064     runnumbersBeam.push_back(stoi(line));
0065   }
0066
0067   const int RUN_MAX = runnumbersBeam.back(); // get the last element in the list
0068   std::cout << "last run in the list of beam runs is:" << RUN_MAX << std::endl;
0069
0070   std::vector<int> runsWritten;   // to collect missing runs
0071
0072   // ---- open anaOut.root ----
0073   TFile* ftemp = TFile::Open("analyze_hits/anaOut.root","READ"); // input file from ana_hits.C macro
0074   TH1F* h_runnumbers_temp = (TH1F*) ftemp->Get("runNumberHist");
0075   TH1F* h_temp_run        = (TH1F*) ftemp->Get("h_temp_run");
0076   int nBins = h_runnumbers_temp->GetNbinsX();
0077
0078   // ---- get temperature for cosmic runs (approx: nearest beam bin) ----
0079   /*
0080   std::vector<float> cosmic_temps;
0081   for (int cos_run : cosmic_runnumbers) {
0082     int minDelta = INT_MAX;
0083     float cos_temp = 0;
0084     int matched_beam = -1;
0085     for (int b=1; b<=nBins; b++) {
0086       int   beam_run = (int)h_runnumbers_temp->GetBinContent(b);
0087       float tbin     =      h_temp_run->GetBinContent(b);
0088       int delta = abs(cos_run - beam_run);
0089       if (delta < minDelta && tbin > 10) {
0090         minDelta = delta;
0091         cos_temp = tbin;
0092     matched_beam = beam_run;
0093       }
0094     }
0095     cosmic_temps.push_back(cos_temp);
0096   // --- print debug info ---
0097   std::cout << "[Cosmic] Run " << cos_run
0098        << " → nearest beam run " << matched_beam
0099        << " with temperature = " << cos_temp << " C"
0100        << std::endl;
0101   }
0102   */
0103
0104   // ---- get temperature for cosmic runs from DB via getCosmicTemps.C ----
0105
0106   std::vector<float> cosmic_temps;
0107   cosmic_temps.reserve(cosmic_runnumbers.size());
0108
0109   for (int cos_run : cosmic_runnumbers) {
0110     float cos_temp = fetch_ohcal_temp_degC(cos_run);   // function from getCosmicTemps.C
0111     cosmic_temps.push_back(cos_temp);
0112
0113     std::cout << "[Cosmic] Run " << cos_run
0114          << " | extracted OHCal temperature = " << cos_temp << " C"
0115          << std::endl;
0116   }
0117
0118
0119   // ---- load half-height corrections ----
0120   std::vector<TowerData> hh_data;
0121   loadCSV("halfHeights_pp_2024_anjaly_dec6.csv", hh_data);
0122
0123   // ---- loop over beam runs ----
0124   for (size_t irb=0; irb<runnumbersBeam.size(); irb++) {
0125     int runnumber = runnumbersBeam[irb];
0126
0127     int runnumber_last;
0128     if (irb+1  < runnumbersBeam.size() ) runnumber_last = runnumbersBeam[irb+1]-1;
0129     else runnumber_last = RUN_MAX;
0130
0131     std::cout << "----------------------------------" << std::endl;
0132     std::cout << "starting run " << runnumber << std::endl;
0133
0134     int asoCosmic = 0;
0135     int minDelta  = INT_MAX;
0136     float asoCosmicTemp = 0;
0137     for (size_t ic=0; ic<cosmic_runnumbers.size(); ic++) {
0138       int cos_run = cosmic_runnumbers[ic];
0139       int delta   =  std::abs(runnumber - cos_run);  // absolute difference
0140       //int delta   =  runnumber - cos_run; //  if (delta < minDelta && delta <= 0)
0141       if (delta < minDelta)  {
0142         minDelta       = delta;
0143         asoCosmic      = cos_run;
0144         asoCosmicTemp  = cosmic_temps[ic];
0145       }
0146     }
0147
0148     // get beam temp directly from histogram
0149     float temp_beam = 0;
0150     for (int it = 1; it <=nBins; it++) {
0151       if ((int)h_runnumbers_temp->GetBinContent(it) == runnumber) {
0152         temp_beam = h_temp_run->GetBinContent(it);
0153         break;
0154       }
0155     }
0156
0157     if (temp_beam < 10) {
0158       std::cout << "Run " << runnumber << " skipped (temp < 10)" << std::endl;
0159       continue;
0160     }
0161
0162     float alpha = 0.045;
0163     float deltaT = temp_beam - asoCosmicTemp;
0164     float correction = 1 + alpha*deltaT;
0165
0166     std::cout << "Beam run " << runnumber
0167          << " uses cosmic " << asoCosmic
0168          << " : T_beam=" << temp_beam
0169          << " T_cosmic=" << asoCosmicTemp
0170          << " corr=" << correction << std::endl;
0171
0172     csv_log << runnumber << "," << asoCosmic << "," << temp_beam << "," << asoCosmicTemp << "," << correction << "\n";
0173
0174
0175     std::string cosmic_CDB_file_oh = std::format("cdbFiles/ohcal_cdb_tsc_cos_calib_{}.root",asoCosmic);  //input file needed
0176     auto hh_towersForRun = getTowersForRun(hh_data, runnumber);
0177     std::string cdbFileName_hh = std::format("final_ohcal_calib_run2pp_hh_abscosmic_nov23/cdb_ohcal_tsc_temp_cosmic_hh_{}_{}.root",
0178                                  runnumber, runnumber_last); //output file name
0179     genCdbCorr(correction, cosmic_CDB_file_oh, cdbFileName_hh, 0);                       // Without HH correction
0180     //genCdbCorr_HH(correction, cosmic_CDB_file_oh, cdbFileName_hh, 0, hh_towersForRun); // With HH correction
0181     runsWritten.push_back(runnumber);
0182   } //beam run loop ends
0183
0184   csv_log << "skipped_run\n";
0185   std::cout << "\n=== Summary of runs with no file written ===\n";
0186   int skippedCount = 0;
0187   for (int r : runnumbersBeam) {
0188     if (std::find(runsWritten.begin(), runsWritten.end(), r) == runsWritten.end()) {
0189         std::cout << r << " ";
0190     csv_log << r << "\n";
0191     skippedCount++;
0192     }
0193 }
0194   std::cout << "\nTotal skipped runs: " << skippedCount << std::endl;
0195   csv_log.close();
0196
0197
0198 }//end macro
0199
0200 // =============== helpers ===============
0201
0202 void loadCSV(const std::string& filename, std::vector<TowerData>& data) {
0203   std::ifstream file(filename);
0204   std::string line;
0205   while (getline(file, line)) {
0206     if (line.empty() || line[0]=='#') continue;
0207     std::stringstream ss(line);
0208     TowerData t{};
0209     char comma;
0210     ss >> t.tower_ieta >> comma
0211        >> t.tower_iphi >> comma
0212        >> t.calib_factor >> comma
0213        >> t.run_range_low >> comma
0214        >> t.run_range_high;
0215     data.push_back(t);
0216   }
0217 }
0218
0219 std::vector<std::vector<float>> getTowersForRun(const std::vector<TowerData>& data, int run_number) {
0220   std::vector<std::vector<float>> towersForRun;
0221   for (const auto& t : data) {
0222     if (run_number >= t.run_range_low && run_number <= t.run_range_high) {
0223       towersForRun.push_back({
0224         static_cast<float>(t.tower_ieta),
0225         static_cast<float>(t.tower_iphi),
0226         t.calib_factor
0227       });
0228     }
0229   }
0230   return towersForRun;
0231 }
0232