File indexing completed on 2022-12-04 23:57:30
0001 #include "Alignment/OfflineValidation/interface/PreparePVTrends.h"
0002
0003 namespace ph = std::placeholders;
0004 namespace pt = boost::property_tree;
0005
0006 PreparePVTrends::PreparePVTrends(const char *outputFileName, int nWorkers, pt::ptree &json)
0007 : outputFileName_(outputFileName), nWorkers_(nWorkers) {
0008 setDirsAndLabels(json);
0009 }
0010
0011 void PreparePVTrends::setDirsAndLabels(pt::ptree &json) {
0012 DirList.clear();
0013 LabelList.clear();
0014 for (const auto &childTree : json) {
0015 DirList.push_back(childTree.first.c_str());
0016 LabelList.push_back(childTree.second.get<std::string>("title"));
0017 }
0018 }
0019
0020 void PreparePVTrends::multiRunPVValidation(bool useRMS, TString lumiInputFile, bool doUnitTest) {
0021 TStopwatch timer;
0022 timer.Start();
0023
0024 gROOT->ProcessLine("gErrorIgnoreLevel = kError;");
0025
0026 ROOT::EnableThreadSafety();
0027 TH1::AddDirectory(kFALSE);
0028
0029 std::ofstream outfile("log.txt");
0030
0031 const Int_t nDirs_ = DirList.size();
0032 TString LegLabels[10];
0033 const char *dirs[10];
0034
0035 std::vector<int> intersection;
0036 std::vector<double> runs;
0037 std::vector<double> x_ticks;
0038 std::vector<double> ex_ticks = {0.};
0039
0040 for (Int_t j = 0; j < nDirs_; j++) {
0041
0042 LegLabels[j] = LabelList[j];
0043 dirs[j] = DirList[j].data();
0044
0045 std::vector<int> currentList = list_files(dirs[j]);
0046 std::vector<int> tempSwap;
0047
0048 std::sort(currentList.begin(), currentList.end());
0049
0050 if (j == 0) {
0051 intersection = currentList;
0052 }
0053
0054 std::sort(intersection.begin(), intersection.end());
0055
0056 std::set_intersection(
0057 currentList.begin(), currentList.end(), intersection.begin(), intersection.end(), std::back_inserter(tempSwap));
0058
0059 intersection.clear();
0060 intersection = tempSwap;
0061 tempSwap.clear();
0062 }
0063
0064 std::ifstream lumifile(lumiInputFile.Data());
0065
0066 std::string line;
0067 while (std::getline(lumifile, line)) {
0068 std::istringstream iss(line);
0069 std::string a, b;
0070 if (!(iss >> a >> b)) {
0071 break;
0072 }
0073 int run = std::stoi(a);
0074
0075
0076 if (std::find(intersection.begin(), intersection.end(), run) == intersection.end()) {
0077 logWarning << " Run: " << run << " is not found in the intersection" << std::endl;
0078 }
0079 }
0080
0081
0082 alignmentTrend dxyPhiMeans_;
0083 alignmentTrend dxyPhiChi2_;
0084 alignmentTrend dxyPhiHiErr_;
0085 alignmentTrend dxyPhiLoErr_;
0086 alignmentTrend dxyPhiKS_;
0087 alignmentTrend dxyPhiHi_;
0088 alignmentTrend dxyPhiLo_;
0089
0090 alignmentTrend dxyEtaMeans_;
0091 alignmentTrend dxyEtaChi2_;
0092 alignmentTrend dxyEtaHiErr_;
0093 alignmentTrend dxyEtaLoErr_;
0094 alignmentTrend dxyEtaKS_;
0095 alignmentTrend dxyEtaHi_;
0096 alignmentTrend dxyEtaLo_;
0097
0098 alignmentTrend dzPhiMeans_;
0099 alignmentTrend dzPhiChi2_;
0100 alignmentTrend dzPhiHiErr_;
0101 alignmentTrend dzPhiLoErr_;
0102 alignmentTrend dzPhiKS_;
0103 alignmentTrend dzPhiHi_;
0104 alignmentTrend dzPhiLo_;
0105
0106 alignmentTrend dzEtaMeans_;
0107 alignmentTrend dzEtaChi2_;
0108 alignmentTrend dzEtaHiErr_;
0109 alignmentTrend dzEtaLoErr_;
0110 alignmentTrend dzEtaKS_;
0111 alignmentTrend dzEtaHi_;
0112 alignmentTrend dzEtaLo_;
0113
0114
0115
0116 std::map<TString, std::vector<unrolledHisto> > dxyVect;
0117 std::map<TString, std::vector<unrolledHisto> > dzVect;
0118
0119 logInfo << " pre do-stuff: " << runs.size() << std::endl;
0120
0121
0122 auto f_processData =
0123 std::bind(processData, ph::_1, intersection, nDirs_, dirs, LegLabels, useRMS, nWorkers_, doUnitTest);
0124
0125
0126
0127
0128 TProcPool procPool(std::min(nWorkers_, intersection.size()));
0129 std::vector<size_t> range(std::min(nWorkers_, intersection.size()));
0130 std::iota(range.begin(), range.end(), 0);
0131
0132 auto extracts = procPool.Map(f_processData, range);
0133
0134
0135 std::sort(extracts.begin(), extracts.end(), [](const outPVtrends &a, const outPVtrends &b) -> bool {
0136 return a.m_index < b.m_index;
0137 });
0138
0139
0140 for (auto extractedTrend : extracts) {
0141 runs.insert(std::end(runs), std::begin(extractedTrend.m_runs), std::end(extractedTrend.m_runs));
0142
0143 for (const auto &label : LegLabels) {
0144
0145 dxyPhiMeans_[label].insert(std::end(dxyPhiMeans_[label]),
0146 std::begin(extractedTrend.m_dxyPhiMeans[label]),
0147 std::end(extractedTrend.m_dxyPhiMeans[label]));
0148 dxyPhiChi2_[label].insert(std::end(dxyPhiChi2_[label]),
0149 std::begin(extractedTrend.m_dxyPhiChi2[label]),
0150 std::end(extractedTrend.m_dxyPhiChi2[label]));
0151 dxyPhiKS_[label].insert(std::end(dxyPhiKS_[label]),
0152 std::begin(extractedTrend.m_dxyPhiKS[label]),
0153 std::end(extractedTrend.m_dxyPhiKS[label]));
0154
0155 dxyPhiHi_[label].insert(std::end(dxyPhiHi_[label]),
0156 std::begin(extractedTrend.m_dxyPhiHi[label]),
0157 std::end(extractedTrend.m_dxyPhiHi[label]));
0158 dxyPhiLo_[label].insert(std::end(dxyPhiLo_[label]),
0159 std::begin(extractedTrend.m_dxyPhiLo[label]),
0160 std::end(extractedTrend.m_dxyPhiLo[label]));
0161
0162
0163 dzPhiMeans_[label].insert(std::end(dzPhiMeans_[label]),
0164 std::begin(extractedTrend.m_dzPhiMeans[label]),
0165 std::end(extractedTrend.m_dzPhiMeans[label]));
0166 dzPhiChi2_[label].insert(std::end(dzPhiChi2_[label]),
0167 std::begin(extractedTrend.m_dzPhiChi2[label]),
0168 std::end(extractedTrend.m_dzPhiChi2[label]));
0169 dzPhiKS_[label].insert(std::end(dzPhiKS_[label]),
0170 std::begin(extractedTrend.m_dzPhiKS[label]),
0171 std::end(extractedTrend.m_dzPhiKS[label]));
0172
0173 dzPhiHi_[label].insert(std::end(dzPhiHi_[label]),
0174 std::begin(extractedTrend.m_dzPhiHi[label]),
0175 std::end(extractedTrend.m_dzPhiHi[label]));
0176 dzPhiLo_[label].insert(std::end(dzPhiLo_[label]),
0177 std::begin(extractedTrend.m_dzPhiLo[label]),
0178 std::end(extractedTrend.m_dzPhiLo[label]));
0179
0180
0181 dxyEtaMeans_[label].insert(std::end(dxyEtaMeans_[label]),
0182 std::begin(extractedTrend.m_dxyEtaMeans[label]),
0183 std::end(extractedTrend.m_dxyEtaMeans[label]));
0184 dxyEtaChi2_[label].insert(std::end(dxyEtaChi2_[label]),
0185 std::begin(extractedTrend.m_dxyEtaChi2[label]),
0186 std::end(extractedTrend.m_dxyEtaChi2[label]));
0187 dxyEtaKS_[label].insert(std::end(dxyEtaKS_[label]),
0188 std::begin(extractedTrend.m_dxyEtaKS[label]),
0189 std::end(extractedTrend.m_dxyEtaKS[label]));
0190
0191 dxyEtaHi_[label].insert(std::end(dxyEtaHi_[label]),
0192 std::begin(extractedTrend.m_dxyEtaHi[label]),
0193 std::end(extractedTrend.m_dxyEtaHi[label]));
0194 dxyEtaLo_[label].insert(std::end(dxyEtaLo_[label]),
0195 std::begin(extractedTrend.m_dxyEtaLo[label]),
0196 std::end(extractedTrend.m_dxyEtaLo[label]));
0197
0198
0199 dzEtaMeans_[label].insert(std::end(dzEtaMeans_[label]),
0200 std::begin(extractedTrend.m_dzEtaMeans[label]),
0201 std::end(extractedTrend.m_dzEtaMeans[label]));
0202 dzEtaChi2_[label].insert(std::end(dzEtaChi2_[label]),
0203 std::begin(extractedTrend.m_dzEtaChi2[label]),
0204 std::end(extractedTrend.m_dzEtaChi2[label]));
0205 dzEtaKS_[label].insert(std::end(dzEtaKS_[label]),
0206 std::begin(extractedTrend.m_dzEtaKS[label]),
0207 std::end(extractedTrend.m_dzEtaKS[label]));
0208
0209 dzEtaHi_[label].insert(std::end(dzEtaHi_[label]),
0210 std::begin(extractedTrend.m_dzEtaHi[label]),
0211 std::end(extractedTrend.m_dzEtaHi[label]));
0212 dzEtaLo_[label].insert(std::end(dzEtaLo_[label]),
0213 std::begin(extractedTrend.m_dzEtaLo[label]),
0214 std::end(extractedTrend.m_dzEtaLo[label]));
0215
0216
0217 dxyVect[label].insert(std::end(dxyVect[label]),
0218 std::begin(extractedTrend.m_dxyVect[label]),
0219 std::end(extractedTrend.m_dxyVect[label]));
0220 dzVect[label].insert(std::end(dzVect[label]),
0221 std::begin(extractedTrend.m_dzVect[label]),
0222 std::end(extractedTrend.m_dzVect[label]));
0223 }
0224 }
0225
0226
0227 for (const auto &label : LegLabels) {
0228 for (unsigned int it = 0; it < dxyPhiMeans_[label].size(); it++) {
0229 dxyPhiHiErr_[label].push_back(std::abs(dxyPhiHi_[label][it] - dxyPhiMeans_[label][it]));
0230 dxyPhiLoErr_[label].push_back(std::abs(dxyPhiLo_[label][it] - dxyPhiMeans_[label][it]));
0231 dxyEtaHiErr_[label].push_back(std::abs(dxyEtaHi_[label][it] - dxyEtaMeans_[label][it]));
0232 dxyEtaLoErr_[label].push_back(std::abs(dxyEtaLo_[label][it] - dxyEtaMeans_[label][it]));
0233
0234 if (VERBOSE) {
0235 logInfo << "label: " << label << " means:" << dxyEtaMeans_[label][it] << " low: " << dxyEtaLo_[label][it]
0236 << " loErr: " << dxyEtaLoErr_[label][it] << std::endl;
0237 }
0238
0239 dzPhiHiErr_[label].push_back(std::abs(dzPhiHi_[label][it] - dzPhiMeans_[label][it]));
0240 dzPhiLoErr_[label].push_back(std::abs(dzPhiLo_[label][it] - dzPhiMeans_[label][it]));
0241 dzEtaHiErr_[label].push_back(std::abs(dzEtaHi_[label][it] - dzEtaMeans_[label][it]));
0242 dzEtaLoErr_[label].push_back(std::abs(dzEtaLo_[label][it] - dzEtaMeans_[label][it]));
0243 }
0244 }
0245
0246
0247 TGraph *g_dxy_phi_vs_run[nDirs_];
0248 TGraphAsymmErrors *gerr_dxy_phi_vs_run[nDirs_];
0249 TGraph *g_chi2_dxy_phi_vs_run[nDirs_];
0250 TGraph *g_KS_dxy_phi_vs_run[nDirs_];
0251
0252 TGraph *g_dxy_phi_hi_vs_run[nDirs_];
0253 TGraph *g_dxy_phi_lo_vs_run[nDirs_];
0254
0255 TGraph *g_dxy_eta_vs_run[nDirs_];
0256 TGraphAsymmErrors *gerr_dxy_eta_vs_run[nDirs_];
0257 TGraph *g_chi2_dxy_eta_vs_run[nDirs_];
0258 TGraph *g_KS_dxy_eta_vs_run[nDirs_];
0259
0260 TGraph *g_dxy_eta_hi_vs_run[nDirs_];
0261 TGraph *g_dxy_eta_lo_vs_run[nDirs_];
0262
0263 TGraph *g_dz_phi_vs_run[nDirs_];
0264 TGraphAsymmErrors *gerr_dz_phi_vs_run[nDirs_];
0265 TGraph *g_chi2_dz_phi_vs_run[nDirs_];
0266 TGraph *g_KS_dz_phi_vs_run[nDirs_];
0267
0268 TGraph *g_dz_phi_hi_vs_run[nDirs_];
0269 TGraph *g_dz_phi_lo_vs_run[nDirs_];
0270
0271 TGraph *g_dz_eta_vs_run[nDirs_];
0272 TGraphAsymmErrors *gerr_dz_eta_vs_run[nDirs_];
0273 TGraph *g_chi2_dz_eta_vs_run[nDirs_];
0274 TGraph *g_KS_dz_eta_vs_run[nDirs_];
0275
0276 TGraph *g_dz_eta_hi_vs_run[nDirs_];
0277 TGraph *g_dz_eta_lo_vs_run[nDirs_];
0278
0279
0280
0281 TH1F *h_RMS_dxy_phi_vs_run[nDirs_];
0282 TH1F *h_RMS_dxy_eta_vs_run[nDirs_];
0283 TH1F *h_RMS_dz_phi_vs_run[nDirs_];
0284 TH1F *h_RMS_dz_eta_vs_run[nDirs_];
0285
0286
0287
0288 TH2F *h2_scatter_dxy_vs_run[nDirs_];
0289 TH2F *h2_scatter_dz_vs_run[nDirs_];
0290
0291
0292 TString theType = "run number";
0293 TString theTypeLabel = "run number";
0294 x_ticks = runs;
0295
0296 pv::bundle theBundle = pv::bundle(nDirs_, theType, theTypeLabel, useRMS);
0297 theBundle.printAll();
0298
0299 TFile *fout = TFile::Open(outputFileName_, "RECREATE");
0300
0301 for (Int_t j = 0; j < nDirs_; j++) {
0302
0303 logInfo << "x_ticks.size()= " << x_ticks.size() << " dxyPhiMeans_[LegLabels[" << j
0304 << "]].size()=" << dxyPhiMeans_[LegLabels[j]].size() << std::endl;
0305
0306
0307 assert(x_ticks.size() == dxyPhiMeans_[LegLabels[j]].size());
0308
0309
0310
0311
0312
0313 auto dxyPhiInputs =
0314 pv::wrappedTrends(dxyPhiMeans_, dxyPhiLo_, dxyPhiHi_, dxyPhiLoErr_, dxyPhiHiErr_, dxyPhiChi2_, dxyPhiKS_);
0315
0316 outputGraphs(dxyPhiInputs,
0317 x_ticks,
0318 ex_ticks,
0319 g_dxy_phi_vs_run[j],
0320 g_chi2_dxy_phi_vs_run[j],
0321 g_KS_dxy_phi_vs_run[j],
0322 g_dxy_phi_lo_vs_run[j],
0323 g_dxy_phi_hi_vs_run[j],
0324 gerr_dxy_phi_vs_run[j],
0325 h_RMS_dxy_phi_vs_run,
0326 theBundle,
0327 pv::dxyphi,
0328 j,
0329 LegLabels[j]);
0330
0331
0332
0333
0334
0335 auto dxyEtaInputs =
0336 pv::wrappedTrends(dxyEtaMeans_, dxyEtaLo_, dxyEtaHi_, dxyEtaLoErr_, dxyEtaHiErr_, dxyEtaChi2_, dxyEtaKS_);
0337
0338 outputGraphs(dxyEtaInputs,
0339 x_ticks,
0340 ex_ticks,
0341 g_dxy_eta_vs_run[j],
0342 g_chi2_dxy_eta_vs_run[j],
0343 g_KS_dxy_eta_vs_run[j],
0344 g_dxy_eta_lo_vs_run[j],
0345 g_dxy_eta_hi_vs_run[j],
0346 gerr_dxy_eta_vs_run[j],
0347 h_RMS_dxy_eta_vs_run,
0348 theBundle,
0349 pv::dxyeta,
0350 j,
0351 LegLabels[j]);
0352
0353
0354
0355
0356
0357 auto dzPhiInputs =
0358 pv::wrappedTrends(dzPhiMeans_, dzPhiLo_, dzPhiHi_, dzPhiLoErr_, dzPhiHiErr_, dzPhiChi2_, dzPhiKS_);
0359
0360 outputGraphs(dzPhiInputs,
0361 x_ticks,
0362 ex_ticks,
0363 g_dz_phi_vs_run[j],
0364 g_chi2_dz_phi_vs_run[j],
0365 g_KS_dz_phi_vs_run[j],
0366 g_dz_phi_lo_vs_run[j],
0367 g_dz_phi_hi_vs_run[j],
0368 gerr_dz_phi_vs_run[j],
0369 h_RMS_dz_phi_vs_run,
0370 theBundle,
0371 pv::dzphi,
0372 j,
0373 LegLabels[j]);
0374
0375
0376
0377
0378
0379 auto dzEtaInputs =
0380 pv::wrappedTrends(dzEtaMeans_, dzEtaLo_, dzEtaHi_, dzEtaLoErr_, dzEtaHiErr_, dzEtaChi2_, dzEtaKS_);
0381
0382 outputGraphs(dzEtaInputs,
0383 x_ticks,
0384 ex_ticks,
0385 g_dz_eta_vs_run[j],
0386 g_chi2_dz_eta_vs_run[j],
0387 g_KS_dz_eta_vs_run[j],
0388 g_dz_eta_lo_vs_run[j],
0389 g_dz_eta_hi_vs_run[j],
0390 gerr_dz_eta_vs_run[j],
0391 h_RMS_dz_eta_vs_run,
0392 theBundle,
0393 pv::dzeta,
0394 j,
0395 LegLabels[j]);
0396
0397
0398
0399
0400
0401 h2_scatter_dxy_vs_run[j] =
0402 new TH2F(Form("h2_scatter_dxy_%s", LegLabels[j].Data()),
0403 Form("scatter of d_{xy} vs %s;%s;d_{xy} [cm]", theType.Data(), theTypeLabel.Data()),
0404 x_ticks.size() - 1,
0405 &(x_ticks[0]),
0406 dxyVect[LegLabels[j]][0].get_n_bins(),
0407 dxyVect[LegLabels[j]][0].get_y_min(),
0408 dxyVect[LegLabels[j]][0].get_y_max());
0409 h2_scatter_dxy_vs_run[j]->SetStats(kFALSE);
0410 h2_scatter_dxy_vs_run[j]->SetTitle(LegLabels[j]);
0411
0412 for (unsigned int runindex = 0; runindex < x_ticks.size(); runindex++) {
0413 for (unsigned int binindex = 0; binindex < dxyVect[LegLabels[j]][runindex].get_n_bins(); binindex++) {
0414 h2_scatter_dxy_vs_run[j]->SetBinContent(runindex + 1,
0415 binindex + 1,
0416 dxyVect[LegLabels[j]][runindex].get_bin_contents().at(binindex) /
0417 dxyVect[LegLabels[j]][runindex].get_integral());
0418 }
0419 }
0420
0421
0422
0423
0424
0425 h2_scatter_dz_vs_run[j] =
0426 new TH2F(Form("h2_scatter_dz_%s", LegLabels[j].Data()),
0427 Form("scatter of d_{z} vs %s;%s;d_{z} [cm]", theType.Data(), theTypeLabel.Data()),
0428 x_ticks.size() - 1,
0429 &(x_ticks[0]),
0430 dzVect[LegLabels[j]][0].get_n_bins(),
0431 dzVect[LegLabels[j]][0].get_y_min(),
0432 dzVect[LegLabels[j]][0].get_y_max());
0433 h2_scatter_dz_vs_run[j]->SetStats(kFALSE);
0434 h2_scatter_dz_vs_run[j]->SetTitle(LegLabels[j]);
0435
0436 for (unsigned int runindex = 0; runindex < x_ticks.size(); runindex++) {
0437 for (unsigned int binindex = 0; binindex < dzVect[LegLabels[j]][runindex].get_n_bins(); binindex++) {
0438 h2_scatter_dz_vs_run[j]->SetBinContent(runindex + 1,
0439 binindex + 1,
0440 dzVect[LegLabels[j]][runindex].get_bin_contents().at(binindex) /
0441 dzVect[LegLabels[j]][runindex].get_integral());
0442 }
0443 }
0444
0445 TString modified_label = (LegLabels[j].ReplaceAll(" ", "_"));
0446 g_dxy_phi_vs_run[j]->Write("mean_" + modified_label + "_dxy_phi_vs_run");
0447 g_chi2_dxy_phi_vs_run[j]->Write("chi2_" + modified_label + "_dxy_phi_vs_run");
0448 g_KS_dxy_phi_vs_run[j]->Write("KS_" + modified_label + "_dxy_phi_vs_run");
0449 g_dxy_phi_hi_vs_run[j]->Write("hi_" + modified_label + "_dxy_phi_hi_vs_run");
0450 g_dxy_phi_lo_vs_run[j]->Write("lo_" + modified_label + "_dxy_phi_lo_vs_run");
0451
0452 g_dxy_eta_vs_run[j]->Write("mean_" + modified_label + "_dxy_eta_vs_run");
0453 g_chi2_dxy_eta_vs_run[j]->Write("chi2_" + modified_label + "_dxy_eta_vs_run");
0454 g_KS_dxy_eta_vs_run[j]->Write("KS_" + modified_label + "_dxy_eta_vs_run");
0455 g_dxy_eta_hi_vs_run[j]->Write("hi_" + modified_label + "_dxy_eta_hi_vs_run");
0456 g_dxy_eta_lo_vs_run[j]->Write("lo_" + modified_label + "_dxy_eta_lo_vs_run");
0457
0458 g_dz_phi_vs_run[j]->Write("mean_" + modified_label + "_dz_phi_vs_run");
0459 g_chi2_dz_phi_vs_run[j]->Write("chi2_" + modified_label + "_dz_phi_vs_run");
0460 g_KS_dz_phi_vs_run[j]->Write("KS_" + modified_label + "_dz_phi_vs_run");
0461 g_dz_phi_hi_vs_run[j]->Write("hi_" + modified_label + "_dz_phi_hi_vs_run");
0462 g_dz_phi_lo_vs_run[j]->Write("lo_" + modified_label + "_dz_phi_lo_vs_run");
0463
0464 g_dz_eta_vs_run[j]->Write("mean_" + modified_label + "_dz_eta_vs_run");
0465 g_chi2_dz_eta_vs_run[j]->Write("chi2_" + modified_label + "_dz_eta_vs_run");
0466 g_KS_dz_eta_vs_run[j]->Write("KS_" + modified_label + "_dz_eta_vs_run");
0467 g_dz_eta_hi_vs_run[j]->Write("hi_" + modified_label + "_dz_eta_hi_vs_run");
0468 g_dz_eta_lo_vs_run[j]->Write("lo_" + modified_label + "_dz_eta_lo_vs_run");
0469
0470 h_RMS_dxy_phi_vs_run[j]->Write("RMS_" + modified_label + "_dxy_phi_vs_run");
0471 h_RMS_dxy_eta_vs_run[j]->Write("RMS_" + modified_label + "_dxy_eta_vs_run");
0472 h_RMS_dz_phi_vs_run[j]->Write("RMS_" + modified_label + "_dz_phi_vs_run");
0473 h_RMS_dz_eta_vs_run[j]->Write("RMS_" + modified_label + "_dz_eta_vs_run");
0474
0475
0476
0477 h2_scatter_dxy_vs_run[j]->Write("Scatter_" + modified_label + "_dxy_vs_run");
0478 h2_scatter_dz_vs_run[j]->Write("Scatter_" + modified_label + "_dz_vs_run");
0479 }
0480
0481
0482 for (int iDir = 0; iDir < nDirs_; iDir++) {
0483 delete g_dxy_phi_vs_run[iDir];
0484 delete g_chi2_dxy_phi_vs_run[iDir];
0485 delete g_KS_dxy_phi_vs_run[iDir];
0486 delete g_dxy_phi_hi_vs_run[iDir];
0487 delete g_dxy_phi_lo_vs_run[iDir];
0488
0489 delete g_dxy_eta_vs_run[iDir];
0490 delete g_chi2_dxy_eta_vs_run[iDir];
0491 delete g_KS_dxy_eta_vs_run[iDir];
0492 delete g_dxy_eta_hi_vs_run[iDir];
0493 delete g_dxy_eta_lo_vs_run[iDir];
0494
0495 delete g_dz_phi_vs_run[iDir];
0496 delete g_chi2_dz_phi_vs_run[iDir];
0497 delete g_KS_dz_phi_vs_run[iDir];
0498 delete g_dz_phi_hi_vs_run[iDir];
0499 delete g_dz_phi_lo_vs_run[iDir];
0500
0501 delete g_dz_eta_vs_run[iDir];
0502 delete g_chi2_dz_eta_vs_run[iDir];
0503 delete g_KS_dz_eta_vs_run[iDir];
0504 delete g_dz_eta_hi_vs_run[iDir];
0505 delete g_dz_eta_lo_vs_run[iDir];
0506
0507 delete h_RMS_dxy_phi_vs_run[iDir];
0508 delete h_RMS_dxy_eta_vs_run[iDir];
0509 delete h_RMS_dz_phi_vs_run[iDir];
0510 delete h_RMS_dz_eta_vs_run[iDir];
0511
0512 delete h2_scatter_dxy_vs_run[iDir];
0513 delete h2_scatter_dz_vs_run[iDir];
0514 }
0515
0516 fout->Close();
0517
0518 timer.Stop();
0519 timer.Print();
0520 }
0521
0522
0523
0524
0525
0526
0527 void PreparePVTrends::outputGraphs(const pv::wrappedTrends &allInputs,
0528 const std::vector<double> &ticks,
0529 const std::vector<double> &ex_ticks,
0530 TGraph *&g_mean,
0531 TGraph *&g_chi2,
0532 TGraph *&g_KS,
0533 TGraph *&g_low,
0534 TGraph *&g_high,
0535 TGraphAsymmErrors *&g_asym,
0536 TH1F *h_RMS[],
0537 const pv::bundle &mybundle,
0538 const pv::view &theView,
0539 const int index,
0540 const TString &label)
0541
0542 {
0543 g_mean = new TGraph(ticks.size(), &(ticks[0]), &((allInputs.getMean()[label])[0]));
0544 g_chi2 = new TGraph(ticks.size(), &(ticks[0]), &((allInputs.getChi2()[label])[0]));
0545 g_KS = new TGraph(ticks.size(), &(ticks[0]), &((allInputs.getKS()[label])[0]));
0546 g_high = new TGraph(ticks.size(), &(ticks[0]), &((allInputs.getHigh()[label])[0]));
0547 g_low = new TGraph(ticks.size(), &(ticks[0]), &((allInputs.getLow()[label])[0]));
0548
0549 g_asym = new TGraphAsymmErrors(ticks.size(),
0550 &(ticks[0]),
0551 &((allInputs.getMean()[label])[0]),
0552 &(ex_ticks[0]),
0553 &(ex_ticks[0]),
0554 &((allInputs.getLowErr()[label])[0]),
0555 &((allInputs.getHighErr()[label])[0]));
0556
0557 g_mean->SetTitle(label);
0558 g_asym->SetTitle(label);
0559
0560
0561 h_RMS[index] = new TH1F(Form("h_RMS_dz_eta_%s", label.Data()), label, ticks.size() - 1, &(ticks[0]));
0562 h_RMS[index]->SetStats(kFALSE);
0563
0564 for (size_t bincounter = 1; bincounter < ticks.size(); bincounter++) {
0565 h_RMS[index]->SetBinContent(
0566 bincounter, std::abs(allInputs.getHigh()[label][bincounter - 1] - allInputs.getLow()[label][bincounter - 1]));
0567 h_RMS[index]->SetBinError(bincounter, 0.01);
0568 }
0569 }
0570
0571
0572
0573
0574
0575
0576 std::vector<int> PreparePVTrends::list_files(const char *dirname, const char *ext)
0577
0578 {
0579 std::vector<int> theRunNumbers;
0580
0581 TSystemDirectory dir(dirname, dirname);
0582 TList *files = dir.GetListOfFiles();
0583 if (files) {
0584 TSystemFile *file;
0585 TString fname;
0586 TIter next(files);
0587 while ((file = (TSystemFile *)next())) {
0588 fname = file->GetName();
0589 if (!file->IsDirectory() && fname.EndsWith(ext) && fname.BeginsWith("PVValidation")) {
0590
0591 TObjArray *bits = fname.Tokenize("_");
0592 TString theRun = bits->At(2)->GetName();
0593
0594 TString formatRun = (theRun.ReplaceAll(".root", "")).ReplaceAll("_", "");
0595
0596 theRunNumbers.push_back(formatRun.Atoi());
0597 }
0598 }
0599 }
0600 return theRunNumbers;
0601 }
0602
0603
0604
0605
0606
0607
0608 TH1F *PreparePVTrends::drawConstantWithErr(TH1F *hist, Int_t iter, Double_t theConst)
0609
0610 {
0611 Int_t nbins = hist->GetNbinsX();
0612 Double_t lowedge = hist->GetBinLowEdge(1);
0613 Double_t highedge = hist->GetBinLowEdge(nbins + 1);
0614
0615 TH1F *hzero = new TH1F(Form("hconst_%s_%i", hist->GetName(), iter),
0616 Form("hconst_%s_%i", hist->GetName(), iter),
0617 nbins,
0618 lowedge,
0619 highedge);
0620 for (Int_t i = 0; i <= hzero->GetNbinsX(); i++) {
0621 hzero->SetBinContent(i, theConst);
0622 hzero->SetBinError(i, hist->GetBinError(i));
0623 }
0624 hzero->SetLineWidth(2);
0625 hzero->SetLineStyle(9);
0626 hzero->SetLineColor(kMagenta);
0627
0628 return hzero;
0629 }
0630
0631
0632
0633
0634
0635
0636 unrolledHisto PreparePVTrends::getUnrolledHisto(TH1F *hist)
0637
0638 {
0639
0640
0641
0642
0643
0644 Double_t y_min = -0.1;
0645 Double_t y_max = 0.1;
0646
0647 std::vector<Double_t> contents;
0648 for (int j = 0; j < hist->GetNbinsX(); j++) {
0649 if (std::abs(hist->GetXaxis()->GetBinCenter(j)) <= 0.1)
0650 contents.push_back(hist->GetBinContent(j + 1));
0651 }
0652
0653 auto ret = unrolledHisto(y_min, y_max, contents.size(), contents);
0654 return ret;
0655 }
0656
0657
0658
0659
0660
0661
0662 pv::biases PreparePVTrends::getBiases(TH1F *hist)
0663
0664 {
0665 int nbins = hist->GetNbinsX();
0666
0667
0668 if (nbins <= 0) {
0669 logError << "No bins in the input histogram";
0670 return pv::biases();
0671 }
0672
0673
0674 double *y = new double[nbins];
0675 double *err = new double[nbins];
0676
0677
0678
0679 for (int j = 0; j < nbins; j++) {
0680 y[j] = hist->GetBinContent(j + 1);
0681 if (hist->GetBinError(j + 1) != 0.) {
0682 err[j] = 1. / (hist->GetBinError(j + 1) * hist->GetBinError(j + 1));
0683 } else {
0684 err[j] = 0.;
0685 }
0686 }
0687
0688 Double_t w_mean = TMath::Mean(nbins, y, err);
0689 Double_t w_rms = TMath::RMS(nbins, y, err);
0690
0691 Double_t mean = TMath::Mean(nbins, y);
0692 Double_t rms = TMath::RMS(nbins, y);
0693
0694 Double_t max = hist->GetMaximum();
0695 Double_t min = hist->GetMinimum();
0696
0697
0698 hist->Fit("pol0", "Q0+");
0699 TF1 *f = (TF1 *)hist->FindObject("pol0");
0700
0701
0702 Double_t chi2 = f->GetChisquare();
0703 Int_t ndf = f->GetNDF();
0704
0705 TH1F *theZero = drawConstantWithErr(hist, 1, 1.);
0706 TH1F *displaced = (TH1F *)hist->Clone("displaced");
0707 displaced->Add(theZero);
0708 Double_t ksScore = std::max(-20., TMath::Log10(displaced->KolmogorovTest(theZero)));
0709
0710
0711
0712
0713
0714
0715
0716
0717 pv::biases result(mean, rms, w_mean, w_rms, min, max, chi2, ndf, ksScore);
0718
0719 delete theZero;
0720 delete displaced;
0721 delete[] y;
0722 delete[] err;
0723 return result;
0724 }
0725
0726
0727
0728
0729
0730
0731 outPVtrends PreparePVTrends::processData(size_t iter,
0732 std::vector<int> intersection,
0733 const Int_t nDirs_,
0734 const char *dirs[10],
0735 TString LegLabels[10],
0736 bool useRMS,
0737 const size_t nWorkers,
0738 bool doUnitTest)
0739
0740 {
0741 outPVtrends ret;
0742
0743 unsigned int effSize = std::min(nWorkers, intersection.size());
0744
0745 unsigned int pitch = std::floor(intersection.size() / effSize);
0746 unsigned int first = iter * pitch;
0747 unsigned int last = (iter == (effSize - 1)) ? intersection.size() : ((iter + 1) * pitch);
0748
0749 logInfo << "iter:" << iter << "| pitch: " << pitch << " [" << first << "-" << last << ")" << std::endl;
0750
0751 ret.m_index = iter;
0752
0753 for (unsigned int n = first; n < last; n++) {
0754
0755
0756
0757
0758
0759
0760 if (VERBOSE) {
0761 logInfo << "iter: " << iter << " " << n << " " << intersection.at(n) << std::endl;
0762 }
0763
0764 TFile *fins[nDirs_];
0765
0766 TH1F *dxyPhiMeanTrend[nDirs_];
0767 TH1F *dxyPhiWidthTrend[nDirs_];
0768 TH1F *dzPhiMeanTrend[nDirs_];
0769 TH1F *dzPhiWidthTrend[nDirs_];
0770
0771
0772
0773
0774
0775
0776
0777
0778
0779
0780
0781 TH1F *dxyEtaMeanTrend[nDirs_];
0782 TH1F *dxyEtaWidthTrend[nDirs_];
0783 TH1F *dzEtaMeanTrend[nDirs_];
0784 TH1F *dzEtaWidthTrend[nDirs_];
0785
0786 TH1F *dxyNormPhiWidthTrend[nDirs_];
0787 TH1F *dxyNormEtaWidthTrend[nDirs_];
0788 TH1F *dzNormPhiWidthTrend[nDirs_];
0789 TH1F *dzNormEtaWidthTrend[nDirs_];
0790
0791 TH1F *dxyNormPtWidthTrend[nDirs_];
0792 TH1F *dzNormPtWidthTrend[nDirs_];
0793 TH1F *dxyPtWidthTrend[nDirs_];
0794 TH1F *dzPtWidthTrend[nDirs_];
0795
0796 TH1F *dxyIntegralTrend[nDirs_];
0797 TH1F *dzIntegralTrend[nDirs_];
0798
0799 bool areAllFilesOK = true;
0800 Int_t lastOpen = 0;
0801
0802
0803 for (Int_t j = 0; j < nDirs_; j++) {
0804
0805 size_t position = std::string(dirs[j]).find('/');
0806 std::string stem = std::string(dirs[j]).substr(position + 1);
0807
0808 fins[j] = new TFile(Form("%s/PVValidation_%s_%i.root", dirs[j], stem.c_str(), intersection[n]));
0809 if (fins[j]->IsZombie()) {
0810 logError << Form("%s/PVValidation_%s_%i.root", dirs[j], stem.c_str(), intersection[n])
0811 << " is a Zombie! cannot combine" << std::endl;
0812 areAllFilesOK = false;
0813 lastOpen = j;
0814 break;
0815 }
0816
0817 if (VERBOSE) {
0818 logInfo << Form("%s/PVValidation_%s_%i.root", dirs[j], stem.c_str(), intersection[n])
0819 << " has size: " << fins[j]->GetSize() << " b ";
0820 }
0821
0822
0823 TH1F *h_tracks = (TH1F *)fins[j]->Get("PVValidation/EventFeatures/h_nTracks");
0824 Double_t numEvents = h_tracks->GetEntries();
0825
0826 if (!doUnitTest) {
0827 if (numEvents < 2500) {
0828 logWarning << "excluding run " << intersection[n] << " because it has less than 2.5k events" << std::endl;
0829 areAllFilesOK = false;
0830 lastOpen = j;
0831 break;
0832 }
0833 } else {
0834 if (numEvents == 0) {
0835 logWarning << "excluding run " << intersection[n] << " because it has 0 events" << std::endl;
0836 areAllFilesOK = false;
0837 lastOpen = j;
0838 break;
0839 }
0840 }
0841
0842 dxyPhiMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dxy_phi");
0843 dxyPhiWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dxy_phi");
0844 dzPhiMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dz_phi");
0845 dzPhiWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dz_phi");
0846
0847
0848
0849
0850
0851
0852 dxyEtaMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dxy_eta");
0853 dxyEtaWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dxy_eta");
0854 dzEtaMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dz_eta");
0855 dzEtaWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dz_eta");
0856
0857
0858
0859
0860
0861
0862 dxyNormPhiWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dxy_phi");
0863 dxyNormEtaWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dxy_eta");
0864 dzNormPhiWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dz_phi");
0865 dzNormEtaWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dz_eta");
0866
0867 dxyNormPtWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dxy_pTCentral");
0868 dzNormPtWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dz_pTCentral");
0869 dxyPtWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dxy_pTCentral");
0870 dzPtWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dz_pTCentral");
0871
0872 dxyIntegralTrend[j] = (TH1F *)fins[j]->Get("PVValidation/ProbeTrackFeatures/h_probedxyRefitV");
0873 dzIntegralTrend[j] = (TH1F *)fins[j]->Get("PVValidation/ProbeTrackFeatures/h_probedzRefitV");
0874
0875
0876
0877 auto dxyPhiBiases = getBiases(dxyPhiMeanTrend[j]);
0878
0879
0880
0881
0882
0883
0884 ret.m_dxyPhiMeans[LegLabels[j]].push_back(dxyPhiBiases.getWeightedMean());
0885 ret.m_dxyPhiChi2[LegLabels[j]].push_back(TMath::Log10(dxyPhiBiases.getNormChi2()));
0886 ret.m_dxyPhiKS[LegLabels[j]].push_back(dxyPhiBiases.getKSScore());
0887
0888
0889
0890 useRMS
0891 ? ret.m_dxyPhiLo[LegLabels[j]].push_back(dxyPhiBiases.getWeightedMean() - 2 * dxyPhiBiases.getWeightedRMS())
0892 : ret.m_dxyPhiLo[LegLabels[j]].push_back(dxyPhiBiases.getMin());
0893 useRMS
0894 ? ret.m_dxyPhiHi[LegLabels[j]].push_back(dxyPhiBiases.getWeightedMean() + 2 * dxyPhiBiases.getWeightedRMS())
0895 : ret.m_dxyPhiHi[LegLabels[j]].push_back(dxyPhiBiases.getMax());
0896
0897 auto dxyEtaBiases = getBiases(dxyEtaMeanTrend[j]);
0898 ret.m_dxyEtaMeans[LegLabels[j]].push_back(dxyEtaBiases.getWeightedMean());
0899 ret.m_dxyEtaChi2[LegLabels[j]].push_back(TMath::Log10(dxyEtaBiases.getNormChi2()));
0900 ret.m_dxyEtaKS[LegLabels[j]].push_back(dxyEtaBiases.getKSScore());
0901 useRMS
0902 ? ret.m_dxyEtaLo[LegLabels[j]].push_back(dxyEtaBiases.getWeightedMean() - 2 * dxyEtaBiases.getWeightedRMS())
0903 : ret.m_dxyEtaLo[LegLabels[j]].push_back(dxyEtaBiases.getMin());
0904 useRMS
0905 ? ret.m_dxyEtaHi[LegLabels[j]].push_back(dxyEtaBiases.getWeightedMean() + 2 * dxyEtaBiases.getWeightedRMS())
0906 : ret.m_dxyEtaHi[LegLabels[j]].push_back(dxyEtaBiases.getMax());
0907
0908 auto dzPhiBiases = getBiases(dzPhiMeanTrend[j]);
0909 ret.m_dzPhiMeans[LegLabels[j]].push_back(dzPhiBiases.getWeightedMean());
0910 ret.m_dzPhiChi2[LegLabels[j]].push_back(TMath::Log10(dzPhiBiases.getNormChi2()));
0911 ret.m_dzPhiKS[LegLabels[j]].push_back(dzPhiBiases.getKSScore());
0912 useRMS ? ret.m_dzPhiLo[LegLabels[j]].push_back(dzPhiBiases.getWeightedMean() - 2 * dzPhiBiases.getWeightedRMS())
0913 : ret.m_dzPhiLo[LegLabels[j]].push_back(dzPhiBiases.getMin());
0914 useRMS ? ret.m_dzPhiHi[LegLabels[j]].push_back(dzPhiBiases.getWeightedMean() + 2 * dzPhiBiases.getWeightedRMS())
0915 : ret.m_dzPhiHi[LegLabels[j]].push_back(dzPhiBiases.getMax());
0916
0917 auto dzEtaBiases = getBiases(dzEtaMeanTrend[j]);
0918 ret.m_dzEtaMeans[LegLabels[j]].push_back(dzEtaBiases.getWeightedMean());
0919 ret.m_dzEtaChi2[LegLabels[j]].push_back(TMath::Log10(dzEtaBiases.getNormChi2()));
0920 ret.m_dzEtaKS[LegLabels[j]].push_back(dzEtaBiases.getKSScore());
0921 useRMS ? ret.m_dzEtaLo[LegLabels[j]].push_back(dzEtaBiases.getWeightedMean() - 2 * dzEtaBiases.getWeightedRMS())
0922 : ret.m_dzEtaLo[LegLabels[j]].push_back(dzEtaBiases.getMin());
0923 useRMS ? ret.m_dzEtaHi[LegLabels[j]].push_back(dzEtaBiases.getWeightedMean() + 2 * dzEtaBiases.getWeightedRMS())
0924 : ret.m_dzEtaHi[LegLabels[j]].push_back(dzEtaBiases.getMax());
0925
0926
0927 ret.m_dxyVect[LegLabels[j]].push_back(getUnrolledHisto(dxyIntegralTrend[j]));
0928 ret.m_dzVect[LegLabels[j]].push_back(getUnrolledHisto(dzIntegralTrend[j]));
0929 }
0930
0931 if (!areAllFilesOK) {
0932
0933 logWarning << "====> not all files are OK" << std::endl;
0934
0935 for (int i = 0; i < lastOpen; i++) {
0936 fins[i]->Close();
0937 }
0938 continue;
0939 } else {
0940 ret.m_runs.push_back(intersection.at(n));
0941 }
0942
0943 if (VERBOSE) {
0944 logInfo << "I am still here - runs.size(): " << ret.m_runs.size() << std::endl;
0945 }
0946
0947
0948
0949 for (int i = 0; i < nDirs_; i++) {
0950 delete dxyPhiMeanTrend[i];
0951 delete dzPhiMeanTrend[i];
0952 delete dxyEtaMeanTrend[i];
0953 delete dzEtaMeanTrend[i];
0954
0955 delete dxyPhiWidthTrend[i];
0956 delete dzPhiWidthTrend[i];
0957 delete dxyEtaWidthTrend[i];
0958 delete dzEtaWidthTrend[i];
0959
0960 delete dxyNormPhiWidthTrend[i];
0961 delete dxyNormEtaWidthTrend[i];
0962 delete dzNormPhiWidthTrend[i];
0963 delete dzNormEtaWidthTrend[i];
0964
0965 delete dxyNormPtWidthTrend[i];
0966 delete dzNormPtWidthTrend[i];
0967 delete dxyPtWidthTrend[i];
0968 delete dzPtWidthTrend[i];
0969
0970 fins[i]->Close();
0971 }
0972
0973 if (VERBOSE) {
0974 logInfo << std::endl;
0975 }
0976 }
0977
0978 return ret;
0979 }