Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-10-29 22:50:44 

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 // double puFactorRho(type, ieta, rho, double eHcal)
0011 //      Returns a multiplicative factor as PU correction evaluated from rho
0012 // double puweight(vtx)
0013 //      Return PU weight for QCD PU sample
0014 // bool fillChain(chain, inputFileList)
0015 //      Prepares a Tchain by chaining several ROOT files specified
0016 // std::vector<std::string> splitString (fLine)
0017 //      Splits a string into several items which are separated by blank in i/p
0018 // CalibCorrFactor(infile, useScale, scale, etamax, debug)
0019 //      A class which reads a file with correction factors and provides
0020 //        bool   doCorr() : flag saying if correction is available
0021 //        double getCorr(id): correction factor for a given cell
0022 // CalibCorr(infile, flag, debug)
0023 //      A second class which reads a file and provide corrections through
0024 //        float getCorr(run, id): Run and ID dependent correction
0025 //        double getCorr(entry): Entry # (in the file) dependent correction
0026 //        bool absent(entry) : if correction factor absent
0027 //        bool present(entry): or present (relevant for ML motivated)
0028 // CalibSelectRBX(rbx, debug)
0029 //      A class for selecting a given Read Out Box and provides
0030 //        bool isItRBX(detId): if it/they is in the chosen RBX
0031 //        bool isItRBX(ieta, iphi): if it is in the chosen RBX
0032 // void CalibCorrTest(infile, flag)
0033 //      Tests a file which contains correction factors used by CalibCorr
0034 //////////////////////////////////////////////////////////////////////////////
0035 
0036 #include <algorithm>
0037 #include <iomanip>
0038 #include <iostream>
0039 #include <fstream>
0040 #include <sstream>
0041 #include <map>
0042 #include <string>
0043 #include <vector>
0044 
0045 #include <TChain.h>
0046 #include <TROOT.h>
0047 
0048 void unpackDetId(unsigned int detId, int& subdet, int& zside, int& ieta, int& iphi, int& depth) {
0049   // The maskings are defined in DataFormats/DetId/interface/DetId.h
0050   //                      and in DataFormats/HcalDetId/interface/HcalDetId.h
0051   // The macro does not invoke the classes there and use them
0052   subdet = ((detId >> 25) & (0x7));
0053   if ((detId & 0x1000000) == 0) {
0054     ieta = ((detId >> 7) & 0x3F);
0055     zside = (detId & 0x2000) ? (1) : (-1);
0056     depth = ((detId >> 14) & 0x1F);
0057     iphi = (detId & 0x3F);
0058   } else {
0059     ieta = ((detId >> 10) & 0x1FF);
0060     zside = (detId & 0x80000) ? (1) : (-1);
0061     depth = ((detId >> 20) & 0xF);
0062     iphi = (detId & 0x3FF);
0063   }
0064 }
0065 
0066 unsigned int truncateId(unsigned int detId, int truncateFlag, bool debug = false) {
0067   //Truncate depth information of DetId's
0068   unsigned int id(detId);
0069   if (debug) {
0070     std::cout << "Truncate 1 " << std::hex << detId << " " << id << std::dec << " Flag " << truncateFlag << std::endl;
0071   }
0072   int subdet, depth, zside, ieta, iphi;
0073   unpackDetId(detId, subdet, zside, ieta, iphi, depth);
0074   if (truncateFlag == 1) {
0075     //Ignore depth index of ieta values of 15 and 16 of HB
0076     if ((subdet == 1) && (ieta > 14))
0077       depth = 1;
0078   } else if (truncateFlag == 2) {
0079     //Ignore depth index of all ieta values
0080     depth = 1;
0081   } else if (truncateFlag == 3) {
0082     //Ignore depth index for depth > 1 in HE
0083     if ((subdet == 2) && (depth > 1))
0084       depth = 2;
0085     else
0086       depth = 1;
0087   } else if (truncateFlag == 4) {
0088     //Ignore depth index for depth > 1 in HB
0089     if ((subdet == 1) && (depth > 1))
0090       depth = 2;
0091     else
0092       depth = 1;
0093   } else if (truncateFlag == 5) {
0094     //Ignore depth index for depth > 1 in HB and HE
0095     if (depth > 1)
0096       depth = 2;
0097   }
0098   id = (subdet << 25) | (0x1000000) | ((depth & 0xF) << 20) | ((zside > 0) ? (0x80000 | (ieta << 10)) : (ieta << 10));
0099   if (debug) {
0100     std::cout << "Truncate 2: " << subdet << " " << zside * ieta << " " << depth << " " << std::hex << id << " input "
0101               << detId << std::dec << std::endl;
0102   }
0103   return id;
0104 }
0105 
0106 unsigned int repackId(const std::string& det, int eta, int depth) {
0107   int subdet = (det == "HE") ? 2 : 1;
0108   int zside = (eta >= 0) ? 1 : -1;
0109   int ieta = (eta >= 0) ? eta : -eta;
0110   unsigned int id =
0111       (subdet << 25) | (0x1000000) | ((depth & 0xF) << 20) | ((zside > 0) ? (0x80000 | (ieta << 10)) : (ieta << 10));
0112   return id;
0113 }
0114 
0115 unsigned int repackId(int eta, int depth) {
0116   int zside = (eta >= 0) ? 1 : -1;
0117   int ieta = (eta >= 0) ? eta : -eta;
0118   int subdet = ((ieta > 16) || ((ieta == 16) && (depth > 3))) ? 2 : 1;
0119   unsigned int id =
0120       (subdet << 25) | (0x1000000) | ((depth & 0xF) << 20) | ((zside > 0) ? (0x80000 | (ieta << 10)) : (ieta << 10));
0121   return id;
0122 }
0123 
0124 double puFactor(int type, int ieta, double pmom, double eHcal, double ediff, bool debug = false) {
0125   double fac(1.0);
0126   if (debug)
0127     std::cout << "Input Type " << type << " ieta " << ieta << " pmon " << pmom << " E " << eHcal << ":" << ediff;
0128   if (type <= 2) {
0129     double frac = (type == 1) ? 0.02 : 0.03;
0130     if (pmom > 0 && ediff > frac * pmom) {
0131       double a1(0), a2(0);
0132       if (type == 1) {
0133         a1 = -0.35;
0134         a2 = -0.65;
0135         if (std::abs(ieta) == 25) {
0136           a2 = -0.30;
0137         } else if (std::abs(ieta) > 25) {
0138           a1 = -0.45;
0139           a2 = -0.10;
0140         }
0141       } else {
0142         a1 = -0.39;
0143         a2 = -0.59;
0144         if (std::abs(ieta) >= 25) {
0145           a1 = -0.283;
0146           a2 = -0.272;
0147         } else if (std::abs(ieta) > 22) {
0148           a1 = -0.238;
0149           a2 = -0.241;
0150         }
0151       }
0152       fac = (1.0 + a1 * (eHcal / pmom) * (ediff / pmom) * (1 + a2 * (ediff / pmom)));
0153       if (debug)
0154         std::cout << " coeff " << a1 << ":" << a2 << " Fac " << fac;
0155     }
0156   } else {
0157     int jeta = std::abs(ieta);
0158     double d2p = (ediff / pmom);
0159     const double DELTA_CUT = 0.03;
0160     if (type == 3) {  // 16pu
0161       const double CONST_COR_COEF[4] = {0.971, 1.008, 0.985, 1.086};
0162       const double LINEAR_COR_COEF[4] = {0, -0.359, -0.251, -0.535};
0163       const double SQUARE_COR_COEF[4] = {0, 0, 0.048, 0.143};
0164       const int PU_IETA_1 = 9;
0165       const int PU_IETA_2 = 16;
0166       const int PU_IETA_3 = 25;
0167       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3));
0168       if (d2p > DELTA_CUT)
0169         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0170       if (debug)
0171         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0172                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0173     } else if (type == 4) {  // 17pu
0174       const double CONST_COR_COEF[4] = {0.974, 1.023, 0.989, 1.077};
0175       const double LINEAR_COR_COEF[4] = {0, -0.524, -0.268, -0.584};
0176       const double SQUARE_COR_COEF[4] = {0, 0, 0.053, 0.170};
0177       const int PU_IETA_1 = 9;
0178       const int PU_IETA_2 = 18;
0179       const int PU_IETA_3 = 25;
0180       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3));
0181       if (d2p > DELTA_CUT)
0182         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0183       if (debug)
0184         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0185                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0186     } else if (type == 5) {  // 18pu
0187       const double CONST_COR_COEF[4] = {0.973, 0.998, 0.992, 0.965};
0188       const double LINEAR_COR_COEF[4] = {0, -0.318, -0.261, -0.406};
0189       const double SQUARE_COR_COEF[4] = {0, 0, 0.047, 0.089};
0190       const int PU_IETA_1 = 7;
0191       const int PU_IETA_2 = 16;
0192       const int PU_IETA_3 = 25;
0193       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3));
0194       if (d2p > DELTA_CUT)
0195         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0196       if (debug)
0197         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0198                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0199     } else if (type == 6) {  // 21pu (old)
0200       const double CONST_COR_COEF[6] = {0.98913, 0.982008, 0.974011, 0.496234, 0.368110, 0.294053};
0201       const double LINEAR_COR_COEF[6] = {-0.0491388, -0.124058, -0.249718, -0.0667390, -0.0770766, -0.0580492};
0202       const double SQUARE_COR_COEF[6] = {0, 0, 0.0368657, 0.00656337, 0.00724508, 0.00568967};
0203       const int PU_IETA_1 = 7;
0204       const int PU_IETA_2 = 16;
0205       const int PU_IETA_3 = 25;
0206       const int PU_IETA_4 = 26;
0207       const int PU_IETA_5 = 27;
0208       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0209                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0210       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0211       if (d2p > deltaCut)
0212         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0213       if (debug)
0214         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0215                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0216     } else if (type == 97) {  // dlphin Try 3
0217       const double CONST_COR_COEF[6] = {0.987617, 0.983421, 0.938622, 0.806662, 0.738354, 0.574195};
0218       const double LINEAR_COR_COEF[6] = {-0.07018610, -0.2494880, -0.1997290, -0.1769320, -0.2427950, -0.1230480};
0219       const double SQUARE_COR_COEF[6] = {0, 0, 0.0263541, 0.0257008, 0.0426584, 0.0200361};
0220       const int PU_IETA_1 = 7;
0221       const int PU_IETA_2 = 16;
0222       const int PU_IETA_3 = 25;
0223       const int PU_IETA_4 = 26;
0224       const int PU_IETA_5 = 27;
0225       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0226                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0227       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0228       if (d2p > deltaCut)
0229         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0230       if (debug)
0231         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0232                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0233     } else if (type == 98) {  // dlphin Try 2
0234       const double CONST_COR_COEF[6] = {0.987665, 0.983468, 0.938628, 0.807241, 0.739132, 0.529059};
0235       const double LINEAR_COR_COEF[6] = {-0.0708906, -0.249995, -0.199683, -0.177692, -0.243436, -0.0668783};
0236       const double SQUARE_COR_COEF[6] = {0, 0, 0.0263163, 0.0260158, 0.0426864, 0.00398778};
0237       const int PU_IETA_1 = 7;
0238       const int PU_IETA_2 = 16;
0239       const int PU_IETA_3 = 25;
0240       const int PU_IETA_4 = 26;
0241       const int PU_IETA_5 = 27;
0242       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0243                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0244       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0245       if (d2p > deltaCut)
0246         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0247       if (debug)
0248         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0249                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0250     } else if (type == 99) {  // dlphin Try 1
0251       const double CONST_COR_COEF[6] = {0.98312, 0.978532, 0.972211, 0.756004, 0.638075, 0.547192};
0252       const double LINEAR_COR_COEF[6] = {-0.0472436, -0.186206, -0.247339, -0.166062, -0.159781, -0.118747};
0253       const double SQUARE_COR_COEF[6] = {0, 0, 0.0356827, 0.0202461, 0.01785078, 0.0123003};
0254       const int PU_IETA_1 = 7;
0255       const int PU_IETA_2 = 16;
0256       const int PU_IETA_3 = 25;
0257       const int PU_IETA_4 = 26;
0258       const int PU_IETA_5 = 27;
0259       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0260                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0261       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0262       if (d2p > deltaCut)
0263         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0264       if (debug)
0265         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0266                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0267     } else {  // 21pu (June, 2021)
0268       const double CONST_COR_COEF[6] = {0.989727, 0.981923, 0.97571, 0.562475, 0.467947, 0.411831};
0269       const double LINEAR_COR_COEF[6] = {-0.0469558, -0.125805, -0.251383, -0.0668994, -0.0964236, -0.0947158};
0270       const double SQUARE_COR_COEF[6] = {0, 0, 0.0399785, 0.00610104, 0.00952528, 0.0100645};
0271       const int PU_IETA_1 = 7;
0272       const int PU_IETA_2 = 16;
0273       const int PU_IETA_3 = 25;
0274       const int PU_IETA_4 = 26;
0275       const int PU_IETA_5 = 27;
0276       unsigned icor = (unsigned(jeta >= PU_IETA_1) + unsigned(jeta >= PU_IETA_2) + unsigned(jeta >= PU_IETA_3) +
0277                        unsigned(jeta >= PU_IETA_4) + unsigned(jeta >= PU_IETA_5));
0278       double deltaCut = (icor > 2) ? 1.0 : DELTA_CUT;
0279       if (d2p > deltaCut)
0280         fac = (CONST_COR_COEF[icor] + LINEAR_COR_COEF[icor] * d2p + SQUARE_COR_COEF[icor] * d2p * d2p);
0281       if (debug)
0282         std::cout << " d2p " << d2p << ":" << DELTA_CUT << " coeff " << icor << ":" << CONST_COR_COEF[icor] << ":"
0283                   << LINEAR_COR_COEF[icor] << ":" << SQUARE_COR_COEF[icor] << " Fac " << fac;
0284     }
0285   }
0286   if (fac < 0 || fac > 1)
0287     fac = 0;
0288   if (debug)
0289     std::cout << " Final factor " << fac << std::endl;
0290   return fac;
0291 }
0292 
0293 double puFactorRho(int type, int ieta, double rho, double eHcal) {
0294   // type = 1: 2017 Data;  2: 2017 MC; 3: 2018 MC; 4: 2018AB; 5: 2018BC
0295   //        6: 2016 MC;
0296   double par[36] = {0.0205395,  -43.0914,   2.67115,    0.239674,    -0.0228009, 0.000476963, 0.137566,  -32.8386,
0297                     3.25886,    0.0863636,  -0.0165639, 0.000413894, 0.206168,   -145.828,    10.3191,   0.531418,
0298                     -0.0578416, 0.00118905, 0.175356,   -175.543,    14.3414,    0.294718,    -0.049836, 0.00106228,
0299                     0.134314,   -175.809,   13.5307,    0.395943,    -0.0539062, 0.00111573,  0.145342,  -98.1904,
0300                     8.14001,    0.205526,   -0.0327818, 0.000726059};
0301   double energy(eHcal);
0302   if (type >= 1 && type <= 6) {
0303     int eta = std::abs(ieta);
0304     int it = 6 * (type - 1);
0305     double ea =
0306         (eta < 20)
0307             ? par[it]
0308             : ((((par[it + 5] * eta + par[it + 4]) * eta + par[it + 3]) * eta + par[it + 2]) * eta + par[it + 1]);
0309     energy -= (rho * ea);
0310   }
0311   return energy;
0312 }
0313 
0314 double puweight(double vtx) {  ///////for QCD PU sample
0315   double a(1.0);
0316   if (vtx < 11)
0317     a = 0.120593;
0318   else if (vtx < 21)
0319     a = 0.58804;
0320   else if (vtx < 31)
0321     a = 1.16306;
0322   else if (vtx < 41)
0323     a = 1.45892;
0324   else if (vtx < 51)
0325     a = 1.35528;
0326   else if (vtx < 61)
0327     a = 1.72032;
0328   else if (vtx < 71)
0329     a = 3.34812;
0330   else if (vtx < 81)
0331     a = 9.70097;
0332   else if (vtx < 91)
0333     a = 9.24839;
0334   else if (vtx < 101)
0335     a = 23.0816;
0336   return a;
0337 }
0338 
0339 bool fillChain(TChain* chain, const char* inputFileList) {
0340   std::string fname(inputFileList);
0341   if (fname.substr(fname.size() - 5, 5) == ".root") {
0342     chain->Add(fname.c_str());
0343   } else {
0344     ifstream infile(inputFileList);
0345     if (!infile.is_open()) {
0346       std::cout << "** ERROR: Can't open '" << inputFileList << "' for input" << std::endl;
0347       return false;
0348     }
0349     while (1) {
0350       infile >> fname;
0351       if (!infile.good())
0352         break;
0353       chain->Add(fname.c_str());
0354     }
0355     infile.close();
0356   }
0357   std::cout << "No. of Entries in this tree : " << chain->GetEntries() << std::endl;
0358   return true;
0359 }
0360 
0361 std::vector<std::string> splitString(const std::string& fLine) {
0362   std::vector<std::string> result;
0363   int start = 0;
0364   bool empty = true;
0365   for (unsigned i = 0; i <= fLine.size(); i++) {
0366     if (fLine[i] == ' ' || i == fLine.size()) {
0367       if (!empty) {
0368         std::string item(fLine, start, i - start);
0369         result.push_back(item);
0370         empty = true;
0371       }
0372       start = i + 1;
0373     } else {
0374       if (empty)
0375         empty = false;
0376     }
0377   }
0378   return result;
0379 }
0380 
0381 class CalibCorrFactor {
0382 public:
0383   CalibCorrFactor(const char* infile, int useScale, double scale, bool etamax, bool marina, bool debug);
0384   ~CalibCorrFactor() {}
0385 
0386   bool doCorr() const { return (corrE_ || (useScale_ != 0)); }
0387   double getCorr(unsigned int id);
0388 
0389 private:
0390   bool readCorrFactor(const char* fName, bool marina);
0391   double getFactor(const int& ieta);
0392 
0393   const int useScale_;
0394   const double scale_;
0395   const bool etaMax_, debug_;
0396   static const int depMax_ = 10;
0397   bool corrE_;
0398   int etamp_, etamn_;
0399   double cfacmp_[depMax_], cfacmn_[depMax_];
0400   std::map<std::pair<int, int>, double> cfactors_;
0401 };
0402 
0403 class CalibCorr {
0404 public:
0405   CalibCorr(const char* infile, int flag, bool debug);
0406   ~CalibCorr() {}
0407 
0408   float getCorr(int run, unsigned int id);
0409   double getCorr(const Long64_t& entry);
0410   double getTrueCorr(const Long64_t& entry);
0411   bool absent(const Long64_t& entry);
0412   bool present(const Long64_t& entry);
0413 
0414 private:
0415   void readCorrRun(const char* infile);
0416   void readCorrDepth(const char* infile);
0417   void readCorrResp(const char* infile);
0418   void readCorrPU(const char* infile);
0419   unsigned int getDetIdHE(int ieta, int iphi, int depth);
0420   unsigned int getDetId(int subdet, int ieta, int iphi, int depth);
0421   unsigned int correctDetId(const unsigned int& detId);
0422 
0423   static const unsigned int nmax_ = 10;
0424   int flag_;
0425   bool debug_;
0426   std::map<unsigned int, float> corrFac_[nmax_], corrFacDepth_, corrFacResp_;
0427   std::map<Long64_t, double> cfactors_;
0428   std::vector<int> runlow_;
0429 };
0430 
0431 class CalibSelectRBX {
0432 public:
0433   CalibSelectRBX(int rbx, bool debug = false);
0434   ~CalibSelectRBX() {}
0435 
0436   bool isItRBX(const unsigned int);
0437   bool isItRBX(const std::vector<unsigned int>*);
0438   bool isItRBX(const int, const int);
0439 
0440 private:
0441   bool debug_;
0442   int subdet_, zside_;
0443   std::vector<int> phis_;
0444 };
0445 
0446 CalibCorrFactor::CalibCorrFactor(const char* infile, int useScale, double scale, bool etamax, bool marina, bool debug)
0447     : useScale_(useScale), scale_(scale), etaMax_(etamax), debug_(debug), etamp_(0), etamn_(0) {
0448   for (int i = 0; i < depMax_; ++i) {
0449     cfacmp_[i] = 1.0;
0450     cfacmn_[i] = 1.0;
0451   }
0452   if (std::string(infile) != "") {
0453     corrE_ = readCorrFactor(infile, marina);
0454     std::cout << "Reads " << cfactors_.size() << " correction factors from " << infile << " with flag " << corrE_
0455               << std::endl
0456               << "Flag for scale " << useScale_ << " with scale " << scale_ << "; flag for etaMax " << etaMax_
0457               << " and flag for Format " << marina << std::endl;
0458   } else {
0459     corrE_ = false;
0460     std::cout << "No correction factors provided; Flag for scale " << useScale_ << " with scale " << scale_
0461               << "; flag for etaMax " << etaMax_ << " and flag for Format " << marina << std::endl;
0462   }
0463 }
0464 
0465 double CalibCorrFactor::getCorr(unsigned int id) {
0466   double cfac(1.0);
0467   if (corrE_) {
0468     int subdet, zside, ieta, iphi, depth;
0469     unpackDetId(id, subdet, zside, ieta, iphi, depth);
0470     std::map<std::pair<int, int>, double>::const_iterator itr =
0471         cfactors_.find(std::pair<int, int>(zside * ieta, depth));
0472     if (itr != cfactors_.end()) {
0473       cfac = itr->second;
0474     } else if (etaMax_) {
0475       if (zside > 0 && ieta > etamp_)
0476         cfac = (depth < depMax_) ? cfacmp_[depth] : cfacmp_[depMax_ - 1];
0477       if (zside < 0 && ieta > -etamn_)
0478         cfac = (depth < depMax_) ? cfacmn_[depth] : cfacmn_[depMax_ - 1];
0479     }
0480   } else if (useScale_ != 0) {
0481     int subdet, zside, ieta, iphi, depth;
0482     unpackDetId(id, subdet, zside, ieta, iphi, depth);
0483     cfac = getFactor(ieta);
0484   }
0485   return cfac;
0486 }
0487 
0488 bool CalibCorrFactor::readCorrFactor(const char* fname, bool marina) {
0489   bool ok(false);
0490   if (std::string(fname) != "") {
0491     std::ifstream fInput(fname);
0492     if (!fInput.good()) {
0493       std::cout << "Cannot open file " << fname << std::endl;
0494     } else {
0495       char buffer[1024];
0496       unsigned int all(0), good(0);
0497       while (fInput.getline(buffer, 1024)) {
0498         ++all;
0499         if (buffer[0] == '#')
0500           continue;  //ignore comment
0501         std::vector<std::string> items = splitString(std::string(buffer));
0502         if (items.size() != 5) {
0503           std::cout << "Ignore  line: " << buffer << std::endl;
0504         } else {
0505           ++good;
0506           int ieta = (marina) ? std::atoi(items[0].c_str()) : std::atoi(items[1].c_str());
0507           int depth = (marina) ? std::atoi(items[1].c_str()) : std::atoi(items[2].c_str());
0508           float corrf = std::atof(items[3].c_str());
0509           double scale = getFactor(std::abs(ieta));
0510           cfactors_[std::pair<int, int>(ieta, depth)] = scale * corrf;
0511           if (ieta > etamp_ && depth == 1)
0512             etamp_ = ieta;
0513           if (ieta == etamp_ && depth < depMax_)
0514             cfacmp_[depth] = scale * corrf;
0515           if (ieta < etamn_ && depth == 1)
0516             etamn_ = ieta;
0517           if (ieta == etamn_ && depth < depMax_)
0518             cfacmn_[depth] = scale * corrf;
0519         }
0520       }
0521       fInput.close();
0522       std::cout << "Reads total of " << all << " and " << good << " good records"
0523                 << " Max eta (z>0) " << etamp_ << " eta (z<0) " << etamn_ << std::endl;
0524       for (int i = 0; i < depMax_; ++i)
0525         std::cout << "[" << i << "] C+ " << cfacmp_[i] << " C- " << cfacmn_[i] << std::endl;
0526       if (good > 0)
0527         ok = true;
0528     }
0529   }
0530   return ok;
0531 }
0532 
0533 double CalibCorrFactor::getFactor(const int& ieta) {
0534   double scale(1.0);
0535   if (ieta < 16) {
0536     if ((useScale_ == 1) || (useScale_ == 3))
0537       scale = scale_;
0538   } else {
0539     if ((useScale_ == 2) || (useScale_ == 3))
0540       scale = scale_;
0541   }
0542   return scale;
0543 }
0544 
0545 CalibCorr::CalibCorr(const char* infile, int flag, bool debug) : flag_(flag), debug_(debug) {
0546   std::cout << "CalibCorr is created with flag " << flag << ":" << flag_ << " for i/p file " << infile << std::endl;
0547   if (flag == 1)
0548     readCorrDepth(infile);
0549   else if (flag == 2)
0550     readCorrResp(infile);
0551   else if (flag == 3)
0552     readCorrPU(infile);
0553   else
0554     readCorrRun(infile);
0555 }
0556 
0557 float CalibCorr::getCorr(int run, unsigned int id) {
0558   float cfac(1.0);
0559   unsigned idx = correctDetId(id);
0560   if (flag_ == 1) {
0561     std::map<unsigned int, float>::iterator itr = corrFacDepth_.find(idx);
0562     if (itr != corrFacDepth_.end())
0563       cfac = itr->second;
0564   } else if (flag_ == 2) {
0565     std::map<unsigned int, float>::iterator itr = corrFacResp_.find(idx);
0566     if (itr != corrFacResp_.end())
0567       cfac = itr->second;
0568   } else {
0569     int ip(-1);
0570     for (unsigned int k = 0; k < runlow_.size(); ++k) {
0571       unsigned int i = runlow_.size() - k - 1;
0572       if (run >= runlow_[i]) {
0573         ip = (int)(i);
0574         break;
0575       }
0576     }
0577     if (debug_) {
0578       std::cout << "Run " << run << " Perdiod " << ip << std::endl;
0579     }
0580     if (ip >= 0) {
0581       std::map<unsigned int, float>::iterator itr = corrFac_[ip].find(idx);
0582       if (itr != corrFac_[ip].end())
0583         cfac = itr->second;
0584     }
0585   }
0586   if (debug_) {
0587     int subdet, zside, ieta, iphi, depth;
0588     unpackDetId(idx, subdet, zside, ieta, iphi, depth);
0589     std::cout << "ID " << std::hex << id << std::dec << " (Sub " << subdet << " eta " << zside * ieta << " phi " << iphi
0590               << " depth " << depth << ")  Factor " << cfac << std::endl;
0591   }
0592   return cfac;
0593 }
0594 
0595 double CalibCorr::getCorr(const Long64_t& entry) {
0596   double cfac(0.0);
0597   std::map<Long64_t, double>::iterator itr = cfactors_.find(entry);
0598   if (itr != cfactors_.end())
0599     cfac = std::min(itr->second, 10.0);
0600   return cfac;
0601 }
0602 
0603 double CalibCorr::getTrueCorr(const Long64_t& entry) {
0604   double cfac(0.0);
0605   std::map<Long64_t, double>::iterator itr = cfactors_.find(entry);
0606   if (itr != cfactors_.end())
0607     cfac = itr->second;
0608   return cfac;
0609 }
0610 
0611 bool CalibCorr::absent(const Long64_t& entry) { return (cfactors_.find(entry) == cfactors_.end()); }
0612 
0613 bool CalibCorr::present(const Long64_t& entry) { return (cfactors_.find(entry) != cfactors_.end()); }
0614 
0615 void CalibCorr::readCorrRun(const char* infile) {
0616   std::cout << "Enters readCorrRun for " << infile << std::endl;
0617   std::ifstream fInput(infile);
0618   unsigned int ncorr(0);
0619   if (!fInput.good()) {
0620     std::cout << "Cannot open file " << infile << std::endl;
0621   } else {
0622     char buffer[1024];
0623     unsigned int all(0), good(0);
0624     while (fInput.getline(buffer, 1024)) {
0625       ++all;
0626       std::string bufferString(buffer);
0627       if (bufferString.substr(0, 5) == "#IOVs") {
0628         std::vector<std::string> items = splitString(bufferString.substr(6));
0629         ncorr = items.size() - 1;
0630         for (unsigned int n = 0; n < ncorr; ++n) {
0631           int run = std::atoi(items[n].c_str());
0632           runlow_.push_back(run);
0633         }
0634         std::cout << ncorr << ":" << runlow_.size() << " Run ranges" << std::endl;
0635         for (unsigned int n = 0; n < runlow_.size(); ++n)
0636           std::cout << " [" << n << "] " << runlow_[n];
0637         std::cout << std::endl;
0638       } else if (buffer[0] == '#') {
0639         continue;  //ignore other comments
0640       } else {
0641         std::vector<std::string> items = splitString(bufferString);
0642         if (items.size() != ncorr + 3) {
0643           std::cout << "Ignore  line: " << buffer << std::endl;
0644         } else {
0645           ++good;
0646           int ieta = std::atoi(items[0].c_str());
0647           int iphi = std::atoi(items[1].c_str());
0648           int depth = std::atoi(items[2].c_str());
0649           unsigned int id = getDetIdHE(ieta, iphi, depth);
0650           for (unsigned int n = 0; n < ncorr; ++n) {
0651             float corrf = std::atof(items[n + 3].c_str());
0652             if (n < nmax_)
0653               corrFac_[n][id] = corrf;
0654           }
0655           if (debug_) {
0656             std::cout << "ID " << std::hex << id << std::dec << ":" << id << " (eta " << ieta << " phi " << iphi
0657                       << " depth " << depth << ")";
0658             for (unsigned int n = 0; n < ncorr; ++n)
0659               std::cout << " " << corrFac_[n][id];
0660             std::cout << std::endl;
0661           }
0662         }
0663       }
0664     }
0665     fInput.close();
0666     std::cout << "Reads total of " << all << " and " << good << " good records of run dependent corrections from "
0667               << infile << std::endl;
0668   }
0669 }
0670 
0671 void CalibCorr::readCorrDepth(const char* infile) {
0672   std::cout << "Enters readCorrDepth for " << infile << std::endl;
0673   std::ifstream fInput(infile);
0674   if (!fInput.good()) {
0675     std::cout << "Cannot open file " << infile << std::endl;
0676   } else {
0677     char buffer[1024];
0678     unsigned int all(0), good(0);
0679     while (fInput.getline(buffer, 1024)) {
0680       ++all;
0681       std::string bufferString(buffer);
0682       if (bufferString.substr(0, 5) == "depth") {
0683         continue;  //ignore other comments
0684       } else {
0685         std::vector<std::string> items = splitString(bufferString);
0686         if (items.size() != 3) {
0687           std::cout << "Ignore  line: " << buffer << " Size " << items.size();
0688           for (unsigned int k = 0; k < items.size(); ++k)
0689             std::cout << " [" << k << "] : " << items[k];
0690           std::cout << std::endl;
0691         } else {
0692           ++good;
0693           int ieta = std::atoi(items[1].c_str());
0694           int depth = std::atoi(items[0].c_str());
0695           float corrf = std::atof(items[2].c_str());
0696           int nphi = (std::abs(ieta) > 20) ? 36 : 72;
0697           for (int i = 1; i <= nphi; ++i) {
0698             int iphi = (nphi > 36) ? i : (2 * i - 1);
0699             unsigned int id = getDetIdHE(ieta, iphi, depth);
0700             corrFacDepth_[id] = corrf;
0701             if (debug_) {
0702               std::cout << "ID " << std::hex << id << std::dec << ":" << id << " (eta " << ieta << " phi " << iphi
0703                         << " depth " << depth << ") " << corrFacDepth_[id] << std::endl;
0704             }
0705           }
0706         }
0707       }
0708     }
0709     fInput.close();
0710     std::cout << "Reads total of " << all << " and " << good << " good records of depth dependent factors from "
0711               << infile << std::endl;
0712   }
0713 }
0714 
0715 void CalibCorr::readCorrResp(const char* infile) {
0716   std::cout << "Enters readCorrResp for " << infile << std::endl;
0717   std::ifstream fInput(infile);
0718   if (!fInput.good()) {
0719     std::cout << "Cannot open file " << infile << std::endl;
0720   } else {
0721     char buffer[1024];
0722     unsigned int all(0), good(0), other(0);
0723     while (fInput.getline(buffer, 1024)) {
0724       ++all;
0725       std::string bufferString(buffer);
0726       if (bufferString.substr(0, 1) == "#") {
0727         continue;  //ignore other comments
0728       } else {
0729         std::vector<std::string> items = splitString(bufferString);
0730         if (items.size() < 5) {
0731           std::cout << "Ignore  line: " << buffer << " Size " << items.size();
0732           for (unsigned int k = 0; k < items.size(); ++k)
0733             std::cout << " [" << k << "] : " << items[k];
0734           std::cout << std::endl;
0735         } else if (items[3] == "HB" || items[3] == "HE") {
0736           ++good;
0737           int ieta = std::atoi(items[0].c_str());
0738           int iphi = std::atoi(items[1].c_str());
0739           int depth = std::atoi(items[2].c_str());
0740           int subdet = (items[3] == "HE") ? 2 : 1;
0741           float corrf = std::atof(items[4].c_str());
0742           unsigned int id = getDetId(subdet, ieta, iphi, depth);
0743           corrFacResp_[id] = corrf;
0744           if (debug_) {
0745             std::cout << "ID " << std::hex << id << std::dec << ":" << id << " (subdet " << items[3] << ":" << subdet
0746                       << " eta " << ieta << " phi " << iphi << " depth " << depth << ") " << corrFacResp_[id]
0747                       << std::endl;
0748           }
0749         } else {
0750           ++other;
0751         }
0752       }
0753     }
0754     fInput.close();
0755     std::cout << "Reads total of " << all << " and " << good << " good and " << other
0756               << " detector records of depth dependent factors from " << infile << std::endl;
0757   }
0758 }
0759 
0760 void CalibCorr::readCorrPU(const char* infile) {
0761   if (std::string(infile) != "") {
0762     std::ifstream fInput(infile);
0763     if (!fInput.good()) {
0764       std::cout << "Cannot open file " << infile << std::endl;
0765     } else {
0766       double val1, val2;
0767       cfactors_.clear();
0768       while (1) {
0769         fInput >> val1 >> val2;
0770         if (!fInput.good())
0771           break;
0772         Long64_t entry = (Long64_t)(val1);
0773         cfactors_[entry] = val2;
0774       }
0775       fInput.close();
0776     }
0777   }
0778   std::cout << "Reads " << cfactors_.size() << " PU correction factors from " << infile << std::endl;
0779 }
0780 
0781 unsigned int CalibCorr::getDetIdHE(int ieta, int iphi, int depth) { return getDetId(2, ieta, iphi, depth); }
0782 
0783 unsigned int CalibCorr::getDetId(int subdet, int ieta, int iphi, int depth) {
0784   // All numbers used here are described as masks/offsets in
0785   // DataFormats/HcalDetId/interface/HcalDetId.h
0786   unsigned int id_ = ((4 << 28) | ((subdet & 0x7) << 25));
0787   id_ |= ((0x1000000) | ((depth & 0xF) << 20) | ((ieta > 0) ? (0x80000 | (ieta << 10)) : ((-ieta) << 10)) |
0788           (iphi & 0x3FF));
0789   return id_;
0790 }
0791 
0792 unsigned int CalibCorr::correctDetId(const unsigned int& detId) {
0793   int subdet, ieta, zside, depth, iphi;
0794   unpackDetId(detId, subdet, zside, ieta, iphi, depth);
0795   if (subdet == 0) {
0796     if (ieta > 16)
0797       subdet = 2;
0798     else if (ieta == 16 && depth > 2)
0799       subdet = 2;
0800     else
0801       subdet = 1;
0802   }
0803   unsigned int id = getDetId(subdet, ieta * zside, iphi, depth);
0804   if ((id != detId) && debug_) {
0805     std::cout << "Correct Id " << std::hex << detId << " to " << id << std::dec << "(Sub " << subdet << " eta "
0806               << ieta * zside << " phi " << iphi << " depth " << depth << ")" << std::endl;
0807   }
0808   return id;
0809 }
0810 
0811 CalibSelectRBX::CalibSelectRBX(int rbx, bool debug) : debug_(debug) {
0812   zside_ = (rbx > 0) ? 1 : -1;
0813   subdet_ = (std::abs(rbx) / 100) % 10;
0814   if (subdet_ != 1)
0815     subdet_ = 2;
0816   int iphis = std::abs(rbx) % 100;
0817   if (iphis > 0 && iphis <= 18) {
0818     for (int i = 0; i < 4; ++i) {
0819       int iphi = (iphis - 2) * 4 + 3 + i;
0820       if (iphi < 1)
0821         iphi += 72;
0822       phis_.push_back(iphi);
0823     }
0824   }
0825   std::cout << "Select RBX " << rbx << " ==> Subdet " << subdet_ << " zside " << zside_ << " with " << phis_.size()
0826             << " iphi values:";
0827   for (unsigned int i = 0; i < phis_.size(); ++i)
0828     std::cout << " " << phis_[i];
0829   std::cout << std::endl;
0830 }
0831 
0832 bool CalibSelectRBX::isItRBX(const unsigned int detId) {
0833   bool ok(true);
0834   if (phis_.size() == 4) {
0835     int subdet, ieta, zside, depth, iphi;
0836     unpackDetId(detId, subdet, zside, ieta, iphi, depth);
0837     ok = ((subdet == subdet_) && (zside == zside_) && (std::find(phis_.begin(), phis_.end(), iphi) != phis_.end()));
0838 
0839     if (debug_) {
0840       std::cout << "isItRBX:subdet|zside|iphi " << subdet << ":" << zside << ":" << iphi << " OK " << ok << std::endl;
0841     }
0842   }
0843   return ok;
0844 }
0845 
0846 bool CalibSelectRBX::isItRBX(const std::vector<unsigned int>* detId) {
0847   bool ok(true);
0848   if (phis_.size() == 4) {
0849     ok = true;
0850     for (unsigned int i = 0; i < detId->size(); ++i) {
0851       int subdet, ieta, zside, depth, iphi;
0852       unpackDetId((*detId)[i], subdet, zside, ieta, iphi, depth);
0853       ok = ((subdet == subdet_) && (zside == zside_) && (std::find(phis_.begin(), phis_.end(), iphi) != phis_.end()));
0854       if (debug_) {
0855         std::cout << "isItRBX: subdet|zside|iphi " << subdet << ":" << zside << ":" << iphi << std::endl;
0856       }
0857       if (ok)
0858         break;
0859     }
0860   }
0861   if (debug_)
0862     std::cout << "isItRBX: size " << detId->size() << " OK " << ok << std::endl;
0863   return ok;
0864 }
0865 
0866 bool CalibSelectRBX::isItRBX(const int ieta, const int iphi) {
0867   bool ok(true);
0868   if (phis_.size() == 4) {
0869     int zside = (ieta > 0) ? 1 : -1;
0870     int subd1 = (std::abs(ieta) <= 16) ? 1 : 2;
0871     int subd2 = (std::abs(ieta) >= 16) ? 2 : 1;
0872     ok = (((subd1 == subdet_) || (subd2 == subdet_)) && (zside == zside_) &&
0873           (std::find(phis_.begin(), phis_.end(), iphi) != phis_.end()));
0874   }
0875   if (debug_) {
0876     std::cout << "isItRBX: ieta " << ieta << " iphi " << iphi << " OK " << ok << std::endl;
0877   }
0878   return ok;
0879 }
0880 
0881 void CalibCorrTest(const char* infile, int flag) {
0882   CalibCorr* c1 = new CalibCorr(infile, flag, true);
0883   for (int ieta = 1; ieta < 29; ++ieta) {
0884     int subdet = (ieta > 16) ? 2 : 1;
0885     int depth = (ieta > 16) ? 2 : 1;
0886     unsigned int id1 = ((4 << 28) | ((subdet & 0x7) << 25));
0887     id1 |= ((0x1000000) | ((depth & 0xF) << 20) | (ieta << 10) | 1);
0888     c1->getCorr(0, id1);
0889     id1 |= (0x80000);
0890     c1->getCorr(0, id1);
0891   }
0892 }
0893 
0894 unsigned int stringTest(const std::string& str) {
0895   std::cout << str << " has " << str.size() << " characters\n";
0896   return str.size();
0897 }
[ 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!