Merge pull request #45758 from bsunanda/Run3-ala248B

Run3-alca248B Changes to the calibration scripts for HCAL Isotrack analysis of 2024 data sets

Author: cmsbuild

Calibration/HcalCalibAlgos/macros/CalibCorr.C

@@ -34,6 +34,17 @@
 //        double getCorr(entry): Entry # (in the file) dependent correction
 //        bool absent(entry) : if correction factor absent
 //        bool present(entry): or present (relevant for ML motivated)
+//        double getPhiCorr(unsigned int id): correction factor from phi
+//                                            symmery studies
+//        double getPhiCorrs(int run, unsigned int id): era-dependent correction
+//                                            factor from phi symmetry studies
+//      *infile* is the name of the input file
+//      *flag* can have a value between 0:5 specifying the type of correction
+//             5 for run-dependent correctons using results from several phi
+//             symmetry studies; 4 for using results from one phi-symmetry
+//             study; 3 for pileup correction using machine learning method;
+//             2 for overall response corrections; 1 for depth dependence
+//             corrections; 0 for raddam corrections)
 // CalibSelectRBX(rbxFile, debug)
 //      A class for selecting a given set of Read Out Box's and provides
 //        bool isItRBX(detId): if it/they is in the chosen RBXs
@@ -42,12 +53,17 @@
 //      A class for either rejecting duplicate entries or giving depth
 //        dependent weight. flag is 0 for keeping a list of duplicate
 //        emtries; 1 is to keep depth dependent weight for each ieta;
-//        2 is to keep a list of ieta, iphi for channels to be selected.
+//        2 is to keep a list of ieta, iphi for channels to be selected;
+//        3 is to read the gain correction file
 //        bool isDuplicate(entry): if it is a duplicate entry
 //        double getWeight(ieta, depth): get the dependent weight
 //        bool select(int ieta, int iphi): channels to be selected
+//        double getCorr(int run, int ieta, int depth): correction factor
+//        to modify the rechit energy due to change in gain fctor
 // void CalibCorrTest(infile, flag)
 //      Tests a file which contains correction factors used by CalibCorr
+//      if *flag* is >= 0; for negative values it initializes CalibDuplicate
+//      with -flag and *debug* = tue
 // void CalibCorrScale(infile, oufile, scale)
 //      Scales all contents of correction factors by "scale" from "infile"
 //      to "outfile"
@@ -591,16 +607,19 @@ public:
   double getCorr(const Long64_t& entry);
   double getTrueCorr(const Long64_t& entry);
   double getPhiCorr(unsigned int id);
+  double getPhiCorrs(int run, unsigned int id);
   bool absent(const Long64_t& entry);
   bool absent() { return (good_ == 0); }
   bool present(const Long64_t& entry);
 
 private:
+  unsigned int runId(const int& run);
   unsigned int readCorrRun(const char* infile);
   unsigned int readCorrDepth(const char* infile);
   unsigned int readCorrResp(const char* infile);
   unsigned int readCorrPU(const char* infile);
   unsigned int readCorrPhi(const char* infile);
+  unsigned int readCorrPhis(const char* infile);
   unsigned int getDetIdHE(int ieta, int iphi, int depth);
   unsigned int getDetId(int subdet, int ieta, int iphi, int depth);
   unsigned int correctDetId(const unsigned int& detId);
@@ -613,6 +632,8 @@ private:
   std::map<Long64_t, double> cfactors_;
   std::vector<int> runlow_;
   std::map<unsigned int, double> corrPhiSym_;
+  std::vector<std::pair<int, int> > runs_, ids_;
+  std::map<std::pair<unsigned int, unsigned int>, double> corrPhiSyms_;
 };
 
 class CalibSelectRBX {
@@ -636,15 +657,22 @@ public:
 
   bool isDuplicate(long entry);
   double getWeight(const unsigned int);
-  bool doCorr() { return ((flag_ == 1) && ok_); }
+  bool doCorr() { return (((flag_ == 1) || (flag_ == 3)) && ok_); }
+  bool doCorr(int flag) { return ((flag_ == flag) && ok_); }
   bool select(int ieta, int iphi);
+  double getCorr(int run, int ieta, int depth);
 
 private:
+  bool readCorrFactor(const char* fName);
+  unsigned int runId(const int& run);
+
   int flag_;
   double debug_, ok_;
   std::vector<Long64_t> entries_;
   std::map<int, std::vector<double> > weights_;
   std::vector<std::pair<int, int> > etaphi_;
+  std::vector<std::pair<int, int> > runs_, ids_;
+  std::map<std::pair<int, int>, std::vector<double> > corrs_;
 };
 
 CalibCorrFactor::CalibCorrFactor(const char* infile, int useScale, double scale, bool etamax, bool marina, bool debug)
@@ -758,6 +786,8 @@ CalibCorr::CalibCorr(const char* infile, int flag, bool debug) : flag_(flag), de
     good_ = readCorrPU(infile);
   else if (flag == 4)
     good_ = readCorrPhi(infile);
