Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-05-10 02:20:18

0001 //////////////////////////////////////////////////////////////////////////////
0002 // Contains several utility functions used by CalibMonitor/CalibSort/CalibTree:
0003 //
0004 // void unpackDetId(detId, subdet, zside, ieta, iphi, depth)
0005 //      Unpacks rawId for Hadron calorimeter
0006 // unsigned int truncateId(detId, truncateFlag, debug)
0007 //      Defines depth of the DetId guided by the value of truncateFlag
0008 // double puFactor(type, ieta, pmom, eHcal, ediff, debug)
0009 //      Returns a multiplicative factor to the energy as PU correction
0010 //        int type = The type of data being considered
0011 //                   (1: Run 1 old; 2: Run1 new; 3: Run2 2016;
0012 //                    4: Run 2 2017; 5: Run2 2018; 6: Run3 Old;
0013 //                    7: Run 3 June 2021; 8: Run3 Mahi Jan2022;
0014 //                    9: Run3 M0 Jan2022l; 97: dlphin Try 3;
0015 //                    98: dlphin Try 2; 99: dlphin Try 1)
0016 // double puFactorRho(type, ieta, rho, double eHcal)
0017 //      Returns a multiplicative factor as PU correction evaluated from rho
0018 //        int type = The type of data being considered
0019 //                   (1: 2017 Data;  2: 2017 MC; 3: 2018 MC; 4: 2018AB;
0020 //                    5: 2018BC; 6: 2016 MC)
0021 // double puweight(vtx)
0022 //      Return PU weight for QCD PU sample
0023 // bool fillChain(chain, inputFileList)
0024 //      Prepares a Tchain by chaining several ROOT files specified
0025 // std::vector<std::string> splitString (fLine)
0026 //      Splits a string into several items which are separated by blank in i/p
0027 // CalibCorrFactor(infile, useScale, scale, etamax, debug)
0028 //      A class which reads a file with correction factors and provides
0029 //        bool   doCorr() : flag saying if correction is available
0030 //        double getCorr(id): correction factor for a given cell
0031 // CalibCorr(infile, flag, debug)
0032 //      A second class which reads a file and provide corrections through
0033 //        float getCorr(run, id): Run and ID dependent correction
0034 //        double getCorr(entry): Entry # (in the file) dependent correction
0035 //        bool absent(entry) : if correction factor absent
0036 //        bool present(entry): or present (relevant for ML motivated)
0037 // CalibSelectRBX(rbxFile, debug)
0038 //      A class for selecting a given set of Read Out Box's and provides
0039 //        bool isItRBX(detId): if it/they is in the chosen RBXs
0040 //        bool isItRBX(ieta, iphi): if it is in the chosen RBXs
0041 // CalibDuplicate(infile, flag, debug)
0042 //      A class for either rejecting duplicate entries or giving depth
0043 //        dependent weight. flag is 0 for keeping a list of duplicate
0044 //        emtries; 1 is to keep depth dependent weight for each ieta;
0045 //        2 is to keep a list of ieta, iphi for channels to be selected.
0046 //        bool isDuplicate(entry): if it is a duplicate entry
0047 //        double getWeight(ieta, depth): get the dependent weight
0048 //        bool select(int ieta, int iphi): channels to be selected
0049 // void CalibCorrTest(infile, flag)
0050 //      Tests a file which contains correction factors used by CalibCorr
0051 // void CalibCorrScale(infile, oufile, scale)
0052 //      Scales all contents of correction factors by "scale" from "infile"
0053 //      to "outfile"
0054 //////////////////////////////////////////////////////////////////////////////
0055 #ifndef CalibrationHcalCalibAlgosCalibCorr_h
0056 #define CalibrationHcalCalibAlgosCalibCorr_h
0057 
0058 #include <algorithm>
0059 #include <iomanip>
0060 #include <iostream>
0061 #include <fstream>
0062 #include <sstream>
0063 #include <map>
0064 #include <string>
0065 #include <vector>
0066 
0067 #include <TChain.h>
0068 #include <TROOT.h>
0069 
0070 void unpackDetId(unsigned int detId, int& subdet, int& zside, int& ieta, int& iphi, int& depth) {
0071   // The maskings are defined in DataFormats/DetId/interface/DetId.h
0072   //                      and in DataFormats/HcalDetId/interface/HcalDetId.h
0073   // The macro does not invoke the classes there and use them
0074   subdet = ((detId >> 25) & (0x7));
0075   if ((detId & 0x1000000) == 0) {
0076     ieta = ((detId >> 7) & 0x3F);
0077     zside = (detId & 0x2000) ? (1) : (-1);
0078     depth = ((detId >> 14) & 0x1F);
0079     iphi = (detId & 0x3F);
0080   } else {
0081     ieta = ((detId >> 10) & 0x1FF);
0082     zside = (detId & 0x80000) ? (1) : (-1);
0083     depth = ((detId >> 20) & 0xF);
0084     iphi = (detId & 0x3FF);
0085   }
0086 }
0087 
0088 unsigned int packDetId(int subdet, int ieta, int iphi, int depth) {
0089   // The maskings are defined in DataFormats/DetId/interface/DetId.h
0090   //                      and in DataFormats/HcalDetId/interface/HcalDetId.h
0091   // The macro does not invoke the classes there and use them
0092   const int det(4);
0093   unsigned int id = (((det & 0xF) << 28) | ((subdet & 0x7) << 25));
0094   id |= ((0x1000000) | ((depth & 0xF) << 20) | ((ieta > 0) ? (0x80000 | (ieta << 10)) : ((-ieta) << 10)) |
0095          (iphi & 0x3FF));
0096   return id;
0097 }
0098 
0099 unsigned int matchDetId(unsigned int detId) {
0100   if ((detId & 0x1000000) == 0) {
0101     int subdet = ((detId >> 25) & (0x7));
0102     int ieta = ((detId >> 7) & 0x3F);
0103     int zside = (detId & 0x2000) ? (1) : (-1);
0104     int depth = ((detId >> 14) & 0x1F);
0105     int iphi = (detId & 0x3F);
0106     return packDetId(subdet, (ieta * zside), iphi, depth);
0107   } else {
0108     return detId;
0109   }
0110 }
0111 
0112 unsigned int truncateId(unsigned int detId, int truncateFlag, bool debug = false) {
0113   //Truncate depth information of DetId's
0114   unsigned int id(detId);
0115   if (debug) {
0116     std::cout << "Truncate 1 " << std::hex << detId << " " << id << std::dec << " Flag " << truncateFlag << std::endl;
0117   }
0118   int truncate0 = ((truncateFlag / 1) % 10);
0119   int truncate1 = ((truncateFlag / 10) % 10);
0120   int subdet, depth, zside, ieta, iphi;
0121   unpackDetId(detId, subdet, zside, ieta, iphi, depth);
0122   if (truncate1 == 1)
0123     zside = 1;
0124   if (truncate0 == 1) {
0125     //Ignore depth index of ieta values of 15 and 16 of HB
0126     if ((subdet == 1) && (ieta > 14))
0127       depth = 1;
0128   } else if (truncate0 == 2) {
0129     //Ignore depth index of all ieta values
0130     depth = 1;
0131   } else if (truncate0 == 3) {
0132     //Ignore depth index for depth > 1 in HE
0133     if ((subdet == 2) && (depth > 1))
0134       depth = 2;
0135     else
0136       depth = 1;
0137   } else if (truncate0 == 4) {
0138     //Ignore depth index for depth > 1 in HB
0139     if ((subdet == 1) && (depth > 1))
0140       depth = 2;
0141     else
0142       depth = 1;
0143   } else if (truncate0 == 5) {
0144     //Ignore depth index for depth > 1 in HB and HE
0145     if (depth > 1)
0146       depth = 2;
0147   }
0148   id = (subdet << 25) | (0x1000000) | ((depth & 0xF) << 20) | ((zside > 0) ? (0x80000 | (ieta << 10)) : (ieta << 10));
0149   if (debug) {
0150     std::cout << "Truncate 2: " << subdet << " " << zside * ieta << " " << depth << " " << std::hex << id << " input "
0151               << detId << std::dec << std::endl;
0152   }
0153   return id;
0154 }
0155 
0156 unsigned int repackId(const std::string& det, int eta, int depth) {
0157   int subdet = (det == "HE") ? 2 : 1;
0158   int zside = (eta >= 0) ? 1 : -1;
0159   int ieta = (eta >= 0) ? eta : -eta;
0160   unsigned int id =
0161       (subdet << 25) | (0x1000000) | ((depth & 0xF) << 20) | ((zside > 0) ? (0x80000 | (ieta << 10)) : (ieta << 10));
0162   return id;
0163 }
0164 
0165 unsigned int repackId(int eta, int depth) {
0166   int zside = (eta >= 0) ? 1 : -1;
0167   int ieta = (eta >= 0) ? eta : -eta;
0168   int subdet = ((ieta > 16) || ((ieta == 16) && (depth > 3))) ? 2 : 1;
0169   unsigned int id =
0170       (subdet << 25) | (0x1000000) | ((depth & 0xF) << 20) | ((zside > 0) ? (0x80000 | (ieta << 10)) : (ieta << 10));
0171   return id;
0172 }
0173 
0174 unsigned int repackId(int subdet, int ieta, int iphi, int depth) {
0175   unsigned int id = ((subdet & 0x7) << 25);
0176   id |= ((0x1000000) | ((depth & 0xF) << 20) | ((ieta > 0) ? (0x80000 | (ieta << 10)) : ((-ieta) << 10)) |
0177          (iphi & 0x3FF));
0178   return id;
0179 }
0180 
0181 bool ifHB(int ieta, int depth) { return ((std::abs(ieta) < 16) || ((std::abs(ieta) == 16) && (depth != 4))); }
0182 
0183 int truncateDepth(int ieta, int depth, int truncateFlag) {
0184   int truncate0 = ((truncateFlag / 1) % 10);
0185   int d(depth);
0186   if (truncate0 == 5) {
0187     d = (depth == 1) ? 1 : 2;
0188   } else if (truncate0 == 4) {
0189     d = ifHB(ieta, depth) ? ((depth == 1) ? 1 : 2) : depth;
0190   } else if (truncate0 == 3) {
0191     d = (!ifHB(ieta, depth)) ? ((depth == 1) ? 1 : 2) : depth;
0192   } else if (truncate0 == 2) {
0193     d = 1;
0194   } else if (truncate0 == 1) {
0195     d = ((std::abs(ieta) == 15) || (std::abs(ieta) == 16)) ? 1 : depth;
0196   }
0197   return d;
0198 }
0199 
0200 double threshold(int subdet, int depth, int form) {
0201   double cutHE[7] = {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2};
0202   double cutHB[3][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}};
0203   double thr(0);
0204   if (form > 0) {
0205     if (subdet == 2)
0206       thr = cutHE[depth - 1];
0207     else
0208       thr = cutHB[form - 1][depth - 1];
0209   }
0210   return thr;
0211 }
0212 
0213 double threshold(unsigned int detId, int form) {
0214   int subdet = ((detId >> 25) & (0x7));
0215   int depth = ((detId & 0x1000000) == 0) ? ((detId >> 14) & 0x1F) : ((detId >> 20) & 0xF);
0216   return threshold(subdet, depth, form);
0217 }
0218 
0219 double puFactor(int type, int ieta, double pmom, double eHcal, double ediff, bool debug = false) {
0220   double fac(1.0);
0221   if (debug)
0222     std::cout << "Input Type " << type << " ieta " << ieta << " pmon " << pmom << " E " << eHcal << ":" << ediff;
0223   if (type <= 2) {
0224     double frac = (type == 1) ? 0.02 : 0.03;
0225     if (pmom > 0 && ediff > frac * pmom) {
0226       double a1(0), a2(0);
0227       if (type == 1) {
0228         a1 = -0.35;
0229         a2 = -0.65;
0230         if (std::abs(ieta) == 25) {
0231           a2 = -0.30;
0232         } else if (std::abs(ieta) > 25) {
0233           a1 = -0.45;
0234           a2 = -0.10;
0235         }
0236       } else {
0237         a1 = -0.39;
0238         a2 = -0.59;
0239         if (std::abs(ieta) >= 25) {
0240           a1 = -0.283;
0241           a2 = -0.272;
0242         } else if (std::abs(ieta) > 22) {
0243           a1 = -0.238;
0244           a2 = -0.241;
0245         }
0246       }
0247       fac = (1.0 + a1 * (eHcal / pmom) * (ediff / pmom) * (1 + a2 * (ediff / pmom)));
0248       if (debug)
0249         std::cout << " coeff " << a1 << ":" << a2 << " Fac " << fac;
0250     }
0251   } else {
0252     int jeta = std::abs(ieta);
0253     double d2p = (ediff / pmom);
0254     const double DELTA_CUT = 0.03;
0255     if (type == 3) {  // 16pu
0256       const double CONST_COR_COEF[4] = {0.971, 1.008, 0.985, 1.086};
0257       const double LINEAR_COR_COEF[4] = {0, -0.359, -0.251, -0.535};
0258       const double SQUARE_COR_COEF[4] = {0, 0, 0.048, 0.143};
0259       const int PU_IETA_1 = 9;
0260       const int PU_IETA_2 = 16;
0261       const int PU_IETA_3 = 25;
0262       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3));
0263       if (d2p > DELTA_CUT)
0264         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0265       if (debug)
0266         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0267                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0268     } else if (type == 4) {  // 17pu
0269       const double CONST_COR_COEF[4] = {0.974, 1.023, 0.989, 1.077};
0270       const double LINEAR_COR_COEF[4] = {0, -0.524, -0.268, -0.584};
0271       const double SQUARE_COR_COEF[4] = {0, 0, 0.053, 0.170};
0272       const int PU_IETA_1 = 9;
0273       const int PU_IETA_2 = 18;
0274       const int PU_IETA_3 = 25;
0275       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3));
0276       if (d2p > DELTA_CUT)
0277         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0278       if (debug)
0279         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0280                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0281     } else if (type == 5) {  // 18pu
0282       const double CONST_COR_COEF[4] = {0.973, 0.998, 0.992, 0.965};
0283       const double LINEAR_COR_COEF[4] = {0, -0.318, -0.261, -0.406};
0284       const double SQUARE_COR_COEF[4] = {0, 0, 0.047, 0.089};
0285       const int PU_IETA_1 = 7;
0286       const int PU_IETA_2 = 16;
0287       const int PU_IETA_3 = 25;
0288       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3));
0289       if (d2p > DELTA_CUT)
0290         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0291       if (debug)
0292         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0293                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0294     } else if (type == 6) {  // 21pu (old)
0295       const double CONST_COR_COEF[6] = {0.98913, 0.982008, 0.974011, 0.496234, 0.368110, 0.294053};
0296       const double LINEAR_COR_COEF[6] = {-0.0491388, -0.124058, -0.249718, -0.0667390, -0.0770766, -0.0580492};
0297       const double SQUARE_COR_COEF[6] = {0, 0, 0.0368657, 0.00656337, 0.00724508, 0.00568967};
0298       const int PU_IETA_1 = 7;
0299       const int PU_IETA_2 = 16;
0300       const int PU_IETA_3 = 25;
0301       const int PU_IETA_4 = 26;
0302       const int PU_IETA_5 = 27;
0303       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0304                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0305       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0306       if (d2p > deltaCut)
0307         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0308       if (debug)
0309         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0310                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0311     } else if (type == 97) {  // dlphin Try 3
0312       const double CONST_COR_COEF[6] = {0.987617, 0.983421, 0.938622, 0.806662, 0.738354, 0.574195};
0313       const double LINEAR_COR_COEF[6] = {-0.07018610, -0.2494880, -0.1997290, -0.1769320, -0.2427950, -0.1230480};
0314       const double SQUARE_COR_COEF[6] = {0, 0, 0.0263541, 0.0257008, 0.0426584, 0.0200361};
0315       const int PU_IETA_1 = 7;
0316       const int PU_IETA_2 = 16;
0317       const int PU_IETA_3 = 25;
0318       const int PU_IETA_4 = 26;
0319       const int PU_IETA_5 = 27;
0320       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0321                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0322       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0323       if (d2p > deltaCut)
0324         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0325       if (debug)
0326         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0327                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0328     } else if (type == 98) {  // dlphin Try 2
0329       const double CONST_COR_COEF[6] = {0.987665, 0.983468, 0.938628, 0.807241, 0.739132, 0.529059};
0330       const double LINEAR_COR_COEF[6] = {-0.0708906, -0.249995, -0.199683, -0.177692, -0.243436, -0.0668783};
0331       const double SQUARE_COR_COEF[6] = {0, 0, 0.0263163, 0.0260158, 0.0426864, 0.00398778};
0332       const int PU_IETA_1 = 7;
0333       const int PU_IETA_2 = 16;
0334       const int PU_IETA_3 = 25;
0335       const int PU_IETA_4 = 26;
0336       const int PU_IETA_5 = 27;
0337       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0338                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0339       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0340       if (d2p > deltaCut)
0341         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0342       if (debug)
0343         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0344                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0345     } else if (type == 99) {  // dlphin Try 1
0346       const double CONST_COR_COEF[6] = {0.98312, 0.978532, 0.972211, 0.756004, 0.638075, 0.547192};
0347       const double LINEAR_COR_COEF[6] = {-0.0472436, -0.186206, -0.247339, -0.166062, -0.159781, -0.118747};
0348       const double SQUARE_COR_COEF[6] = {0, 0, 0.0356827, 0.0202461, 0.01785078, 0.0123003};
0349       const int PU_IETA_1 = 7;
0350       const int PU_IETA_2 = 16;
0351       const int PU_IETA_3 = 25;
0352       const int PU_IETA_4 = 26;
0353       const int PU_IETA_5 = 27;
0354       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0355                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0356       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0357       if (d2p > deltaCut)
0358         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0359       if (debug)
0360         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0361                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0362     } else if (type == 7) {  // 21pu (June, 2021)
0363       const double CONST_COR_COEF[6] = {0.989727, 0.981923, 0.97571, 0.562475, 0.467947, 0.411831};
0364       const double LINEAR_COR_COEF[6] = {-0.0469558, -0.125805, -0.251383, -0.0668994, -0.0964236, -0.0947158};
0365       const double SQUARE_COR_COEF[6] = {0, 0, 0.0399785, 0.00610104, 0.00952528, 0.0100645};
0366       const int PU_IETA_1 = 7;
0367       const int PU_IETA_2 = 16;
0368       const int PU_IETA_3 = 25;
0369       const int PU_IETA_4 = 26;
0370       const int PU_IETA_5 = 27;
0371       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0372                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0373       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0374       if (d2p > deltaCut)
0375         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0376       if (debug)
0377         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0378                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0379     } else if (type == 9) {  // M0 22pu (Jan, 2022)
0380       const double CONST_COR_COEF[6] = {0.980941, 0.973156, 0.970749, 0.726582, 0.532628, 0.473727};
0381       const double LINEAR_COR_COEF[6] = {-0.0770642, -0.178295, -0.241338, -0.122956, -0.122346, -0.112574};
0382       const double SQUARE_COR_COEF[6] = {0, 0, 0.0401732, 0.00989908, 0.0108291, 0.0100508};
0383       const int PU_IETA_1 = 7;
0384       const int PU_IETA_2 = 16;
0385       const int PU_IETA_3 = 25;
0386       const int PU_IETA_4 = 26;
0387       const int PU_IETA_5 = 27;
0388       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0389                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0390       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0391       if (d2p > deltaCut)
0392         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0393       if (debug)
0394         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0395                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0396     } else if (type == 10) {  // MAHI (Jan, 2023)
0397       const double CONST_COR_COEF[6] = {0.987967, 0.983376, 0.954840, 0.676950, 0.513111, 0.430349};
0398       const double LINEAR_COR_COEF[6] = {-0.0399269, -0.101755, -0.156848, -0.0969012 - 0.107831, -0.0911755};
0399       const double SQUARE_COR_COEF[6] = {0, 0, 0.0133473, 0.00727513, 0.00863409, 0.00727055};
0400       const int PU_IETA_1 = 7;
0401       const int PU_IETA_2 = 16;
0402       const int PU_IETA_3 = 25;
0403       const int PU_IETA_4 = 26;
0404       const int PU_IETA_5 = 27;
0405       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0406                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0407       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0408       if (d2p > deltaCut)
0409         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0410       if (debug)
0411         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0412                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0413     } else {  // Mahi 22pu (Jan, 2022)
0414       const double CONST_COR_COEF[6] = {0.995902, 0.991240, 0.981019, 0.788052, 0.597956, 0.538731};
0415       const double LINEAR_COR_COEF[6] = {-0.0540563, -0.104361, -0.215936, -0.147801, -0.160845, -0.154359};
0416       const double SQUARE_COR_COEF[6] = {0, 0, 0.0365911, 0.0161266, 0.0180053, 0.0184295};
0417       const int PU_IETA_1 = 7;
0418       const int PU_IETA_2 = 16;
0419       const int PU_IETA_3 = 25;
0420       const int PU_IETA_4 = 26;
0421       const int PU_IETA_5 = 27;
0422       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0423                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0424       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0425       if (d2p > deltaCut)
0426         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0427       if (debug)
0428         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0429                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0430     }
0431   }
0432   if (fac < 0 || fac > 1)
0433     fac = 0;
0434   if (debug)
0435     std::cout << " Final factor " << fac << std::endl;
0436   return fac;
0437 }
0438 
0439 double puFactorRho(int type, int ieta, double rho, double eHcal) {
0440   // type = 1: 2017 Data;  2: 2017 MC; 3: 2018 MC; 4: 2018AB; 5: 2018BC
0441   //        6: 2016 MC;
0442   double par[36] = {0.0205395,  -43.0914,   2.67115,    0.239674,    -0.0228009, 0.000476963, 0.137566,  -32.8386,
0443                     3.25886,    0.0863636,  -0.0165639, 0.000413894, 0.206168,   -145.828,    10.3191,   0.531418,
0444                     -0.0578416, 0.00118905, 0.175356,   -175.543,    14.3414,    0.294718,    -0.049836, 0.00106228,
0445                     0.134314,   -175.809,   13.5307,    0.395943,    -0.0539062, 0.00111573,  0.145342,  -98.1904,
0446                     8.14001,    0.205526,   -0.0327818, 0.000726059};
0447   double energy(eHcal);
0448   if (type >= 1 && type <= 6) {
0449     int eta = std::abs(ieta);
0450     int it = 6 * (type - 1);
0451     double ea =
0452         (eta < 20)
0453             ? par[it]
0454             : ((((par[it + 5] * eta + par[it + 4]) * eta + par[it + 3]) * eta + par[it + 2]) * eta + par[it + 1]);
0455     energy -= (rho * ea);
0456   }
0457   return energy;
0458 }
0459 
0460 double puweight(double vtx) {  ///////for QCD PU sample
0461   double a(1.0);
0462   if (vtx < 11)
0463     a = 0.120593;
0464   else if (vtx < 21)
0465     a = 0.58804;
0466   else if (vtx < 31)
0467     a = 1.16306;
0468   else if (vtx < 41)
0469     a = 1.45892;
0470   else if (vtx < 51)
0471     a = 1.35528;
0472   else if (vtx < 61)
0473     a = 1.72032;
0474   else if (vtx < 71)
0475     a = 3.34812;
0476   else if (vtx < 81)
0477     a = 9.70097;
0478   else if (vtx < 91)
0479     a = 9.24839;
0480   else if (vtx < 101)
0481     a = 23.0816;
0482   return a;
0483 }
0484 
0485 bool fillChain(TChain* chain, const char* inputFileList) {
0486   std::string fname(inputFileList);
0487   if (fname.substr(fname.size() - 5, 5) == ".root") {
0488     chain->Add(fname.c_str());
0489   } else {
0490     std::ifstream infile(inputFileList);
0491     if (!infile.is_open()) {
0492       std::cout << "** ERROR: Can't open '" << inputFileList << "' for input" << std::endl;
0493       return false;
0494     }
0495     while (1) {
0496       infile >> fname;
0497       if (!infile.good())
0498         break;
0499       chain->Add(fname.c_str());
0500     }
0501     infile.close();
0502   }
0503   std::cout << "No. of Entries in this tree : " << chain->GetEntries() << std::endl;
0504   return true;
0505 }
0506 
0507 std::vector<std::string> splitString(const std::string& fLine) {
0508   std::vector<std::string> result;
0509   int start = 0;
0510   bool empty = true;
0511   for (unsigned i = 0; i <= fLine.size(); i++) {
0512     if (fLine[i] == ' ' || i == fLine.size()) {
0513       if (!empty) {
0514         std::string item(fLine, start, i - start);
0515         result.push_back(item);
0516         empty = true;
0517       }
0518       start = i + 1;
0519     } else {
0520       if (empty)
0521         empty = false;
0522     }
0523   }
0524   return result;
0525 }
0526 
0527 class CalibCorrFactor {
0528 public:
0529   CalibCorrFactor(const char* infile, int useScale, double scale, bool etamax, bool marina, bool debug);
0530   ~CalibCorrFactor() {}
0531 
0532   bool doCorr() const { return (corrE_ || (useScale_ != 0)); }
0533   double getCorr(unsigned int id);
0534 
0535 private:
0536   bool readCorrFactor(const char* fName, bool marina);
0537   double getFactor(const int& ieta);
0538 
0539   const int useScale_;
0540   const double scale_;
0541   const bool etaMax_, debug_;
0542   static const int depMax_ = 10;
0543   bool corrE_;
0544   int etamp_, etamn_;
0545   double cfacmp_[depMax_], cfacmn_[depMax_];
0546   std::map<std::pair<int, int>, double> cfactors_;
0547 };
0548 
0549 class CalibCorr {
0550 public:
0551   CalibCorr(const char* infile, int flag, bool debug);
0552   ~CalibCorr() {}
0553 
0554   float getCorr(int run, unsigned int id);
0555   double getCorr(const Long64_t& entry);
0556   double getTrueCorr(const Long64_t& entry);
0557   double getPhiCorr(unsigned int id);
0558   bool absent(const Long64_t& entry);
0559   bool absent() { return (good_ == 0); }
0560   bool present(const Long64_t& entry);
0561 
0562 private:
0563   unsigned int readCorrRun(const char* infile);
0564   unsigned int readCorrDepth(const char* infile);
0565   unsigned int readCorrResp(const char* infile);
0566   unsigned int readCorrPU(const char* infile);
0567   unsigned int readCorrPhi(const char* infile);
0568   unsigned int getDetIdHE(int ieta, int iphi, int depth);
0569   unsigned int getDetId(int subdet, int ieta, int iphi, int depth);
0570   unsigned int correctDetId(const unsigned int& detId);
0571 
0572   static const unsigned int nmax_ = 10;
0573   int flag_;
0574   bool debug_;
0575   unsigned int good_;
0576   std::map<unsigned int, float> corrFac_[nmax_], corrFacDepth_, corrFacResp_;
0577   std::map<Long64_t, double> cfactors_;
0578   std::vector<int> runlow_;
0579   std::map<unsigned int, double> corrPhiSym_;
0580 };
0581 
0582 class CalibSelectRBX {
0583 public:
0584   CalibSelectRBX(const char* rbxFile, bool debug = false);
0585   ~CalibSelectRBX() {}
0586 
0587   bool isItRBX(const unsigned int);
0588   bool isItRBX(const std::vector<unsigned int>*);
0589   bool isItRBX(const int, const int);
0590 
0591 private:
0592   bool debug_;
0593   std::vector<int> zsphis_;
0594 };
0595 
0596 class CalibDuplicate {
0597 public:
0598   CalibDuplicate(const char* infile, int flag, bool debug);
0599   ~CalibDuplicate() {}
0600 
0601   bool isDuplicate(long entry);
0602   double getWeight(const unsigned int);
0603   bool doCorr() { return ((flag_ == 1) && ok_); }
0604   bool select(int ieta, int iphi);
0605 
0606 private:
0607   int flag_;
0608   double debug_, ok_;
0609   std::vector<Long64_t> entries_;
0610   std::map<int, std::vector<double> > weights_;
0611   std::vector<std::pair<int, int> > etaphi_;
0612 };
0613 
0614 CalibCorrFactor::CalibCorrFactor(const char* infile, int useScale, double scale, bool etamax, bool marina, bool debug)
0615     : useScale_(useScale), scale_(scale), etaMax_(etamax), debug_(debug), etamp_(0), etamn_(0) {
0616   for (int i = 0; i < depMax_; ++i) {
0617     cfacmp_[i] = 1.0;
0618     cfacmn_[i] = 1.0;
0619   }
0620   if (std::string(infile) != "") {
0621     corrE_ = readCorrFactor(infile, marina);
0622     std::cout << "Reads " << cfactors_.size() << " correction factors from " << infile << " with flag " << corrE_
0623               << std::endl
0624               << "Flag for scale " << useScale_ << " with scale " << scale_ << "; flag for etaMax " << etaMax_
0625               << " and flag for Format " << marina << std::endl;
0626   } else {
0627     corrE_ = false;
0628     std::cout << "No correction factors provided; Flag for scale " << useScale_ << " with scale " << scale_
0629               << "; flag for etaMax " << etaMax_ << " and flag for Format " << marina << std::endl;
0630   }
0631 }
0632 
0633 double CalibCorrFactor::getCorr(unsigned int id) {
0634   double cfac(1.0);
0635   if (corrE_) {
0636     if (cfactors_.size() > 0) {
0637       int subdet, zside, ieta, iphi, depth;
0638       unpackDetId(id, subdet, zside, ieta, iphi, depth);
0639       std::map<std::pair<int, int>, double>::const_iterator itr =
0640           cfactors_.find(std::pair<int, int>(zside * ieta, depth));
0641       if (itr != cfactors_.end()) {
0642         cfac = itr->second;
0643       } else if (etaMax_) {
0644         if (zside > 0 && ieta > etamp_)
0645           cfac = (depth < depMax_) ? cfacmp_[depth] : cfacmp_[depMax_ - 1];
0646         if (zside < 0 && ieta > -etamn_)
0647           cfac = (depth < depMax_) ? cfacmn_[depth] : cfacmn_[depMax_ - 1];
0648       }
0649     }
0650   } else if (useScale_ != 0) {
0651     int subdet, zside, ieta, iphi, depth;
0652     unpackDetId(id, subdet, zside, ieta, iphi, depth);
0653     cfac = getFactor(ieta);
0654   }
0655   return cfac;
0656 }
0657 
0658 bool CalibCorrFactor::readCorrFactor(const char* fname, bool marina) {
0659   bool ok(false);
0660   if (std::string(fname) != "") {
0661     std::ifstream fInput(fname);
0662     if (!fInput.good()) {
0663       std::cout << "Cannot open file " << fname << std::endl;
0664     } else {
0665       char buffer[1024];
0666       unsigned int all(0), good(0);
0667       while (fInput.getline(buffer, 1024)) {
0668         ++all;
0669         if (buffer[0] == '#')
0670           continue;  //ignore comment
0671         std::vector<std::string> items = splitString(std::string(buffer));
0672         if (items.size() != 5) {
0673           std::cout << "Ignore  line: " << buffer << std::endl;
0674         } else {
0675           ++good;
0676           int ieta = (marina) ? std::atoi(items[0].c_str()) : std::atoi(items[1].c_str());
0677           int depth = (marina) ? std::atoi(items[1].c_str()) : std::atoi(items[2].c_str());
0678           float corrf = std::atof(items[3].c_str());
0679           double scale = getFactor(std::abs(ieta));
0680           cfactors_[std::pair<int, int>(ieta, depth)] = scale * corrf;
0681           if (ieta > etamp_ && depth == 1)
0682             etamp_ = ieta;
0683           if (ieta == etamp_ && depth < depMax_)
0684             cfacmp_[depth] = scale * corrf;
0685           if (ieta < etamn_ && depth == 1)
0686             etamn_ = ieta;
0687           if (ieta == etamn_ && depth < depMax_)
0688             cfacmn_[depth] = scale * corrf;
0689         }
0690       }
0691       fInput.close();
0692       std::cout << "Reads total of " << all << " and " << good << " good records"
0693                 << " Max eta (z>0) " << etamp_ << " eta (z<0) " << etamn_ << std::endl;
0694       for (int i = 0; i < depMax_; ++i)
0695         std::cout << "[" << i << "] C+ " << cfacmp_[i] << " C- " << cfacmn_[i] << std::endl;
0696       if (good > 0)
0697         ok = true;
0698     }
0699   }
0700   return ok;
0701 }
0702 
0703 double CalibCorrFactor::getFactor(const int& ieta) {
0704   double scale(1.0);
0705   if (ieta < 16) {
0706     if ((useScale_ == 1) || (useScale_ == 3))
0707       scale = scale_;
0708   } else {
0709     if ((useScale_ == 2) || (useScale_ == 3))
0710       scale = scale_;
0711   }
0712   return scale;
0713 }
0714 
0715 CalibCorr::CalibCorr(const char* infile, int flag, bool debug) : flag_(flag), debug_(debug) {
0716   std::cout << "CalibCorr is created with flag " << flag << ":" << flag_ << " for i/p file " << infile << std::endl;
0717   if (flag == 1)
0718     good_ = readCorrDepth(infile);
0719   else if (flag == 2)
0720     good_ = readCorrResp(infile);
0721   else if (flag == 3)
0722     good_ = readCorrPU(infile);
0723   else if (flag == 4)
0724     good_ = readCorrPhi(infile);
0725   else
0726     good_ = readCorrRun(infile);
0727 }
0728 
0729 float CalibCorr::getCorr(int run, unsigned int id) {
0730   float cfac(1.0);
0731   if (good_ == 0)
0732     return cfac;
0733   unsigned idx = correctDetId(id);
0734   if (flag_ == 1) {
0735     std::map<unsigned int, float>::iterator itr = corrFacDepth_.find(idx);
0736     if (itr != corrFacDepth_.end())
0737       cfac = itr->second;
0738   } else if (flag_ == 2) {
0739     std::map<unsigned int, float>::iterator itr = corrFacResp_.find(idx);
0740     if (itr != corrFacResp_.end())
0741       cfac = itr->second;
0742   } else if (flag_ == 4) {
0743     std::map<unsigned int, double>::iterator itr = corrPhiSym_.find(idx);
0744     if (itr != corrPhiSym_.end())
0745       cfac = itr->second;
0746   } else {
0747     int ip(-1);
0748     for (unsigned int k = 0; k < runlow_.size(); ++k) {
0749       unsigned int i = runlow_.size() - k - 1;
0750       if (run >= runlow_[i]) {
0751         ip = (int)(i);
0752         break;
0753       }
0754     }
0755     if (debug_) {
0756       std::cout << "Run " << run << " Perdiod " << ip << std::endl;
0757     }
0758     if (ip >= 0) {
0759       std::map<unsigned int, float>::iterator itr = corrFac_[ip].find(idx);
0760       if (itr != corrFac_[ip].end())
0761         cfac = itr->second;
0762     }
0763   }
0764   if (debug_) {
0765     int subdet, zside, ieta, iphi, depth;
0766     unpackDetId(idx, subdet, zside, ieta, iphi, depth);
0767     std::cout << "ID " << std::hex << id << std::dec << " (Sub " << subdet << " eta " << zside * ieta << " phi " << iphi
0768               << " depth " << depth << ")  Factor " << cfac << std::endl;
0769   }
0770   return cfac;
0771 }
0772 
0773 double CalibCorr::getCorr(const Long64_t& entry) {
0774   if (good_ == 0)
0775     return 1.0;
0776   double cfac(0.0);
0777   std::map<Long64_t, double>::iterator itr = cfactors_.find(entry);
0778   if (itr != cfactors_.end())
0779     cfac = std::min(itr->second, 10.0);
0780   return cfac;
0781 }
0782 
0783 double CalibCorr::getTrueCorr(const Long64_t& entry) {
0784   if (good_ == 0)
0785     return 1.0;
0786   double cfac(0.0);
0787   std::map<Long64_t, double>::iterator itr = cfactors_.find(entry);
0788   if (itr != cfactors_.end())
0789     cfac = itr->second;
0790   return cfac;
0791 }
0792 
0793 double CalibCorr::getPhiCorr(unsigned int idx) {
0794   double cfac(1.0);
0795   if (good_ == 0)
0796     return cfac;
0797   std::map<unsigned int, double>::iterator itr = corrPhiSym_.find(idx);
0798   if (itr != corrPhiSym_.end())
0799     cfac = itr->second;
0800   if (debug_) {
0801     int subdet, zside, ieta, iphi, depth;
0802     unpackDetId(idx, subdet, zside, ieta, iphi, depth);
0803     std::cout << "ID " << std::hex << idx << std::dec << " (Sub " << subdet << " eta " << zside * ieta << " phi "
0804               << iphi << " depth " << depth << ")  Factor " << cfac << std::endl;
0805   }
0806   return cfac;
0807 }
0808 
0809 bool CalibCorr::absent(const Long64_t& entry) { return (cfactors_.find(entry) == cfactors_.end()); }
0810 
0811 bool CalibCorr::present(const Long64_t& entry) { return (cfactors_.find(entry) != cfactors_.end()); }
0812 
0813 unsigned int CalibCorr::readCorrRun(const char* infile) {
0814   std::cout << "Enters readCorrRun for " << infile << std::endl;
0815   std::ifstream fInput(infile);
0816   unsigned int all(0), good(0), ncorr(0);
0817   if (!fInput.good()) {
0818     std::cout << "Cannot open file " << infile << std::endl;
0819   } else {
0820     char buffer[1024];
0821     while (fInput.getline(buffer, 1024)) {
0822       ++all;
0823       std::string bufferString(buffer);
0824       if (bufferString.substr(0, 5) == "#IOVs") {
0825         std::vector<std::string> items = splitString(bufferString.substr(6));
0826         ncorr = items.size() - 1;
0827         for (unsigned int n = 0; n < ncorr; ++n) {
0828           int run = std::atoi(items[n].c_str());
0829           runlow_.push_back(run);
0830         }
0831         std::cout << ncorr << ":" << runlow_.size() << " Run ranges" << std::endl;
0832         for (unsigned int n = 0; n < runlow_.size(); ++n)
0833           std::cout << " [" << n << "] " << runlow_[n];
0834         std::cout << std::endl;
0835       } else if (buffer[0] == '#') {
0836         continue;  //ignore other comments
0837       } else {
0838         std::vector<std::string> items = splitString(bufferString);
0839         if (items.size() != ncorr + 3) {
0840           std::cout << "Ignore  line: " << buffer << std::endl;
0841         } else {
0842           ++good;
0843           int ieta = std::atoi(items[0].c_str());
0844           int iphi = std::atoi(items[1].c_str());
0845           int depth = std::atoi(items[2].c_str());
0846           unsigned int id = getDetIdHE(ieta, iphi, depth);
0847           for (unsigned int n = 0; n < ncorr; ++n) {
0848             float corrf = std::atof(items[n + 3].c_str());
0849             if (n < nmax_)
0850               corrFac_[n][id] = corrf;
0851           }
0852           if (debug_) {
0853             std::cout << "ID " << std::hex << id << std::dec << ":" << id << " (eta " << ieta << " phi " << iphi
0854                       << " depth " << depth << ")";
0855             for (unsigned int n = 0; n < ncorr; ++n)
0856               std::cout << " " << corrFac_[n][id];
0857             std::cout << std::endl;
0858           }
0859         }
0860       }
0861     }
0862     fInput.close();
0863     std::cout << "Reads total of " << all << " and " << good << " good records of run dependent corrections from "
0864               << infile << std::endl;
0865   }
0866   return good;
0867 }
0868 
0869 unsigned int CalibCorr::readCorrDepth(const char* infile) {
0870   std::cout << "Enters readCorrDepth for " << infile << std::endl;
0871   unsigned int all(0), good(0);
0872   std::ifstream fInput(infile);
0873   if (!fInput.good()) {
0874     std::cout << "Cannot open file " << infile << std::endl;
0875   } else {
0876     char buffer[1024];
0877     while (fInput.getline(buffer, 1024)) {
0878       ++all;
0879       std::string bufferString(buffer);
0880       if (bufferString.substr(0, 5) == "depth") {
0881         continue;  //ignore other comments
0882       } else {
0883         std::vector<std::string> items = splitString(bufferString);
0884         if (items.size() != 3) {
0885           std::cout << "Ignore  line: " << buffer << " Size " << items.size();
0886           for (unsigned int k = 0; k < items.size(); ++k)
0887             std::cout << " [" << k << "] : " << items[k];
0888           std::cout << std::endl;
0889         } else {
0890           ++good;
0891           int ieta = std::atoi(items[1].c_str());
0892           int depth = std::atoi(items[0].c_str());
0893           float corrf = std::atof(items[2].c_str());
0894           int nphi = (std::abs(ieta) > 20) ? 36 : 72;
0895           for (int i = 1; i <= nphi; ++i) {
0896             int iphi = (nphi > 36) ? i : (2 * i - 1);
0897             unsigned int id = getDetIdHE(ieta, iphi, depth);
0898             corrFacDepth_[id] = corrf;
0899             if (debug_) {
0900               std::cout << "ID " << std::hex << id << std::dec << ":" << id << " (eta " << ieta << " phi " << iphi
0901                         << " depth " << depth << ") " << corrFacDepth_[id] << std::endl;
0902             }
0903           }
0904         }
0905       }
0906     }
0907     fInput.close();
0908     std::cout << "Reads total of " << all << " and " << good << " good records of depth dependent factors from "
0909               << infile << std::endl;
0910   }
0911   return good;
0912 }
0913 
0914 unsigned int CalibCorr::readCorrResp(const char* infile) {
0915   std::cout << "Enters readCorrResp for " << infile << std::endl;
0916   unsigned int all(0), good(0), other(0);
0917   std::ifstream fInput(infile);
0918   if (!fInput.good()) {
0919     std::cout << "Cannot open file " << infile << std::endl;
0920   } else {
0921     char buffer[1024];
0922     while (fInput.getline(buffer, 1024)) {
0923       ++all;
0924       std::string bufferString(buffer);
0925       if (bufferString.substr(0, 1) == "#") {
0926         continue;  //ignore other comments
0927       } else {
0928         std::vector<std::string> items = splitString(bufferString);
0929         if (items.size() < 5) {
0930           std::cout << "Ignore  line: " << buffer << " Size " << items.size();
0931           for (unsigned int k = 0; k < items.size(); ++k)
0932             std::cout << " [" << k << "] : " << items[k];
0933           std::cout << std::endl;
0934         } else if (items[3] == "HB" || items[3] == "HE") {
0935           ++good;
0936           int ieta = std::atoi(items[0].c_str());
0937           int iphi = std::atoi(items[1].c_str());
0938           int depth = std::atoi(items[2].c_str());
0939           int subdet = (items[3] == "HE") ? 2 : 1;
0940           float corrf = std::atof(items[4].c_str());
0941           unsigned int id = getDetId(subdet, ieta, iphi, depth);
0942           corrFacResp_[id] = corrf;
0943           if (debug_) {
0944             std::cout << "ID " << std::hex << id << std::dec << ":" << id << " (subdet " << items[3] << ":" << subdet
0945                       << " eta " << ieta << " phi " << iphi << " depth " << depth << ") " << corrFacResp_[id]
0946                       << std::endl;
0947           }
0948         } else {
0949           ++other;
0950         }
0951       }
0952     }
0953     fInput.close();
0954     std::cout << "Reads total of " << all << " and " << good << " good and " << other
0955               << " detector records of depth dependent factors from " << infile << std::endl;
0956   }
0957   return good;
0958 }
0959 
0960 unsigned int CalibCorr::readCorrPU(const char* infile) {
0961   if (std::string(infile) != "") {
0962     std::ifstream fInput(infile);
0963     if (!fInput.good()) {
0964       std::cout << "Cannot open file " << infile << std::endl;
0965     } else {
0966       double val1, val2;
0967       cfactors_.clear();
0968       while (1) {
0969         fInput >> val1 >> val2;
0970         if (!fInput.good())
0971           break;
0972         Long64_t entry = (Long64_t)(val1);
0973         cfactors_[entry] = val2;
0974       }
0975       fInput.close();
0976     }
0977   }
0978   std::cout << "Reads " << cfactors_.size() << " PU correction factors from " << infile << std::endl;
0979   return cfactors_.size();
0980 }
0981 
0982 unsigned int CalibCorr::readCorrPhi(const char* infile) {
0983   std::cout << "Enters readCorrPhi for " << infile << std::endl;
0984   unsigned int all(0), good(0);
0985   std::ifstream fInput(infile);
0986   if (!fInput.good()) {
0987     std::cout << "Cannot open file " << infile << std::endl;
0988   } else {
0989     char buffer[1024];
0990     while (fInput.getline(buffer, 1024)) {
0991       ++all;
0992       std::string bufferString(buffer);
0993       if (bufferString.substr(0, 1) == "#") {
0994         continue;  //ignore other comments
0995       } else {
0996         std::vector<std::string> items = splitString(bufferString);
0997         if (items.size() < 5) {
0998           std::cout << "Ignore  line: " << buffer << " Size " << items.size();
0999           for (unsigned int k = 0; k < items.size(); ++k)
1000             std::cout << " [" << k << "] : " << items[k];
1001           std::cout << std::endl;
1002         } else {
1003           ++good;
1004           int subdet = std::atoi(items[0].c_str());
1005           int ieta = std::atoi(items[1].c_str());
1006           int iphi = std::atoi(items[2].c_str());
1007           int depth = std::atoi(items[3].c_str());
1008           double corrf = std::atof(items[4].c_str());
1009           unsigned int id = packDetId(subdet, ieta, iphi, depth);
1010           corrPhiSym_[id] = corrf;
1011           if (debug_)
1012             std::cout << "ID " << std::hex << id << std::dec << ":" << id << " (subdet " << subdet << " eta " << ieta
1013                       << " phi " << iphi << " depth " << depth << ") " << corrPhiSym_[id] << std::endl;
1014         }
1015       }
1016     }
1017     fInput.close();
1018     std::cout << "Reads total of " << all << " and " << good << " good records of phi-symmetry factors from " << infile
1019               << std::endl;
1020   }
1021   return good;
1022 }
1023 
1024 unsigned int CalibCorr::getDetIdHE(int ieta, int iphi, int depth) { return getDetId(2, ieta, iphi, depth); }
1025 
1026 unsigned int CalibCorr::getDetId(int subdet, int ieta, int iphi, int depth) {
1027   // All numbers used here are described as masks/offsets in
1028   // DataFormats/HcalDetId/interface/HcalDetId.h
1029   unsigned int id_ = ((4 << 28) | ((subdet & 0x7) << 25));
1030   id_ |= ((0x1000000) | ((depth & 0xF) << 20) | ((ieta > 0) ? (0x80000 | (ieta << 10)) : ((-ieta) << 10)) |
1031           (iphi & 0x3FF));
1032   return id_;
1033 }
1034 
1035 unsigned int CalibCorr::correctDetId(const unsigned int& detId) {
1036   int subdet, ieta, zside, depth, iphi;
1037   unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1038   if (subdet == 0) {
1039     if (ieta > 16)
1040       subdet = 2;
1041     else if (ieta == 16 && depth > 2)
1042       subdet = 2;
1043     else
1044       subdet = 1;
1045   }
1046   unsigned int id = getDetId(subdet, ieta * zside, iphi, depth);
1047   if ((id != detId) && debug_) {
1048     std::cout << "Correct Id " << std::hex << detId << " to " << id << std::dec << "(Sub " << subdet << " eta "
1049               << ieta * zside << " phi " << iphi << " depth " << depth << ")" << std::endl;
1050   }
1051   return id;
1052 }
1053 
1054 CalibSelectRBX::CalibSelectRBX(const char* rbxFile, bool debug) : debug_(debug) {
1055   std::cout << "Enters CalibSelectRBX for " << rbxFile << std::endl;
1056   unsigned int all(0), good(0);
1057   std::vector<int> rbxs;
1058   std::ifstream fInput(rbxFile);
1059   if (!fInput.good()) {
1060     std::cout << "Cannot open file " << rbxFile << std::endl;
1061   } else {
1062     char buffer[1024];
1063     while (fInput.getline(buffer, 1024)) {
1064       ++all;
1065       std::string bufferString(buffer);
1066       if (bufferString.substr(0, 1) == "#") {
1067         continue;  //ignore other comments
1068       } else {
1069         std::vector<std::string> items = splitString(bufferString);
1070         if (items.size() != 1) {
1071           std::cout << "Ignore  line: " << buffer << " Size " << items.size() << std::endl;
1072         } else {
1073           ++good;
1074           int rbx = std::atoi(items[0].c_str());
1075           rbxs.push_back(rbx);
1076           int zside = (rbx > 0) ? 1 : -1;
1077           int subdet = (std::abs(rbx) / 100) % 10;
1078           if (subdet != 1)
1079             subdet = 2;
1080           int iphis = std::abs(rbx) % 100;
1081           if (iphis > 0 && iphis <= 18) {
1082             for (int i = 0; i < 4; ++i) {
1083               int iphi = (iphis - 2) * 4 + 3 + i;
1084               if (iphi < 1)
1085                 iphi += 72;
1086               int zsphi = zside * (subdet * 100 + iphi);
1087               zsphis_.push_back(zsphi);
1088             }
1089           }
1090         }
1091       }
1092     }
1093     fInput.close();
1094   }
1095   std::cout << "Select a set of RBXs " << rbxs.size() << " by reading " << all << ":" << good << " records from "
1096             << rbxFile << " with " << zsphis_.size() << " iphi values:";
1097   for (unsigned int i = 0; i < zsphis_.size(); ++i)
1098     std::cout << " " << zsphis_[i];
1099   std::cout << std::endl;
1100 }
1101 
1102 bool CalibSelectRBX::isItRBX(const unsigned int detId) {
1103   bool ok(true);
1104   if (zsphis_.size() > 0) {
1105     int subdet, ieta, zside, depth, iphi;
1106     unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1107     int zsphi = zside * (subdet * 100 + iphi);
1108     ok = (std::find(zsphis_.begin(), zsphis_.end(), zsphi) != zsphis_.end());
1109 
1110     if (debug_) {
1111       std::cout << "isItRBX:subdet|zside|iphi " << subdet << ":" << zside << ":" << iphi << " OK " << ok << std::endl;
1112     }
1113   }
1114   return ok;
1115 }
1116 
1117 bool CalibSelectRBX::isItRBX(const std::vector<unsigned int>* detId) {
1118   bool ok(true);
1119   if (zsphis_.size() > 0) {
1120     ok = true;
1121     for (unsigned int i = 0; i < detId->size(); ++i) {
1122       int subdet, ieta, zside, depth, iphi;
1123       unpackDetId((*detId)[i], subdet, zside, ieta, iphi, depth);
1124       int zsphi = zside * (subdet * 100 + iphi);
1125       ok = (std::find(zsphis_.begin(), zsphis_.end(), zsphi) != zsphis_.end());
1126       if (debug_) {
1127         std::cout << "isItRBX: subdet|zside|iphi " << subdet << ":" << zside << ":" << iphi << " OK " << ok
1128                   << std::endl;
1129       }
1130       if (ok)
1131         break;
1132     }
1133   }
1134   if (debug_)
1135     std::cout << "isItRBX: size " << detId->size() << " OK " << ok << std::endl;
1136   return ok;
1137 }
1138 
1139 bool CalibSelectRBX::isItRBX(const int ieta, const int iphi) {
1140   bool ok(true);
1141   if (zsphis_.size() == 4) {
1142     int zside = (ieta > 0) ? 1 : -1;
1143     int subd1 = (std::abs(ieta) <= 16) ? 1 : 2;
1144     int subd2 = (std::abs(ieta) >= 16) ? 2 : 1;
1145     int zsphi1 = zside * (subd1 * 100 + iphi);
1146     int zsphi2 = zside * (subd2 * 100 + iphi);
1147     ok = ((std::find(zsphis_.begin(), zsphis_.end(), zsphi1) != zsphis_.end()) ||
1148           (std::find(zsphis_.begin(), zsphis_.end(), zsphi2) != zsphis_.end()));
1149   }
1150   if (debug_) {
1151     std::cout << "isItRBX: ieta " << ieta << " iphi " << iphi << " OK " << ok << std::endl;
1152   }
1153   return ok;
1154 }
1155 
1156 CalibDuplicate::CalibDuplicate(const char* fname, int flag, bool debug) : flag_(flag), debug_(debug), ok_(false) {
1157   if (flag_ == 0) {
1158     if (strcmp(fname, "") != 0) {
1159       std::ifstream infile(fname);
1160       if (!infile.is_open()) {
1161         std::cout << "Cannot open duplicate file " << fname << std::endl;
1162       } else {
1163         while (1) {
1164           Long64_t jentry;
1165           infile >> jentry;
1166           if (!infile.good())
1167             break;
1168           entries_.push_back(jentry);
1169         }
1170         infile.close();
1171         std::cout << "Reads a list of " << entries_.size() << " events from " << fname << std::endl;
1172         if (entries_.size() > 0)
1173           ok_ = true;
1174       }
1175     } else {
1176       std::cout << "No duplicate events in the input file" << std::endl;
1177     }
1178   } else if (flag_ == 1) {
1179     if (strcmp(fname, "") != 0) {
1180       std::ifstream infile(fname);
1181       if (!infile.is_open()) {
1182         std::cout << "Cannot open depth dependence file " << fname << std::endl;
1183       } else {
1184         unsigned int all(0), good(0);
1185         char buffer[1024];
1186         while (infile.getline(buffer, 1024)) {
1187           ++all;
1188           std::string bufferString(buffer);
1189           if (bufferString.substr(0, 1) == "#") {
1190             continue;  //ignore other comments
1191           } else {
1192             std::vector<std::string> items = splitString(bufferString);
1193             if (items.size() < 3) {
1194               std::cout << "Ignore  line: " << buffer << " Size " << items.size();
1195               for (unsigned int k = 0; k < items.size(); ++k)
1196                 std::cout << " [" << k << "] : " << items[k];
1197               std::cout << std::endl;
1198             } else {
1199               ++good;
1200               int ieta = std::atoi(items[0].c_str());
1201               std::vector<double> weights;
1202               for (unsigned int k = 1; k < items.size(); ++k) {
1203                 double corrf = std::atof(items[k].c_str());
1204                 weights.push_back(corrf);
1205               }
1206               weights_[ieta] = weights;
1207               if (debug_) {
1208                 std::cout << "Eta " << ieta << " with " << weights.size() << " depths having weights:";
1209                 for (unsigned int k = 0; k < weights.size(); ++k)
1210                   std::cout << " " << weights[k];
1211                 std::cout << std::endl;
1212               }
1213             }
1214           }
1215         }
1216         infile.close();
1217         std::cout << "Reads total of " << all << " and " << good << " good records of depth dependent factors from "
1218                   << fname << std::endl;
1219         if (good > 0)
1220           ok_ = true;
1221       }
1222     }
1223   } else {
1224     flag_ = 2;
1225     if (strcmp(fname, "") != 0) {
1226       std::ifstream infile(fname);
1227       if (!infile.is_open()) {
1228         std::cout << "Cannot open rejection file " << fname << std::endl;
1229       } else {
1230         unsigned int all(0), good(0);
1231         char buffer[1024];
1232         while (infile.getline(buffer, 1024)) {
1233           ++all;
1234           std::string bufferString(buffer);
1235           if (bufferString.substr(0, 1) == "#") {
1236             continue;  //ignore other comments
1237           } else {
1238             std::vector<std::string> items = splitString(bufferString);
1239             if (items.size() != 2) {
1240               std::cout << "Ignore  line: " << buffer << " Size " << items.size();
1241               for (unsigned int k = 0; k < items.size(); ++k)
1242                 std::cout << " [" << k << "] : " << items[k];
1243               std::cout << std::endl;
1244             } else {
1245               ++good;
1246               int ieta = std::atoi(items[0].c_str());
1247               int iphi = std::atoi(items[1].c_str());
1248               etaphi_.push_back(std::pair<int, int>(ieta, iphi));
1249               if (debug_)
1250                 std::cout << "Select channels with iEta " << ieta << " iPhi " << iphi << std::endl;
1251             }
1252           }
1253         }
1254         infile.close();
1255         std::cout << "Reads total of " << all << " and " << good << " good records of rejection candidates from "
1256                   << fname << std::endl;
1257         if (good > 0)
1258           ok_ = true;
1259       }
1260     }
1261   }
1262 }
1263 
1264 bool CalibDuplicate::isDuplicate(long entry) {
1265   if (ok_)
1266     return (std::find(entries_.begin(), entries_.end(), entry) != entries_.end());
1267   else
1268     return false;
1269 }
1270 
1271 double CalibDuplicate::getWeight(unsigned int detId) {
1272   double wt(1.0);
1273   if (ok_) {
1274     int subdet, ieta, zside, depth, iphi;
1275     unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1276     std::map<int, std::vector<double> >::const_iterator itr = weights_.find(ieta);
1277     --depth;
1278     if (depth < 0)
1279       std::cout << "Strange Depth value " << depth << " in " << std::hex << detId << std::dec << std::endl;
1280     if (itr != weights_.end()) {
1281       if (depth < static_cast<int>(itr->second.size()))
1282         wt = itr->second[depth];
1283     }
1284   }
1285   return wt;
1286 }
1287 
1288 bool CalibDuplicate::select(int ieta, int iphi) {
1289   bool flag(false);
1290   if ((ok_) && (flag_ == 2))
1291     flag = (std::find(etaphi_.begin(), etaphi_.end(), std::pair<int, int>(ieta, iphi)) != etaphi_.end());
1292   if (debug_)
1293     std::cout << "Input " << ieta << ":" << iphi << " Flags " << ok_ << ":" << flag_ << ":" << flag << std::endl;
1294   return flag;
1295 }
1296 
1297 void CalibCorrTest(const char* infile, int flag) {
1298   CalibCorr* c1 = new CalibCorr(infile, flag, true);
1299   for (int ieta = 1; ieta < 29; ++ieta) {
1300     int subdet = (ieta > 16) ? 2 : 1;
1301     int depth = (ieta > 16) ? 2 : 1;
1302     unsigned int id1 = ((4 << 28) | ((subdet & 0x7) << 25));
1303     id1 |= ((0x1000000) | ((depth & 0xF) << 20) | (ieta << 10) | 1);
1304     c1->getCorr(0, id1);
1305     id1 |= (0x80000);
1306     c1->getCorr(0, id1);
1307   }
1308 }
1309 
1310 unsigned int stringTest(const std::string& str) {
1311   std::cout << str << " has " << str.size() << " characters\n";
1312   return str.size();
1313 }
1314 
1315 void CalibCorrScale(const char* infile, const char* outfile, double scale) {
1316   std::ofstream myfile;
1317   myfile.open(outfile);
1318   if (!myfile.is_open()) {
1319     std::cout << "** ERROR: Can't open '" << outfile << std::endl;
1320   } else {
1321     if (std::string(infile) != "") {
1322       std::ifstream fInput(infile);
1323       if (!fInput.good()) {
1324         std::cout << "Cannot open file " << infile << std::endl;
1325       } else {
1326         char buffer[1024];
1327         unsigned int all(0), good(0), comment(0);
1328         while (fInput.getline(buffer, 1024)) {
1329           ++all;
1330           if (buffer[0] == '#') {
1331             myfile << buffer << std::endl;
1332             ++comment;
1333             continue;  //ignore comment
1334           }
1335           std::vector<std::string> items = splitString(std::string(buffer));
1336           if (items.size() != 5) {
1337             std::cout << "Ignore  line: " << buffer << std::endl;
1338           } else {
1339             ++good;
1340             int ieta = std::atoi(items[1].c_str());
1341             int depth = std::atoi(items[2].c_str());
1342             float corrf = scale * std::atof(items[3].c_str());
1343             float dcorr = scale * std::atof(items[4].c_str());
1344             myfile << std::setw(10) << items[0] << std::setw(10) << std::dec << ieta << std::setw(10) << depth
1345                    << std::setw(10) << corrf << " " << std::setw(10) << dcorr << std::endl;
1346           }
1347         }
1348         fInput.close();
1349         std::cout << "Reads total of " << all << ", " << comment << " and " << good << " good records from " << infile
1350                   << " and copied to " << outfile << std::endl;
1351       }
1352     }
1353     myfile.close();
1354   }
1355 }
1356 #endif