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