File indexing completed on 2023-03-17 10:39:57
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
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
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
0064 TFile *file_out = TFile::Open(outputdir + "GCPTrends.root", "RECREATE");
0065 file_out->cd();
0066
0067 gErrorIgnoreLevel = kError;
0068
0069
0070 std::map<Int_t, TString> Sublevel_Subdetector_Map = {
0071 {1, "PXB"}, {2, "PXF"}, {3, "TIB"}, {4, "TID"}, {5, "TOB"}, {6, "TEC"}};
0072
0073
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
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
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
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
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
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
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
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
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
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 }