OfflineValidation/src/PrepareDMRTrends.cc

0001 #include "Alignment/OfflineValidation/interface/PrepareDMRTrends.h"
0002
0003 using namespace std;
0004 namespace fs = std::experimental::filesystem;
0005 namespace pt = boost::property_tree;
0006
0007 PrepareDMRTrends::PrepareDMRTrends(const char *outputFileName, pt::ptree &json) : outputFileName_(outputFileName) {
0008   geometries.clear();
0009   for (const auto &childTree : json) {
0010     geometries.push_back(childTree.second.get<std::string>("title"));
0011   }
0012 }
0013
0014 /*! \fn getName
0015  *  \brief Function used to get a string containing information on the high level structure, the layer/disc and the geometry.
0016  */
0017
0018 TString PrepareDMRTrends::getName(TString structure, int layer, TString geometry) {
0019   geometry.ReplaceAll(" ", "_");
0020   TString name = geometry + "_" + structure;
0021   if (layer != 0) {
0022     if (structure == "TID" || structure == "TEC")
0023       name += "_disc";
0024     else
0025       name += "_layer";
0026     name += layer;
0027   }
0028
0029   return name;
0030 };
0031
0032 /*! \fn compileDMRTrends
0033  *  \brief  Create a file where the DMR trends are stored in the form of TGraph.
0034  */
0035
0036 void PrepareDMRTrends::compileDMRTrends(vector<int> IOVlist,
0037                                         TString Variable,
0038                                         std::vector<std::string> inputFiles,
0039                                         vector<TString> structures,
0040                                         const map<TString, int> nlayers,
0041                                         bool FORCE) {
0042   gROOT->SetBatch();
0043
0044   float ScaleFactor = DMRFactor;
0045   if (Variable == "DrmsNR")
0046     ScaleFactor = 1;
0047
0048   map<pair<pair<TString, int>, TString>, Geometry> mappoints;  // pair = (structure, layer), geometry
0049   Point *point = nullptr;
0050   TFile *f = nullptr;
0051
0052   for (unsigned int i = 0; i < inputFiles.size(); ++i) {
0053     if (fs::is_empty(inputFiles.at(i).c_str())) {
0054       cout << "ERROR: Empty file " << inputFiles.at(i).c_str() << endl;
0055       continue;
0056     }
0057
0058     int runN = IOVlist.at(i);
0059
0060     f = new TFile(inputFiles.at(i).c_str(), "READ");
0061     std::cout << inputFiles.at(i) << std::endl;
0062
0063     for (TString &structure : structures) {
0064       TString structname = structure;
0065       structname.ReplaceAll("_y", "");
0066       size_t layersnumber = nlayers.at(structname);
0067       for (size_t layer = 0; layer <= layersnumber; layer++) {
0068         for (const string &geometry : geometries) {
0069           TString name = Variable + "_" + getName(structure, layer, geometry);
0070           TH1F *histo = dynamic_cast<TH1F *>(f->Get(name));
0071
0072           // Three possibilities:
0073           //  - All histograms are produced correctly
0074           //  - Only the non-split histograms are produced
0075           //  - No histogram is produced correctly
0076           //  FORCE means that the Point is not added to the points collection in the chosen geometry for that structure
0077           //  If FORCE is not enabled a default value for the Point is used (-9999) which will appear in the plots
0078           if (!histo) {
0079             //cout << "Run" << runN << " Histogram: " << name << " not found" << endl;
0080             if (FORCE)
0081               continue;
0082             point = new Point(runN, ScaleFactor);
0083           } else if (structure != "TID" && structure != "TEC") {
0084             TH1F *histoplus = dynamic_cast<TH1F *>(f->Get((name + "_plus")));
0085             TH1F *histominus = dynamic_cast<TH1F *>(f->Get((name + "_minus")));
0086             if (!histoplus || !histominus) {
0087               //cout << "Run" << runN << " Histogram: " << name << " plus or minus not found" << endl;
0088               if (FORCE)
0089                 continue;
0090               point = new Point(runN, ScaleFactor, histo);
0091             } else
0092               point = new Point(runN, ScaleFactor, histo, histoplus, histominus);
0093
0094           } else
0095             point = new Point(runN, ScaleFactor, histo);
0096           mappoints[make_pair(make_pair(structure, layer), geometry)].points.push_back(*point);
0097         }
0098       }
0099     }
0100     f->Close();
0101   }
0102
0103   TFile *fout = TFile::Open(outputFileName_, "RECREATE");
0104   TGraphErrors *g = nullptr;
0105   for (TString &structure : structures) {
0106     TString structname = structure;
0107     structname.ReplaceAll("_y", "");
0108     size_t layersnumber = nlayers.at(structname);
0109     for (size_t layer = 0; layer <= layersnumber; layer++) {
0110       for (const string &geometry : geometries) {
0111         TString name = Variable + "_" + getName(structure, layer, geometry);
0112         Geometry geom = mappoints[make_pair(make_pair(structure, layer), geometry)];
0113         using Trend = vector<float> (Geometry::*)() const;
0114         vector<Trend> trends{&Geometry::Mu,
0115                              &Geometry::Sigma,
0116                              &Geometry::MuPlus,
0117                              &Geometry::SigmaPlus,
0118                              &Geometry::MuMinus,
0119                              &Geometry::SigmaMinus,
0120                              &Geometry::DeltaMu,
0121                              &Geometry::SigmaDeltaMu};
0122         vector<TString> variables{
0123             "mu", "sigma", "muplus", "sigmaplus", "muminus", "sigmaminus", "deltamu", "sigmadeltamu"};
0124         vector<float> runs = geom.Run();
0125         size_t n = runs.size();
0126         vector<float> emptyvec;
0127         for (size_t i = 0; i < runs.size(); i++)
0128           emptyvec.push_back(0.);
0129         for (size_t iVar = 0; iVar < variables.size(); iVar++) {
0130           Trend trend = trends.at(iVar);
0131           g = new TGraphErrors(n, runs.data(), (geom.*trend)().data(), emptyvec.data(), emptyvec.data());
0132           g->SetTitle(geometry.c_str());
0133           g->Write(name + "_" + variables.at(iVar));
0134         }
0135         vector<pair<Trend, Trend>> trendspair{make_pair(&Geometry::Mu, &Geometry::Sigma),
0136                                               make_pair(&Geometry::MuPlus, &Geometry::SigmaPlus),
0137                                               make_pair(&Geometry::MuMinus, &Geometry::SigmaMinus),
0138                                               make_pair(&Geometry::DeltaMu, &Geometry::SigmaDeltaMu)};
0139         vector<pair<TString, TString>> variablepairs{make_pair("mu", "sigma"),
0140                                                      make_pair("muplus", "sigmaplus"),
0141                                                      make_pair("muminus", "sigmaminus"),
0142                                                      make_pair("deltamu", "sigmadeltamu")};
0143         for (size_t iVar = 0; iVar < variablepairs.size(); iVar++) {
0144           Trend meantrend = trendspair.at(iVar).first;
0145           Trend sigmatrend = trendspair.at(iVar).second;
0146           g = new TGraphErrors(
0147               n, runs.data(), (geom.*meantrend)().data(), emptyvec.data(), (geom.*sigmatrend)().data());
0148           g->SetTitle(geometry.c_str());
0149           TString graphname = name + "_" + variablepairs.at(iVar).first;
0150           graphname += variablepairs.at(iVar).second;
0151           g->Write(graphname);
0152         }
0153       }
0154     }
0155   }
0156   fout->Close();
0157
0158   delete point;
0159   delete f;
0160   delete g;
0161 }