Merge pull request cms-sw#33999 from bsunanda/Run3-alca186

Run3-alca186 Update the calibration scripts for HCAL

Author: cmsbuild

Calibration/HcalCalibAlgos/macros/CalibCorr.C

@@ -178,7 +178,7 @@ double puFactor(int type, int ieta, double pmom, double eHcal, double ediff, boo
       if (debug)
         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
-    } else {  // 21pu
+    } else if (type == 6) {  // 21pu (old)
       const double CONST_COR_COEF[6] = {0.98913, 0.982008, 0.974011, 0.496234, 0.368110, 0.294053};
       const double LINEAR_COR_COEF[6] = {-0.0491388, -0.124058, -0.249718, -0.0667390, -0.0770766, -0.0580492};
       const double SQUARE_COR_COEF[6] = {0, 0, 0.0368657, 0.00656337, 0.00724508, 0.00568967};
@@ -195,6 +195,40 @@ double puFactor(int type, int ieta, double pmom, double eHcal, double ediff, boo
       if (debug)
         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
+    } else if (type == 99) {  // dlphin
+      const double CONST_COR_COEF[6] = {0.98312, 0.978532, 0.972211, 0.756004, 0.638075, 0.547192};
+      const double LINEAR_COR_COEF[6] = {-0.0472436, -0.186206, -0.247339, -0.166062, -0.159781, -0.118747};
+      const double SQUARE_COR_COEF[6] = {0, 0, 0.0356827, 0.0202461, 0.01785078, 0.0123003};
+      const int PU_IETA_1 = 7;
+      const int PU_IETA_2 = 16;
+      const int PU_IETA_3 = 25;
+      const int PU_IETA_4 = 26;
+      const int PU_IETA_5 = 27;
+      unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
+                       unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
+      double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
+      if (d2p > deltaCut)
+        fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
+      if (debug)
+        std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
+                  << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
+    } else {  // 21pu (June, 2021)
+      const double CONST_COR_COEF[6] = {0.989727, 0.981923, 0.97571, 0.562475, 0.467947, 0.411831};
+      const double LINEAR_COR_COEF[6] = {-0.0469558, -0.125805, -0.251383, -0.0668994, -0.0964236, -0.0947158};
+      const double SQUARE_COR_COEF[6] = {0, 0, 0.0399785, 0.00610104, 0.00952528, 0.0100645};
+      const int PU_IETA_1 = 7;
+      const int PU_IETA_2 = 16;
+      const int PU_IETA_3 = 25;
+      const int PU_IETA_4 = 26;
+      const int PU_IETA_5 = 27;
+      unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
+                       unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
+      double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
+      if (d2p > deltaCut)
+        fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
+      if (debug)
+        std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
+                  << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
     }
   }
   if (fac < 0 || fac > 1)
@@ -307,9 +341,10 @@ private:
   const int useScale_;
   const double scale_;
   const bool etaMax_, debug_;
+  static const int depMax_ = 10;
   bool corrE_;
   int etamp_, etamn_;
-  double cfacmp_, cfacmn_;
+  double cfacmp_[depMax_], cfacmn_[depMax_];
   std::map<std::pair<int, int>, double> cfactors_;
 };
 
@@ -357,7 +392,11 @@ private:
 };
 
 CalibCorrFactor::CalibCorrFactor(const char* infile, int useScale, double scale, bool etamax, bool debug)
