CMS 3D CMS Logo

Functions
GCPtrends.cc File Reference
#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"

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 22 of file GCPtrends.cc.

References simBeamSpotPI::alpha, dir2webdir::argc, GCPpyPlots::argv, cms::cuda::assert(), HLT_2024v12_cff::beta, gather_cfg::cout, Calorimetry_cff::da, dqmiodatasetharvest::db, l1ctLayer1_cff::dr, PVValHelper::dx, PVValHelper::dy, PVValHelper::dz, PVValHelper::eta, beamvalidation::exit(), writedatasetfile::firstrun, CustomPhysics_cfi::gamma, utils::gErrorIgnoreLevel, plotFactory::Histos, fwlog::kError, writedatasetfile::lastrun, personalPlayback::level, Trend::lgd, SiStripPI::mean, dtNoiseDBValidation_cfg::outputdir, castor_dqm_sourceclient_file_cfg::path, PVValHelper::phi, alignCSCRings::r, contentValuesCheck::run_number, AlCaHLTBitMon_QueryRunRegistry::string, to_string(), units(), HPSPFTaus_cff::variable, x, detailsBasic3DVector::y, and detailsBasic3DVector::z.

22  {
23  if (argc < 2) {
24  std::cout << "Usage of the program : GCPtrends [gcp_trend_config.json]" << std::endl;
25  std::cout << "gcp_trend_config.json : file with configuration in json format" << std::endl;
26  exit(1);
27  }
28 
29  // parsing input file
30  std::string lumi_file(argv[1]);
31  boost::property_tree::ptree config_root;
32  boost::property_tree::read_json(lumi_file, config_root);
33 
34  // configuration
35  TString outputdir = config_root.get<std::string>("outputdir", ".");
36  outputdir = TString(outputdir) + TString(outputdir.EndsWith("/") ? "" : "/");
37 
38  bool usebadmodules = config_root.get<bool>("usebadmodules", true);
39 
40  bool splitpixel = config_root.get<bool>("splitpixel", false);
41 
42  Int_t trackerphase = config_root.get<int>("trackerphase");
43 
44  TString lumitype = config_root.get<std::string>("trends.lumitype");
45 
46  Int_t firstrun = config_root.get<int>("trends.firstrun");
47 
48  Int_t lastrun = config_root.get<int>("trends.lastrun");
49 
50  std::string lumibyrunfile = config_root.get<std::string>("lumibyrunfile", "");
51  std::ifstream lumifile(lumibyrunfile);
52  assert(lumifile.good());
53  const Run2Lumi lumi_per_run(lumibyrunfile, firstrun, lastrun, 1000);
54 
55  std::map<TString, Int_t> comparison_map;
56  for (boost::property_tree::ptree::value_type &comparison : config_root.get_child("comparisons")) {
57  TString path = comparison.first;
58  Int_t run_label = comparison.second.get_value<int>();
59 
60  comparison_map[path] = run_label;
61  }
62 
63  // create output file
64  TFile *file_out = TFile::Open(outputdir + "GCPTrends.root", "RECREATE");
65  file_out->cd();
66 
68 
69  // sub-detector mapping
70  std::map<Int_t, TString> Sublevel_Subdetector_Map = {
71  {1, "PXB"}, {2, "PXF"}, {3, "TIB"}, {4, "TID"}, {5, "TOB"}, {6, "TEC"}};
72 
73  // wheels/layers per tracker phase: <pahse,<sublevel,number of layers/wheels>>
74  const std::map<Int_t, std::map<Int_t, Int_t>> Phase_Subdetector_Layers_Map = {
75  {0, {{1, 3}, {2, 2}}}, {1, {{1, 4}, {2, 3}}}, {2, {{1, 4}, {2, 12}}}};
76 
77  // adding layers/wheels in case Pixel is requested further split
78  if (splitpixel) {
79  assert(trackerphase < 3);
80 
81  for (auto &sub_layer : Phase_Subdetector_Layers_Map.at(trackerphase)) {
82  for (int iLayer = 1; iLayer <= sub_layer.second; iLayer++) {
83  // subid*100+subsubid
84  if (sub_layer.first % 2 != 0) {
85  Sublevel_Subdetector_Map[100 * sub_layer.first + iLayer] =
86  Sublevel_Subdetector_Map[sub_layer.first] + "Layer" + TString(std::to_string(iLayer));
87 
88  } else {
89  Sublevel_Subdetector_Map[100 * sub_layer.first + (1 - 1) * sub_layer.second + iLayer] =
90  Sublevel_Subdetector_Map[sub_layer.first] + "Wheel" + TString(std::to_string(iLayer)) + "Side1";
91  Sublevel_Subdetector_Map[100 * sub_layer.first + (2 - 1) * sub_layer.second + iLayer] =
92  Sublevel_Subdetector_Map[sub_layer.first] + "Wheel" + TString(std::to_string(iLayer)) + "Side2";
93  }
94  }
95  }
96  }
97 
98  // variable mapping
99  const std::map<Int_t, TString> Index_Variable_Map = {
100  {0, "DX"}, {1, "DY"}, {2, "DZ"}, {3, "DAlpha"}, {4, "DBeta"}, {5, "DGamma"}};
101  const std::map<TString, TString> Variable_Name_Map = {{"DX", "#Delta x"},
102  {"DY", "#Delta y"},
103  {"DZ", "#Delta z"},
104  {"DAlpha", "#Delta #alpha"},
105  {"DBeta", "#Delta #beta"},
106  {"DGamma", "#Delta #gamma"}};
107  // estimator mapping
108  const std::map<Int_t, TString> Index_Estimator_Map = {{0, "Mean"}, {1, "Sigma"}};
109  const std::map<TString, TString> Estimator_Name_Map = {{"Mean", "#mu"}, {"Sigma", "#sigma"}};
110 
111  // constant unit conversion
112  const int convert_cm_to_mum = 10000;
113  const int convert_rad_to_mrad = 1000;
114 
115  std::cout << std::endl;
116  std::cout << " ==> " << comparison_map.size() << " GCPs to be processed in configuration file ... " << std::endl;
117  std::cout << std::endl;
118 
119  // booking TGraphs and TH1D
120  TH1::StatOverflows(kTRUE);
121  std::map<Int_t, std::map<Int_t, std::map<Int_t, TGraphErrors *>>> Graphs;
122  std::map<Int_t, std::map<Int_t, TH1D *>> Histos;
123  for (auto &level : Sublevel_Subdetector_Map) {
124  for (auto &variable : Index_Variable_Map) {
125  Histos[level.first][variable.first] = new TH1D(
126  "Histo_" + level.second + "_" + variable.second, "Histo_" + level.second + "_" + variable.second, 1, -1, 1);
127 
128  for (auto &estimator : Index_Estimator_Map) {
129  Graphs[level.first][variable.first][estimator.first] = new TGraphErrors();
130  }
131  }
132  }
133 
134  // loop over the comparisons (GCPs)
135  for (auto &path_run_map : comparison_map) {
136  TChain Events("alignTree", "alignTree");
137  Events.Add(path_run_map.first);
138 
139  Int_t run_number = path_run_map.second;
140 
141  std::cout << " --> processing GCPtree file: " << path_run_map.first << std::endl;
142 
143  TTreeReader Reader(&Events);
144  Reader.Restart();
145 
146  TTreeReaderValue<Int_t> id = {Reader, "id"};
147  TTreeReaderValue<Int_t> badModuleQuality = {Reader, "badModuleQuality"};
148  TTreeReaderValue<Int_t> inModuleList = {Reader, "inModuleList"};
149  TTreeReaderValue<Int_t> level = {Reader, "level"};
150  TTreeReaderValue<Int_t> mid = {Reader, "mid"};
151  TTreeReaderValue<Int_t> mlevel = {Reader, "mlevel"};
152  TTreeReaderValue<Int_t> sublevel = {Reader, "sublevel"};
153  TTreeReaderValue<Float_t> x = {Reader, "x"};
154  TTreeReaderValue<Float_t> y = {Reader, "y"};
155  TTreeReaderValue<Float_t> z = {Reader, "z"};
156  TTreeReaderValue<Float_t> r = {Reader, "r"};
157  TTreeReaderValue<Float_t> phi = {Reader, "phi"};
158  TTreeReaderValue<Float_t> eta = {Reader, "eta"};
159  TTreeReaderValue<Float_t> alpha = {Reader, "alpha"};
160  TTreeReaderValue<Float_t> beta = {Reader, "beta"};
161  TTreeReaderValue<Float_t> gamma = {Reader, "gamma"};
162  TTreeReaderValue<Float_t> dx = {Reader, "dx"};
163  TTreeReaderValue<Float_t> dy = {Reader, "dy"};
164  TTreeReaderValue<Float_t> dz = {Reader, "dz"};
165  TTreeReaderValue<Float_t> dr = {Reader, "dr"};
166  TTreeReaderValue<Float_t> dphi = {Reader, "dphi"};
167  TTreeReaderValue<Float_t> dalpha = {Reader, "dalpha"};
168  TTreeReaderValue<Float_t> dbeta = {Reader, "dbeta"};
169  TTreeReaderValue<Float_t> dgamma = {Reader, "dgamma"};
170  TTreeReaderValue<Float_t> du = {Reader, "du"};
171  TTreeReaderValue<Float_t> dv = {Reader, "dv"};
172  TTreeReaderValue<Float_t> dw = {Reader, "dw"};
173  TTreeReaderValue<Float_t> da = {Reader, "da"};
174  TTreeReaderValue<Float_t> db = {Reader, "db"};
175  TTreeReaderValue<Float_t> dg = {Reader, "dg"};
176  TTreeReaderValue<Int_t> useDetId = {Reader, "useDetId"};
177  TTreeReaderValue<Int_t> detDim = {Reader, "detDim"};
178  TTreeReaderArray<Int_t> identifiers = {Reader, "identifiers"};
179 
180  // loop over modules
181  while (Reader.Next()) {
182  if (!usebadmodules)
183  if (*badModuleQuality.Get())
184  continue;
185 
186  Int_t sublevel_idx = *sublevel.Get();
187 
188  Histos[sublevel_idx][0]->Fill(*dx.Get() * convert_cm_to_mum);
189  Histos[sublevel_idx][1]->Fill(*dy.Get() * convert_cm_to_mum);
190  Histos[sublevel_idx][2]->Fill(*dz.Get() * convert_cm_to_mum);
191  Histos[sublevel_idx][3]->Fill(*dalpha.Get() * convert_rad_to_mrad);
192  Histos[sublevel_idx][4]->Fill(*dbeta.Get() * convert_rad_to_mrad);
193  Histos[sublevel_idx][5]->Fill(*dgamma.Get() * convert_rad_to_mrad);
194 
195  if (splitpixel && sublevel_idx <= 2) {
196  Int_t layer_index;
197 
198  if (sublevel_idx % 2 != 0)
199  layer_index = 100 * sublevel_idx + identifiers[2];
200  else
201  layer_index = 100 * sublevel_idx +
202  (identifiers[4] - 1) * Phase_Subdetector_Layers_Map.at(trackerphase).at(sublevel_idx) +
203  identifiers[3];
204 
205  Histos[layer_index][0]->Fill(*dx.Get() * convert_cm_to_mum);
206  Histos[layer_index][1]->Fill(*dy.Get() * convert_cm_to_mum);
207  Histos[layer_index][2]->Fill(*dz.Get() * convert_cm_to_mum);
208  Histos[layer_index][3]->Fill(*dalpha.Get() * convert_rad_to_mrad);
209  Histos[layer_index][4]->Fill(*dbeta.Get() * convert_rad_to_mrad);
210  Histos[layer_index][5]->Fill(*dgamma.Get() * convert_rad_to_mrad);
211  }
212  }
213 
214  for (auto &level : Sublevel_Subdetector_Map) {
215  for (auto &variable : Index_Variable_Map) {
216  Double_t mean = Histos[level.first][variable.first]->GetMean();
217  Double_t sigma = Histos[level.first][variable.first]->GetStdDev();
218 
219  Graphs[level.first][variable.first][0]->AddPoint(run_number, mean);
220  if (std::fabs(mean) > Graphs[level.first][variable.first][0]->GetMaximum() && std::fabs(mean) > 0.) {
221  Graphs[level.first][variable.first][0]->SetMaximum(std::fabs(mean));
222  Graphs[level.first][variable.first][0]->SetMinimum(-std::fabs(mean));
223  }
224 
225  Graphs[level.first][variable.first][1]->AddPoint(run_number, sigma);
226  if (sigma > Graphs[level.first][variable.first][1]->GetMaximum() && sigma > 0.) {
227  Graphs[level.first][variable.first][1]->SetMaximum(sigma);
228  Graphs[level.first][variable.first][1]->SetMinimum(0.);
229  }
230 
231  Histos[level.first][variable.first]->Reset("ICESM");
232  }
233  }
234  }
235 
236  // saving TGraphs
237  for (auto &level : Sublevel_Subdetector_Map) {
238  for (auto &variable : Index_Variable_Map) {
239  for (auto &estimator : Index_Estimator_Map) {
240  Graphs[level.first][variable.first][estimator.first]->Write("Graph_" + level.second + "_" + variable.second +
241  "_" + estimator.second);
242 
243  TString units = "mrad";
244  if (variable.second.Contains("DX") || variable.second.Contains("DY") || variable.second.Contains("DZ"))
245  units = "#mum";
246 
247  Trend trend("Graph_" + level.second + "_" + variable.second + "_" + estimator.second,
248  "output",
249  "Trend",
250  Estimator_Name_Map.at(estimator.second) + "_{" + Variable_Name_Map.at(variable.second) + "} [" +
251  units + "]",
252  -2 * std::fabs(Graphs[level.first][variable.first][estimator.first]->GetMinimum()),
253  2 * std::fabs(Graphs[level.first][variable.first][estimator.first]->GetMaximum()),
254  config_root,
255  lumi_per_run,
256  lumitype);
257 
258  Graphs[level.first][variable.first][estimator.first]->SetFillColor(4);
259 
260  trend.lgd.SetHeader(level.second, "center");
261  trend.lgd.AddEntry(Graphs[level.first][variable.first][estimator.first], "Average over all modules", "pl");
262 
263  trend(Graphs[level.first][variable.first][estimator.first], "p", "pf", false);
264  }
265  }
266  }
267 
268  file_out->Close();
269 
270  std::cout << std::endl;
271  std::cout << " ==> done." << std::endl;
272  std::cout << std::endl;
273 
274  return 0;
275 }
Definition: Trend.h:22
assert(be >=bs)
float float float z
static std::string to_string(const XMLCh *ch)
Definition: Trend.h:78
gErrorIgnoreLevel
Definition: utils.py:27
TString units(TString variable, Char_t axis)
float x
Definition: Histos.h:19
def exit(msg="")