File indexing completed on 2021-10-29 22:50:44
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 #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
0050
0051
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
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
0076 if ((subdet == 1) && (ieta > 14))
0077 depth = 1;
0078 } else if (truncateFlag == 2) {
0079
0080 depth = 1;
0081 } else if (truncateFlag == 3) {
0082
0083 if ((subdet == 2) && (depth > 1))
0084 depth = 2;
0085 else
0086 depth = 1;
0087 } else if (truncateFlag == 4) {
0088
0089 if ((subdet == 1) && (depth > 1))
0090 depth = 2;
0091 else
0092 depth = 1;
0093 } else if (truncateFlag == 5) {
0094
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) {
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) {
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) {
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) {
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) {
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) {
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) {
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 {
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
0295
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) {
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;
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;
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;
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;
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
0785
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 }