Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-06-29 22:57:36 

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 //        double getPhiCorr(unsigned int id): correction factor from phi
0038 //                                            symmery studies
0039 //        double getPhiCorrs(int run, unsigned int id): era-dependent correction
0040 //                                            factor from phi symmetry studies
0041 //      *infile* is the name of the input file
0042 //      *flag* can have a value between 0:5 specifying the type of correction
0043 //             5 for run-dependent correctons using results from several phi
0044 //             symmetry studies; 4 for using results from one phi-symmetry
0045 //             study; 3 for pileup correction using machine learning method;
0046 //             2 for overall response corrections; 1 for depth dependence
0047 //             corrections; 0 for raddam corrections)
0048 // CalibSelectRBX(rbxFile, debug)
0049 //      A class for selecting a given set of Read Out Box's and provides
0050 //        bool isItRBX(detId): if it/they is in the chosen RBXs
0051 //        bool isItRBX(ieta, iphi): if it is in the chosen RBXs
0052 // CalibDuplicate(infile, flag, debug)
0053 //      A class for either rejecting duplicate entries or giving depth
0054 //        dependent weight. flag is 0 for keeping a list of duplicate
0055 //        emtries; 1 is to keep depth dependent weight for each ieta;
0056 //        2 is to keep a list of ieta, iphi for channels to be selected;
0057 //        3 is to read the gain correction file
0058 //        bool isDuplicate(entry): if it is a duplicate entry
0059 //        double getWeight(ieta, depth): get the dependent weight
0060 //        bool select(int ieta, int iphi): channels to be selected
0061 //        double getCorr(int run, int ieta, int depth): correction factor
0062 //        to modify the rechit energy due to change in gain fctor
0063 // void CalibCorrTest(infile, flag)
0064 //      Tests a file which contains correction factors used by CalibCorr
0065 //      if *flag* is >= 0; for negative values it initializes CalibDuplicate
0066 //      with -flag and *debug* = tue
0067 // void CalibCorrScale(infile, oufile, scale)
0068 //      Scales all contents of correction factors by "scale" from "infile"
0069 //      to "outfile"
0070 // void CalibCorrScale2(infile, oufile, scaleB, scaleT, scaleE)
0071 //      Scales all contents of correction factors in the barrel, transition and
0072 //      endcap regions by "scaleB", "scaleT", "scaleE" from "infile" and writes
0073 //      them to "outfile"
0074 //////////////////////////////////////////////////////////////////////////////
0075 #ifndef CalibrationHcalCalibAlgosCalibCorr_h
0076 #define CalibrationHcalCalibAlgosCalibCorr_h
0077 
0078 #include <algorithm>
0079 #include <iomanip>
0080 #include <iostream>
0081 #include <fstream>
0082 #include <sstream>
0083 #include <map>
0084 #include <string>
0085 #include <vector>
0086 
0087 #include <TChain.h>
0088 #include <TROOT.h>
0089 
0090 void unpackDetId(unsigned int detId, int& subdet, int& zside, int& ieta, int& iphi, int& depth) {
0091   // The maskings are defined in DataFormats/DetId/interface/DetId.h
0092   //                      and in DataFormats/HcalDetId/interface/HcalDetId.h
0093   // The macro does not invoke the classes there and use them
0094   subdet = ((detId >> 25) & (0x7));
0095   if ((detId & 0x1000000) == 0) {
0096     ieta = ((detId >> 7) & 0x3F);
0097     zside = (detId & 0x2000) ? (1) : (-1);
0098     depth = ((detId >> 14) & 0x1F);
0099     iphi = (detId & 0x3F);
0100   } else {
0101     ieta = ((detId >> 10) & 0x1FF);
0102     zside = (detId & 0x80000) ? (1) : (-1);
0103     depth = ((detId >> 20) & 0xF);
0104     iphi = (detId & 0x3FF);
0105   }
0106 }
0107 
0108 unsigned int packDetId(int subdet, int ieta, int iphi, int depth) {
0109   // The maskings are defined in DataFormats/DetId/interface/DetId.h
0110   //                      and in DataFormats/HcalDetId/interface/HcalDetId.h
0111   // The macro does not invoke the classes there and use them
0112   const int det(4);
0113   unsigned int id = (((det & 0xF) << 28) | ((subdet & 0x7) << 25));
0114   id |= ((0x1000000) | ((depth & 0xF) << 20) | ((ieta > 0) ? (0x80000 | (ieta << 10)) : ((-ieta) << 10)) |
0115          (iphi & 0x3FF));
0116   return id;
0117 }
0118 
0119 unsigned int matchDetId(unsigned int detId) {
0120   if ((detId & 0x1000000) == 0) {
0121     int subdet = ((detId >> 25) & (0x7));
0122     int ieta = ((detId >> 7) & 0x3F);
0123     int zside = (detId & 0x2000) ? (1) : (-1);
0124     int depth = ((detId >> 14) & 0x1F);
0125     int iphi = (detId & 0x3F);
0126     return packDetId(subdet, (ieta * zside), iphi, depth);
0127   } else {
0128     return detId;
0129   }
0130 }
0131 
0132 unsigned int truncateId(unsigned int detId, int truncateFlag, bool debug = false) {
0133   //Truncate depth information of DetId's
0134   unsigned int id(detId);
0135   if (debug) {
0136     std::cout << "Truncate 1 " << std::hex << detId << " " << id << std::dec << " Flag " << truncateFlag << std::endl;
0137   }
0138   int truncate0 = ((truncateFlag / 1) % 10);
0139   int truncate1 = ((truncateFlag / 10) % 10);
0140   int subdet, depth, zside, ieta, iphi;
0141   unpackDetId(detId, subdet, zside, ieta, iphi, depth);
0142   if (truncate1 == 1)
0143     zside = 1;
0144   if (truncate0 == 1) {
0145     //Ignore depth index of ieta values of 15 and 16 of HB
0146     if ((subdet == 1) && (ieta > 14))
0147       depth = 1;
0148   } else if (truncate0 == 2) {
0149     //Ignore depth index of all ieta values
0150     depth = 1;
0151   } else if (truncate0 == 3) {
0152     //Ignore depth index for HE
0153     if (subdet == 2)
0154       depth = 1;
0155   } else if (truncate0 == 4) {
0156     //Ignore depth index for HB
0157     if (subdet == 1)
0158       depth = 1;
0159   } else if (truncate0 == 5) {
0160     //Ignore depth index for depth > 1 in HB and HE
0161     if (depth > 1)
0162       depth = 2;
0163   } else if (truncate0 == 6) {
0164     //Ignore depth index for depth > 2 in HB and HE; depths 1, 2 considered as 1
0165     if (depth <= 2)
0166       depth = 1;
0167     else
0168       depth = 2;
0169   } else if (truncate0 == 7) {
0170     //Ignore depth index for HB; depth 1, 2 considered as 1; deepth > 2 as 2
0171     //Ignore depth index for depth > 2 in HE; all depths considered as 1
0172     if (subdet == 1) {
0173       if (depth <= 2)
0174         depth = 1;
0175       else
0176         depth = 2;
0177     } else {
0178       depth = 1;
0179     }
0180   } else if (truncate0 == 8) {
0181     //Assign all depths > 4 as depth = 5
0182     if (depth >= 5)
0183       depth = 5;
0184   } else if (truncate0 == 9) {
0185     //Depths 1 and 2 considered as depth = 2; rest of the depths are kept as it is
0186     if (depth <= 2)
0187       depth = 2;
0188   }
0189   id = (subdet << 25) | (0x1000000) | ((depth & 0xF) << 20) | ((zside > 0) ? (0x80000 | (ieta << 10)) : (ieta << 10));
0190   if (debug) {
0191     std::cout << "Truncate 2: " << subdet << " " << zside * ieta << " " << depth << " " << std::hex << id << " input "
0192               << detId << std::dec << std::endl;
0193   }
0194   return id;
0195 }
0196 
0197 unsigned int repackId(const std::string& det, int eta, int depth) {
0198   int subdet = (det == "HE") ? 2 : 1;
0199   int zside = (eta >= 0) ? 1 : -1;
0200   int ieta = (eta >= 0) ? eta : -eta;
0201   unsigned int id =
0202       (subdet << 25) | (0x1000000) | ((depth & 0xF) << 20) | ((zside > 0) ? (0x80000 | (ieta << 10)) : (ieta << 10));
0203   return id;
0204 }
0205 
0206 unsigned int repackId(int eta, int depth) {
0207   int zside = (eta >= 0) ? 1 : -1;
0208   int ieta = (eta >= 0) ? eta : -eta;
0209   int subdet = ((ieta > 16) || ((ieta == 16) && (depth > 3))) ? 2 : 1;
0210   unsigned int id =
0211       (subdet << 25) | (0x1000000) | ((depth & 0xF) << 20) | ((zside > 0) ? (0x80000 | (ieta << 10)) : (ieta << 10));
0212   return id;
0213 }
0214 
0215 unsigned int repackId(int subdet, int ieta, int iphi, int depth) {
0216   unsigned int id = ((subdet & 0x7) << 25);
0217   id |= ((0x1000000) | ((depth & 0xF) << 20) | ((ieta > 0) ? (0x80000 | (ieta << 10)) : ((-ieta) << 10)) |
0218          (iphi & 0x3FF));
0219   return id;
0220 }
0221 
0222 bool ifHB(int ieta, int depth) { return ((std::abs(ieta) < 16) || ((std::abs(ieta) == 16) && (depth != 4))); }
0223 
0224 int truncateDepth(int ieta, int depth, int truncateFlag) {
0225   int truncate0 = ((truncateFlag / 1) % 10);
0226   int d(depth);
0227   if (truncate0 == 5) {
0228     d = (depth == 1) ? 1 : 2;
0229   } else if (truncate0 == 4) {
0230     d = ifHB(ieta, depth) ? ((depth == 1) ? 1 : 2) : depth;
0231   } else if (truncate0 == 3) {
0232     d = (!ifHB(ieta, depth)) ? ((depth == 1) ? 1 : 2) : depth;
0233   } else if (truncate0 == 2) {
0234     d = 1;
0235   } else if (truncate0 == 1) {
0236     d = ((std::abs(ieta) == 15) || (std::abs(ieta) == 16)) ? 1 : depth;
0237   }
0238   return d;
0239 }
0240 
0241 double puFactor(int type, int ieta, double pmom, double eHcal, double ediff, bool debug = false) {
0242   double fac(1.0);
0243   if (debug)
0244     std::cout << "Input Type " << type << " ieta " << ieta << " pmon " << pmom << " E " << eHcal << ":" << ediff;
0245   if (type <= 2) {
0246     double frac = (type == 1) ? 0.02 : 0.03;
0247     if (pmom > 0 && ediff > frac * pmom) {
0248       double a1(0), a2(0);
0249       if (type == 1) {
0250         a1 = -0.35;
0251         a2 = -0.65;
0252         if (std::abs(ieta) == 25) {
0253           a2 = -0.30;
0254         } else if (std::abs(ieta) > 25) {
0255           a1 = -0.45;
0256           a2 = -0.10;
0257         }
0258       } else {
0259         a1 = -0.39;
0260         a2 = -0.59;
0261         if (std::abs(ieta) >= 25) {
0262           a1 = -0.283;
0263           a2 = -0.272;
0264         } else if (std::abs(ieta) > 22) {
0265           a1 = -0.238;
0266           a2 = -0.241;
0267         }
0268       }
0269       fac = (1.0 + a1 * (eHcal / pmom) * (ediff / pmom) * (1 + a2 * (ediff / pmom)));
0270       if (debug)
0271         std::cout << " coeff " << a1 << ":" << a2 << " Fac " << fac;
0272     }
0273   } else {
0274     int jeta = std::abs(ieta);
0275     double d2p = (ediff / pmom);
0276     const double DELTA_CUT = 0.03;
0277     if (type == 3) {  // 16pu
0278       const double CONST_COR_COEF[4] = {0.971, 1.008, 0.985, 1.086};
0279       const double LINEAR_COR_COEF[4] = {0, -0.359, -0.251, -0.535};
0280       const double SQUARE_COR_COEF[4] = {0, 0, 0.048, 0.143};
0281       const int PU_IETA_1 = 9;
0282       const int PU_IETA_2 = 16;
0283       const int PU_IETA_3 = 25;
0284       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3));
0285       if (d2p > DELTA_CUT)
0286         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0287       if (debug)
0288         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0289                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0290     } else if (type == 4) {  // 17pu
0291       const double CONST_COR_COEF[4] = {0.974, 1.023, 0.989, 1.077};
0292       const double LINEAR_COR_COEF[4] = {0, -0.524, -0.268, -0.584};
0293       const double SQUARE_COR_COEF[4] = {0, 0, 0.053, 0.170};
0294       const int PU_IETA_1 = 9;
0295       const int PU_IETA_2 = 18;
0296       const int PU_IETA_3 = 25;
0297       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3));
0298       if (d2p > DELTA_CUT)
0299         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0300       if (debug)
0301         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0302                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0303     } else if (type == 5) {  // 18pu
0304       const double CONST_COR_COEF[4] = {0.973, 0.998, 0.992, 0.965};
0305       const double LINEAR_COR_COEF[4] = {0, -0.318, -0.261, -0.406};
0306       const double SQUARE_COR_COEF[4] = {0, 0, 0.047, 0.089};
0307       const int PU_IETA_1 = 7;
0308       const int PU_IETA_2 = 16;
0309       const int PU_IETA_3 = 25;
0310       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3));
0311       if (d2p > DELTA_CUT)
0312         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0313       if (debug)
0314         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0315                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0316     } else if (type == 6) {  // 21pu (old)
0317       const double CONST_COR_COEF[6] = {0.98913, 0.982008, 0.974011, 0.496234, 0.368110, 0.294053};
0318       const double LINEAR_COR_COEF[6] = {-0.0491388, -0.124058, -0.249718, -0.0667390, -0.0770766, -0.0580492};
0319       const double SQUARE_COR_COEF[6] = {0, 0, 0.0368657, 0.00656337, 0.00724508, 0.00568967};
0320       const int PU_IETA_1 = 7;
0321       const int PU_IETA_2 = 16;
0322       const int PU_IETA_3 = 25;
0323       const int PU_IETA_4 = 26;
0324       const int PU_IETA_5 = 27;
0325       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0326                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0327       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0328       if (d2p > deltaCut)
0329         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0330       if (debug)
0331         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0332                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0333     } else if (type == 97) {  // dlphin Try 3
0334       const double CONST_COR_COEF[6] = {0.987617, 0.983421, 0.938622, 0.806662, 0.738354, 0.574195};
0335       const double LINEAR_COR_COEF[6] = {-0.07018610, -0.2494880, -0.1997290, -0.1769320, -0.2427950, -0.1230480};
0336       const double SQUARE_COR_COEF[6] = {0, 0, 0.0263541, 0.0257008, 0.0426584, 0.0200361};
0337       const int PU_IETA_1 = 7;
0338       const int PU_IETA_2 = 16;
0339       const int PU_IETA_3 = 25;
0340       const int PU_IETA_4 = 26;
0341       const int PU_IETA_5 = 27;
0342       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0343                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0344       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0345       if (d2p > deltaCut)
0346         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0347       if (debug)
0348         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0349                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0350     } else if (type == 98) {  // dlphin Try 2
0351       const double CONST_COR_COEF[6] = {0.987665, 0.983468, 0.938628, 0.807241, 0.739132, 0.529059};
0352       const double LINEAR_COR_COEF[6] = {-0.0708906, -0.249995, -0.199683, -0.177692, -0.243436, -0.0668783};
0353       const double SQUARE_COR_COEF[6] = {0, 0, 0.0263163, 0.0260158, 0.0426864, 0.00398778};
0354       const int PU_IETA_1 = 7;
0355       const int PU_IETA_2 = 16;
0356       const int PU_IETA_3 = 25;
0357       const int PU_IETA_4 = 26;
0358       const int PU_IETA_5 = 27;
0359       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0360                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0361       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0362       if (d2p > deltaCut)
0363         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0364       if (debug)
0365         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0366                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0367     } else if (type == 99) {  // dlphin Try 1
0368       const double CONST_COR_COEF[6] = {0.98312, 0.978532, 0.972211, 0.756004, 0.638075, 0.547192};
0369       const double LINEAR_COR_COEF[6] = {-0.0472436, -0.186206, -0.247339, -0.166062, -0.159781, -0.118747};
0370       const double SQUARE_COR_COEF[6] = {0, 0, 0.0356827, 0.0202461, 0.01785078, 0.0123003};
0371       const int PU_IETA_1 = 7;
0372       const int PU_IETA_2 = 16;
0373       const int PU_IETA_3 = 25;
0374       const int PU_IETA_4 = 26;
0375       const int PU_IETA_5 = 27;
0376       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0377                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0378       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0379       if (d2p > deltaCut)
0380         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0381       if (debug)
0382         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0383                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0384     } else if (type == 7) {  // 21pu (June, 2021)
0385       const double CONST_COR_COEF[6] = {0.989727, 0.981923, 0.97571, 0.562475, 0.467947, 0.411831};
0386       const double LINEAR_COR_COEF[6] = {-0.0469558, -0.125805, -0.251383, -0.0668994, -0.0964236, -0.0947158};
0387       const double SQUARE_COR_COEF[6] = {0, 0, 0.0399785, 0.00610104, 0.00952528, 0.0100645};
0388       const int PU_IETA_1 = 7;
0389       const int PU_IETA_2 = 16;
0390       const int PU_IETA_3 = 25;
0391       const int PU_IETA_4 = 26;
0392       const int PU_IETA_5 = 27;
0393       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0394                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0395       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0396       if (d2p > deltaCut)
0397         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0398       if (debug)
0399         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0400                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0401     } else if (type == 9) {  // M0 22pu (Jan, 2022)
0402       const double CONST_COR_COEF[6] = {0.980941, 0.973156, 0.970749, 0.726582, 0.532628, 0.473727};
0403       const double LINEAR_COR_COEF[6] = {-0.0770642, -0.178295, -0.241338, -0.122956, -0.122346, -0.112574};
0404       const double SQUARE_COR_COEF[6] = {0, 0, 0.0401732, 0.00989908, 0.0108291, 0.0100508};
0405       const int PU_IETA_1 = 7;
0406       const int PU_IETA_2 = 16;
0407       const int PU_IETA_3 = 25;
0408       const int PU_IETA_4 = 26;
0409       const int PU_IETA_5 = 27;
0410       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0411                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0412       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0413       if (d2p > deltaCut)
0414         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0415       if (debug)
0416         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0417                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0418     } else if (type == 10) {  // MAHI (Jan, 2023)
0419       const double CONST_COR_COEF[6] = {0.987967, 0.983376, 0.954840, 0.676950, 0.513111, 0.430349};
0420       const double LINEAR_COR_COEF[6] = {-0.0399269, -0.101755, -0.156848, -0.0969012 - 0.107831, -0.0911755};
0421       const double SQUARE_COR_COEF[6] = {0, 0, 0.0133473, 0.00727513, 0.00863409, 0.00727055};
0422       const int PU_IETA_1 = 7;
0423       const int PU_IETA_2 = 16;
0424       const int PU_IETA_3 = 25;
0425       const int PU_IETA_4 = 26;
0426       const int PU_IETA_5 = 27;
0427       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0428                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0429       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0430       if (d2p > deltaCut)
0431         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0432       if (debug)
0433         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0434                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0435     } else {  // Mahi 22pu (Jan, 2022)
0436       const double CONST_COR_COEF[6] = {0.995902, 0.991240, 0.981019, 0.788052, 0.597956, 0.538731};
0437       const double LINEAR_COR_COEF[6] = {-0.0540563, -0.104361, -0.215936, -0.147801, -0.160845, -0.154359};
0438       const double SQUARE_COR_COEF[6] = {0, 0, 0.0365911, 0.0161266, 0.0180053, 0.0184295};
0439       const int PU_IETA_1 = 7;
0440       const int PU_IETA_2 = 16;
0441       const int PU_IETA_3 = 25;
0442       const int PU_IETA_4 = 26;
0443       const int PU_IETA_5 = 27;
0444       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0445                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0446       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0447       if (d2p > deltaCut)
0448         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0449       if (debug)
0450         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0451                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0452     }
0453   }
0454   if (fac < 0 || fac > 1)
0455     fac = 0;
0456   if (debug)
0457     std::cout << " Final factor " << fac << std::endl;
0458   return fac;
0459 }
0460 
0461 double puFactorRho(int type, int ieta, double rho, double eHcal) {
0462   // type = 1: 2017 Data;  2: 2017 MC; 3: 2018 MC; 4: 2018AB; 5: 2018BC
0463   //        6: 2016 MC;
0464   double par[36] = {0.0205395,  -43.0914,   2.67115,    0.239674,    -0.0228009, 0.000476963, 0.137566,  -32.8386,
0465                     3.25886,    0.0863636,  -0.0165639, 0.000413894, 0.206168,   -145.828,    10.3191,   0.531418,
0466                     -0.0578416, 0.00118905, 0.175356,   -175.543,    14.3414,    0.294718,    -0.049836, 0.00106228,
0467                     0.134314,   -175.809,   13.5307,    0.395943,    -0.0539062, 0.00111573,  0.145342,  -98.1904,
0468                     8.14001,    0.205526,   -0.0327818, 0.000726059};
0469   double energy(eHcal);
0470   if (type >= 1 && type <= 6) {
0471     int eta = std::abs(ieta);
0472     int it = 6 * (type - 1);
0473     double ea =
0474         (eta < 20)
0475             ? par[it]
0476             : ((((par[it + 5] * eta + par[it + 4]) * eta + par[it + 3]) * eta + par[it + 2]) * eta + par[it + 1]);
0477     energy -= (rho * ea);
0478   }
0479   return energy;
0480 }
0481 
0482 double puweight(double vtx) {  ///////for QCD PU sample
0483   double a(1.0);
0484   if (vtx < 11)
0485     a = 0.120593;
0486   else if (vtx < 21)
0487     a = 0.58804;
0488   else if (vtx < 31)
0489     a = 1.16306;
0490   else if (vtx < 41)
0491     a = 1.45892;
0492   else if (vtx < 51)
0493     a = 1.35528;
0494   else if (vtx < 61)
0495     a = 1.72032;
0496   else if (vtx < 71)
0497     a = 3.34812;
0498   else if (vtx < 81)
0499     a = 9.70097;
0500   else if (vtx < 91)
0501     a = 9.24839;
0502   else if (vtx < 101)
0503     a = 23.0816;
0504   return a;
0505 }
0506 
0507 bool fillChain(TChain* chain, const char* inputFileList) {
0508   std::string fname(inputFileList);
0509   if (fname.substr(fname.size() - 5, 5) == ".root") {
0510     chain->Add(fname.c_str());
0511   } else {
0512     std::ifstream infile(inputFileList);
0513     if (!infile.is_open()) {
0514       std::cout << "** ERROR: Can't open '" << inputFileList << "' for input" << std::endl;
0515       return false;
0516     }
0517     while (1) {
0518       infile >> fname;
0519       if (!infile.good())
0520         break;
0521       chain->Add(fname.c_str());
0522     }
0523     infile.close();
0524   }
0525   std::cout << "No. of Entries in this tree : " << chain->GetEntries() << std::endl;
0526   return true;
0527 }
0528 
0529 std::vector<std::string> splitString(const std::string& fLine) {
0530   std::vector<std::string> result;
0531   int start = 0;
0532   bool empty = true;
0533   for (unsigned i = 0; i <= fLine.size(); i++) {
0534     if (fLine[i] == ' ' || i == fLine.size()) {
0535       if (!empty) {
0536         std::string item(fLine, start, i - start);
0537         result.push_back(item);
0538         empty = true;
0539       }
0540       start = i + 1;
0541     } else {
0542       if (empty)
0543         empty = false;
0544     }
0545   }
0546   return result;
0547 }
0548 
0549 class CalibThreshold {
0550 public:
0551   CalibThreshold(int form);
0552   ~CalibThreshold() {}
0553 
0554   double threshold(unsigned int detId);
0555 
0556 private:
0557   double threshold(int subdet, int ieta, int depth);
0558   bool fileThreshold(const char* fname);
0559 
0560   int form_;
0561   std::map<std::pair<int, int>, double> thresh_;
0562   bool ok_;
0563 };
0564 
0565 CalibThreshold::CalibThreshold(int form) : form_(form) {
0566   if (form_ == 5)
0567     ok_ = fileThreshold("PFCuts362975.txt");
0568   else if ((form_ < 1) || (form_ > 5))
0569     ok_ = false;
0570   else
0571     ok_ = true;
0572   std::cout << "CalibThreshold initialized with flag " << ok_ << " for form " << form_ << std::endl;
0573 }
0574 
0575 bool CalibThreshold::fileThreshold(const char* fname) {
0576   bool ok(false);
0577   if (std::string(fname) != "") {
0578     std::ifstream fInput(fname);
0579     if (!fInput.good()) {
0580       std::cout << "Cannot open file " << fname << std::endl;
0581     } else {
0582       char buffer[1024];
0583       unsigned int all(0), good(0), bad(0);
0584       while (fInput.getline(buffer, 1024)) {
0585         ++all;
0586         if (buffer[0] == '#')
0587           continue;  //ignore comment
0588         std::vector<std::string> items = splitString(std::string(buffer));
0589         if (items.size() != 7) {
0590           ++bad;
0591           std::cout << "Ignore  line: " << buffer << std::endl;
0592         } else {
0593           ++good;
0594           int ieta = std::atoi(items[0].c_str());
0595           int depth = std::atoi(items[2].c_str());
0596           double thr = std::atof(items[4].c_str());
0597           thresh_[std::pair<int, int>(ieta, depth)] = thr;
0598         }
0599       }
0600       fInput.close();
0601       std::cout << "Reads " << all << " entries from " << fname << " with " << good << " Good and " << bad
0602                 << " bad records" << std::endl;
0603       if (good > 0)
0604         ok = true;
0605     }
0606   }
0607   return ok;
0608 }
0609 
0610 double CalibThreshold::threshold(int subdet, int ieta, int depth) {
0611   double cutHE[4][7] = {{0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2},
0612                         {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2},
0613                         {0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2},
0614                         {0.2, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3}};
0615   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}};
0616 
0617   double thr(0);
0618   if (ok_) {
0619     if ((form_ > 0) && (form_ <= 4)) {
0620       if (subdet == 2)
0621         thr = cutHE[form_ - 1][depth - 1];
0622       else
0623         thr = cutHB[form_ - 1][depth - 1];
0624     } else {
0625       std::map<std::pair<int, int>, double>::const_iterator itr = thresh_.find(std::pair<int, int>(ieta, depth));
0626       if (itr != thresh_.end())
0627         thr = itr->second;
0628     }
0629   }
0630   return thr;
0631 }
0632 
0633 double CalibThreshold::threshold(unsigned int detId) {
0634   int subdet, zside, ieta, iphi, depth;
0635   unpackDetId(detId, subdet, zside, ieta, iphi, depth);
0636   return threshold(subdet, (zside * ieta), depth);
0637 }
0638 
0639 class CalibCorrFactor {
0640 public:
0641   CalibCorrFactor(const char* infile, int useScale, double scale, bool etamax, bool marina, bool debug);
0642   ~CalibCorrFactor() {}
0643 
0644   bool doCorr() const { return (corrE_ || (useScale_ != 0)); }
0645   double getCorr(unsigned int id);
0646 
0647 private:
0648   bool readCorrFactor(const char* fName, bool marina);
0649   double getFactor(const int& ieta);
0650 
0651   const int useScale_;
0652   const double scale_;
0653   const bool etaMax_, debug_;
0654   static const int depMax_ = 10;
0655   bool corrE_;
0656   int etamp_, etamn_;
0657   double cfacmp_[depMax_], cfacmn_[depMax_];
0658   std::map<std::pair<int, int>, double> cfactors_;
0659 };
0660 
0661 class CalibCorr {
0662 public:
0663   CalibCorr(const char* infile, int flag, bool debug);
0664   ~CalibCorr() {}
0665 
0666   float getCorr(int run, unsigned int id);
0667   double getCorr(const Long64_t& entry);
0668   double getTrueCorr(const Long64_t& entry);
0669   double getPhiCorr(unsigned int id);
0670   double getPhiCorrs(int run, unsigned int id);
0671   bool absent(const Long64_t& entry);
0672   bool absent() { return (good_ == 0); }
0673   bool present(const Long64_t& entry);
0674 
0675 private:
0676   unsigned int runId(const int& run);
0677   unsigned int readCorrRun(const char* infile);
0678   unsigned int readCorrDepth(const char* infile);
0679   unsigned int readCorrResp(const char* infile);
0680   unsigned int readCorrPU(const char* infile);
0681   unsigned int readCorrPhi(const char* infile);
0682   unsigned int readCorrPhis(const char* infile);
0683   unsigned int getDetIdHE(int ieta, int iphi, int depth);
0684   unsigned int getDetId(int subdet, int ieta, int iphi, int depth);
0685   unsigned int correctDetId(const unsigned int& detId);
0686 
0687   static const unsigned int nmax_ = 10;
0688   int flag_;
0689   bool debug_;
0690   unsigned int good_;
0691   std::map<unsigned int, float> corrFac_[nmax_], corrFacDepth_, corrFacResp_;
0692   std::map<Long64_t, double> cfactors_;
0693   std::vector<int> runlow_;
0694   std::map<unsigned int, double> corrPhiSym_;
0695   std::vector<std::pair<int, int> > runs_, ids_;
0696   std::map<std::pair<unsigned int, unsigned int>, double> corrPhiSyms_;
0697 };
0698 
0699 class CalibSelectRBX {
0700 public:
0701   CalibSelectRBX(const char* rbxFile, bool debug = false);
0702   ~CalibSelectRBX() {}
0703 
0704   bool isItRBX(const unsigned int);
0705   bool isItRBX(const std::vector<unsigned int>*);
0706   bool isItRBX(const int, const int);
0707 
0708 private:
0709   bool debug_;
0710   std::vector<int> zsphis_;
0711 };
0712 
0713 class CalibDuplicate {
0714 public:
0715   CalibDuplicate(const char* infile, int flag, bool debug);
0716   ~CalibDuplicate() {}
0717 
0718   bool isDuplicate(long entry);
0719   double getWeight(const unsigned int);
0720   bool doCorr() { return (((flag_ == 1) || (flag_ == 3)) && ok_); }
0721   bool doCorr(int flag) { return ((flag_ == flag) && ok_); }
0722   bool select(int ieta, int iphi);
0723   double getCorr(int run, int ieta, int depth);
0724 
0725 private:
0726   bool readCorrFactor(const char* fName);
0727   unsigned int runId(const int& run);
0728 
0729   int flag_, nocc_;
0730   double debug_, ok_;
0731   std::vector<Long64_t> entries_;
0732   std::map<int, std::vector<double> > weights_;
0733   std::vector<std::pair<int, int> > etaphi_;
0734   std::vector<std::pair<int, int> > runs_, ids_;
0735   std::map<std::pair<int, int>, std::vector<double> > corrs_;
0736 };
0737 
0738 CalibCorrFactor::CalibCorrFactor(const char* infile, int useScale, double scale, bool etamax, bool marina, bool debug)
0739     : useScale_(useScale), scale_(scale), etaMax_(etamax), debug_(debug), etamp_(0), etamn_(0) {
0740   for (int i = 0; i < depMax_; ++i) {
0741     cfacmp_[i] = 1.0;
0742     cfacmn_[i] = 1.0;
0743   }
0744   if (std::string(infile) != "") {
0745     corrE_ = readCorrFactor(infile, marina);
0746     std::cout << "Reads " << cfactors_.size() << " correction factors from " << infile << " with flag " << corrE_
0747               << std::endl
0748               << "Flag for scale " << useScale_ << " with scale " << scale_ << "; flag for etaMax " << etaMax_
0749               << " and flag for Format " << marina << std::endl;
0750   } else {
0751     corrE_ = false;
0752     std::cout << "No correction factors provided; Flag for scale " << useScale_ << " with scale " << scale_
0753               << "; flag for etaMax " << etaMax_ << " and flag for Format " << marina << std::endl;
0754   }
0755 }
0756 
0757 double CalibCorrFactor::getCorr(unsigned int id) {
0758   double cfac(1.0);
0759   if (corrE_) {
0760     if (cfactors_.size() > 0) {
0761       int subdet, zside, ieta, iphi, depth;
0762       unpackDetId(id, subdet, zside, ieta, iphi, depth);
0763       std::map<std::pair<int, int>, double>::const_iterator itr =
0764           cfactors_.find(std::pair<int, int>(zside * ieta, depth));
0765       if (itr != cfactors_.end()) {
0766         cfac = itr->second;
0767       } else if (etaMax_) {
0768         if (zside > 0 && ieta > etamp_)
0769           cfac = (depth < depMax_) ? cfacmp_[depth] : cfacmp_[depMax_ - 1];
0770         if (zside < 0 && ieta > -etamn_)
0771           cfac = (depth < depMax_) ? cfacmn_[depth] : cfacmn_[depMax_ - 1];
0772       }
0773     }
0774   } else if (useScale_ != 0) {
0775     int subdet, zside, ieta, iphi, depth;
0776     unpackDetId(id, subdet, zside, ieta, iphi, depth);
0777     cfac = getFactor(ieta);
0778   }
0779   return cfac;
0780 }
0781 
0782 bool CalibCorrFactor::readCorrFactor(const char* fname, bool marina) {
0783   bool ok(false);
0784   if (std::string(fname) != "") {
0785     std::ifstream fInput(fname);
0786     if (!fInput.good()) {
0787       std::cout << "Cannot open file " << fname << std::endl;
0788     } else {
0789       char buffer[1024];
0790       unsigned int all(0), good(0);
0791       while (fInput.getline(buffer, 1024)) {
0792         ++all;
0793         if (buffer[0] == '#')
0794           continue;  //ignore comment
0795         std::vector<std::string> items = splitString(std::string(buffer));
0796         if (items.size() != 5) {
0797           std::cout << "Ignore  line: " << buffer << std::endl;
0798         } else {
0799           ++good;
0800           int ieta = (marina) ? std::atoi(items[0].c_str()) : std::atoi(items[1].c_str());
0801           int depth = (marina) ? std::atoi(items[1].c_str()) : std::atoi(items[2].c_str());
0802           float corrf = std::atof(items[3].c_str());
0803           double scale = getFactor(std::abs(ieta));
0804           cfactors_[std::pair<int, int>(ieta, depth)] = scale * corrf;
0805           if (ieta > etamp_ && depth == 1)
0806             etamp_ = ieta;
0807           if (ieta == etamp_ && depth < depMax_)
0808             cfacmp_[depth] = scale * corrf;
0809           if (ieta < etamn_ && depth == 1)
0810             etamn_ = ieta;
0811           if (ieta == etamn_ && depth < depMax_)
0812             cfacmn_[depth] = scale * corrf;
0813         }
0814       }
0815       fInput.close();
0816       std::cout << "Reads total of " << all << " and " << good << " good records"
0817                 << " Max eta (z>0) " << etamp_ << " eta (z<0) " << etamn_ << std::endl;
0818       for (int i = 0; i < depMax_; ++i)
0819         std::cout << "[" << i << "] C+ " << cfacmp_[i] << " C- " << cfacmn_[i] << std::endl;
0820       if (good > 0)
0821         ok = true;
0822     }
0823   }
0824   return ok;
0825 }
0826 
0827 double CalibCorrFactor::getFactor(const int& ieta) {
0828   double scale(1.0);
0829   if (ieta < 16) {
0830     if ((useScale_ == 1) || (useScale_ == 3))
0831       scale = scale_;
0832   } else {
0833     if ((useScale_ == 2) || (useScale_ == 3))
0834       scale = scale_;
0835   }
0836   return scale;
0837 }
0838 
0839 CalibCorr::CalibCorr(const char* infile, int flag, bool debug) : flag_(flag), debug_(debug) {
0840   std::cout << "CalibCorr is created with flag " << flag << ":" << flag_ << " for i/p file " << infile << std::endl;
0841   if (flag == 1)
0842     good_ = readCorrDepth(infile);
0843   else if (flag == 2)
0844     good_ = readCorrResp(infile);
0845   else if (flag == 3)
0846     good_ = readCorrPU(infile);
0847   else if (flag == 4)
0848     good_ = readCorrPhi(infile);
0849   else if (flag == 5)
0850     good_ = readCorrPhis(infile);
0851   else
0852     good_ = readCorrRun(infile);
0853 }
0854 
0855 float CalibCorr::getCorr(int run, unsigned int id) {
0856   float cfac(1.0);
0857   if (good_ == 0)
0858     return cfac;
0859   unsigned idx = correctDetId(id);
0860   if (flag_ == 1) {
0861     std::map<unsigned int, float>::iterator itr = corrFacDepth_.find(idx);
0862     if (itr != corrFacDepth_.end())
0863       cfac = itr->second;
0864   } else if (flag_ == 2) {
0865     std::map<unsigned int, float>::iterator itr = corrFacResp_.find(idx);
0866     if (itr != corrFacResp_.end())
0867       cfac = itr->second;
0868   } else if (flag_ == 4) {
0869     cfac = getPhiCorr(idx);
0870   } else if (flag_ == 5) {
0871     cfac = getPhiCorrs(run, idx);
0872   } else {
0873     int ip(-1);
0874     for (unsigned int k = 0; k < runlow_.size(); ++k) {
0875       unsigned int i = runlow_.size() - k - 1;
0876       if (run >= runlow_[i]) {
0877         ip = (int)(i);
0878         break;
0879       }
0880     }
0881     if (debug_) {
0882       std::cout << "Run " << run << " Perdiod " << ip << std::endl;
0883     }
0884     if (ip >= 0) {
0885       std::map<unsigned int, float>::iterator itr = corrFac_[ip].find(idx);
0886       if (itr != corrFac_[ip].end())
0887         cfac = itr->second;
0888     }
0889   }
0890   if (debug_) {
0891     int subdet, zside, ieta, iphi, depth;
0892     unpackDetId(idx, subdet, zside, ieta, iphi, depth);
0893     std::cout << "ID " << std::hex << id << std::dec << " (Sub " << subdet << " eta " << zside * ieta << " phi " << iphi
0894               << " depth " << depth << ")  Factor " << cfac << std::endl;
0895   }
0896   return cfac;
0897 }
0898 
0899 double CalibCorr::getCorr(const Long64_t& entry) {
0900   if (good_ == 0)
0901     return 1.0;
0902   double cfac(0.0);
0903   std::map<Long64_t, double>::iterator itr = cfactors_.find(entry);
0904   if (itr != cfactors_.end())
0905     cfac = std::min(itr->second, 10.0);
0906   return cfac;
0907 }
0908 
0909 double CalibCorr::getTrueCorr(const Long64_t& entry) {
0910   if (good_ == 0)
0911     return 1.0;
0912   double cfac(0.0);
0913   std::map<Long64_t, double>::iterator itr = cfactors_.find(entry);
0914   if (itr != cfactors_.end())
0915     cfac = itr->second;
0916   return cfac;
0917 }
0918 
0919 double CalibCorr::getPhiCorr(unsigned int idx) {
0920   double cfac(1.0);
0921   if (good_ == 0)
0922     return cfac;
0923   std::map<unsigned int, double>::iterator itr = corrPhiSym_.find(idx);
0924   if (itr != corrPhiSym_.end())
0925     cfac = itr->second;
0926   if (debug_) {
0927     int subdet, zside, ieta, iphi, depth;
0928     unpackDetId(idx, subdet, zside, ieta, iphi, depth);
0929     std::cout << "ID " << std::hex << idx << std::dec << " (Sub " << subdet << " eta " << zside * ieta << " phi "
0930               << iphi << " depth " << depth << ")  Factor " << cfac << std::endl;
0931   }
0932   return cfac;
0933 }
0934 
0935 double CalibCorr::getPhiCorrs(int run, unsigned int idx) {
0936   double cfac(1.0);
0937   if (good_ == 0)
0938     return cfac;
0939   unsigned int runid = runId(run);
0940   std::map<std::pair<unsigned int, unsigned int>, double>::iterator itr =
0941       corrPhiSyms_.find(std::pair<unsigned int, unsigned int>(runid, idx));
0942   if (itr != corrPhiSyms_.end())
0943     cfac = itr->second;
0944   if (debug_) {
0945     int subdet, zside, ieta, iphi, depth;
0946     unpackDetId(idx, subdet, zside, ieta, iphi, depth);
0947     std::cout << "Run " << run << " ID " << runid << ":" << std::hex << idx << std::dec << " (Sub " << subdet << " eta "
0948               << zside * ieta << " phi " << iphi << " depth " << depth << ")  Factor " << cfac << std::endl;
0949   }
0950   return cfac;
0951 }
0952 
0953 bool CalibCorr::absent(const Long64_t& entry) { return (cfactors_.find(entry) == cfactors_.end()); }
0954 
0955 bool CalibCorr::present(const Long64_t& entry) { return (cfactors_.find(entry) != cfactors_.end()); }
0956 
0957 unsigned int CalibCorr::readCorrRun(const char* infile) {
0958   std::cout << "Enters readCorrRun for " << infile << std::endl;
0959   std::ifstream fInput(infile);
0960   unsigned int all(0), good(0), ncorr(0);
0961   if (!fInput.good()) {
0962     std::cout << "Cannot open file " << infile << std::endl;
0963   } else {
0964     char buffer[1024];
0965     while (fInput.getline(buffer, 1024)) {
0966       ++all;
0967       std::string bufferString(buffer);
0968       if (bufferString.substr(0, 5) == "#IOVs") {
0969         std::vector<std::string> items = splitString(bufferString.substr(6));
0970         ncorr = items.size() - 1;
0971         for (unsigned int n = 0; n < ncorr; ++n) {
0972           int run = std::atoi(items[n].c_str());
0973           runlow_.push_back(run);
0974         }
0975         std::cout << ncorr << ":" << runlow_.size() << " Run ranges" << std::endl;
0976         for (unsigned int n = 0; n < runlow_.size(); ++n)
0977           std::cout << " [" << n << "] " << runlow_[n];
0978         std::cout << std::endl;
0979       } else if (buffer[0] == '#') {
0980         continue;  //ignore other comments
0981       } else {
0982         std::vector<std::string> items = splitString(bufferString);
0983         if (items.size() != ncorr + 3) {
0984           std::cout << "Ignore  line: " << buffer << std::endl;
0985         } else {
0986           ++good;
0987           int ieta = std::atoi(items[0].c_str());
0988           int iphi = std::atoi(items[1].c_str());
0989           int depth = std::atoi(items[2].c_str());
0990           unsigned int id = getDetIdHE(ieta, iphi, depth);
0991           for (unsigned int n = 0; n < ncorr; ++n) {
0992             float corrf = std::atof(items[n + 3].c_str());
0993             if (n < nmax_)
0994               corrFac_[n][id] = corrf;
0995           }
0996           if (debug_) {
0997             std::cout << "ID " << std::hex << id << std::dec << ":" << id << " (eta " << ieta << " phi " << iphi
0998                       << " depth " << depth << ")";
0999             for (unsigned int n = 0; n < ncorr; ++n)
1000               std::cout << " " << corrFac_[n][id];
1001             std::cout << std::endl;
1002           }
1003         }
1004       }
1005     }
1006     fInput.close();
1007     std::cout << "Reads total of " << all << " and " << good << " good records of run dependent corrections from "
1008               << infile << std::endl;
1009   }
1010   return good;
1011 }
1012 
1013 unsigned int CalibCorr::readCorrDepth(const char* infile) {
1014   std::cout << "Enters readCorrDepth for " << infile << std::endl;
1015   unsigned int all(0), good(0);
1016   std::ifstream fInput(infile);
1017   if (!fInput.good()) {
1018     std::cout << "Cannot open file " << infile << std::endl;
1019   } else {
1020     char buffer[1024];
1021     while (fInput.getline(buffer, 1024)) {
1022       ++all;
1023       std::string bufferString(buffer);
1024       if (bufferString.substr(0, 5) == "depth") {
1025         continue;  //ignore other comments
1026       } else {
1027         std::vector<std::string> items = splitString(bufferString);
1028         if (items.size() != 3) {
1029           std::cout << "Ignore  line: " << buffer << " Size " << items.size();
1030           for (unsigned int k = 0; k < items.size(); ++k)
1031             std::cout << " [" << k << "] : " << items[k];
1032           std::cout << std::endl;
1033         } else {
1034           ++good;
1035           int ieta = std::atoi(items[1].c_str());
1036           int depth = std::atoi(items[0].c_str());
1037           float corrf = std::atof(items[2].c_str());
1038           int nphi = (std::abs(ieta) > 20) ? 36 : 72;
1039           for (int i = 1; i <= nphi; ++i) {
1040             int iphi = (nphi > 36) ? i : (2 * i - 1);
1041             unsigned int id = getDetIdHE(ieta, iphi, depth);
1042             corrFacDepth_[id] = corrf;
1043             if (debug_) {
1044               std::cout << "ID " << std::hex << id << std::dec << ":" << id << " (eta " << ieta << " phi " << iphi
1045                         << " depth " << depth << ") " << corrFacDepth_[id] << std::endl;
1046             }
1047           }
1048         }
1049       }
1050     }
1051     fInput.close();
1052     std::cout << "Reads total of " << all << " and " << good << " good records of depth dependent factors from "
1053               << infile << std::endl;
1054   }
1055   return good;
1056 }
1057 
1058 unsigned int CalibCorr::readCorrResp(const char* infile) {
1059   std::cout << "Enters readCorrResp for " << infile << std::endl;
1060   unsigned int all(0), good(0), other(0);
1061   std::ifstream fInput(infile);
1062   if (!fInput.good()) {
1063     std::cout << "Cannot open file " << infile << std::endl;
1064   } else {
1065     char buffer[1024];
1066     while (fInput.getline(buffer, 1024)) {
1067       ++all;
1068       std::string bufferString(buffer);
1069       if (bufferString.substr(0, 1) == "#") {
1070         continue;  //ignore other comments
1071       } else {
1072         std::vector<std::string> items = splitString(bufferString);
1073         if (items.size() < 5) {
1074           std::cout << "Ignore  line: " << buffer << " Size " << items.size();
1075           for (unsigned int k = 0; k < items.size(); ++k)
1076             std::cout << " [" << k << "] : " << items[k];
1077           std::cout << std::endl;
1078         } else if (items[3] == "HB" || items[3] == "HE") {
1079           ++good;
1080           int ieta = std::atoi(items[0].c_str());
1081           int iphi = std::atoi(items[1].c_str());
1082           int depth = std::atoi(items[2].c_str());
1083           int subdet = (items[3] == "HE") ? 2 : 1;
1084           float corrf = std::atof(items[4].c_str());
1085           unsigned int id = getDetId(subdet, ieta, iphi, depth);
1086           corrFacResp_[id] = corrf;
1087           if (debug_) {
1088             std::cout << "ID " << std::hex << id << std::dec << ":" << id << " (subdet " << items[3] << ":" << subdet
1089                       << " eta " << ieta << " phi " << iphi << " depth " << depth << ") " << corrFacResp_[id]
1090                       << std::endl;
1091           }
1092         } else {
1093           ++other;
1094         }
1095       }
1096     }
1097     fInput.close();
1098     std::cout << "Reads total of " << all << " and " << good << " good and " << other
1099               << " detector records of depth dependent factors from " << infile << std::endl;
1100   }
1101   return good;
1102 }
1103 
1104 unsigned int CalibCorr::readCorrPU(const char* infile) {
1105   if (std::string(infile) != "") {
1106     std::ifstream fInput(infile);
1107     if (!fInput.good()) {
1108       std::cout << "Cannot open file " << infile << std::endl;
1109     } else {
1110       double val1, val2;
1111       cfactors_.clear();
1112       while (1) {
1113         fInput >> val1 >> val2;
1114         if (!fInput.good())
1115           break;
1116         Long64_t entry = (Long64_t)(val1);
1117         cfactors_[entry] = val2;
1118       }
1119       fInput.close();
1120     }
1121   }
1122   std::cout << "Reads " << cfactors_.size() << " PU correction factors from " << infile << std::endl;
1123   return cfactors_.size();
1124 }
1125 
1126 unsigned int CalibCorr::readCorrPhi(const char* infile) {
1127   std::cout << "Enters readCorrPhi for " << infile << std::endl;
1128   unsigned int all(0), good(0);
1129   std::ifstream fInput(infile);
1130   if (!fInput.good()) {
1131     std::cout << "Cannot open file " << infile << std::endl;
1132   } else {
1133     char buffer[1024];
1134     while (fInput.getline(buffer, 1024)) {
1135       ++all;
1136       std::string bufferString(buffer);
1137       if (bufferString.substr(0, 1) == "#") {
1138         continue;  //ignore other comments
1139       } else {
1140         std::vector<std::string> items = splitString(bufferString);
1141         if (items.size() < 5) {
1142           std::cout << "Ignore  line: " << buffer << " Size " << items.size();
1143           for (unsigned int k = 0; k < items.size(); ++k)
1144             std::cout << " [" << k << "] : " << items[k];
1145           std::cout << std::endl;
1146         } else {
1147           ++good;
1148           int subdet = std::atoi(items[0].c_str());
1149           int ieta = std::atoi(items[1].c_str());
1150           int iphi = std::atoi(items[2].c_str());
1151           int depth = std::atoi(items[3].c_str());
1152           double corrf = std::atof(items[4].c_str());
1153           unsigned int id = packDetId(subdet, ieta, iphi, depth);
1154           corrPhiSym_[id] = corrf;
1155           if (debug_)
1156             std::cout << "ID " << std::hex << id << std::dec << ":" << id << " (subdet " << subdet << " eta " << ieta
1157                       << " phi " << iphi << " depth " << depth << ") " << corrPhiSym_[id] << std::endl;
1158         }
1159       }
1160     }
1161     fInput.close();
1162     std::cout << "Reads total of " << all << " and " << good << " good records of phi-symmetry factors from " << infile
1163               << std::endl;
1164   }
1165   return good;
1166 }
1167 
1168 unsigned int CalibCorr::readCorrPhis(const char* infile0) {
1169   std::cout << "Enters readCorrPhis for " << infile0 << std::endl;
1170   unsigned int all(0), good(0);
1171   std::ifstream fIn(infile0);
1172   if (!fIn.good()) {
1173     std::cout << "Cannot open file " << infile0 << std::endl;
1174   } else {
1175     char buffer0[1024];
1176     while (fIn.getline(buffer0, 1024)) {
1177       std::string bufferString0(buffer0);
1178       if (bufferString0.substr(0, 1) == "#") {
1179         continue;  //ignore other comments
1180       } else {
1181         std::vector<std::string> item = splitString(bufferString0);
1182         if (item.size() > 2) {
1183           int run1 = std::atoi(item[0].c_str());
1184           int run2 = std::atoi(item[1].c_str());
1185           std::string infile(item[2]);
1186           unsigned int runid = runs_.size();
1187           runs_.push_back(std::pair<int, int>(run1, run2));
1188           std::ifstream fInput(infile.c_str());
1189           if (!fInput.good()) {
1190             std::cout << "Cannot open file " << infile << std::endl;
1191           } else {
1192             char buffer[1024];
1193             while (fInput.getline(buffer, 1024)) {
1194               ++all;
1195               std::string bufferString(buffer);
1196               if (bufferString.substr(0, 1) == "#") {
1197                 continue;  //ignore other comments
1198               } else {
1199                 std::vector<std::string> items = splitString(bufferString);
1200                 if (items.size() < 5) {
1201                   std::cout << "Ignore  line: " << buffer << " Size " << items.size();
1202                   for (unsigned int k = 0; k < items.size(); ++k)
1203                     std::cout << " [" << k << "] : " << items[k];
1204                   std::cout << std::endl;
1205                 } else {
1206                   ++good;
1207                   int subdet = std::atoi(items[0].c_str());
1208                   int ieta = std::atoi(items[1].c_str());
1209                   int iphi = std::atoi(items[2].c_str());
1210                   int depth = std::atoi(items[3].c_str());
1211                   double corrf = std::atof(items[4].c_str());
1212                   unsigned int id = packDetId(subdet, ieta, iphi, depth);
1213                   corrPhiSyms_[std::pair<int, int>(runid, id)] = corrf;
1214                   if (debug_)
1215                     std::cout << "ID " << runid << ":" << std::hex << id << std::dec << ":" << id << " (subdet "
1216                               << subdet << " eta " << ieta << " phi " << iphi << " depth " << depth << ") "
1217                               << corrPhiSym_[id] << std::endl;
1218                 }
1219               }
1220             }
1221             fInput.close();
1222           }
1223         }
1224       }
1225     }
1226     fIn.close();
1227     std::cout << "Reads total of " << all << " and " << good << " good records of phi-symmetry factors from "
1228               << runs_.size() << " files as given in " << infile0 << std::endl;
1229   }
1230   return good;
1231 }
1232 
1233 unsigned int CalibCorr::runId(const int& run) {
1234   unsigned int id(runs_.size());
1235   for (unsigned int k = 0; k < runs_.size(); ++k) {
1236     if ((run >= runs_[k].first) && (run <= runs_[k].second)) {
1237       id = k;
1238       break;
1239     }
1240   }
1241   return id;
1242 }
1243 
1244 unsigned int CalibCorr::getDetIdHE(int ieta, int iphi, int depth) { return getDetId(2, ieta, iphi, depth); }
1245 
1246 unsigned int CalibCorr::getDetId(int subdet, int ieta, int iphi, int depth) {
1247   // All numbers used here are described as masks/offsets in
1248   // DataFormats/HcalDetId/interface/HcalDetId.h
1249   unsigned int id_ = ((4 << 28) | ((subdet & 0x7) << 25));
1250   id_ |= ((0x1000000) | ((depth & 0xF) << 20) | ((ieta > 0) ? (0x80000 | (ieta << 10)) : ((-ieta) << 10)) |
1251           (iphi & 0x3FF));
1252   return id_;
1253 }
1254 
1255 unsigned int CalibCorr::correctDetId(const unsigned int& detId) {
1256   int subdet, ieta, zside, depth, iphi;
1257   unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1258   if (subdet == 0) {
1259     if (ieta > 16)
1260       subdet = 2;
1261     else if (ieta == 16 && depth > 2)
1262       subdet = 2;
1263     else
1264       subdet = 1;
1265   }
1266   unsigned int id = getDetId(subdet, ieta * zside, iphi, depth);
1267   if ((id != detId) && debug_) {
1268     std::cout << "Correct Id " << std::hex << detId << " to " << id << std::dec << "(Sub " << subdet << " eta "
1269               << ieta * zside << " phi " << iphi << " depth " << depth << ")" << std::endl;
1270   }
1271   return id;
1272 }
1273 
1274 CalibSelectRBX::CalibSelectRBX(const char* rbxFile, bool debug) : debug_(debug) {
1275   std::cout << "Enters CalibSelectRBX for " << rbxFile << std::endl;
1276   unsigned int all(0), good(0);
1277   std::vector<int> rbxs;
1278   std::ifstream fInput(rbxFile);
1279   if (!fInput.good()) {
1280     std::cout << "Cannot open file " << rbxFile << std::endl;
1281   } else {
1282     char buffer[1024];
1283     while (fInput.getline(buffer, 1024)) {
1284       ++all;
1285       std::string bufferString(buffer);
1286       if (bufferString.substr(0, 1) == "#") {
1287         continue;  //ignore other comments
1288       } else {
1289         std::vector<std::string> items = splitString(bufferString);
1290         if (items.size() != 1) {
1291           std::cout << "Ignore  line: " << buffer << " Size " << items.size() << std::endl;
1292         } else {
1293           ++good;
1294           int rbx = std::atoi(items[0].c_str());
1295           rbxs.push_back(rbx);
1296           int zside = (rbx > 0) ? 1 : -1;
1297           int subdet = (std::abs(rbx) / 100) % 10;
1298           if (subdet != 1)
1299             subdet = 2;
1300           int iphis = std::abs(rbx) % 100;
1301           if (iphis > 0 && iphis <= 18) {
1302             for (int i = 0; i < 4; ++i) {
1303               int iphi = (iphis - 2) * 4 + 3 + i;
1304               if (iphi < 1)
1305                 iphi += 72;
1306               int zsphi = zside * (subdet * 100 + iphi);
1307               zsphis_.push_back(zsphi);
1308             }
1309           }
1310         }
1311       }
1312     }
1313     fInput.close();
1314   }
1315   std::cout << "Select a set of RBXs " << rbxs.size() << " by reading " << all << ":" << good << " records from "
1316             << rbxFile << " with " << zsphis_.size() << " iphi values:";
1317   for (unsigned int i = 0; i < zsphis_.size(); ++i)
1318     std::cout << " " << zsphis_[i];
1319   std::cout << std::endl;
1320 }
1321 
1322 bool CalibSelectRBX::isItRBX(const unsigned int detId) {
1323   bool ok(true);
1324   if (zsphis_.size() > 0) {
1325     int subdet, ieta, zside, depth, iphi;
1326     unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1327     int zsphi = zside * (subdet * 100 + iphi);
1328     ok = (std::find(zsphis_.begin(), zsphis_.end(), zsphi) != zsphis_.end());
1329 
1330     if (debug_) {
1331       std::cout << "isItRBX:subdet|zside|iphi " << subdet << ":" << zside << ":" << iphi << " OK " << ok << std::endl;
1332     }
1333   }
1334   return ok;
1335 }
1336 
1337 bool CalibSelectRBX::isItRBX(const std::vector<unsigned int>* detId) {
1338   bool ok(true);
1339   if (zsphis_.size() > 0) {
1340     ok = true;
1341     for (unsigned int i = 0; i < detId->size(); ++i) {
1342       int subdet, ieta, zside, depth, iphi;
1343       unpackDetId((*detId)[i], subdet, zside, ieta, iphi, depth);
1344       int zsphi = zside * (subdet * 100 + iphi);
1345       ok = (std::find(zsphis_.begin(), zsphis_.end(), zsphi) != zsphis_.end());
1346       if (debug_) {
1347         std::cout << "isItRBX: subdet|zside|iphi " << subdet << ":" << zside << ":" << iphi << " OK " << ok
1348                   << std::endl;
1349       }
1350       if (ok)
1351         break;
1352     }
1353   }
1354   if (debug_)
1355     std::cout << "isItRBX: size " << detId->size() << " OK " << ok << std::endl;
1356   return ok;
1357 }
1358 
1359 bool CalibSelectRBX::isItRBX(const int ieta, const int iphi) {
1360   bool ok(true);
1361   if (zsphis_.size() == 4) {
1362     int zside = (ieta > 0) ? 1 : -1;
1363     int subd1 = (std::abs(ieta) <= 16) ? 1 : 2;
1364     int subd2 = (std::abs(ieta) >= 16) ? 2 : 1;
1365     int zsphi1 = zside * (subd1 * 100 + iphi);
1366     int zsphi2 = zside * (subd2 * 100 + iphi);
1367     ok = ((std::find(zsphis_.begin(), zsphis_.end(), zsphi1) != zsphis_.end()) ||
1368           (std::find(zsphis_.begin(), zsphis_.end(), zsphi2) != zsphis_.end()));
1369   }
1370   if (debug_) {
1371     std::cout << "isItRBX: ieta " << ieta << " iphi " << iphi << " OK " << ok << std::endl;
1372   }
1373   return ok;
1374 }
1375 
1376 CalibDuplicate::CalibDuplicate(const char* fname, int flag, bool debug)
1377     : flag_(flag), nocc_(0), debug_(debug), ok_(false) {
1378   std::cout << "CalibDuplicate called with Flag " << flag_ << " i/p file " << fname << std::endl;
1379   if (flag_ == 0) {
1380     if (strcmp(fname, "") != 0) {
1381       std::ifstream infile(fname);
1382       if (!infile.is_open()) {
1383         std::cout << "Cannot open duplicate file " << fname << std::endl;
1384       } else {
1385         while (1) {
1386           Long64_t jentry;
1387           infile >> jentry;
1388           if (!infile.good())
1389             break;
1390           entries_.push_back(jentry);
1391         }
1392         infile.close();
1393         std::cout << "Reads a list of " << entries_.size() << " events from " << fname << std::endl;
1394         if (entries_.size() > 0)
1395           ok_ = true;
1396       }
1397     } else {
1398       std::cout << "No duplicate events in the input file" << std::endl;
1399     }
1400   } else if (flag_ == 1) {
1401     if (strcmp(fname, "") != 0) {
1402       std::ifstream infile(fname);
1403       if (!infile.is_open()) {
1404         std::cout << "Cannot open depth dependence file " << fname << std::endl;
1405       } else {
1406         unsigned int all(0), good(0);
1407         char buffer[1024];
1408         while (infile.getline(buffer, 1024)) {
1409           ++all;
1410           std::string bufferString(buffer);
1411           if (bufferString.substr(0, 1) == "#") {
1412             continue;  //ignore other comments
1413           } else {
1414             std::vector<std::string> items = splitString(bufferString);
1415             if (items.size() < 3) {
1416               std::cout << "Ignore  line: " << buffer << " Size " << items.size();
1417               for (unsigned int k = 0; k < items.size(); ++k)
1418                 std::cout << " [" << k << "] : " << items[k];
1419               std::cout << std::endl;
1420             } else {
1421               ++good;
1422               int ieta = std::atoi(items[0].c_str());
1423               std::vector<double> weights;
1424               for (unsigned int k = 1; k < items.size(); ++k) {
1425                 double corrf = std::atof(items[k].c_str());
1426                 weights.push_back(corrf);
1427               }
1428               weights_[ieta] = weights;
1429               if (debug_) {
1430                 std::cout << "Eta " << ieta << " with " << weights.size() << " depths having weights:";
1431                 for (unsigned int k = 0; k < weights.size(); ++k)
1432                   std::cout << " " << weights[k];
1433                 std::cout << std::endl;
1434               }
1435             }
1436           }
1437         }
1438         infile.close();
1439         std::cout << "Reads total of " << all << " and " << good << " good records of depth dependent factors from "
1440                   << fname << std::endl;
1441         if (good > 0)
1442           ok_ = true;
1443       }
1444     }
1445   } else if (flag_ == 3) {
1446     ok_ = readCorrFactor(fname);
1447   } else {
1448     flag_ = 2;
1449     if (strcmp(fname, "") != 0) {
1450       std::ifstream infile(fname);
1451       if (!infile.is_open()) {
1452         std::cout << "Cannot open rejection file " << fname << std::endl;
1453       } else {
1454         unsigned int all(0), good(0);
1455         char buffer[1024];
1456         while (infile.getline(buffer, 1024)) {
1457           ++all;
1458           std::string bufferString(buffer);
1459           if (bufferString.substr(0, 1) == "#") {
1460             continue;  //ignore other comments
1461           } else {
1462             std::vector<std::string> items = splitString(bufferString);
1463             if (items.size() != 2) {
1464               std::cout << "Ignore  line: " << buffer << " Size " << items.size();
1465               for (unsigned int k = 0; k < items.size(); ++k)
1466                 std::cout << " [" << k << "] : " << items[k];
1467               std::cout << std::endl;
1468             } else {
1469               ++good;
1470               int ieta = std::atoi(items[0].c_str());
1471               int iphi = std::atoi(items[1].c_str());
1472               etaphi_.push_back(std::pair<int, int>(ieta, iphi));
1473               if (debug_)
1474                 std::cout << "Select channels with iEta " << ieta << " iPhi " << iphi << std::endl;
1475             }
1476           }
1477         }
1478         infile.close();
1479         std::cout << "Reads total of " << all << " and " << good << " good records of rejection candidates from "
1480                   << fname << std::endl;
1481         if (good > 0)
1482           ok_ = true;
1483       }
1484     }
1485   }
1486 }
1487 
1488 bool CalibDuplicate::isDuplicate(long entry) {
1489   if (ok_)
1490     return (std::find(entries_.begin(), entries_.end(), entry) != entries_.end());
1491   else
1492     return false;
1493 }
1494 
1495 double CalibDuplicate::getWeight(unsigned int detId) {
1496   double wt(1.0);
1497   if (ok_) {
1498     int subdet, ieta, zside, depth, iphi;
1499     unpackDetId(detId, subdet, zside, ieta, iphi, depth);
1500     std::map<int, std::vector<double> >::const_iterator itr = weights_.find(ieta);
1501     --depth;
1502     if (depth < 0)
1503       std::cout << "Strange Depth value " << depth << " in " << std::hex << detId << std::dec << std::endl;
1504     if (itr != weights_.end()) {
1505       if (depth < static_cast<int>(itr->second.size()))
1506         wt = itr->second[depth];
1507     }
1508   }
1509   return wt;
1510 }
1511 
1512 bool CalibDuplicate::select(int ieta, int iphi) {
1513   bool flag(false);
1514   if ((ok_) && (flag_ == 2))
1515     flag = (std::find(etaphi_.begin(), etaphi_.end(), std::pair<int, int>(ieta, iphi)) != etaphi_.end());
1516   if (debug_)
1517     std::cout << "Input " << ieta << ":" << iphi << " Flags " << ok_ << ":" << flag_ << ":" << flag << std::endl;
1518   return flag;
1519 }
1520 
1521 double CalibDuplicate::getCorr(int run, int ieta, int depth) {
1522   ++nocc_;
1523   std::map<std::pair<int, int>, std::vector<double> >::const_iterator itr =
1524       corrs_.find(std::pair<int, int>(ieta, depth));
1525   double corr(1.0);
1526   if (itr != corrs_.end()) {
1527     unsigned int irun = runId(run);
1528     corr = (itr->second)[irun];
1529   }
1530   if (debug_ && (nocc_ < 1000))
1531     std::cout << "Run " << run << " ieta:depth " << ieta << ":" << depth << " found " << (itr == corrs_.end())
1532               << " Corr " << corr << std::endl;
1533   return corr;
1534 }
1535 
1536 bool CalibDuplicate::readCorrFactor(const char* infile) {
1537   bool flag(false);
1538   std::ifstream fInput(infile);
1539   if (!fInput.good()) {
1540     std::cout << "CalibDuplicate::readCorrFactor:Cannot open file " << infile << std::endl;
1541   } else {
1542     int nrun, neta;
1543     fInput >> nrun >> neta;
1544     int indx, run1, run2;
1545     for (int k = 0; k < nrun; ++k) {
1546       fInput >> indx >> run1 >> run2;
1547       runs_.push_back(std::pair<int, int>(run1, run2));
1548     }
1549     int eta, depth;
1550     double corr;
1551     for (int k = 0; k < neta; ++k) {
1552       fInput >> eta >> depth;
1553       std::vector<double> corrs;
1554       for (int i = 0; i < nrun; ++i) {
1555         fInput >> corr;
1556         corrs.push_back(corr);
1557       }
1558       corrs_[std::pair<int, int>(eta, depth)] = corrs;
1559     }
1560     fInput.close();
1561     std::cout << "CalibDuplicate::readCorrFactor:Reads information of " << runs_.size() << " runs and " << corrs_.size()
1562               << " channels from " << infile << std::endl;
1563     flag = true;
1564     if (debug_) {
1565       for (unsigned int k = 0; k < runs_.size(); ++k)
1566         std::cout << "Run range[" << k << "] " << runs_[k].first << ":" << runs_[k].second << std::endl;
1567       std::map<std::pair<int, int>, std::vector<double> >::const_iterator itr;
1568       for (itr = corrs_.begin(); itr != corrs_.end(); ++itr) {
1569         std::cout << "eta:depth [" << (itr->first).first << ":" << (itr->first).second << "]";
1570         for (unsigned int i = 0; i < (itr->second).size(); ++i)
1571           std::cout << " " << (itr->second)[i];
1572         std::cout << std::endl;
1573       }
1574     }
1575   }
1576   return flag;
1577 }
1578 
1579 unsigned int CalibDuplicate::runId(const int& run) {
1580   unsigned int id(runs_.size());
1581   for (unsigned int k = 0; k < runs_.size(); ++k) {
1582     if ((run >= runs_[k].first) && (run <= runs_[k].second)) {
1583       id = k;
1584       break;
1585     }
1586   }
1587   return id;
1588 }
1589 
1590 void CalibCorrTest(const char* infile, int flag) {
1591   if (flag >= 0) {
1592     CalibCorr* c1 = new CalibCorr(infile, flag, true);
1593     for (int ieta = 1; ieta < 29; ++ieta) {
1594       int subdet = (ieta > 16) ? 2 : 1;
1595       int depth = (ieta > 16) ? 2 : 1;
1596       unsigned int id1 = ((4 << 28) | ((subdet & 0x7) << 25));
1597       id1 |= ((0x1000000) | ((depth & 0xF) << 20) | (ieta << 10) | 1);
1598       c1->getCorr(0, id1);
1599       id1 |= (0x80000);
1600       c1->getCorr(0, id1);
1601     }
1602     delete c1;
1603   } else {
1604     CalibDuplicate* c1 = new CalibDuplicate(infile, -flag, true);
1605     delete c1;
1606   }
1607 }
1608 
1609 unsigned int stringTest(const std::string& str) {
1610   std::cout << str << " has " << str.size() << " characters\n";
1611   return str.size();
1612 }
1613 
1614 void CalibCorrScale(const char* infile, const char* outfile, double scale) {
1615   std::ofstream myfile;
1616   myfile.open(outfile);
1617   if (!myfile.is_open()) {
1618     std::cout << "** ERROR: Can't open '" << outfile << std::endl;
1619   } else {
1620     if (std::string(infile) != "") {
1621       std::ifstream fInput(infile);
1622       if (!fInput.good()) {
1623         std::cout << "Cannot open file " << infile << std::endl;
1624       } else {
1625         char buffer[1024];
1626         unsigned int all(0), good(0), comment(0);
1627         while (fInput.getline(buffer, 1024)) {
1628           ++all;
1629           if (buffer[0] == '#') {
1630             myfile << buffer << std::endl;
1631             ++comment;
1632             continue;  //ignore comment
1633           }
1634           std::vector<std::string> items = splitString(std::string(buffer));
1635           if (items.size() != 5) {
1636             std::cout << "Ignore  line: " << buffer << std::endl;
1637           } else {
1638             ++good;
1639             int ieta = std::atoi(items[1].c_str());
1640             int depth = std::atoi(items[2].c_str());
1641             float corrf = scale * std::atof(items[3].c_str());
1642             float dcorr = scale * std::atof(items[4].c_str());
1643             myfile << std::setw(10) << items[0] << std::setw(10) << std::dec << ieta << std::setw(10) << depth
1644                    << std::setw(10) << corrf << " " << std::setw(10) << dcorr << std::endl;
1645           }
1646         }
1647         fInput.close();
1648         std::cout << "Reads total of " << all << ", " << comment << " and " << good << " good records from " << infile
1649                   << " and copied to " << outfile << std::endl;
1650       }
1651     }
1652     myfile.close();
1653   }
1654 }
1655 
1656 void CalibCorrScale2(const char* infile, const char* outfile, double scaleB, double scaleT, double scaleE) {
1657   int ietasL[3] = {0, 13, 17};
1658   int ietasH[3] = {14, 18, 29};
1659   double scale[3] = {scaleB, scaleT, scaleE};
1660   std::ofstream myfile;
1661   myfile.open(outfile);
1662   if (!myfile.is_open()) {
1663     std::cout << "** ERROR: Can't open '" << outfile << std::endl;
1664   } else {
1665     if (std::string(infile) != "") {
1666       std::ifstream fInput(infile);
1667       if (!fInput.good()) {
1668         std::cout << "Cannot open file " << infile << std::endl;
1669       } else {
1670         char buffer[1024];
1671         unsigned int all(0), good(0), comment(0);
1672         while (fInput.getline(buffer, 1024)) {
1673           ++all;
1674           if (buffer[0] == '#') {
1675             myfile << buffer << std::endl;
1676             ++comment;
1677             continue;  //ignore comment
1678           }
1679           std::vector<std::string> items = splitString(std::string(buffer));
1680           if (items.size() != 5) {
1681             std::cout << "Ignore  line: " << buffer << std::endl;
1682           } else {
1683             ++good;
1684             int ieta = std::atoi(items[1].c_str());
1685             int depth = std::atoi(items[2].c_str());
1686             int jp(-1);
1687             for (int j = 0; j < 3; ++j) {
1688               if (std::abs(ieta) > ietasL[j] && std::abs(ieta) <= ietasH[j]) {
1689                 if (jp < 0)
1690                   jp = j;
1691               }
1692             }
1693             if (jp < 0)
1694               jp = 2;
1695             float corrf = scale[jp] * std::atof(items[3].c_str());
1696             float dcorr = scale[jp] * std::atof(items[4].c_str());
1697             myfile << std::setw(10) << items[0] << std::setw(10) << std::dec << ieta << std::setw(10) << depth
1698                    << std::setw(10) << corrf << " " << std::setw(10) << dcorr << std::endl;
1699           }
1700         }
1701         fInput.close();
1702         std::cout << "Reads total of " << all << ", " << comment << " and " << good << " good records from " << infile
1703                   << " and copied to " << outfile << std::endl;
1704       }
1705     }
1706     myfile.close();
1707   }
1708 }
1709 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!