Correct and extend the macros for HCAL isoTrac calibration and these …

…are now used in geeting response correction factors for 2023 data

Code check

Small fixes

Calibration/HcalCalibAlgos/macros/CalibCorr.C

@@ -648,7 +648,7 @@ CalibCorr::CalibCorr(const char* infile, int flag, bool debug) : flag_(flag), de
     good_ = readCorrResp(infile);
   else if (flag == 3)
     good_ = readCorrPU(infile);
-  else if (flag == 3)
+  else if (flag == 4)
     good_ = readCorrPhi(infile);
   else
     good_ = readCorrRun(infile);
@@ -1132,6 +1132,9 @@ double CalibDuplicate::getWeight(unsigned int detId) {
     int subdet, ieta, zside, depth, iphi;
     unpackDetId(detId, subdet, zside, ieta, iphi, depth);
     std::map<int, std::vector<double> >::const_iterator itr = weights_.find(ieta);
+    --depth;
+    if (depth < 0)
+      std::cout << "Strange Depth value " << depth << " in " << std::hex << detId << std::dec << std::endl;
     if (itr != weights_.end()) {
       if (depth < static_cast<int>(itr->second.size()))
         wt = itr->second[depth];

Calibration/HcalCalibAlgos/macros/CalibMonitor.C

@@ -296,9 +296,9 @@ private:
   const int etalo_, etahi_;
   int runlo_, runhi_;
   const int phimin_, phimax_, zside_, nvxlo_, nvxhi_, rbx_;
-  bool exclude_, corrE_, cutL1T_, selRBX_;
+  bool exclude_, cutL1T_, selRBX_;
   bool includeRun_;
-  int coarseBin_, etamp_, etamn_, plotType_;
+  int coarseBin_, plotType_;
   int flexibleSelect_, ifDepth_, duplicate_;
   double log2by18_;
   std::ofstream fileout_;
@@ -607,7 +607,7 @@ void CalibMonitor::Init(TChain *tree, const char *comFileName, const char *outFi
       nxbin = 1000;
   }
 
-  char name[20], title[200];
+  char name[100], title[500];
   std::string titl[5] = {
       "All tracks", "Good quality tracks", "Selected tracks", "Tracks with charge isolation", "Tracks MIP in ECAL"};
   for (int i = 0; i < 9; ++i)

Calibration/HcalCalibAlgos/macros/CalibPlotProperties.C

@@ -2036,59 +2036,52 @@ void CalibSplit::copyTree() {
   tout_qltyMissFlag = t_qltyMissFlag;
   tout_qltyPVFlag = t_qltyPVFlag;
   tout_gentrackP = t_gentrackP;
+  tout_DetIds->clear();
   if (t_DetIds != nullptr) {
     tout_DetIds->reserve(t_DetIds->size());
     for (unsigned int i = 0; i < t_DetIds->size(); ++i)
       tout_DetIds->push_back((*t_DetIds)[i]);
-  } else {
-    tout_DetIds->clear();
   }
+  tout_HitEnergies->clear();
   if (t_HitEnergies != nullptr) {
     tout_HitEnergies->reserve(t_HitEnergies->size());
     for (unsigned int i = 0; i < t_HitEnergies->size(); ++i)
       tout_HitEnergies->push_back((*t_HitEnergies)[i]);
-  } else {
-    tout_HitEnergies->clear();
   }
+  tout_trgbits->clear();
   if (t_trgbits != nullptr) {
     tout_trgbits->reserve(t_trgbits->size());
     for (unsigned int i = 0; i < t_trgbits->size(); ++i)
       tout_trgbits->push_back((*t_trgbits)[i]);
-  } else {
-    tout_trgbits->clear();
   }
+  tout_DetIds1->clear();
   if (t_DetIds1 != nullptr) {
     tout_DetIds1->reserve(t_DetIds1->size());
     for (unsigned int i = 0; i < t_DetIds1->size(); ++i)
       tout_DetIds1->push_back((*t_DetIds1)[i]);
-  } else {
-    tout_DetIds1->clear();
   }
+  tout_DetIds3->clear();
   if (t_DetIds3 != nullptr) {
     tout_DetIds3->reserve(t_DetIds3->size());
     for (unsigned int i = 0; i < t_DetIds3->size(); ++i)
       tout_DetIds3->push_back((*t_DetIds3)[i]);
-  } else {
-    tout_DetIds3->clear();
   }
+  tout_HitEnergies1->clear();
   if (t_HitEnergies1 != nullptr) {
     tout_HitEnergies1->reserve(t_HitEnergies1->size());
     for (unsigned int i = 0; i < t_HitEnergies1->size(); ++i)
       tout_HitEnergies1->push_back((*t_HitEnergies1)[i]);
-  } else {
-    tout_HitEnergies1->clear();
   }
+  tout_HitEnergies3->clear();
   if (t_HitEnergies1 != nullptr) {
     tout_HitEnergies3->reserve(t_HitEnergies3->size());
     for (unsigned int i = 0; i < t_HitEnergies3->size(); ++i)
       tout_HitEnergies3->push_back((*t_HitEnergies3)[i]);
-  } else {
-    tout_HitEnergies3->clear();
   }
 }
 
 void CalibSplit::close() {
-  if (!outputFile_) {
+  if (outputFile_) {
     outputDir_->cd();
     std::cout << "file yet to be Written" << std::endl;
     outputTree_->Write();

Calibration/HcalCalibAlgos/macros/CalibSort.C

@@ -109,6 +109,27 @@
 //   dirnm        (std::string) = name of the directory where Tree resides
 //                               (default "HcalIsoTrkAnalyzer")
 //
+// .L CalibSort.C+g
+//  combineML(const char* inputFileList, const char* outfile)
+//
+//  Combines the ML values otained from the analysis of muon analysis to
+//  determine depth dependent correction factors
+//
+//   inputFileList (const char*) = file containing filenames having the ML
+//                                 values for a given depth
+//   outfile       (const char*) = name of the output file where the
+//                                 depth dependent correction factors
+//                                 will be stored
+//   Example of a inputFileList:
+//      depth Name of the file
+//          1 ml_values_depth1.txt
+//          2 ml_values_depth2.txt
+//          3 ml_values_depth3.txt
+//          4 ml_values_depth4.txt
+//   where each file conatins
+//   (4 quantities per line with no tab separating each item):
+//   #ieta depth  ml    uncertainity-in-ml
+//
 //////////////////////////////////////////////////////////////////////////////
 
 #include <TCanvas.h>
@@ -166,6 +187,13 @@ struct recordEventLess {
   }
 };
 
+struct listML {
+  listML(int dep = 0, double ml = 0, double dml = 0) : depth_(dep), ml_(ml), dml_(dml) {}
+
+  int depth_;
+  double ml_, dml_;
+};
+
 class CalibSort {
 public:
   CalibSort(const char *fname,
@@ -2304,3 +2332,87 @@ void CalibMerge::close() {
   outputFile_ = nullptr;
   std::cout << "now doing return" << std::endl;
 }
+
+void combineML(const char *inputFileList, const char *outfile) {
+  std::map<int, std::string> files;
+  std::ifstream infile(inputFileList);
+  if (!infile.is_open()) {
+    std::cout << "** ERROR: Can't open '" << inputFileList << "' for input" << std::endl;
+  } else {
+    while (1) {
+      int depth;
+      std::string fname;
+      infile >> depth >> fname;
+      if (!infile.good())
+        break;
+      files[depth] = fname;
+    }
+    infile.close();
+  }
+  std::cout << "Gets a list of " << files.size() << " file names from " << inputFileList << std::endl;
+  if (files.size() > 0) {
+    std::map<int, std::vector<listML> > mlList;
+    for (std::map<int, std::string>::const_iterator itr = files.begin(); itr != files.end(); ++itr) {
+      int depth = itr->first;
+      std::string fname = itr->second;
+      std::ifstream fInput(fname.c_str());
+      if (!fInput.good()) {
+        std::cout << "Cannot open file " << fname << std::endl;
+      } else {
+        char buffer[1024];
+        unsigned int all(0), good1(0), good2(0);
+        while (fInput.getline(buffer, 1024)) {
+          ++all;
+          if (buffer[0] == '#')
+            continue;  //ignore comment
+          std::vector<std::string> items = splitString(std::string(buffer));
+          if (items.size() != 4) {
+            std::cout << "Ignore  line: " << buffer << std::endl;
+          } else {
+            ++good1;
+            int depth0 = std::atoi(items[1].c_str());
+            if (depth0 == depth) {
+              ++good2;
+              int ieta = std::atoi(items[0].c_str());
+              double ml = std::atoi(items[2].c_str());
+              double dml = std::atoi(items[3].c_str());
+              listML l0(depth, ml, dml);
+              if (mlList.find(ieta) == mlList.end()) {
+                std::vector<listML> l0v;
+                mlList[ieta] = l0v;
+              }
+              (mlList[ieta]).push_back(l0);
+            }
+          }
+        }
+        fInput.close();
+        std::cout << "Reads total of " << all << " and " << good1 << ":" << good2 << " good records for depth " << depth
+                  << std::endl;
+      }
+    }
+    if (mlList.size() > 0) {
+      std::ofstream fout(outfile);
+      for (std::map<int, std::vector<listML> >::const_iterator itr = mlList.begin(); itr != mlList.end(); ++itr) {
+        int ieta = itr->first;
+        std::vector<listML> l0v = itr->second;
+        double den(0), dden(0);
+        for (unsigned int k = 0; k < l0v.size(); ++k) {
+          if (l0v[k].depth_ == 2) {
+            den = l0v[k].ml_;
+            dden = l0v[k].dml_;
+          }
+        }
+        if (den > 0) {
+          fout << std::setw(4) << ieta << "   " << l0v.size();
+          for (unsigned int k = 0; k < l0v.size(); ++k) {
+            double ml = den / l0v[k].ml_;
+            double dml = l0v[k].dml_ * den / (l0v[k].ml_ * l0v[k].ml_);
+            fout << "  " << l0v[k].depth_ << "  " << std::setw(6) << ml << "  " << std::setw(6) << dml;
+          }
+          fout << std::endl;
+        }
+      }
+      fout.close();
+    }
+  }
+}