ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Mat_map_surface_plot.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Mat_map_surface_plot.C
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2020 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include <TROOT.h>
10 
11 #include "materialPlotHelper.cpp"
12 
13 #include <fstream>
14 #include <iostream>
15 #include <sstream>
16 
18 
19 void plot(std::vector<TH2F*> Map, const sinfo& surface_info, const std::string& name){
20 
21  std::string out_name = name+"/"+surface_info.name+"/"+surface_info.name+"_"+surface_info.idname;
22  gSystem->Exec( Form("mkdir %s", (name+"/"+surface_info.name).c_str()) );
23 
24  // Disk
25  if(surface_info.type == 2){
26 
27  TText *vol = new TText(.1,.95,surface_info.name.c_str());
28  vol->SetNDC();
29  TText *surface = new TText(.1,.9,surface_info.id.c_str());
30  surface->SetNDC();
31  TText *surface_z = new TText(.1,.85,("Z = " + to_string(surface_info.pos)).c_str() );
32  surface_z->SetNDC();
33 
34  TCanvas *c1 = new TCanvas("c1","mat_X0",1200,1200);
35  c1->SetRightMargin(0.14);
36  c1->SetTopMargin(0.14);
37  c1->SetLeftMargin(0.14);
38  c1->SetBottomMargin(0.14);
39  Map[0]->Draw("COLZ");
40  vol->Draw();
41  surface->Draw();
42  surface_z->Draw();
43  c1->Print( (out_name+"_X0.pdf").c_str());
44  //c1->Print( (out_name+"_X0.root").c_str());
45 
46  delete c1;
47 
48  delete vol;
49  delete surface;
50  delete surface_z;
51  }
52 
53  // Cylinder
54  if(surface_info.type == 1){
55 
56  TText *vol = new TText(.1,.95,surface_info.name.c_str());
57  vol->SetNDC();
58  TText *surface = new TText(.1,.9,surface_info.id.c_str());
59  surface->SetNDC();
60  TText *surface_r = new TText(.1,.85,("R = " + to_string(surface_info.pos)).c_str() );
61  surface_r->SetNDC();
62  TCanvas *c1 = new TCanvas("c1","mat_X0",1200,1200);
63  c1->SetRightMargin(0.14);
64  c1->SetTopMargin(0.14);
65  c1->SetLeftMargin(0.14);
66  c1->SetBottomMargin(0.14);
67  Map[0]->Draw("COLZ");
68  vol->Draw();
69  surface->Draw();
70  surface_r->Draw();
71  c1->Print( (out_name+"_X0.pdf").c_str());
72  //c1->Print( (out_name+"_X0.root").c_str());
73 
74  delete c1;
75 
76  delete vol;
77  delete surface;
78  delete surface_r;
79  }
80  return;
81 }
82 
84 
85 void Initialise_hist(std::vector<TH2F*>& surface_hist,
86  const sinfo& surface_info){
87 
88  TH2F * Map_X0;
89  TH2F * Map_L0;
90 
91  TH2F * Map_scale;
92 
93  if(surface_info.type == 1){
94  Map_X0 = new TH2F(("Map_X0_"+surface_info.idname).c_str(),("Map_X0_"+surface_info.idname).c_str(),
95  50,-6,6,50,-3.2,3.2);
96  Map_L0 = new TH2F(("Map_L0_"+surface_info.idname).c_str(),("Map_L0_"+surface_info.idname).c_str(),
97  50,-6,6,50,-3.2,3.2);
98  Map_scale = new TH2F(("Map_scale_"+surface_info.idname).c_str(),("Map_scale_"+surface_info.idname).c_str(),
99  50,-6,6,50,-3.2,3.2);
100  Map_X0->GetXaxis()->SetTitle("Eta");
101  Map_X0->GetYaxis()->SetTitle("Phi");
102  Map_X0->GetZaxis()->SetTitle("X0");
103  Map_L0->GetXaxis()->SetTitle("Eta");
104  Map_L0->GetYaxis()->SetTitle("Phi");
105  Map_L0->GetZaxis()->SetTitle("L0");
106  }
107 
108  if(surface_info.type == 2){
109  Map_X0 = new TH2F(("Map_X0_"+surface_info.idname).c_str(),("Map_X0_"+surface_info.idname).c_str(),
110  50,-1*surface_info.range_max,surface_info.range_max,50,-1*surface_info.range_max, surface_info.range_max);
111  Map_L0 = new TH2F(("Map_L0_"+surface_info.idname).c_str(),("Map_L0_"+surface_info.idname).c_str(),
112  50,-1*surface_info.range_max,surface_info.range_max,50,-1*surface_info.range_max, surface_info.range_max);
113  Map_scale = new TH2F(("Map_scale_"+surface_info.idname).c_str(),("Map_scale_"+surface_info.idname).c_str(),
114  50,-1*surface_info.range_max,surface_info.range_max,50,-1*surface_info.range_max, surface_info.range_max);
115  Map_X0->GetXaxis()->SetTitle("X [mm]");
116  Map_X0->GetYaxis()->SetTitle("Y [mm]");
117  Map_X0->GetZaxis()->SetTitle("X0");
118  Map_L0->GetXaxis()->SetTitle("X [mm]");
119  Map_L0->GetYaxis()->SetTitle("Y [mm]");
120  Map_L0->GetZaxis()->SetTitle("L0");
121  }
122  std::vector<TH2F*> v_hist;
123  v_hist.push_back(Map_X0);
124  v_hist.push_back(Map_L0);
125  v_hist.push_back(Map_scale);
126  surface_hist = v_hist;
127 }
128 
130 
131 void Fill(std::map<uint64_t,std::vector<TH2F*>>& surface_hist, std::map<uint64_t,sinfo>& surface_info,
132  const std::string& input_file, const std::string& json_surface_file, const int& nbprocess){
133 
134  nlohmann::json Det;
135  std::map<std::string,std::string> surface_name;
136 
137  if(json_surface_file != ""){
138  std::ifstream lfile(json_surface_file.c_str());
139  lfile >> Det;
140 
141  Parse_Json(Det, surface_name);
142  }
143 
144  std::map<uint64_t,float> surface_weight;
145 
146  //Get old file, old tree and set top branch address
147  TFile *tfile = new TFile(input_file.c_str());
148  TTree *tree = (TTree*)tfile->Get("material-tracks");
149 
150  float v_phi = 0;
151  float v_eta = 0;
152  std::vector<float> *mat_X0 = 0;
153  std::vector<float> *mat_L0 = 0;
154  std::vector<float> *mat_step_length = 0;
155 
156  std::vector<uint64_t> *sur_id = 0;
157  std::vector<int32_t> *sur_type = 0;
158  std::vector<float> *sur_x = 0;
159  std::vector<float> *sur_y = 0;
160  std::vector<float> *sur_z = 0;
161  std::vector<float> *sur_range_min = 0;
162  std::vector<float> *sur_range_max = 0;
163 
164  tree->SetBranchAddress("v_phi",&v_phi);
165  tree->SetBranchAddress("v_eta",&v_eta);
166  tree->SetBranchAddress("mat_X0",&mat_X0);
167  tree->SetBranchAddress("mat_L0",&mat_L0);
168  tree->SetBranchAddress("mat_step_length",&mat_step_length);
169 
170  tree->SetBranchAddress("sur_id",&sur_id);
171  tree->SetBranchAddress("sur_type",&sur_type);
172  tree->SetBranchAddress("sur_x",&sur_x);
173  tree->SetBranchAddress("sur_y",&sur_y);
174  tree->SetBranchAddress("sur_z",&sur_z);
175  tree->SetBranchAddress("sur_range_min",&sur_range_min);
176  tree->SetBranchAddress("sur_range_max",&sur_range_max);
177 
178  int nentries = tree->GetEntries();
179  if(nentries > nbprocess && nbprocess != -1) nentries = nbprocess;
180  // Loop over all the material tracks.
181  for (Long64_t i=0;i<nentries; i++) {
182  if(i%10000==0) std::cout << "processed " << i << " events out of " << nentries << std::endl;
183  tree->GetEntry(i);
184 
185  // Reset the weight
186  for (auto weight_it = surface_weight.begin(); weight_it != surface_weight.end(); weight_it++){
187  weight_it->second = 0;
188  }
189  // loop over all the material hit to do initialisation and compute weight
190  for(int j=0; j<mat_X0->size(); j++ ){
191 
192  // Ignore surface of incorrect type
193  if(sur_type->at(j) == -1) continue;
194  // If a surface was never encountered initialise the hist, info and weight
195  if(surface_hist.find(sur_id->at(j))==surface_hist.end()){
196  float pos;
197  if(sur_type->at(j) == 1){
198  pos = sqrt(sur_x->at(j)*sur_x->at(j)+sur_y->at(j)*sur_y->at(j));
199  }
200  if(sur_type->at(j) == 2){
201  pos = sur_z->at(j);
202  }
203  surface_weight[sur_id->at(j)] = 0;
204  Initialise_info(surface_info[sur_id->at(j)], surface_name, sur_id->at(j), sur_type->at(j), pos, sur_range_min->at(j), sur_range_max->at(j));
205  Initialise_hist(surface_hist[sur_id->at(j)], surface_info[sur_id->at(j)]);
206  }
207  // Weight for each surface = number of hit associated to it.
208  surface_weight[sur_id->at(j)]++;
209  }
210 
211  // loop over all the material hit to fill the histogram
212  for(int j=0; j<mat_X0->size(); j++ ){
213 
214  // Ignore surface of incorrect type
215  if(sur_type->at(j) == -1) continue;
216 
217  if(sur_type->at(j) == 1){
218  surface_hist[sur_id->at(j)][0]->Fill(v_eta, v_phi, (mat_step_length->at(j)/mat_X0->at(j)));
219  surface_hist[sur_id->at(j)][1]->Fill(v_eta, v_phi, (mat_step_length->at(j)/mat_L0->at(j)));
220  surface_hist[sur_id->at(j)][2]->Fill(v_eta, v_phi, (1/surface_weight[sur_id->at(j)]));
221  }
222  if(sur_type->at(j) == 2){
223  surface_hist[sur_id->at(j)][0]->Fill(sur_x->at(j), sur_y->at(j), (mat_step_length->at(j)/mat_X0->at(j)));
224  surface_hist[sur_id->at(j)][1]->Fill(sur_x->at(j), sur_y->at(j), (mat_step_length->at(j)/mat_L0->at(j)));
225  surface_hist[sur_id->at(j)][2]->Fill(sur_x->at(j), sur_y->at(j), (1/surface_weight[sur_id->at(j)]));
226  }
227  }
228  }
229  // Normalise the histograms
230  for (auto hist_it = surface_hist.begin(); hist_it != surface_hist.end(); hist_it++){
231  hist_it->second[0]->Divide(hist_it->second[2]);
232  hist_it->second[1]->Divide(hist_it->second[2]);
233  }
234 }
235 
242 
243 void Mat_map_surface_plot(std::string input_file = "", std::string json_surface_file = "", int nbprocess = -1, std::string name = ""){
244 
245  gStyle->SetOptStat(0);
246  gStyle->SetOptTitle(0);
247 
248  std::map<uint64_t,std::vector<TH2F*>> surface_hist;
249  std::map<uint64_t,sinfo> surface_info;
250 
251  Fill(surface_hist, surface_info, input_file, json_surface_file, nbprocess);
252  for (auto hist_it = surface_hist.begin(); hist_it != surface_hist.end(); hist_it++){
253  plot(hist_it->second, surface_info[hist_it->first], name);
254  for (auto hist : hist_it->second){
255  delete hist;
256  }
257  hist_it->second.clear();
258  }
259 }