EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Mat_map_surface_plot_1D.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Mat_map_surface_plot_1D.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 
17 // Draw the two 1D histograms for each surface.
18 
19 void plot(std::vector<TH1F*> 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_R",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("P HIST");
40  vol->Draw();
41  surface->Draw();
42  surface_z->Draw();
43  c1->Print( (out_name+"_R.pdf").c_str());
44  //c1->Print( (out_name+"_X0.root").c_str());
45 
46  TCanvas *c2 = new TCanvas("c2","mat_Phi",1200,1200);
47  c2->SetRightMargin(0.14);
48  c2->SetTopMargin(0.14);
49  c2->SetLeftMargin(0.14);
50  c2->SetBottomMargin(0.14);
51  Map[1]->Draw("P HIST");
52  vol->Draw();
53  surface->Draw();
54  surface_z->Draw();
55  c2->Print( (out_name+"_Phi.pdf").c_str());
56  //c2->Print( (out_name+"_Phi.root").c_str());
57 
58  delete c1;
59  delete c2;
60 
61  delete vol;
62  delete surface;
63  delete surface_z;
64  }
65 
66  // Cylinder
67  if(surface_info.type == 1){
68 
69  TText *vol = new TText(.1,.95,surface_info.name.c_str());
70  vol->SetNDC();
71  TText *surface = new TText(.1,.9,surface_info.id.c_str());
72  surface->SetNDC();
73  TText *surface_r = new TText(.1,.85,("R = " + to_string(surface_info.pos)).c_str() );
74  surface_r->SetNDC();
75  TCanvas *c1 = new TCanvas("c1","mat_Z",1200,1200);
76  c1->SetRightMargin(0.14);
77  c1->SetTopMargin(0.14);
78  c1->SetLeftMargin(0.14);
79  c1->SetBottomMargin(0.14);
80  Map[0]->Draw("P HIST");
81  vol->Draw();
82  surface->Draw();
83  surface_r->Draw();
84  c1->Print( (out_name+"_Z.pdf").c_str());
85  //c1->Print( (out_name+"_X0.root").c_str());
86 
87  TCanvas *c2 = new TCanvas("c2","mat_Phi",1200,1200);
88  c2->SetRightMargin(0.14);
89  c2->SetTopMargin(0.14);
90  c2->SetLeftMargin(0.14);
91  c2->SetBottomMargin(0.14);
92  Map[1]->Draw("P HIST");
93  vol->Draw();
94  surface->Draw();
95  surface_r->Draw();
96  c2->Print( (out_name+"_Phi.pdf").c_str());
97  //c2->Print( (out_name+"_Phi.root").c_str());
98 
99  delete c1;
100  delete c2;
101 
102  delete vol;
103  delete surface;
104  delete surface_r;
105  }
106  return;
107 }
108 
110 
111 void Initialise_hist(std::vector<TH1F*>& surface_hist,
112  const sinfo& surface_info){
113 
114  TH1F * Map;
115  TH1F * Map_Phi;
116  TH1F * Map_scale;
117  TH1F * Map_scale_Phi;
118 
119  if(surface_info.type == 1){
120  Map = new TH1F(("Map_"+surface_info.idname).c_str(),("Map_"+surface_info.idname).c_str(),
121  50,-1*surface_info.range_max, surface_info.range_max);
122  Map_Phi = new TH1F(("Map_Phi_"+surface_info.idname).c_str(),("Map_Phi_"+surface_info.idname).c_str(),
123  50,-3.2,3.2);
124  Map_scale = new TH1F(("Map_scale_"+surface_info.idname).c_str(),("Map_scale_"+surface_info.idname).c_str(),
125  50,-1*surface_info.range_max, surface_info.range_max);
126  Map_scale_Phi = new TH1F(("Map_scale_Phi_"+surface_info.idname).c_str(),("Map_scale_Phi_"+surface_info.idname).c_str(),
127  50,-3.2,3.2);
128  Map->GetXaxis()->SetTitle("Z [mm]");
129  Map->GetYaxis()->SetTitle("X0");
130  Map_Phi->GetXaxis()->SetTitle("Phi");
131  Map_Phi->GetYaxis()->SetTitle("X0");
132  }
133 
134  if(surface_info.type == 2){
135  Map = new TH1F(("Map_"+surface_info.idname).c_str(),("Map_"+surface_info.idname).c_str(),
136  50,surface_info.range_min, surface_info.range_max);
137  Map_Phi = new TH1F(("Map_Phi_"+surface_info.idname).c_str(),("Map_Phi_"+surface_info.idname).c_str(),
138  50,-3.2,3.2);
139  Map_scale = new TH1F(("Map_scale_"+surface_info.idname).c_str(),("Map_scale_"+surface_info.idname).c_str(),
140  50,surface_info.range_min, surface_info.range_max);
141  Map_scale_Phi = new TH1F(("Map_scale_Phi_"+surface_info.idname).c_str(),("Map_scale_Phi_"+surface_info.idname).c_str(),
142  50,-3.2,3.2);
143  Map->GetXaxis()->SetTitle("R [mm]");
144  Map->GetYaxis()->SetTitle("X0");
145  Map_Phi->GetXaxis()->SetTitle("Phi");
146  Map_Phi->GetYaxis()->SetTitle("X0");
147  }
148 
149  Map->SetMarkerStyle(20);
150  Map_Phi->SetMarkerStyle(20);
151 
152  std::vector<TH1F*> v_hist;
153  v_hist.push_back(Map);
154  v_hist.push_back(Map_Phi);
155  v_hist.push_back(Map_scale);
156  v_hist.push_back(Map_scale_Phi);
157  surface_hist = v_hist;
158 }
159 
161 
162 void Fill(std::map<uint64_t,std::vector<TH1F*>>& surface_hist, std::map<uint64_t,sinfo>& surface_info,
163  const std::string& input_file, const std::string& json_surface_file, const int& nbprocess){
164 
165  nlohmann::json Det;
166  std::map<std::string,std::string> surface_name;
167 
168  if(json_surface_file != ""){
169  std::ifstream lfile(json_surface_file.c_str());
170  lfile >> Det;
171 
172  Parse_Json(Det, surface_name);
173  }
174 
175  std::map<uint64_t,float> surface_weight;
176 
177  //Get file, tree and set top branch address
178  TFile *tfile = new TFile(input_file.c_str());
179  TTree *tree = (TTree*)tfile->Get("material-tracks");
180 
181  float v_phi = 0;
182  float v_eta = 0;
183  std::vector<float> *mat_X0 = 0;
184  std::vector<float> *mat_L0 = 0;
185  std::vector<float> *mat_step_length = 0;
186 
187  std::vector<uint64_t> *sur_id = 0;
188  std::vector<int32_t> *sur_type = 0;
189  std::vector<float> *sur_x = 0;
190  std::vector<float> *sur_y = 0;
191  std::vector<float> *sur_z = 0;
192  std::vector<float> *sur_range_min = 0;
193  std::vector<float> *sur_range_max = 0;
194 
195  tree->SetBranchAddress("v_phi",&v_phi);
196  tree->SetBranchAddress("v_eta",&v_eta);
197  tree->SetBranchAddress("mat_X0",&mat_X0);
198  tree->SetBranchAddress("mat_L0",&mat_L0);
199  tree->SetBranchAddress("mat_step_length",&mat_step_length);
200 
201  tree->SetBranchAddress("sur_id",&sur_id);
202  tree->SetBranchAddress("sur_type",&sur_type);
203  tree->SetBranchAddress("sur_x",&sur_x);
204  tree->SetBranchAddress("sur_y",&sur_y);
205  tree->SetBranchAddress("sur_z",&sur_z);
206  tree->SetBranchAddress("sur_range_min",&sur_range_min);
207  tree->SetBranchAddress("sur_range_max",&sur_range_max);
208 
209  int nentries = tree->GetEntries();
210  if(nentries > nbprocess && nbprocess != -1) nentries = nbprocess;
211  // Loop over all the material tracks.
212  for (Long64_t i=0;i<nentries; i++) {
213  if(i%10000==0) std::cout << "processed " << i << " events out of " << nentries << std::endl;
214  tree->GetEntry(i);
215 
216  // Reset the weight
217  for (auto weight_it = surface_weight.begin(); weight_it != surface_weight.end(); weight_it++){
218  weight_it->second = 0;
219  }
220  // loop over all the material hits to do initialisation and compute weight
221  for(int j=0; j<mat_X0->size(); j++ ){
222 
223  // Ignore surface of incorrect type
224  if(sur_type->at(j) == -1) continue;
225  // If a surface was never encountered initialise the hist, info and weight
226  if(surface_hist.find(sur_id->at(j))==surface_hist.end()){
227 
228  float pos;
229  if(sur_type->at(j) == 1){
230  pos = sqrt(sur_x->at(j)*sur_x->at(j)+sur_y->at(j)*sur_y->at(j));
231  }
232  if(sur_type->at(j) == 2){
233  pos = sur_z->at(j);
234  }
235  // Weight for each surface = number of hit associated to it.
236  surface_weight[sur_id->at(j)] = 0;
237  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));
238  Initialise_hist(surface_hist[sur_id->at(j)], surface_info[sur_id->at(j)]);
239  }
240  // Weight for each surface = number of hit associated to it.
241  surface_weight[sur_id->at(j)]++;
242  }
243 
244  // loop over all the material hit to fill the histogram
245  for(int j=0; j<mat_X0->size(); j++ ){
246 
247  // Ignore surface of incorrect type
248  if(sur_type->at(j) == -1) continue;
249 
250  if(sur_type->at(j) == 1){
251  surface_hist[sur_id->at(j)][0]->Fill(sur_z->at(j), (mat_step_length->at(j)/mat_X0->at(j)));
252  surface_hist[sur_id->at(j)][1]->Fill(v_phi, (mat_step_length->at(j)/mat_L0->at(j)));
253  surface_hist[sur_id->at(j)][2]->Fill(sur_z->at(j), (1/surface_weight[sur_id->at(j)]));
254  surface_hist[sur_id->at(j)][3]->Fill(v_phi, (1/surface_weight[sur_id->at(j)]));
255  }
256  if(sur_type->at(j) == 2){
257  surface_hist[sur_id->at(j)][0]->Fill(sqrt(sur_x->at(j)*sur_x->at(j)+sur_y->at(j)*sur_y->at(j)), (mat_step_length->at(j)/mat_X0->at(j)));
258  surface_hist[sur_id->at(j)][1]->Fill(v_phi, (mat_step_length->at(j)/mat_L0->at(j)));
259  surface_hist[sur_id->at(j)][2]->Fill(sqrt(sur_x->at(j)*sur_x->at(j)+sur_y->at(j)*sur_y->at(j)), (1/surface_weight[sur_id->at(j)]));
260  surface_hist[sur_id->at(j)][3]->Fill(v_phi, (1/surface_weight[sur_id->at(j)]));
261  }
262  }
263  }
264  // Normalise the histograms
265  for (auto hist_it = surface_hist.begin(); hist_it != surface_hist.end(); hist_it++){
266  hist_it->second[0]->Divide(hist_it->second[2]);
267  hist_it->second[1]->Divide(hist_it->second[3]);
268  }
269 }
270 
277 
278 void Mat_map_surface_plot_1D(std::string input_file = "", std::string json_surface_file = "", int nbprocess = -1, std::string name = ""){
279 
280  gStyle->SetOptStat(0);
281  gStyle->SetOptTitle(0);
282 
283  std::map<uint64_t,std::vector<TH1F*>> surface_hist;
284  std::map<uint64_t,sinfo> surface_info;
285 
286  Fill(surface_hist, surface_info, input_file, json_surface_file, nbprocess);
287  for (auto hist_it = surface_hist.begin(); hist_it != surface_hist.end(); hist_it++){
288  plot(hist_it->second, surface_info[hist_it->first], name);
289  for (auto hist : hist_it->second){
290  delete hist;
291  }
292  hist_it->second.clear();
293  }
294 }