File indexing completed on 2024-05-24 04:07:57
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 std::vector<double> newTicks = ticks;
0562 newTicks.insert(newTicks.end(), ticks.back() + 1.);
0563 h_RMS[index] = new TH1F(Form("h_RMS_dz_eta_%s", label.Data()), label, newTicks.size() - 1, &(newTicks[0]));
0564 h_RMS[index]->SetStats(kFALSE);
0565
0566 for (size_t bincounter = 1; bincounter < ticks.size() + 1; bincounter++) {
0567 h_RMS[index]->SetBinContent(
0568 bincounter, std::abs(allInputs.getHigh()[label][bincounter - 1] - allInputs.getLow()[label][bincounter - 1]));
0569 h_RMS[index]->SetBinError(bincounter, 0.01);
0570 }
0571 }
0572
0573
0574
0575
0576
0577
0578 std::vector<int> PreparePVTrends::list_files(const char *dirname, const char *ext)
0579
0580 {
0581 std::vector<int> theRunNumbers;
0582
0583 TSystemDirectory dir(dirname, dirname);
0584 TList *files = dir.GetListOfFiles();
0585 if (files) {
0586 TSystemFile *file;
0587 TString fname;
0588 TIter next(files);
0589 while ((file = (TSystemFile *)next())) {
0590 fname = file->GetName();
0591 if (!file->IsDirectory() && fname.EndsWith(ext) && fname.BeginsWith("PVValidation")) {
0592
0593 TObjArray *bits = fname.Tokenize("_");
0594 TString theRun = bits->At(2)->GetName();
0595
0596 TString formatRun = (theRun.ReplaceAll(".root", "")).ReplaceAll("_", "");
0597
0598 theRunNumbers.push_back(formatRun.Atoi());
0599 }
0600 }
0601 }
0602 return theRunNumbers;
0603 }
0604
0605
0606
0607
0608
0609
0610 TH1F *PreparePVTrends::drawConstantWithErr(TH1F *hist, Int_t iter, Double_t theConst)
0611
0612 {
0613 Int_t nbins = hist->GetNbinsX();
0614 Double_t lowedge = hist->GetBinLowEdge(1);
0615 Double_t highedge = hist->GetBinLowEdge(nbins + 1);
0616
0617 TH1F *hzero = new TH1F(Form("hconst_%s_%i", hist->GetName(), iter),
0618 Form("hconst_%s_%i", hist->GetName(), iter),
0619 nbins,
0620 lowedge,
0621 highedge);
0622 for (Int_t i = 0; i <= hzero->GetNbinsX(); i++) {
0623 hzero->SetBinContent(i, theConst);
0624 hzero->SetBinError(i, hist->GetBinError(i));
0625 }
0626 hzero->SetLineWidth(2);
0627 hzero->SetLineStyle(9);
0628 hzero->SetLineColor(kMagenta);
0629
0630 return hzero;
0631 }
0632
0633
0634
0635
0636
0637
0638 unrolledHisto PreparePVTrends::getUnrolledHisto(TH1F *hist)
0639
0640 {
0641
0642
0643
0644
0645
0646 Double_t y_min = -0.1;
0647 Double_t y_max = 0.1;
0648
0649 std::vector<Double_t> contents;
0650 for (int j = 0; j < hist->GetNbinsX(); j++) {
0651 if (std::abs(hist->GetXaxis()->GetBinCenter(j)) <= 0.1)
0652 contents.push_back(hist->GetBinContent(j + 1));
0653 }
0654
0655 auto ret = unrolledHisto(y_min, y_max, contents.size(), contents);
0656 return ret;
0657 }
0658
0659
0660
0661
0662
0663
0664 pv::biases PreparePVTrends::getBiases(TH1F *hist)
0665
0666 {
0667 int nbins = hist->GetNbinsX();
0668
0669
0670 if (nbins <= 0) {
0671 logError << "No bins in the input histogram";
0672 return pv::biases();
0673 }
0674
0675
0676 double *y = new double[nbins];
0677 double *err = new double[nbins];
0678
0679
0680
0681 for (int j = 0; j < nbins; j++) {
0682 y[j] = hist->GetBinContent(j + 1);
0683 if (hist->GetBinError(j + 1) != 0.) {
0684 err[j] = 1. / (hist->GetBinError(j + 1) * hist->GetBinError(j + 1));
0685 } else {
0686 err[j] = 0.;
0687 }
0688 }
0689
0690 Double_t w_mean = TMath::Mean(nbins, y, err);
0691 Double_t w_rms = TMath::RMS(nbins, y, err);
0692
0693 Double_t mean = TMath::Mean(nbins, y);
0694 Double_t rms = TMath::RMS(nbins, y);
0695
0696 Double_t max = hist->GetMaximum();
0697 Double_t min = hist->GetMinimum();
0698
0699
0700 hist->Fit("pol0", "Q0+");
0701 TF1 *f = (TF1 *)hist->FindObject("pol0");
0702
0703
0704 Double_t chi2 = f->GetChisquare();
0705 Int_t ndf = f->GetNDF();
0706
0707 TH1F *theZero = drawConstantWithErr(hist, 1, 1.);
0708 TH1F *displaced = (TH1F *)hist->Clone("displaced");
0709 displaced->Add(theZero);
0710 Double_t ksScore = std::max(-20., TMath::Log10(displaced->KolmogorovTest(theZero)));
0711
0712
0713
0714
0715
0716
0717
0718
0719 pv::biases result(mean, rms, w_mean, w_rms, min, max, chi2, ndf, ksScore);
0720
0721 delete theZero;
0722 delete displaced;
0723 delete[] y;
0724 delete[] err;
0725 return result;
0726 }
0727
0728
0729
0730
0731
0732
0733 outPVtrends PreparePVTrends::processData(size_t iter,
0734 std::vector<int> intersection,
0735 const Int_t nDirs_,
0736 const char *dirs[10],
0737 TString LegLabels[10],
0738 bool useRMS,
0739 const size_t nWorkers,
0740 bool doUnitTest)
0741
0742 {
0743 outPVtrends ret;
0744
0745 unsigned int effSize = std::min(nWorkers, intersection.size());
0746
0747 unsigned int pitch = std::floor(intersection.size() / effSize);
0748 unsigned int first = iter * pitch;
0749 unsigned int last = (iter == (effSize - 1)) ? intersection.size() : ((iter + 1) * pitch);
0750
0751 logInfo << "iter:" << iter << "| pitch: " << pitch << " [" << first << "-" << last << ")" << std::endl;
0752
0753 ret.m_index = iter;
0754
0755 for (unsigned int n = first; n < last; n++) {
0756
0757
0758
0759
0760
0761
0762 if (VERBOSE) {
0763 logInfo << "iter: " << iter << " " << n << " " << intersection.at(n) << std::endl;
0764 }
0765
0766 TFile *fins[nDirs_];
0767
0768 TH1F *dxyPhiMeanTrend[nDirs_];
0769 TH1F *dxyPhiWidthTrend[nDirs_];
0770 TH1F *dzPhiMeanTrend[nDirs_];
0771 TH1F *dzPhiWidthTrend[nDirs_];
0772
0773
0774
0775
0776
0777
0778
0779
0780
0781
0782
0783 TH1F *dxyEtaMeanTrend[nDirs_];
0784 TH1F *dxyEtaWidthTrend[nDirs_];
0785 TH1F *dzEtaMeanTrend[nDirs_];
0786 TH1F *dzEtaWidthTrend[nDirs_];
0787
0788 TH1F *dxyNormPhiWidthTrend[nDirs_];
0789 TH1F *dxyNormEtaWidthTrend[nDirs_];
0790 TH1F *dzNormPhiWidthTrend[nDirs_];
0791 TH1F *dzNormEtaWidthTrend[nDirs_];
0792
0793 TH1F *dxyNormPtWidthTrend[nDirs_];
0794 TH1F *dzNormPtWidthTrend[nDirs_];
0795 TH1F *dxyPtWidthTrend[nDirs_];
0796 TH1F *dzPtWidthTrend[nDirs_];
0797
0798 TH1F *dxyIntegralTrend[nDirs_];
0799 TH1F *dzIntegralTrend[nDirs_];
0800
0801 bool areAllFilesOK = true;
0802 Int_t lastOpen = 0;
0803
0804
0805 for (Int_t j = 0; j < nDirs_; j++) {
0806
0807 size_t position = std::string(dirs[j]).find('/');
0808 std::string stem = std::string(dirs[j]).substr(position + 1);
0809
0810 fins[j] = new TFile(Form("%s/PVValidation_%s_%i.root", dirs[j], stem.c_str(), intersection[n]));
0811 if (fins[j]->IsZombie()) {
0812 logError << Form("%s/PVValidation_%s_%i.root", dirs[j], stem.c_str(), intersection[n])
0813 << " is a Zombie! cannot combine" << std::endl;
0814 areAllFilesOK = false;
0815 lastOpen = j;
0816 break;
0817 }
0818
0819 if (VERBOSE) {
0820 logInfo << Form("%s/PVValidation_%s_%i.root", dirs[j], stem.c_str(), intersection[n])
0821 << " has size: " << fins[j]->GetSize() << " b ";
0822 }
0823
0824
0825 TH1F *h_tracks = (TH1F *)fins[j]->Get("PVValidation/EventFeatures/h_nTracks");
0826 Double_t numEvents = h_tracks->GetEntries();
0827
0828 if (!doUnitTest) {
0829 if (numEvents < 2500) {
0830 logWarning << "excluding run " << intersection[n] << " because it has less than 2.5k events" << std::endl;
0831 areAllFilesOK = false;
0832 lastOpen = j;
0833 break;
0834 }
0835 } else {
0836 if (numEvents == 0) {
0837 logWarning << "excluding run " << intersection[n] << " because it has 0 events" << std::endl;
0838 areAllFilesOK = false;
0839 lastOpen = j;
0840 break;
0841 }
0842 }
0843
0844 dxyPhiMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dxy_phi");
0845 dxyPhiWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dxy_phi");
0846 dzPhiMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dz_phi");
0847 dzPhiWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dz_phi");
0848
0849
0850
0851
0852
0853
0854 dxyEtaMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dxy_eta");
0855 dxyEtaWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dxy_eta");
0856 dzEtaMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dz_eta");
0857 dzEtaWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dz_eta");
0858
0859
0860
0861
0862
0863
0864 dxyNormPhiWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dxy_phi");
0865 dxyNormEtaWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dxy_eta");
0866 dzNormPhiWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dz_phi");
0867 dzNormEtaWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dz_eta");
0868
0869 dxyNormPtWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dxy_pTCentral");
0870 dzNormPtWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/norm_widths_dz_pTCentral");
0871 dxyPtWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dxy_pTCentral");
0872 dzPtWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dz_pTCentral");
0873
0874 dxyIntegralTrend[j] = (TH1F *)fins[j]->Get("PVValidation/ProbeTrackFeatures/h_probedxyRefitV");
0875 dzIntegralTrend[j] = (TH1F *)fins[j]->Get("PVValidation/ProbeTrackFeatures/h_probedzRefitV");
0876
0877
0878
0879 auto dxyPhiBiases = getBiases(dxyPhiMeanTrend[j]);
0880
0881
0882
0883
0884
0885
0886 ret.m_dxyPhiMeans[LegLabels[j]].push_back(dxyPhiBiases.getWeightedMean());
0887 ret.m_dxyPhiChi2[LegLabels[j]].push_back(TMath::Log10(dxyPhiBiases.getNormChi2()));
0888 ret.m_dxyPhiKS[LegLabels[j]].push_back(dxyPhiBiases.getKSScore());
0889
0890
0891
0892 useRMS
0893 ? ret.m_dxyPhiLo[LegLabels[j]].push_back(dxyPhiBiases.getWeightedMean() - 2 * dxyPhiBiases.getWeightedRMS())
0894 : ret.m_dxyPhiLo[LegLabels[j]].push_back(dxyPhiBiases.getMin());
0895 useRMS
0896 ? ret.m_dxyPhiHi[LegLabels[j]].push_back(dxyPhiBiases.getWeightedMean() + 2 * dxyPhiBiases.getWeightedRMS())
0897 : ret.m_dxyPhiHi[LegLabels[j]].push_back(dxyPhiBiases.getMax());
0898
0899 auto dxyEtaBiases = getBiases(dxyEtaMeanTrend[j]);
0900 ret.m_dxyEtaMeans[LegLabels[j]].push_back(dxyEtaBiases.getWeightedMean());
0901 ret.m_dxyEtaChi2[LegLabels[j]].push_back(TMath::Log10(dxyEtaBiases.getNormChi2()));
0902 ret.m_dxyEtaKS[LegLabels[j]].push_back(dxyEtaBiases.getKSScore());
0903 useRMS
0904 ? ret.m_dxyEtaLo[LegLabels[j]].push_back(dxyEtaBiases.getWeightedMean() - 2 * dxyEtaBiases.getWeightedRMS())
0905 : ret.m_dxyEtaLo[LegLabels[j]].push_back(dxyEtaBiases.getMin());
0906 useRMS
0907 ? ret.m_dxyEtaHi[LegLabels[j]].push_back(dxyEtaBiases.getWeightedMean() + 2 * dxyEtaBiases.getWeightedRMS())
0908 : ret.m_dxyEtaHi[LegLabels[j]].push_back(dxyEtaBiases.getMax());
0909
0910 auto dzPhiBiases = getBiases(dzPhiMeanTrend[j]);
0911 ret.m_dzPhiMeans[LegLabels[j]].push_back(dzPhiBiases.getWeightedMean());
0912 ret.m_dzPhiChi2[LegLabels[j]].push_back(TMath::Log10(dzPhiBiases.getNormChi2()));
0913 ret.m_dzPhiKS[LegLabels[j]].push_back(dzPhiBiases.getKSScore());
0914 useRMS ? ret.m_dzPhiLo[LegLabels[j]].push_back(dzPhiBiases.getWeightedMean() - 2 * dzPhiBiases.getWeightedRMS())
0915 : ret.m_dzPhiLo[LegLabels[j]].push_back(dzPhiBiases.getMin());
0916 useRMS ? ret.m_dzPhiHi[LegLabels[j]].push_back(dzPhiBiases.getWeightedMean() + 2 * dzPhiBiases.getWeightedRMS())
0917 : ret.m_dzPhiHi[LegLabels[j]].push_back(dzPhiBiases.getMax());
0918
0919 auto dzEtaBiases = getBiases(dzEtaMeanTrend[j]);
0920 ret.m_dzEtaMeans[LegLabels[j]].push_back(dzEtaBiases.getWeightedMean());
0921 ret.m_dzEtaChi2[LegLabels[j]].push_back(TMath::Log10(dzEtaBiases.getNormChi2()));
0922 ret.m_dzEtaKS[LegLabels[j]].push_back(dzEtaBiases.getKSScore());
0923 useRMS ? ret.m_dzEtaLo[LegLabels[j]].push_back(dzEtaBiases.getWeightedMean() - 2 * dzEtaBiases.getWeightedRMS())
0924 : ret.m_dzEtaLo[LegLabels[j]].push_back(dzEtaBiases.getMin());
0925 useRMS ? ret.m_dzEtaHi[LegLabels[j]].push_back(dzEtaBiases.getWeightedMean() + 2 * dzEtaBiases.getWeightedRMS())
0926 : ret.m_dzEtaHi[LegLabels[j]].push_back(dzEtaBiases.getMax());
0927
0928
0929 ret.m_dxyVect[LegLabels[j]].push_back(getUnrolledHisto(dxyIntegralTrend[j]));
0930 ret.m_dzVect[LegLabels[j]].push_back(getUnrolledHisto(dzIntegralTrend[j]));
0931 }
0932
0933 if (!areAllFilesOK) {
0934
0935 logWarning << "====> not all files are OK" << std::endl;
0936
0937 for (int i = 0; i < lastOpen; i++) {
0938 fins[i]->Close();
0939 }
0940 continue;
0941 } else {
0942 ret.m_runs.push_back(intersection.at(n));
0943 }
0944
0945 if (VERBOSE) {
0946 logInfo << "I am still here - runs.size(): " << ret.m_runs.size() << std::endl;
0947 }
0948
0949
0950
0951 for (int i = 0; i < nDirs_; i++) {
0952 delete dxyPhiMeanTrend[i];
0953 delete dzPhiMeanTrend[i];
0954 delete dxyEtaMeanTrend[i];
0955 delete dzEtaMeanTrend[i];
0956
0957 delete dxyPhiWidthTrend[i];
0958 delete dzPhiWidthTrend[i];
0959 delete dxyEtaWidthTrend[i];
0960 delete dzEtaWidthTrend[i];
0961
0962 delete dxyNormPhiWidthTrend[i];
0963 delete dxyNormEtaWidthTrend[i];
0964 delete dzNormPhiWidthTrend[i];
0965 delete dzNormEtaWidthTrend[i];
0966
0967 delete dxyNormPtWidthTrend[i];
0968 delete dzNormPtWidthTrend[i];
0969 delete dxyPtWidthTrend[i];
0970 delete dzPtWidthTrend[i];
0971
0972 fins[i]->Close();
0973 }
0974
0975 if (VERBOSE) {
0976 logInfo << std::endl;
0977 }
0978 }
0979
0980 return ret;
0981 }