Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-06-07 02:29:04 

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