Adding new GCP trends tool and commenting unnecessary print out in Trend #40336
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,275 @@ | ||
| #include <iostream> | ||
| #include <fstream> | ||
| #include <sstream> | ||
|
|
||
| #include <boost/property_tree/ptree.hpp> | ||
| #include <boost/property_tree/json_parser.hpp> | ||
|
|
||
| #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<std::string>("outputdir", "."); | ||
| outputdir = TString(outputdir) + TString(outputdir.EndsWith("/") ? "" : "/"); | ||
|
|
||
| bool usebadmodules = config_root.get<bool>("usebadmodules", true); | ||
|
|
||
| bool splitpixel = config_root.get<bool>("splitpixel", false); | ||
|
|
||
| Int_t trackerphase = config_root.get<int>("trackerphase"); | ||
|
|
||
| TString lumitype = config_root.get<std::string>("trends.lumitype"); | ||
|
|
||
| Int_t firstrun = config_root.get<int>("trends.firstrun"); | ||
|
|
||
| Int_t lastrun = config_root.get<int>("trends.lastrun"); | ||
|
|
||
| std::string lumibyrunfile = config_root.get<std::string>("lumibyrunfile", ""); | ||
| std::ifstream lumifile(lumibyrunfile); | ||
| assert(lumifile.good()); | ||
| const Run2Lumi lumi_per_run(lumibyrunfile, firstrun, lastrun, 1000); | ||
|
|
||
| std::map<TString, Int_t> 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<int>(); | ||
|
|
||
| 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<Int_t, TString> Sublevel_Subdetector_Map = { | ||
| {1, "PXB"}, {2, "PXF"}, {3, "TIB"}, {4, "TID"}, {5, "TOB"}, {6, "TEC"}}; | ||
|
|
||
| // wheels/layers per tracker phase: <pahse,<sublevel,number of layers/wheels>> | ||
| const std::map<Int_t, std::map<Int_t, Int_t>> Phase_Subdetector_Layers_Map = { | ||
| {0, {{1, 3}, {2, 2}}}, {1, {{1, 4}, {2, 3}}}, {2, {{1, 4}, {2, 12}}}}; | ||
|
There was a problem hiding this comment. so in the end you dont plan to have this through the topology, right? There was a problem hiding this comment. this is essentially a plotting macro, and it doesn't have access to There was a problem hiding this comment. No, I think it would be over-killing it for a simple 2D mapping. Moreover, this macros is currently in an almost-cmssw-free standalone version, only requiring ROOT & boost, it'd be preferable to leave it simple as that considering it works on standard ROOT TTrees |
||
|
|
||
| // 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<Int_t, TString> Index_Variable_Map = { | ||
| {0, "DX"}, {1, "DY"}, {2, "DZ"}, {3, "DAlpha"}, {4, "DBeta"}, {5, "DGamma"}}; | ||
| const std::map<TString, TString> 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<Int_t, TString> Index_Estimator_Map = {{0, "Mean"}, {1, "Sigma"}}; | ||
| const std::map<TString, TString> 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<Int_t, std::map<Int_t, std::map<Int_t, TGraphErrors *>>> Graphs; | ||
| std::map<Int_t, std::map<Int_t, TH1D *>> 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<Int_t> id = {Reader, "id"}; | ||
| TTreeReaderValue<Int_t> badModuleQuality = {Reader, "badModuleQuality"}; | ||
| TTreeReaderValue<Int_t> inModuleList = {Reader, "inModuleList"}; | ||
| TTreeReaderValue<Int_t> level = {Reader, "level"}; | ||
| TTreeReaderValue<Int_t> mid = {Reader, "mid"}; | ||
| TTreeReaderValue<Int_t> mlevel = {Reader, "mlevel"}; | ||
| TTreeReaderValue<Int_t> sublevel = {Reader, "sublevel"}; | ||
| TTreeReaderValue<Float_t> x = {Reader, "x"}; | ||
| TTreeReaderValue<Float_t> y = {Reader, "y"}; | ||
| TTreeReaderValue<Float_t> z = {Reader, "z"}; | ||
| TTreeReaderValue<Float_t> r = {Reader, "r"}; | ||
| TTreeReaderValue<Float_t> phi = {Reader, "phi"}; | ||
| TTreeReaderValue<Float_t> eta = {Reader, "eta"}; | ||
| TTreeReaderValue<Float_t> alpha = {Reader, "alpha"}; | ||
| TTreeReaderValue<Float_t> beta = {Reader, "beta"}; | ||
| TTreeReaderValue<Float_t> gamma = {Reader, "gamma"}; | ||
| TTreeReaderValue<Float_t> dx = {Reader, "dx"}; | ||
| TTreeReaderValue<Float_t> dy = {Reader, "dy"}; | ||
| TTreeReaderValue<Float_t> dz = {Reader, "dz"}; | ||
| TTreeReaderValue<Float_t> dr = {Reader, "dr"}; | ||
| TTreeReaderValue<Float_t> dphi = {Reader, "dphi"}; | ||
| TTreeReaderValue<Float_t> dalpha = {Reader, "dalpha"}; | ||
| TTreeReaderValue<Float_t> dbeta = {Reader, "dbeta"}; | ||
| TTreeReaderValue<Float_t> dgamma = {Reader, "dgamma"}; | ||
| TTreeReaderValue<Float_t> du = {Reader, "du"}; | ||
| TTreeReaderValue<Float_t> dv = {Reader, "dv"}; | ||
| TTreeReaderValue<Float_t> dw = {Reader, "dw"}; | ||
| TTreeReaderValue<Float_t> da = {Reader, "da"}; | ||
| TTreeReaderValue<Float_t> db = {Reader, "db"}; | ||
| TTreeReaderValue<Float_t> dg = {Reader, "dg"}; | ||
| TTreeReaderValue<Int_t> useDetId = {Reader, "useDetId"}; | ||
| TTreeReaderValue<Int_t> detDim = {Reader, "detDim"}; | ||
| TTreeReaderArray<Int_t> 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; | ||
| } | ||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think C++ native types are generally preferred to ROOT types, but perhaps it's not mandatory here, as this is essentially a plotting macro.