+  else if (flag == 5)
+    good_ = readCorrPhis(infile);
   else
     good_ = readCorrRun(infile);
 }
@@ -776,9 +806,9 @@ float CalibCorr::getCorr(int run, unsigned int id) {
     if (itr != corrFacResp_.end())
       cfac = itr->second;
   } else if (flag_ == 4) {
-    std::map<unsigned int, double>::iterator itr = corrPhiSym_.find(idx);
-    if (itr != corrPhiSym_.end())
-      cfac = itr->second;
+    cfac = getPhiCorr(idx);
+  } else if (flag_ == 5) {
+    cfac = getPhiCorrs(run, idx);
   } else {
     int ip(-1);
     for (unsigned int k = 0; k < runlow_.size(); ++k) {
@@ -842,6 +872,24 @@ double CalibCorr::getPhiCorr(unsigned int idx) {
   return cfac;
 }
 
+double CalibCorr::getPhiCorrs(int run, unsigned int idx) {
+  double cfac(1.0);
+  if (good_ == 0)
+    return cfac;
+  unsigned int runid = runId(run);
+  std::map<std::pair<unsigned int, unsigned int>, double>::iterator itr =
+      corrPhiSyms_.find(std::pair<unsigned int, unsigned int>(runid, idx));
+  if (itr != corrPhiSyms_.end())
+    cfac = itr->second;
+  if (debug_) {
+    int subdet, zside, ieta, iphi, depth;
+    unpackDetId(idx, subdet, zside, ieta, iphi, depth);
+    std::cout << "Run " << run << " ID " << runid << ":" << std::hex << idx << std::dec << " (Sub " << subdet << " eta "
+              << zside * ieta << " phi " << iphi << " depth " << depth << ")  Factor " << cfac << std::endl;
+  }
+  return cfac;
+}
+
 bool CalibCorr::absent(const Long64_t& entry) { return (cfactors_.find(entry) == cfactors_.end()); }
 
 bool CalibCorr::present(const Long64_t& entry) { return (cfactors_.find(entry) != cfactors_.end()); }
@@ -1057,6 +1105,82 @@ unsigned int CalibCorr::readCorrPhi(const char* infile) {
   return good;
 }
 
+unsigned int CalibCorr::readCorrPhis(const char* infile0) {
+  std::cout << "Enters readCorrPhis for " << infile0 << std::endl;
+  unsigned int all(0), good(0);
+  std::ifstream fIn(infile0);
+  if (!fIn.good()) {
+    std::cout << "Cannot open file " << infile0 << std::endl;
+  } else {
+    char buffer0[1024];
+    while (fIn.getline(buffer0, 1024)) {
+      std::string bufferString0(buffer0);
+      if (bufferString0.substr(0, 1) == "#") {
+        continue;  //ignore other comments
+      } else {
+        std::vector<std::string> item = splitString(bufferString0);
+        if (item.size() > 2) {
+          int run1 = std::atoi(item[0].c_str());
+          int run2 = std::atoi(item[1].c_str());
+          std::string infile(item[2]);
+          unsigned int runid = runs_.size();
+          runs_.push_back(std::pair<int, int>(run1, run2));
+          std::ifstream fInput(infile.c_str());
+          if (!fInput.good()) {
+            std::cout << "Cannot open file " << infile << std::endl;
+          } else {
+            char buffer[1024];
+            while (fInput.getline(buffer, 1024)) {
+              ++all;
+              std::string bufferString(buffer);
+              if (bufferString.substr(0, 1) == "#") {
+                continue;  //ignore other comments
+              } else {
+                std::vector<std::string> items = splitString(bufferString);
+                if (items.size() < 5) {
+                  std::cout << "Ignore  line: " << buffer << " Size " << items.size();
+                  for (unsigned int k = 0; k < items.size(); ++k)
+                    std::cout << " [" << k << "] : " << items[k];
+                  std::cout << std::endl;
+                } else {
+                  ++good;
+                  int subdet = std::atoi(items[0].c_str());
+                  int ieta = std::atoi(items[1].c_str());
+                  int iphi = std::atoi(items[2].c_str());
+                  int depth = std::atoi(items[3].c_str());
+                  double corrf = std::atof(items[4].c_str());
+                  unsigned int id = packDetId(subdet, ieta, iphi, depth);
+                  corrPhiSyms_[std::pair<int, int>(runid, id)] = corrf;
+                  if (debug_)
+                    std::cout << "ID " << runid << ":" << std::hex << id << std::dec << ":" << id << " (subdet "
+                              << subdet << " eta " << ieta << " phi " << iphi << " depth " << depth << ") "
+                              << corrPhiSym_[id] << std::endl;
+                }
+              }
+            }
+            fInput.close();
+          }
+        }
+      }
+    }
+    fIn.close();
+    std::cout << "Reads total of " << all << " and " << good << " good records of phi-symmetry factors from "
+              << runs_.size() << " files as given in " << infile0 << std::endl;
+  }
+  return good;
+}
+
+unsigned int CalibCorr::runId(const int& run) {
+  unsigned int id(runs_.size());
+  for (unsigned int k = 0; k < runs_.size(); ++k) {
+    if ((run >= runs_[k].first) && (run <= runs_[k].second)) {
+      id = k;
+      break;
+    }
+  }
+  return id;
+}
+
 unsigned int CalibCorr::getDetIdHE(int ieta, int iphi, int depth) { return getDetId(2, ieta, iphi, depth); }
 
 unsigned int CalibCorr::getDetId(int subdet, int ieta, int iphi, int depth) {
@@ -1256,6 +1380,8 @@ CalibDuplicate::CalibDuplicate(const char* fname, int flag, bool debug) : flag_(
           ok_ = true;
       }
     }
+  } else if (flag_ == 3) {
+    ok_ = readCorrFactor(fname);
   } else {
     flag_ = 2;
     if (strcmp(fname, "") != 0) {
@@ -1330,16 +1456,87 @@ bool CalibDuplicate::select(int ieta, int iphi) {
   return flag;
 }
 
+double CalibDuplicate::getCorr(int run, int ieta, int depth) {
+  std::map<std::pair<int, int>, std::vector<double> >::const_iterator itr =
+      corrs_.find(std::pair<int, int>(ieta, depth));
+  double corr(1.0);
+  if (itr != corrs_.end()) {
+    unsigned int irun = runId(run);
+    corr = (itr->second)[irun];
+  }
+  return corr;
+}
+
+bool CalibDuplicate::readCorrFactor(const char* infile) {
+  bool flag(false);
+  std::ifstream fInput(infile);
+  if (!fInput.good()) {
+    std::cout << "CalibDuplicate::readCorrFactor:Cannot open file " << infile << std::endl;
+  } else {
+    int nrun, neta;
+    fInput >> nrun >> neta;
+    int indx, run1, run2;
+    for (int k = 0; k < nrun; ++k) {
+      fInput >> indx >> run1 >> run2;
+      runs_.push_back(std::pair<int, int>(run1, run2));
+    }
+    int eta, depth;
+    double corr;
+    for (int k = 0; k < neta; ++k) {
+      fInput >> eta >> depth;
+      std::vector<double> corrs;
+      for (int i = 0; i < nrun; ++i) {
+        fInput >> corr;
+        corrs.push_back(corr);
+      }
+      corrs_[std::pair<int, int>(eta, depth)] = corrs;
+    }
+    fInput.close();
+    std::cout << "CalibDuplicate::readCorrFactor:Reads information of " << runs_.size() << " runs and " << corrs_.size()
+              << " channels from " << infile << std::endl;
+    flag = true;
+    if (debug_) {
+      for (unsigned int k = 0; k < runs_.size(); ++k)
+        std::cout << "Run range[" << k << "] " << runs_[k].first << ":" << runs_[k].second << std::endl;
+      std::map<std::pair<int, int>, std::vector<double> >::const_iterator itr;
+      for (itr = corrs_.begin(); itr != corrs_.end(); ++itr) {
+        std::cout << "eta:depth [" << (itr->first).first << ":" << (itr->first).second << "]";
+        for (unsigned int i = 0; i < (itr->second).size(); ++i)
+          std::cout << " " << (itr->second)[i];
+        std::cout << std::endl;
+      }
+    }
+  }
+  return flag;
+}
+
+unsigned int CalibDuplicate::runId(const int& run) {
+  unsigned int id(runs_.size());
+  for (unsigned int k = 0; k < runs_.size(); ++k) {
+    if ((run >= runs_[k].first) && (run <= runs_[k].second)) {
+      id = k;
+      break;
+    }
+  }
+  return id;
+}
+
 void CalibCorrTest(const char* infile, int flag) {
-  CalibCorr* c1 = new CalibCorr(infile, flag, true);
-  for (int ieta = 1; ieta < 29; ++ieta) {
-    int subdet = (ieta > 16) ? 2 : 1;
-    int depth = (ieta > 16) ? 2 : 1;
-    unsigned int id1 = ((4 << 28) | ((subdet & 0x7) << 25));
-    id1 |= ((0x1000000) | ((depth & 0xF) << 20) | (ieta << 10) | 1);
-    c1->getCorr(0, id1);
-    id1 |= (0x80000);
-    c1->getCorr(0, id1);
+  if (flag >= 0) {
+    CalibCorr* c1 = new CalibCorr(infile, flag, true);
+    for (int ieta = 1; ieta < 29; ++ieta) {
+      int subdet = (ieta > 16) ? 2 : 1;
+      int depth = (ieta > 16) ? 2 : 1;
+      unsigned int id1 = ((4 << 28) | ((subdet & 0x7) << 25));
+      id1 |= ((0x1000000) | ((depth & 0xF) << 20) | (ieta << 10) | 1);
+      c1->getCorr(0, id1);
+      id1 |= (0x80000);
+      c1->getCorr(0, id1);
+    }
+    delete c1;
+  } else {
+    CalibDuplicate* c1 = new CalibDuplicate(infile, -flag, true);
+    delete c1;
   }
 }