Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-12-04 23:57:30

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   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 /*! \fn list_files
0572  *  \brief utility function to list of filles in a directory
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         //logInfo << fname.Data() << std::endl;
0591         TObjArray *bits = fname.Tokenize("_");
0592         TString theRun = bits->At(2)->GetName();
0593         //logInfo << theRun << std::endl;
0594         TString formatRun = (theRun.ReplaceAll(".root", "")).ReplaceAll("_", "");
0595         //logInfo << dirname << " "<< formatRun.Atoi() << std::endl;
0596         theRunNumbers.push_back(formatRun.Atoi());
0597       }
0598     }
0599   }
0600   return theRunNumbers;
0601 }
0602 
0603 /*! \fn DrawConstant
0604  *  \brief utility function to draw a constant histogram with erros !=0
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 /*! \fn getUnrolledHisto
0632  *  \brief utility function to tranform a TH1 into a vector of floats
0633  */
0634 
0635 /*--------------------------------------------------------------------*/
0636 unrolledHisto PreparePVTrends::getUnrolledHisto(TH1F *hist)
0637 /*--------------------------------------------------------------------*/
0638 {
0639   /*
0640     Double_t y_min = hist->GetBinLowEdge(1);
0641     Double_t y_max = hist->GetBinLowEdge(hist->GetNbinsX()+1);
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 /*! \fn getBiases
0658  *  \brief utility function to extract characterization of the PV bias plot
0659  */
0660 
0661 /*--------------------------------------------------------------------*/
0662 pv::biases PreparePVTrends::getBiases(TH1F *hist)
0663 /*--------------------------------------------------------------------*/
0664 {
0665   int nbins = hist->GetNbinsX();
0666   // if there are no bins in the histogram then return default constructed object
0667   // shouldn't really ever happen
0668   if (nbins <= 0) {
0669     logError << "No bins in the input histogram";
0670     return pv::biases();
0671   }
0672 
0673   //extract median from histogram
0674   double *y = new double[nbins];
0675   double *err = new double[nbins];
0676 
0677   // remember for weight means <x> = sum_i (x_i* w_i) / sum_i w_i ; where w_i = 1/sigma^2_i
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   // in case one would like to use a pol0 fit
0698   hist->Fit("pol0", "Q0+");
0699   TF1 *f = (TF1 *)hist->FindObject("pol0");
0700   //f->SetLineColor(hist->GetLineColor());
0701   //f->SetLineStyle(hist->GetLineStyle());
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     std::pair<std::pair<Double_t,Double_t>, Double_t> result;
0712     std::pair<Double_t,Double_t> resultBounds;
0713     resultBounds = useRMS_ ? std::make_pair(mean-rms,mean+rms) :  std::make_pair(min,max)  ;
0714     result = make_pair(resultBounds,mean);
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 /*! \fn processData
0727  *  \brief function where the magic happens, take the raw inputs and creates the output Trends
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     //in case of debug, use only 50
0755     //for(unsigned int n=0; n<50;n++){
0756 
0757     //if(intersection.at(n)!=283946)
0758     //  continue;
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     //TH1F *dxyLadderMeanTrend[nDirs_];
0772     //TH1F *dxyLadderWidthTrend[nDirs_];
0773     //TH1F *dzLadderWidthTrend[nDirs_];
0774     //TH1F *dzLadderMeanTrend[nDirs_];
0775 
0776     //TH1F *dxyModZMeanTrend[nDirs_];
0777     //TH1F *dxyModZWidthTrend[nDirs_];
0778     //TH1F *dzModZMeanTrend[nDirs_];
0779     //TH1F *dzModZWidthTrend[nDirs_];
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     // loop over the objects
0803     for (Int_t j = 0; j < nDirs_; j++) {
0804       //fins[j] = TFile::Open(Form("%s/PVValidation_%s_%i.root",dirs[j],dirs[j],intersection[n]));
0805       size_t position = std::string(dirs[j]).find('/');
0806       std::string stem = std::string(dirs[j]).substr(position + 1);  // get from position to the end
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       // sanity check
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       //dxyLadderMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dxy_ladder");
0848       //dxyLadderWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dxy_ladder");
0849       //dzLadderMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dz_ladder");
0850       //dzLadderWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dz_ladder");
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       //dxyModZMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dxy_modZ");
0858       //dxyModZWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dxy_modZ");
0859       //dzModZMeanTrend[j] = (TH1F *)fins[j]->Get("PVValidation/MeanTrends/means_dz_modZ");
0860       //dzModZWidthTrend[j] = (TH1F *)fins[j]->Get("PVValidation/WidthTrends/widths_dz_modZ");
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       // fill the vectors of biases
0876 
0877       auto dxyPhiBiases = getBiases(dxyPhiMeanTrend[j]);
0878 
0879       //logInfo<<"\n" <<j<<" "<< LegLabels[j] << " dxy(phi) mean: "<< dxyPhiBiases.getWeightedMean()
0880       //       <<" dxy(phi) max: "<< dxyPhiBiases.getMax()
0881       //       <<" dxy(phi) min: "<< dxyPhiBiases.getMin()
0882       //       << std::endl;
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       //logInfo<<"\n" <<j<<" "<< LegLabels[j] << " dxy(phi) ks score: "<< dxyPhiBiases.getKSScore() << std::endl;
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       // unrolled histograms
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       // do all the necessary deletions
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     // do all the necessary deletions
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 }