OfflineValidation/src/PreparePVTrends.cc

0001 #include "Alignment/OfflineValidation/interface/PreparePVTrends.h"
0002
0003 namespace ph = std::placeholders;  // for _1, _2, _3...
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     // Retrieve labels
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     }  // error
0073     int run = std::stoi(a);
0074
0075     // check if the run is in the list
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   // book the vectors of values
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   // unrolled histos
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   //we should use std::bind to create a functor and then pass it to the procPool
0122   auto f_processData =
0123       std::bind(processData, ph::_1, intersection, nDirs_, dirs, LegLabels, useRMS, nWorkers_, doUnitTest);
0124
0125   //f_processData(0);
0126   //logInfo<<" post do-stuff: " <<  runs.size() << std::endl;
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   //procPool.Map([&f_processData](size_t a) { f_processData(a); },{1,2,3});
0132   auto extracts = procPool.Map(f_processData, range);
0133
0134   // sort the extracts according to the global index
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   // re-assemble everything together
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   // extra vectors for low and high boundaries
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   // bias on the mean
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   //TGraph *gprime_dxy_phi_vs_run[nDirs_];
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   //TGraph *gprime_dxy_eta_vs_run[nDirs_];
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   //TGraph *gprime_dz_phi_vs_run[nDirs_];
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   //TGraph *gprime_dz_eta_vs_run[nDirs_];
0276   TGraph *g_dz_eta_hi_vs_run[nDirs_];
0277   TGraph *g_dz_eta_lo_vs_run[nDirs_];
0278
0279   // resolutions
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   // scatters of integrated bias
0287
0288   TH2F *h2_scatter_dxy_vs_run[nDirs_];
0289   TH2F *h2_scatter_dz_vs_run[nDirs_];
0290
0291   // decide the type
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     // check on the sanity
0303     logInfo << "x_ticks.size()= " << x_ticks.size() << " dxyPhiMeans_[LegLabels[" << j
0304             << "]].size()=" << dxyPhiMeans_[LegLabels[j]].size() << std::endl;
0305
0306     // otherwise something very bad has happened
0307     assert(x_ticks.size() == dxyPhiMeans_[LegLabels[j]].size());
0308
0309     // *************************************
0310     // dxy vs phi
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     // dxy vs eta
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     // dz vs phi
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     // dz vs eta
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     // Integrated bias dxy scatter plots
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     // Integrated bias dz scatter plots
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     // scatter
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   // do all the deletes
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 /*! \fn outputGraphs
0523  *  \brief function to build the output graphs
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   // scatter or RMS TH1
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 /*! \fn list_files
0574  *  \brief utility function to list of filles in a directory
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         //logInfo << fname.Data() << std::endl;
0593         TObjArray *bits = fname.Tokenize("_");
0594         TString theRun = bits->At(2)->GetName();
0595         //logInfo << theRun << std::endl;
0596         TString formatRun = (theRun.ReplaceAll(".root", "")).ReplaceAll("_", "");
0597         //logInfo << dirname << " "<< formatRun.Atoi() << std::endl;
0598         theRunNumbers.push_back(formatRun.Atoi());
0599       }
0600     }
0601   }
0602   return theRunNumbers;
0603 }
0604
0605 /*! \fn DrawConstant
0606  *  \brief utility function to draw a constant histogram with erros !=0
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 /*! \fn getUnrolledHisto
0634  *  \brief utility function to tranform a TH1 into a vector of floats
0635  */
0636
0637 /*--------------------------------------------------------------------*/
0638 unrolledHisto PreparePVTrends::getUnrolledHisto(TH1F *hist)
0639 /*--------------------------------------------------------------------*/
0640 {
0641   /*
0642     Double_t y_min = hist->GetBinLowEdge(1);
0643     Double_t y_max = hist->GetBinLowEdge(hist->GetNbinsX()+1);
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 /*! \fn getBiases
0660  *  \brief utility function to extract characterization of the PV bias plot
0661  */
0662
0663 /*--------------------------------------------------------------------*/
0664 pv::biases PreparePVTrends::getBiases(TH1F *hist)
0665 /*--------------------------------------------------------------------*/
0666 {
0667   int nbins = hist->GetNbinsX();
0668   // if there are no bins in the histogram then return default constructed object
0669   // shouldn't really ever happen
0670   if (nbins <= 0) {
0671     logError << "No bins in the input histogram";
0672     return pv::biases();
0673   }
0674
0675   //extract median from histogram
0676   double *y = new double[nbins];
0677   double *err = new double[nbins];
0678
0679   // remember for weight means <x> = sum_i (x_i* w_i) / sum_i w_i ; where w_i = 1/sigma^2_i
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   // in case one would like to use a pol0 fit
0700   hist->Fit("pol0", "Q0+");
0701   TF1 *f = (TF1 *)hist->FindObject("pol0");
0702   //f->SetLineColor(hist->GetLineColor());
0703   //f->SetLineStyle(hist->GetLineStyle());
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     std::pair<std::pair<Double_t,Double_t>, Double_t> result;
0714     std::pair<Double_t,Double_t> resultBounds;
0715     resultBounds = useRMS_ ? std::make_pair(mean-rms,mean+rms) :  std::make_pair(min,max)  ;
0716     result = make_pair(resultBounds,mean);
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 /*! \fn processData
0729  *  \brief function where the magic happens, take the raw inputs and creates the output Trends
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     //in case of debug, use only 50
0757     //for(unsigned int n=0; n<50;n++){
0758
0759     //if(intersection.at(n)!=283946)
0760     //  continue;
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     //TH1F *dxyLadderMeanTrend[nDirs_];
0774     //TH1F *dxyLadderWidthTrend[nDirs_];
0775     //TH1F *dzLadderWidthTrend[nDirs_];
0776     //TH1F *dzLadderMeanTrend[nDirs_];
0777
0778     //TH1F *dxyModZMeanTrend[nDirs_];
0779     //TH1F *dxyModZWidthTrend[nDirs_];
0780     //TH1F *dzModZMeanTrend[nDirs_];
0781     //TH1F *dzModZWidthTrend[nDirs_];
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     // loop over the objects
0805     for (Int_t j = 0; j < nDirs_; j++) {
0806       //fins[j] = TFile::Open(Form("%s/PVValidation_%s_%i.root",dirs[j],dirs[j],intersection[n]));
0807       size_t position = std::string(dirs[j]).find('/');
0808       std::string stem = std::string(dirs[j]).substr(position + 1);  // get from position to the end
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       // sanity check
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       //dxyLadderMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dxy_ladder");
0850       //dxyLadderWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dxy_ladder");
0851       //dzLadderMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dz_ladder");
0852       //dzLadderWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dz_ladder");
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       //dxyModZMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dxy_modZ");
0860       //dxyModZWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dxy_modZ");
0861       //dzModZMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dz_modZ");
0862       //dzModZWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dz_modZ");
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       // fill the vectors of biases
0878
0879       auto dxyPhiBiases = getBiases(dxyPhiMeanTrend[j]);
0880
0881       //logInfo<<"\n" <<j<<" "<< LegLabels[j] << " dxy(phi) mean: "<< dxyPhiBiases.getWeightedMean()
0882       //       <<" dxy(phi) max: "<< dxyPhiBiases.getMax()
0883       //       <<" dxy(phi) min: "<< dxyPhiBiases.getMin()
0884       //       << std::endl;
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       //logInfo<<"\n" <<j<<" "<< LegLabels[j] << " dxy(phi) ks score: "<< dxyPhiBiases.getKSScore() << std::endl;
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       // unrolled histograms
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       // do all the necessary deletions
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     // do all the necessary deletions
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 }