-    : useScale_(useScale), scale_(scale), etaMax_(etamax), debug_(debug), etamp_(0), etamn_(0), cfacmp_(1), cfacmn_(1) {
+    : useScale_(useScale), scale_(scale), etaMax_(etamax), debug_(debug), etamp_(0), etamn_(0) {
+  for (int i = 0; i < depMax_; ++i) {
+    cfacmp_[i] = 1.0;
+    cfacmn_[i] = 1.0;
+  }
   if (std::string(infile) != "") {
     corrE_ = readCorrFactor(infile);
     std::cout << "Reads " << cfactors_.size() << " correction factors from " << infile << " with flag " << corrE_
@@ -382,9 +421,9 @@ double CalibCorrFactor::getCorr(unsigned int id) {
       cfac = itr->second;
     } else if (etaMax_) {
       if (zside > 0 && ieta > etamp_)
-        cfac = cfacmp_;
+        cfac = (depth < depMax_) ? cfacmp_[depth] : cfacmp_[depMax_ - 1];
       if (zside < 0 && ieta > -etamn_)
-        cfac = cfacmn_;
+        cfac = (depth < depMax_) ? cfacmn_[depth] : cfacmn_[depMax_ - 1];
     }
   } else if (useScale_ != 0) {
     int subdet, zside, ieta, iphi, depth;
@@ -417,20 +456,21 @@ bool CalibCorrFactor::readCorrFactor(const char* fname) {
           float corrf = std::atof(items[3].c_str());
           double scale = getFactor(std::abs(ieta));
           cfactors_[std::pair<int, int>(ieta, depth)] = scale * corrf;
-          if (ieta > etamp_ && depth == 1) {
+          if (ieta > etamp_ && depth == 1)
             etamp_ = ieta;
-            cfacmp_ = scale * corrf;
-          }
-          if (ieta < etamn_ && depth == 1) {
+          if (ieta == etamp_ && depth < depMax_)
+            cfacmp_[depth] = scale * corrf;
+          if (ieta < etamn_ && depth == 1)
             etamn_ = ieta;
-            cfacmn_ = scale * corrf;
-          }
+          if (ieta == etamn_ && depth < depMax_)
+            cfacmn_[depth] = scale * corrf;
         }
       }
       fInput.close();
       std::cout << "Reads total of " << all << " and " << good << " good records"
-                << " Max eta (z>0) " << etamp_ << ":" << cfacmp_ << " eta (z<0) " << etamn_ << ":" << cfacmn_
-                << std::endl;
+                << " Max eta (z>0) " << etamp_ << " eta (z<0) " << etamn_ << std::endl;
+      for (int i = 0; i < depMax_; ++i)
+        std::cout << "[" << i << "] C+ " << cfacmp_[i] << " C- " << cfacmn_[i] << std::endl;
       if (good > 0)
         ok = true;
     }

Calibration/HcalCalibAlgos/macros/CalibFitPlots.C

@@ -36,8 +36,8 @@
 //
 //             For plotting on the same canvas plots with different
 //             prefixes residing in the same file with approrprate text
-//   PlotTwoHists(infile, prefix1, text1, prefix2, text2, text0, type, iname,
-//                lumi, ener, drawStatBox, save);
+//  PlotTwoHists(infile, prefix1, text1, prefix2, text2, text0, type, iname,
+//               lumi, ener, drawStatBox, save);
 //      Defaults: type=0; iname=2; lumi=0; ener=13; drawStatBox=true;
 //                save=0;
 //      Note prefixN, textN have the same meaning as prefix and text for set N
@@ -46,8 +46,8 @@
 //               =1 plots MPV of response vs RBX #
 //
 //             For plotting stored histograms from CalibTree
-//   PlotFiveHists(infile, text0, prefix0, type, iname, drawStatBox, normalize,
-//                 save, prefix1, text1, prefix2, text2, prefix3, text3,
+//  PlotFiveHists(infile, text0, prefix0, type, iname, drawStatBox, normalize,
+//                save, prefix1, text1, prefix2, text2, prefix3, text3,
 //                 prefix4, text4, prefix5, text5);
 //      Defaults: type=0; iname=0; drawStatBox=true; normalize=false;
 //                save=0; prefixN=""; textN=""; (for N > 0)
@@ -1222,6 +1222,7 @@ void PlotHistEta(const char* infile,
     int eta = (k > numb2) ? (k - numb2) : (k - numb2 - 1);
     sprintf(name, "%s%s%d%d", prefix.c_str(), name0.c_str(), iene, k);
     TH1D* hist1 = (TH1D*)file->FindObjectAny(name);
+    std::cout << name << " at " << hist1 << std::endl;
     if (hist1 != nullptr) {
       TH1D* hist = (TH1D*)(hist1->Clone());
       double ymin(0.90);
@@ -1928,8 +1929,8 @@ void PlotHistCorrFactor(char* infile,
     gStyle->SetOptStat(0);
     gStyle->SetOptFit(0);
   }
-  int colors[6] = {1, 6, 4, 7, 2, 9};
-  int mtype[6] = {20, 21, 22, 23, 24, 33};
+  int colors[7] = {1, 6, 4, 7, 2, 9, 3};
+  int mtype[7] = {20, 21, 22, 23, 24, 33, 25};
   int nbin = etamax - etamin + 1;
   std::vector<TH1D*> hists;
   std::vector<int> entries;

Calibration/HcalCalibAlgos/macros/CalibMonitor.C

@@ -935,7 +935,7 @@ void CalibMonitor::Loop() {
       break;
     nb = fChain->GetEntry(jentry);
     nbytes += nb;
-    if (jentry % 100000 == 0)
+    if (jentry % 1000000 == 0)
       std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
     double pmom = (useGen_ && (t_gentrackP > 0)) ? t_gentrackP : t_p;
     bool p4060 = ((pmom >= 40.0) && (pmom <= 60.0));
@@ -2623,7 +2623,7 @@ void CalibPlotProperties::Loop() {
       break;
     nb = fChain->GetEntry(jentry);
     nbytes += nb;
-    if (jentry % 100000 == 0)
+    if (jentry % 1000000 == 0)
       std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
     bool select = (std::find(entries_.begin(), entries_.end(), jentry) == entries_.end());
     if (!select) {

Calibration/HcalCalibAlgos/macros/CalibPlotProperties.C

@@ -778,7 +778,7 @@ void CalibPlotProperties::Loop(Long64_t nentries) {
       break;
     nb = fChain->GetEntry(jentry);
     nbytes += nb;
-    if (jentry % 100000 == 0)
+    if (jentry % 1000000 == 0)
       std::cout << "Entry " << jentry << " Run " << t_Run << " Event " << t_Event << std::endl;
     bool select = (std::find(entries_.begin(), entries_.end(), jentry) == entries_.end());
     if (!select) {