Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-01-04 00:29:09 

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