diff --git a/Alignment/OfflineValidation/bin/BuildFile.xml b/Alignment/OfflineValidation/bin/BuildFile.xml index b8d68c8419dcf..4034992bc0c58 100644 --- a/Alignment/OfflineValidation/bin/BuildFile.xml +++ b/Alignment/OfflineValidation/bin/BuildFile.xml @@ -12,6 +12,7 @@ + diff --git a/Alignment/OfflineValidation/bin/GCPtrends.cc b/Alignment/OfflineValidation/bin/GCPtrends.cc new file mode 100644 index 0000000000000..8aa0f657b2b77 --- /dev/null +++ b/Alignment/OfflineValidation/bin/GCPtrends.cc @@ -0,0 +1,275 @@ +#include +#include +#include + +#include +#include + +#include "TFile.h" +#include "TChain.h" +#include "TSelector.h" +#include "TTreeReader.h" +#include "TTreeReaderValue.h" +#include "TTreeReaderArray.h" +#include "TGraphErrors.h" +#include "TH1D.h" +#include "TObject.h" +#include "TMath.h" +#include "TString.h" + +#include "Alignment/OfflineValidation/interface/Trend.h" + +int main(int argc, char *argv[]) { + if (argc < 2) { + std::cout << "Usage of the program : GCPtrends [gcp_trend_config.json]" << std::endl; + std::cout << "gcp_trend_config.json : file with configuration in json format" << std::endl; + exit(1); + } + + // parsing input file + std::string lumi_file(argv[1]); + boost::property_tree::ptree config_root; + boost::property_tree::read_json(lumi_file, config_root); + + // configuration + TString outputdir = config_root.get("outputdir", "."); + outputdir = TString(outputdir) + TString(outputdir.EndsWith("/") ? "" : "/"); + + bool usebadmodules = config_root.get("usebadmodules", true); + + bool splitpixel = config_root.get("splitpixel", false); + + Int_t trackerphase = config_root.get("trackerphase"); + + TString lumitype = config_root.get("trends.lumitype"); + + Int_t firstrun = config_root.get("trends.firstrun"); + + Int_t lastrun = config_root.get("trends.lastrun"); + + std::string lumibyrunfile = config_root.get("lumibyrunfile", ""); + std::ifstream lumifile(lumibyrunfile); + assert(lumifile.good()); + const Run2Lumi lumi_per_run(lumibyrunfile, firstrun, lastrun, 1000); + + std::map comparison_map; + for (boost::property_tree::ptree::value_type &comparison : config_root.get_child("comparisons")) { + TString path = comparison.first; + Int_t run_label = comparison.second.get_value(); + + comparison_map[path] = run_label; + } + + // create output file + TFile *file_out = TFile::Open(outputdir + "GCPTrends.root", "RECREATE"); + file_out->cd(); + + gErrorIgnoreLevel = kError; + + // sub-detector mapping + std::map Sublevel_Subdetector_Map = { + {1, "PXB"}, {2, "PXF"}, {3, "TIB"}, {4, "TID"}, {5, "TOB"}, {6, "TEC"}}; + + // wheels/layers per tracker phase: > + const std::map> Phase_Subdetector_Layers_Map = { + {0, {{1, 3}, {2, 2}}}, {1, {{1, 4}, {2, 3}}}, {2, {{1, 4}, {2, 12}}}}; + + // adding layers/wheels in case Pixel is requested further split + if (splitpixel) { + assert(trackerphase < 3); + + for (auto &sub_layer : Phase_Subdetector_Layers_Map.at(trackerphase)) { + for (int iLayer = 1; iLayer <= sub_layer.second; iLayer++) { + // subid*100+subsubid + if (sub_layer.first % 2 != 0) { + Sublevel_Subdetector_Map[100 * sub_layer.first + iLayer] = + Sublevel_Subdetector_Map[sub_layer.first] + "Layer" + TString(std::to_string(iLayer)); + + } else { + Sublevel_Subdetector_Map[100 * sub_layer.first + (1 - 1) * sub_layer.second + iLayer] = + Sublevel_Subdetector_Map[sub_layer.first] + "Wheel" + TString(std::to_string(iLayer)) + "Side1"; + Sublevel_Subdetector_Map[100 * sub_layer.first + (2 - 1) * sub_layer.second + iLayer] = + Sublevel_Subdetector_Map[sub_layer.first] + "Wheel" + TString(std::to_string(iLayer)) + "Side2"; + } + } + } + } + + // variable mapping + const std::map Index_Variable_Map = { + {0, "DX"}, {1, "DY"}, {2, "DZ"}, {3, "DAlpha"}, {4, "DBeta"}, {5, "DGamma"}}; + const std::map Variable_Name_Map = {{"DX", "#Delta x"}, + {"DY", "#Delta y"}, + {"DZ", "#Delta z"}, + {"DAlpha", "#Delta #alpha"}, + {"DBeta", "#Delta #beta"}, + {"DGamma", "#Delta #gamma"}}; + // estimator mapping + const std::map Index_Estimator_Map = {{0, "Mean"}, {1, "Sigma"}}; + const std::map Estimator_Name_Map = {{"Mean", "#mu"}, {"Sigma", "#sigma"}}; + + // constant unit conversion + const int convert_cm_to_mum = 10000; + const int convert_rad_to_mrad = 1000; + + std::cout << std::endl; + std::cout << " ==> " << comparison_map.size() << " GCPs to be processed in configuration file ... " << std::endl; + std::cout << std::endl; + + // booking TGraphs and TH1D + TH1::StatOverflows(kTRUE); + std::map>> Graphs; + std::map> Histos; + for (auto &level : Sublevel_Subdetector_Map) { + for (auto &variable : Index_Variable_Map) { + Histos[level.first][variable.first] = new TH1D( + "Histo_" + level.second + "_" + variable.second, "Histo_" + level.second + "_" + variable.second, 1, -1, 1); + + for (auto &estimator : Index_Estimator_Map) { + Graphs[level.first][variable.first][estimator.first] = new TGraphErrors(); + } + } + } + + // loop over the comparisons (GCPs) + for (auto &path_run_map : comparison_map) { + TChain Events("alignTree", "alignTree"); + Events.Add(path_run_map.first); + + Int_t run_number = path_run_map.second; + + std::cout << " --> processing GCPtree file: " << path_run_map.first << std::endl; + + TTreeReader Reader(&Events); + Reader.Restart(); + + TTreeReaderValue id = {Reader, "id"}; + TTreeReaderValue badModuleQuality = {Reader, "badModuleQuality"}; + TTreeReaderValue inModuleList = {Reader, "inModuleList"}; + TTreeReaderValue level = {Reader, "level"}; + TTreeReaderValue mid = {Reader, "mid"}; + TTreeReaderValue mlevel = {Reader, "mlevel"}; + TTreeReaderValue sublevel = {Reader, "sublevel"}; + TTreeReaderValue x = {Reader, "x"}; + TTreeReaderValue y = {Reader, "y"}; + TTreeReaderValue z = {Reader, "z"}; + TTreeReaderValue r = {Reader, "r"}; + TTreeReaderValue phi = {Reader, "phi"}; + TTreeReaderValue eta = {Reader, "eta"}; + TTreeReaderValue alpha = {Reader, "alpha"}; + TTreeReaderValue beta = {Reader, "beta"}; + TTreeReaderValue gamma = {Reader, "gamma"}; + TTreeReaderValue dx = {Reader, "dx"}; + TTreeReaderValue dy = {Reader, "dy"}; + TTreeReaderValue dz = {Reader, "dz"}; + TTreeReaderValue dr = {Reader, "dr"}; + TTreeReaderValue dphi = {Reader, "dphi"}; + TTreeReaderValue dalpha = {Reader, "dalpha"}; + TTreeReaderValue dbeta = {Reader, "dbeta"}; + TTreeReaderValue dgamma = {Reader, "dgamma"}; + TTreeReaderValue du = {Reader, "du"}; + TTreeReaderValue dv = {Reader, "dv"}; + TTreeReaderValue dw = {Reader, "dw"}; + TTreeReaderValue da = {Reader, "da"}; + TTreeReaderValue db = {Reader, "db"}; + TTreeReaderValue dg = {Reader, "dg"}; + TTreeReaderValue useDetId = {Reader, "useDetId"}; + TTreeReaderValue detDim = {Reader, "detDim"}; + TTreeReaderArray identifiers = {Reader, "identifiers"}; + + // loop over modules + while (Reader.Next()) { + if (!usebadmodules) + if (*badModuleQuality.Get()) + continue; + + Int_t sublevel_idx = *sublevel.Get(); + + Histos[sublevel_idx][0]->Fill(*dx.Get() * convert_cm_to_mum); + Histos[sublevel_idx][1]->Fill(*dy.Get() * convert_cm_to_mum); + Histos[sublevel_idx][2]->Fill(*dz.Get() * convert_cm_to_mum); + Histos[sublevel_idx][3]->Fill(*dalpha.Get() * convert_rad_to_mrad); + Histos[sublevel_idx][4]->Fill(*dbeta.Get() * convert_rad_to_mrad); + Histos[sublevel_idx][5]->Fill(*dgamma.Get() * convert_rad_to_mrad); + + if (splitpixel && sublevel_idx <= 2) { + Int_t layer_index; + + if (sublevel_idx % 2 != 0) + layer_index = 100 * sublevel_idx + identifiers[2]; + else + layer_index = 100 * sublevel_idx + + (identifiers[4] - 1) * Phase_Subdetector_Layers_Map.at(trackerphase).at(sublevel_idx) + + identifiers[3]; + + Histos[layer_index][0]->Fill(*dx.Get() * convert_cm_to_mum); + Histos[layer_index][1]->Fill(*dy.Get() * convert_cm_to_mum); + Histos[layer_index][2]->Fill(*dz.Get() * convert_cm_to_mum); + Histos[layer_index][3]->Fill(*dalpha.Get() * convert_rad_to_mrad); + Histos[layer_index][4]->Fill(*dbeta.Get() * convert_rad_to_mrad); + Histos[layer_index][5]->Fill(*dgamma.Get() * convert_rad_to_mrad); + } + } + + for (auto &level : Sublevel_Subdetector_Map) { + for (auto &variable : Index_Variable_Map) { + Double_t mean = Histos[level.first][variable.first]->GetMean(); + Double_t sigma = Histos[level.first][variable.first]->GetStdDev(); + + Graphs[level.first][variable.first][0]->AddPoint(run_number, mean); + if (std::fabs(mean) > Graphs[level.first][variable.first][0]->GetMaximum() && std::fabs(mean) > 0.) { + Graphs[level.first][variable.first][0]->SetMaximum(std::fabs(mean)); + Graphs[level.first][variable.first][0]->SetMinimum(-std::fabs(mean)); + } + + Graphs[level.first][variable.first][1]->AddPoint(run_number, sigma); + if (sigma > Graphs[level.first][variable.first][1]->GetMaximum() && sigma > 0.) { + Graphs[level.first][variable.first][1]->SetMaximum(sigma); + Graphs[level.first][variable.first][1]->SetMinimum(0.); + } + + Histos[level.first][variable.first]->Reset("ICESM"); + } + } + } + + // saving TGraphs + for (auto &level : Sublevel_Subdetector_Map) { + for (auto &variable : Index_Variable_Map) { + for (auto &estimator : Index_Estimator_Map) { + Graphs[level.first][variable.first][estimator.first]->Write("Graph_" + level.second + "_" + variable.second + + "_" + estimator.second); + + TString units = "mrad"; + if (variable.second.Contains("DX") || variable.second.Contains("DY") || variable.second.Contains("DZ")) + units = "#mum"; + + Trend trend("Graph_" + level.second + "_" + variable.second + "_" + estimator.second, + "output", + "Trend", + Estimator_Name_Map.at(estimator.second) + "_{" + Variable_Name_Map.at(variable.second) + "} [" + + units + "]", + -2 * std::fabs(Graphs[level.first][variable.first][estimator.first]->GetMinimum()), + 2 * std::fabs(Graphs[level.first][variable.first][estimator.first]->GetMaximum()), + config_root, + lumi_per_run, + lumitype); + + Graphs[level.first][variable.first][estimator.first]->SetFillColor(4); + + trend.lgd.SetHeader(level.second, "center"); + trend.lgd.AddEntry(Graphs[level.first][variable.first][estimator.first], "Average over all modules", "pl"); + + trend(Graphs[level.first][variable.first][estimator.first], "p", "pf", false); + } + } + } + + file_out->Close(); + + std::cout << std::endl; + std::cout << " ==> done." << std::endl; + std::cout << std::endl; + + return 0; +} diff --git a/Alignment/OfflineValidation/src/Trend.cc b/Alignment/OfflineValidation/src/Trend.cc index bafeb4bf71800..30c85d38d0a6f 100644 --- a/Alignment/OfflineValidation/src/Trend.cc +++ b/Alignment/OfflineValidation/src/Trend.cc @@ -22,10 +22,8 @@ namespace pt = boost::property_tree; Run2Lumi::Run2Lumi(fs::path file, int first, int last, float convertUnit) : firstRun(first), lastRun(last), convertUnit(convertUnit) { - cout << __func__ << endl; assert(first < last); - cout << file << endl; assert(fs::exists(file)); ifstream f(file); @@ -143,8 +141,6 @@ Trend::Trend(const char* name, JSON(json), GetLumi(GetLumiFunctor), lumiType(lumiAxisType) { - cout << __func__ << endl; - if (JSON.count("CMSlabel")) CMS = Form("#scale[1.1]{#bf{CMS}} #it{%s}", JSON.get("CMSlabel").data()); @@ -164,11 +160,8 @@ Trend::Trend(const char* name, frame->GetYaxis()->SetLabelSize(fontsize); frame->GetYaxis()->SetTitleSize(fontsize); lgd.SetTextSize(fontsize); - cout << "frame->GetXaxis()->GetLabelSize() = " << frame->GetXaxis()->GetLabelSize() << endl; - cout << "frame->GetXaxis()->GetTitleSize() = " << frame->GetXaxis()->GetTitleSize() << endl; if (ymax > 0 && ymin < 0) { - cout << "Plotting horizontal line at zero" << endl; TLine l; l.SetLineColor(kBlack); l.SetLineStyle(kDashed); @@ -210,7 +203,7 @@ Trend::Trend(const char* name, auto currentRun = run.second.get_value(); auto lumi = GetLumi(GetLumi.firstRun, currentRun); - //cout << currentRun << '\t' << lumi << endl; + if (lumi > 0) v->DrawLine(lumi, ymin, lumi, ymax); } @@ -218,20 +211,17 @@ Trend::Trend(const char* name, } void Trend::operator()(TObject* obj, TString drawOpt, TString lgdOpt, bool fullRange) { - cout << __func__ << endl; c.cd(); TString classname = obj->ClassName(); if (classname.Contains("TGraph")) { auto g = dynamic_cast(obj); int n = g->GetN(); - cout << g->GetPointX(n - 1) << ' ' << GetLumi.lastRun; + if (fullRange) { - cout << " -> adding one point" << endl; g->Set(n); g->SetPoint(n, GetLumi.lastRun, 0); - } else - cout << " -> hole between end of graph and right edge" << endl; + } g = GetLumi(g); g->Draw("same" + drawOpt); } else if (classname.Contains("TH1")) { @@ -255,8 +245,6 @@ void Trend::operator()(TObject* obj, TString drawOpt, TString lgdOpt, bool fullR } Trend::~Trend() { - cout << __func__ << endl; - c.cd(); lgd.Draw(); @@ -292,7 +280,7 @@ Trend::~Trend() { auto lumi = max(GetLumi(currentRun), (float)0.01); auto posX = l + (lumi / totLumi) / (l + 1 + r) + 0.02; - cout << currentRun << setw(20) << label << setw(20) << lumi << setw(20) << posX << endl; + label = "#scale[0.8]{" + label + "}"; latex.DrawLatex(posX, posY, label.c_str());