Add the new ways to introduce PF cuts frm a file derivd from DB

Sunanda · Sunanda · commit 648ddf4157eb · 2025-06-28T11:53:51.000+02:00
diff --git a/Calibration/HcalCalibAlgos/macros/CalibCorr.C b/Calibration/HcalCalibAlgos/macros/CalibCorr.C
@@ -238,28 +238,6 @@ int truncateDepth(int ieta, int depth, int truncateFlag) {
   return d;
 }
 
-double threshold(int subdet, int depth, int form) {
-  double cutHE[4][7] = {{0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2},
-                        {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2},
-                        {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2},
-                        {0.2, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3}};
-  double cutHB[4][4] = {{0.1, 0.2, 0.3, 0.3}, {0.25, 0.25, 0.3, 0.3}, {0.4, 0.3, 0.3, 0.3}, {0.6, 0.4, 0.4, 0.5}};
-  double thr(0);
-  if (form > 0) {
-    if (subdet == 2)
-      thr = cutHE[form - 1][depth - 1];
-    else
-      thr = cutHB[form - 1][depth - 1];
-  }
-  return thr;
-}
-
-double threshold(unsigned int detId, int form) {
-  int subdet = ((detId >> 25) & (0x7));
-  int depth = ((detId & 0x1000000) == 0) ? ((detId >> 14) & 0x1F) : ((detId >> 20) & 0xF);
-  return threshold(subdet, depth, form);
-}
-
 double puFactor(int type, int ieta, double pmom, double eHcal, double ediff, bool debug = false) {
   double fac(1.0);
   if (debug)
@@ -568,6 +546,98 @@ std::vector<std::string> splitString(const std::string& fLine) {
   return result;
 }
 
+class CalibThreshold {
+public:
+  CalibThreshold(int form);
+  ~CalibThreshold() {}
+
+  double threshold(unsigned int detId);
+private:
+  double threshold(int subdet, int ieta, int depth);
+  bool fileThreshold(const char* fname);
+
+  int form_;
+  std::map<std::pair<int, int>, double> thresh_;
+  bool ok_;
+};
+
+CalibThreshold::CalibThreshold(int form) : form_(form) {
+  if (form_ == 5)
+    ok_ = fileThreshold("PFCuts362975.txt");
+  else if ((form_ < 1) || (form_ > 5))
+    ok_ = false;
+  else
+    ok_ = true;
+  std::cout << "CalibThreshold initialized with flag " << ok_ << " for form " << form_ << std::endl;
+}
+
+bool CalibThreshold::fileThreshold(const char* fname) {
+  bool ok(false);
+  if (std::string(fname) != "") {
+    std::ifstream fInput(fname);
+    if (!fInput.good()) {
+      std::cout << "Cannot open file " << fname << std::endl;
+    } else {
+      char buffer[1024];
+      unsigned int all(0), good(0), bad(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() != 7) {
+	  ++bad;
+          std::cout << "Ignore  line: " << buffer << std::endl;
+        } else {
+          ++good;
+	  int ieta = std::atoi(items[0].c_str());
+	  int depth = std::atoi(items[2].c_str());
+	  double thr = std::atof(items[4].c_str());
+	  thresh_[std::pair<int, int>(ieta, depth)] = thr;
+	}
+      }
+      fInput.close();
+      std::cout << "Reads " << all << " entries from " << fname << " with " << good << " Good and " << bad << " bad records" << std::endl;
+      if (good > 0)
+	ok = true;
+    }
+  }
+  return ok;
+}
+
+double CalibThreshold::threshold(int subdet, int ieta, int depth) {
+  double cutHE[4][7] = {{0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2},
+                        {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2},
+                        {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2},
+                        {0.2, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3}};
+  double cutHB[4][4] = {{0.1, 0.2, 0.3, 0.3}, 
+			{0.25, 0.25, 0.3, 0.3}, 
+			{0.4, 0.3, 0.3, 0.3}, 
+			{0.6, 0.4, 0.4, 0.5}};
+
+  double thr(0);
+  if (ok_) {
+    if ((form_ > 0) && (form_ <= 4)) {
+      if (subdet == 2)
+	thr = cutHE[form_ - 1][depth - 1];
+      else
+	thr = cutHB[form_ - 1][depth - 1];
+    } else {
+      std::map<std::pair<int, int>, double>::const_iterator itr =
+	thresh_.find(std::pair<int, int>(ieta, depth));
+      if (itr != thresh_.end())
+	thr = itr->second;
+    }
+  }
+  return thr;
+}
+
+double CalibThreshold::threshold(unsigned int detId) {
+  int subdet, zside, ieta, iphi, depth;
+  unpackDetId(detId, subdet, zside, ieta, iphi, depth);
+  return threshold(subdet, (zside * ieta), depth);
+}
+
 class CalibCorrFactor {
 public:
   CalibCorrFactor(const char* infile, int useScale, double scale, bool etamax, bool marina, bool debug);
diff --git a/Calibration/HcalCalibAlgos/macros/CalibFitPlots.C b/Calibration/HcalCalibAlgos/macros/CalibFitPlots.C
@@ -265,6 +265,7 @@
 #include "CalibCorr.C"
 
 const double fitrangeFactor = 1.5;
+const double fitrangeFactor1 = 1.2;
 
 struct cfactors {
   int ieta, depth;
@@ -445,7 +446,7 @@ results fitTwoGauss(TH1D* hist, bool debug) {
   std::pair<double, double> mrms = GetMean(hist, 0.2, 2.0, rms);
   double mean = mrms.first;
   double LowEdge = mean - fitrangeFactor * rms;
-  double HighEdge = mean + fitrangeFactor * rms;
+  double HighEdge = mean + fitrangeFactor1 * rms;
   if (LowEdge < 0.15)
     LowEdge = 0.15;
   std::string option = (hist->GetEntries() > 100) ? "QRS" : "QRWLS";
@@ -482,7 +483,7 @@ results fitTwoGauss(TH1D* hist, bool debug) {
   highValue[5] = 100. * startvalues[5];
   //fitrange[0] = mean - 2.0*rms; fitrange[1] = mean + 2.0*rms;
   fitrange[0] = Fit->Value(1) - fitrangeFactor * Fit->Value(2);
-  fitrange[1] = Fit->Value(1) + fitrangeFactor * Fit->Value(2);
+  fitrange[1] = Fit->Value(1) + fitrangeFactor1 * Fit->Value(2);
   TFitResultPtr Fitfun = functionFit(hist, fitrange, startvalues, lowValue, highValue);
   double wt1 = (Fitfun->Value(0)) * (Fitfun->Value(2));
   double value1 = Fitfun->Value(1);
@@ -514,7 +515,8 @@ results fitOneGauss(TH1D* hist, bool fitTwice, bool debug) {
   std::pair<double, double> mrms = GetMean(hist, 0.2, 2.0, rms);
   double mean = mrms.first;
   double LowEdge = ((mean - fitrangeFactor * rms) < 0.5) ? 0.5 : (mean - fitrangeFactor * rms);
-  double HighEdge = (mean + fitrangeFactor * rms);
+  double diff = mean - LowEdge;
+  double HighEdge = (diff > fitrangeFactor1 * rms) ? (mean + fitrangeFactor1 * rms) : (mean + diff);
   if (debug)
     std::cout << hist->GetName() << " Mean " << mean << " RMS " << rms << " Range " << LowEdge << ":" << HighEdge
               << "\n";
@@ -531,9 +533,10 @@ results fitOneGauss(TH1D* hist, bool fitTwice, bool debug) {
   double werror = Fit1->FitResult::Error(2);
   if (fitTwice) {
     LowEdge = Fit1->Value(1) - fitrangeFactor * Fit1->Value(2);
-    HighEdge = Fit1->Value(1) + fitrangeFactor * Fit1->Value(2);
     if (LowEdge < 0.5)
       LowEdge = 0.5;
+    diff = Fit1->Value(1) - LowEdge;
+    HighEdge = (diff > fitrangeFactor1 * rms) ? (Fit1->Value(1) + fitrangeFactor1 * Fit1->Value(2)) : (Fit1->Value(1) + diff);
     if (HighEdge > 5.0)
       HighEdge = 5.0;
     if (debug)
@@ -2220,6 +2223,7 @@ void PlotHistCorrFactor(char* infile,
   int nbin = etamax - etamin + 1;
   std::vector<TH1D*> hists;
   std::vector<int> entries;
+  std::vector<int> depths;
   char name[100];
   double dy(0);
   int fits(0);
@@ -2254,17 +2258,20 @@ void PlotHistCorrFactor(char* infile,
       TF1* func = new TF1(name, "pol0", etamin, etamax);
       h->Fit(func, "+QWLR", "");
     }
-    h->SetLineColor(colors[j]);
-    h->SetMarkerColor(colors[j]);
-    h->SetMarkerStyle(mtype[j]);
-    h->GetXaxis()->SetTitle("i#eta");
-    h->GetYaxis()->SetTitle("Correction Factor");
-    h->GetYaxis()->SetLabelOffset(0.005);
-    h->GetYaxis()->SetTitleOffset(1.20);
-    h->GetYaxis()->SetRangeUser(0.0, 2.0);
-    hists.push_back(h);
-    entries.push_back(nent);
-    dy += 0.025;
+    if (nent > 0) {
+      h->SetLineColor(colors[j]);
+      h->SetMarkerColor(colors[j]);
+      h->SetMarkerStyle(mtype[j]);
+      h->GetXaxis()->SetTitle("i#eta");
+      h->GetYaxis()->SetTitle("Correction Factor");
+      h->GetYaxis()->SetLabelOffset(0.005);
+      h->GetYaxis()->SetTitleOffset(1.20);
+      h->GetYaxis()->SetRangeUser(0.0, 2.5);
+      hists.push_back(h);
+      entries.push_back(nent);
+      depths.push_back(j + 1);
+      dy += 0.025;
+    }
   }
   sprintf(name, "c_%sCorrFactor", prefixF.c_str());
   TCanvas* pad = new TCanvas(name, name, 700, 500);
@@ -2292,7 +2299,7 @@ void PlotHistCorrFactor(char* infile,
       st1->SetX2NDC(0.90);
       yh -= dy;
     }
-    sprintf(name, "Depth %d (%s)", k + 1, text.c_str());
+    sprintf(name, "Depth %d (%s)", depths[k], text.c_str());
     legend->AddEntry(hists[k], name, "lp");
   }
   legend->Draw("same");
@@ -2493,7 +2500,7 @@ void PlotHistCorrAsymmetry(
   std::vector<int> entries;
   char name[100];
   double dy(0);
-  int maxd = (depth < 0) ? maxdepth : 1;
+  int maxd = (depth <= 0) ? maxdepth : 1;
   for (int j = 0; j < maxd; ++j) {
     int dep = (depth <= 0) ? (j + 1) : depth;
     sprintf(name, "hd%d", dep);
@@ -4330,7 +4337,7 @@ void PlotMeanError(const std::string infilest, int reg = 3, bool resol = false,
   } else {
     ok = true;
     fInput >> nEner >> nType >> nPts;
-    int nmax = nEner * nType;
+    const int nmax = nEner * nType;
     int type, elow, ehigh;
     double v1[4], e1[4], v2[4], e2[4];
     for (int n = 0; n < nmax; ++n) {
@@ -4406,7 +4413,7 @@ void PlotMeanError(const std::string infilest, int reg = 3, bool resol = false,
     legend->SetFillColor(kWhite);
     std::string nameg[ntypmx] = {"MAHI", "M0", "M2"};
     for (int n = 0; n < nType; ++n) {
-      unsigned int nmax0 = energy[n].size();
+      const unsigned int nmax0 = energy[n].size();
       double mom[nmax0], dmom[nmax0], mean[nmax0], dmean[nmax0];
       for (unsigned int k = 0; k < nmax0; ++k) {
         mom[k] = energy[n][k];
diff --git a/Calibration/HcalCalibAlgos/macros/CalibMonitor.C b/Calibration/HcalCalibAlgos/macros/CalibMonitor.C
@@ -314,6 +314,7 @@ private:
   CalibCorr *cFactor_;
   CalibSelectRBX *cSelect_;
   CalibDuplicate *cDuplicate_;
+  CalibThreshold *cThr_;
   const std::string fname_, dirnm_, prefix_, outFileName_;
   const int corrPU_, flag_, numb_;
   const bool isRealData_, useGen_;
@@ -370,6 +371,7 @@ CalibMonitor::CalibMonitor(const char *fname,
       cFactor_(nullptr),
       cSelect_(nullptr),
       cDuplicate_(nullptr),
+      cThr_(nullptr),
       fname_(std::string(fname)),
       dirnm_(std::string(dirnm)),
       prefix_(prefix),
@@ -440,6 +442,8 @@ CalibMonitor::CalibMonitor(const char *fname,
     }
     if (std::string(rbxFile) != "")
       cSelect_ = new CalibSelectRBX(rbxFile, false);
+    if (thrForm_ > 0)
+      cThr_ = new CalibThreshold(thrForm_);
   }
 }
 
@@ -1150,7 +1154,7 @@ void CalibMonitor::Loop(Long64_t nmax, bool debug) {
       eHcal = 0;
       for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
         // Apply thresholds if necessary
-        bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[k] > threshold((*t_DetIds)[k], thrForm_));
+        bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[k] > (cThr_->threshold((*t_DetIds)[k])));
         // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
         if (okcell) {
           double cfac(1.0);
@@ -1572,7 +1576,7 @@ bool CalibMonitor::selectPhi(bool debug) {
     // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
     for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
       // Apply thresholds if necessary
-      bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[k] > threshold((*t_DetIds)[k], thrForm_));
+      bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[k] > (cThr_->threshold((*t_DetIds)[k])));
       if (okcell) {
         int iphi = ((*t_DetIds)[k]) & (0x3FF);
         int zside = ((*t_DetIds)[k] & 0x80000) ? (1) : (-1);
@@ -1760,7 +1764,7 @@ void CalibMonitor::correctEnergy(double &eHcal, const Long64_t &entry) {
       // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
       for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
         // Apply thresholds if necessary
-        bool okcell = (thrForm_ == 0) || ((*t_HitEnergies1)[idet] > threshold((*t_DetIds1)[idet], thrForm_));
+        bool okcell = (thrForm_ == 0) || ((*t_HitEnergies1)[idet] > (cThr_->threshold((*t_DetIds1)[idet])));
         if (okcell) {
           unsigned int id = truncateId((*t_DetIds1)[idet], truncateFlag_, false);
           double cfac = corrFactor_->getCorr(id);
@@ -1779,7 +1783,7 @@ void CalibMonitor::correctEnergy(double &eHcal, const Long64_t &entry) {
       }
       for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
         // Apply thresholds if necessary
-        bool okcell = (thrForm_ == 0) || ((*t_HitEnergies3)[idet] > threshold((*t_DetIds3)[idet], thrForm_));
+        bool okcell = (thrForm_ == 0) || ((*t_HitEnergies3)[idet] > (cThr_->threshold((*t_DetIds3)[idet])));
         if (okcell) {
           unsigned int id = truncateId((*t_DetIds3)[idet], truncateFlag_, false);
           double cfac = corrFactor_->getCorr(id);
diff --git a/Calibration/HcalCalibAlgos/macros/CalibPlotProperties.C b/Calibration/HcalCalibAlgos/macros/CalibPlotProperties.C
@@ -329,6 +329,7 @@ private:
   CalibCorr *cFactor_;
   CalibSelectRBX *cSelect_;
   CalibDuplicate *cDuplicate_;
+  CalibThreshold *cThr_;
   const std::string fname_, dirnm_, prefix_, outFileName_;
   const int corrPU_, flag_;
   const bool isRealData_, useGen_;
@@ -387,6 +388,7 @@ CalibPlotProperties::CalibPlotProperties(const char *fname,
       cFactor_(nullptr),
       cSelect_(nullptr),
       cDuplicate_(nullptr),
+      cThr_(nullptr),
       fname_(fname),
       dirnm_(dirnm),
       prefix_(prefix),
@@ -452,6 +454,8 @@ CalibPlotProperties::CalibPlotProperties(const char *fname,
       cDuplicate_ = new CalibDuplicate(dupFileName, duplicate_, false);
     if (std::string(rbxFile) != "")
       cSelect_ = new CalibSelectRBX(rbxFile, false);
+    if (thrForm_ > 0)
+      cThr_ = new CalibThreshold(thrForm_);
   }
 }
 
@@ -923,7 +927,7 @@ void CalibPlotProperties::Loop(Long64_t nentries) {
       eHcal = 0;
       for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
         // Apply thresholds if necessary
-        bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[k] > threshold((*t_DetIds)[k], thrForm_));
+        bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[k] > (cThr_->threshold((*t_DetIds)[k])));
         if (okcell) {
           // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
           unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
@@ -1024,7 +1028,7 @@ void CalibPlotProperties::Loop(Long64_t nentries) {
             double eb(0), ee(0);
             for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
               // Apply thresholds if necessary
-              bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[k] > threshold((*t_DetIds)[k], thrForm_));
+              bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[k] > (cThr_->threshold((*t_DetIds)[k])));
               if (okcell) {
                 unsigned int id = truncateId((*t_DetIds)[k], truncateFlag_, false);
                 double cfac = corrFactor_->getCorr(id);
@@ -1110,7 +1114,7 @@ bool CalibPlotProperties::selectPhi(bool debug) {
     // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
     for (unsigned int k = 0; k < t_HitEnergies->size(); ++k) {
       // Apply thresholds if necessary
-      bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[k] > threshold((*t_DetIds)[k], thrForm_));
+      bool okcell = (thrForm_ == 0) || ((*t_HitEnergies)[k] > (cThr_->threshold((*t_DetIds)[k])));
       if (okcell) {
         int iphi = ((*t_DetIds)[k]) & (0x3FF);
         int zside = ((*t_DetIds)[k] & 0x80000) ? (1) : (-1);
@@ -1258,7 +1262,7 @@ void CalibPlotProperties::correctEnergy(double &eHcal) {
       // The masks are defined in DataFormats/HcalDetId/interface/HcalDetId.h
       for (unsigned int idet = 0; idet < (*t_DetIds1).size(); idet++) {
         // Apply thresholds if necessary
-        bool okcell = (thrForm_ == 0) || ((*t_HitEnergies1)[idet] > threshold((*t_DetIds1)[idet], thrForm_));
+        bool okcell = (thrForm_ == 0) || ((*t_HitEnergies1)[idet] > (cThr_->threshold((*t_DetIds1)[idet])));
         if (okcell) {
           unsigned int id = truncateId((*t_DetIds1)[idet], truncateFlag_, false);
           double cfac = corrFactor_->getCorr(id);
@@ -1277,7 +1281,7 @@ void CalibPlotProperties::correctEnergy(double &eHcal) {
       }
       for (unsigned int idet = 0; idet < (*t_DetIds3).size(); idet++) {
         // Apply thresholds if necessary
-        bool okcell = (thrForm_ == 0) || ((*t_HitEnergies3)[idet] > threshold((*t_DetIds3)[idet], thrForm_));
+        bool okcell = (thrForm_ == 0) || ((*t_HitEnergies3)[idet] > (cThr_->threshold((*t_DetIds3)[idet])));
         if (okcell) {
           unsigned int id = truncateId((*t_DetIds3)[idet], truncateFlag_, false);
           double cfac = corrFactor_->getCorr(id);
diff --git a/Calibration/HcalCalibAlgos/macros/CalibTree.C b/Calibration/HcalCalibAlgos/macros/CalibTree.C
diff --git a/Calibration/HcalCalibAlgos/macros/PFCuts362975.txt b/Calibration/HcalCalibAlgos/macros/PFCuts362975.txt
