OfflineValidation/bin/GCPtrends.cc

0001 #include <iostream>
0002 #include <fstream>
0003 #include <sstream>
0004
0005 #include <boost/property_tree/ptree.hpp>
0006 #include <boost/property_tree/json_parser.hpp>
0007
0008 #include "TFile.h"
0009 #include "TChain.h"
0010 #include "TSelector.h"
0011 #include "TTreeReader.h"
0012 #include "TTreeReaderValue.h"
0013 #include "TTreeReaderArray.h"
0014 #include "TGraphErrors.h"
0015 #include "TH1D.h"
0016 #include "TObject.h"
0017 #include "TMath.h"
0018 #include "TString.h"
0019
0020 #include "Alignment/OfflineValidation/interface/Trend.h"
0021
0022 int main(int argc, char *argv[]) {
0023   if (argc < 2) {
0024     std::cout << "Usage of the program : GCPtrends [gcp_trend_config.json]" << std::endl;
0025     std::cout << "gcp_trend_config.json : file with configuration in json format" << std::endl;
0026     exit(1);
0027   }
0028
0029   // parsing input file
0030   std::string lumi_file(argv[1]);
0031   boost::property_tree::ptree config_root;
0032   boost::property_tree::read_json(lumi_file, config_root);
0033
0034   // configuration
0035   TString outputdir = config_root.get<std::string>("outputdir", ".");
0036   outputdir = TString(outputdir) + TString(outputdir.EndsWith("/") ? "" : "/");
0037
0038   bool usebadmodules = config_root.get<bool>("usebadmodules", true);
0039
0040   bool splitpixel = config_root.get<bool>("splitpixel", false);
0041
0042   Int_t trackerphase = config_root.get<int>("trackerphase");
0043
0044   TString lumitype = config_root.get<std::string>("trends.lumitype");
0045
0046   Int_t firstrun = config_root.get<int>("trends.firstrun");
0047
0048   Int_t lastrun = config_root.get<int>("trends.lastrun");
0049
0050   std::string lumibyrunfile = config_root.get<std::string>("lumibyrunfile", "");
0051   std::ifstream lumifile(lumibyrunfile);
0052   assert(lumifile.good());
0053   const Run2Lumi lumi_per_run(lumibyrunfile, firstrun, lastrun, 1000);
0054
0055   std::map<TString, Int_t> comparison_map;
0056   for (boost::property_tree::ptree::value_type &comparison : config_root.get_child("comparisons")) {
0057     TString path = comparison.first;
0058     Int_t run_label = comparison.second.get_value<int>();
0059
0060     comparison_map[path] = run_label;
0061   }
0062
0063   // create output file
0064   TFile *file_out = TFile::Open(outputdir + "GCPTrends.root", "RECREATE");
0065   file_out->cd();
0066
0067   gErrorIgnoreLevel = kError;
0068
0069   // sub-detector mapping
0070   std::map<Int_t, TString> Sublevel_Subdetector_Map = {
0071       {1, "PXB"}, {2, "PXF"}, {3, "TIB"}, {4, "TID"}, {5, "TOB"}, {6, "TEC"}};
0072
0073   // wheels/layers per tracker phase: <pahse,<sublevel,number of layers/wheels>>
0074   const std::map<Int_t, std::map<Int_t, Int_t>> Phase_Subdetector_Layers_Map = {
0075       {0, {{1, 3}, {2, 2}}}, {1, {{1, 4}, {2, 3}}}, {2, {{1, 4}, {2, 12}}}};
0076
0077   // adding layers/wheels in case Pixel is requested further split
0078   if (splitpixel) {
0079     assert(trackerphase < 3);
0080
0081     for (auto &sub_layer : Phase_Subdetector_Layers_Map.at(trackerphase)) {
0082       for (int iLayer = 1; iLayer <= sub_layer.second; iLayer++) {
0083         // subid*100+subsubid
0084         if (sub_layer.first % 2 != 0) {
0085           Sublevel_Subdetector_Map[100 * sub_layer.first + iLayer] =
0086               Sublevel_Subdetector_Map[sub_layer.first] + "Layer" + TString(std::to_string(iLayer));
0087
0088         } else {
0089           Sublevel_Subdetector_Map[100 * sub_layer.first + (1 - 1) * sub_layer.second + iLayer] =
0090               Sublevel_Subdetector_Map[sub_layer.first] + "Wheel" + TString(std::to_string(iLayer)) + "Side1";
0091           Sublevel_Subdetector_Map[100 * sub_layer.first + (2 - 1) * sub_layer.second + iLayer] =
0092               Sublevel_Subdetector_Map[sub_layer.first] + "Wheel" + TString(std::to_string(iLayer)) + "Side2";
0093         }
0094       }
0095     }
0096   }
0097
0098   // variable mapping
0099   const std::map<Int_t, TString> Index_Variable_Map = {
0100       {0, "DX"}, {1, "DY"}, {2, "DZ"}, {3, "DAlpha"}, {4, "DBeta"}, {5, "DGamma"}};
0101   const std::map<TString, TString> Variable_Name_Map = {{"DX", "#Delta x"},
0102                                                         {"DY", "#Delta y"},
0103                                                         {"DZ", "#Delta z"},
0104                                                         {"DAlpha", "#Delta #alpha"},
0105                                                         {"DBeta", "#Delta #beta"},
0106                                                         {"DGamma", "#Delta #gamma"}};
0107   // estimator mapping
0108   const std::map<Int_t, TString> Index_Estimator_Map = {{0, "Mean"}, {1, "Sigma"}};
0109   const std::map<TString, TString> Estimator_Name_Map = {{"Mean", "#mu"}, {"Sigma", "#sigma"}};
0110
0111   // constant unit conversion
0112   const int convert_cm_to_mum = 10000;
0113   const int convert_rad_to_mrad = 1000;
0114
0115   std::cout << std::endl;
0116   std::cout << "   ==> " << comparison_map.size() << " GCPs to be processed in configuration file ... " << std::endl;
0117   std::cout << std::endl;
0118
0119   // booking TGraphs and TH1D
0120   TH1::StatOverflows(kTRUE);
0121   std::map<Int_t, std::map<Int_t, std::map<Int_t, TGraphErrors *>>> Graphs;
0122   std::map<Int_t, std::map<Int_t, TH1D *>> Histos;
0123   for (auto &level : Sublevel_Subdetector_Map) {
0124     for (auto &variable : Index_Variable_Map) {
0125       Histos[level.first][variable.first] = new TH1D(
0126           "Histo_" + level.second + "_" + variable.second, "Histo_" + level.second + "_" + variable.second, 1, -1, 1);
0127
0128       for (auto &estimator : Index_Estimator_Map) {
0129         Graphs[level.first][variable.first][estimator.first] = new TGraphErrors();
0130       }
0131     }
0132   }
0133
0134   // loop over the comparisons (GCPs)
0135   for (auto &path_run_map : comparison_map) {
0136     TChain Events("alignTree", "alignTree");
0137     Events.Add(path_run_map.first);
0138
0139     Int_t run_number = path_run_map.second;
0140
0141     std::cout << "       --> processing GCPtree file: " << path_run_map.first << std::endl;
0142
0143     TTreeReader Reader(&Events);
0144     Reader.Restart();
0145
0146     TTreeReaderValue<Int_t> id = {Reader, "id"};
0147     TTreeReaderValue<Int_t> badModuleQuality = {Reader, "badModuleQuality"};
0148     TTreeReaderValue<Int_t> inModuleList = {Reader, "inModuleList"};
0149     TTreeReaderValue<Int_t> level = {Reader, "level"};
0150     TTreeReaderValue<Int_t> mid = {Reader, "mid"};
0151     TTreeReaderValue<Int_t> mlevel = {Reader, "mlevel"};
0152     TTreeReaderValue<Int_t> sublevel = {Reader, "sublevel"};
0153     TTreeReaderValue<Float_t> x = {Reader, "x"};
0154     TTreeReaderValue<Float_t> y = {Reader, "y"};
0155     TTreeReaderValue<Float_t> z = {Reader, "z"};
0156     TTreeReaderValue<Float_t> r = {Reader, "r"};
0157     TTreeReaderValue<Float_t> phi = {Reader, "phi"};
0158     TTreeReaderValue<Float_t> eta = {Reader, "eta"};
0159     TTreeReaderValue<Float_t> alpha = {Reader, "alpha"};
0160     TTreeReaderValue<Float_t> beta = {Reader, "beta"};
0161     TTreeReaderValue<Float_t> gamma = {Reader, "gamma"};
0162     TTreeReaderValue<Float_t> dx = {Reader, "dx"};
0163     TTreeReaderValue<Float_t> dy = {Reader, "dy"};
0164     TTreeReaderValue<Float_t> dz = {Reader, "dz"};
0165     TTreeReaderValue<Float_t> dr = {Reader, "dr"};
0166     TTreeReaderValue<Float_t> dphi = {Reader, "dphi"};
0167     TTreeReaderValue<Float_t> dalpha = {Reader, "dalpha"};
0168     TTreeReaderValue<Float_t> dbeta = {Reader, "dbeta"};
0169     TTreeReaderValue<Float_t> dgamma = {Reader, "dgamma"};
0170     TTreeReaderValue<Float_t> du = {Reader, "du"};
0171     TTreeReaderValue<Float_t> dv = {Reader, "dv"};
0172     TTreeReaderValue<Float_t> dw = {Reader, "dw"};
0173     TTreeReaderValue<Float_t> da = {Reader, "da"};
0174     TTreeReaderValue<Float_t> db = {Reader, "db"};
0175     TTreeReaderValue<Float_t> dg = {Reader, "dg"};
0176     TTreeReaderValue<Int_t> useDetId = {Reader, "useDetId"};
0177     TTreeReaderValue<Int_t> detDim = {Reader, "detDim"};
0178     TTreeReaderArray<Int_t> identifiers = {Reader, "identifiers"};
0179
0180     // loop over modules
0181     while (Reader.Next()) {
0182       if (!usebadmodules)
0183         if (*badModuleQuality.Get())
0184           continue;
0185
0186       Int_t sublevel_idx = *sublevel.Get();
0187
0188       Histos[sublevel_idx][0]->Fill(*dx.Get() * convert_cm_to_mum);
0189       Histos[sublevel_idx][1]->Fill(*dy.Get() * convert_cm_to_mum);
0190       Histos[sublevel_idx][2]->Fill(*dz.Get() * convert_cm_to_mum);
0191       Histos[sublevel_idx][3]->Fill(*dalpha.Get() * convert_rad_to_mrad);
0192       Histos[sublevel_idx][4]->Fill(*dbeta.Get() * convert_rad_to_mrad);
0193       Histos[sublevel_idx][5]->Fill(*dgamma.Get() * convert_rad_to_mrad);
0194
0195       if (splitpixel && sublevel_idx <= 2) {
0196         Int_t layer_index;
0197
0198         if (sublevel_idx % 2 != 0)
0199           layer_index = 100 * sublevel_idx + identifiers[2];
0200         else
0201           layer_index = 100 * sublevel_idx +
0202                         (identifiers[4] - 1) * Phase_Subdetector_Layers_Map.at(trackerphase).at(sublevel_idx) +
0203                         identifiers[3];
0204
0205         Histos[layer_index][0]->Fill(*dx.Get() * convert_cm_to_mum);
0206         Histos[layer_index][1]->Fill(*dy.Get() * convert_cm_to_mum);
0207         Histos[layer_index][2]->Fill(*dz.Get() * convert_cm_to_mum);
0208         Histos[layer_index][3]->Fill(*dalpha.Get() * convert_rad_to_mrad);
0209         Histos[layer_index][4]->Fill(*dbeta.Get() * convert_rad_to_mrad);
0210         Histos[layer_index][5]->Fill(*dgamma.Get() * convert_rad_to_mrad);
0211       }
0212     }
0213
0214     for (auto &level : Sublevel_Subdetector_Map) {
0215       for (auto &variable : Index_Variable_Map) {
0216         Double_t mean = Histos[level.first][variable.first]->GetMean();
0217         Double_t sigma = Histos[level.first][variable.first]->GetStdDev();
0218
0219         Graphs[level.first][variable.first][0]->AddPoint(run_number, mean);
0220         if (std::fabs(mean) > Graphs[level.first][variable.first][0]->GetMaximum() && std::fabs(mean) > 0.) {
0221           Graphs[level.first][variable.first][0]->SetMaximum(std::fabs(mean));
0222           Graphs[level.first][variable.first][0]->SetMinimum(-std::fabs(mean));
0223         }
0224
0225         Graphs[level.first][variable.first][1]->AddPoint(run_number, sigma);
0226         if (sigma > Graphs[level.first][variable.first][1]->GetMaximum() && sigma > 0.) {
0227           Graphs[level.first][variable.first][1]->SetMaximum(sigma);
0228           Graphs[level.first][variable.first][1]->SetMinimum(0.);
0229         }
0230
0231         Histos[level.first][variable.first]->Reset("ICESM");
0232       }
0233     }
0234   }
0235
0236   // saving TGraphs
0237   for (auto &level : Sublevel_Subdetector_Map) {
0238     for (auto &variable : Index_Variable_Map) {
0239       for (auto &estimator : Index_Estimator_Map) {
0240         Graphs[level.first][variable.first][estimator.first]->Write("Graph_" + level.second + "_" + variable.second +
0241                                                                     "_" + estimator.second);
0242
0243         TString units = "mrad";
0244         if (variable.second.Contains("DX") || variable.second.Contains("DY") || variable.second.Contains("DZ"))
0245           units = "#mum";
0246
0247         Trend trend("Graph_" + level.second + "_" + variable.second + "_" + estimator.second,
0248                     "output",
0249                     "Trend",
0250                     Estimator_Name_Map.at(estimator.second) + "_{" + Variable_Name_Map.at(variable.second) + "} [" +
0251                         units + "]",
0252                     -2 * std::fabs(Graphs[level.first][variable.first][estimator.first]->GetMinimum()),
0253                     2 * std::fabs(Graphs[level.first][variable.first][estimator.first]->GetMaximum()),
0254                     config_root,
0255                     lumi_per_run,
0256                     lumitype);
0257
0258         Graphs[level.first][variable.first][estimator.first]->SetFillColor(4);
0259
0260         trend.lgd.SetHeader(level.second, "center");
0261         trend.lgd.AddEntry(Graphs[level.first][variable.first][estimator.first], "Average over all modules", "pl");
0262
0263         trend(Graphs[level.first][variable.first][estimator.first], "p", "pf", false);
0264       }
0265     }
0266   }
0267
0268   file_out->Close();
0269
0270   std::cout << std::endl;
0271   std::cout << "   ==> done." << std::endl;
0272   std::cout << std::endl;
0273
0274   return 0;
0275 }