EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Mat_map_surface_plot_dist.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Mat_map_surface_plot_dist.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 plot for each surface.
18 void plot(TGraph* Dist, const sinfo& surface_info, const std::string& name){
19 
20  std::string out_name = name+"/"+surface_info.name+"/"+surface_info.name+"_"+surface_info.idname;
21  gSystem->Exec( Form("mkdir %s", (name+"/"+surface_info.name).c_str()) );
22 
23  TCanvas *c = new TCanvas("c","dist",1200,1200);
24  c->SetRightMargin(0.14);
25  c->SetTopMargin(0.14);
26  c->SetLeftMargin(0.14);
27  c->SetBottomMargin(0.14);
28 
29  Dist->Draw("AP");
30  TLine *line_pos;
31 
32  TText *vol = new TText(.1,.95,surface_info.name.c_str());
33  TText *surface = new TText(.1,.9,surface_info.id.c_str());
34  TText *surface_z = new TText(.1,.85,("Z = " + to_string(surface_info.pos)).c_str() );
35  TText *surface_r = new TText(.1,.85,("R = " + to_string(surface_info.pos)).c_str() );
36  vol->SetNDC();
37  surface->SetNDC();
38  surface_z->SetNDC();
39  surface_r->SetNDC();
40 
41  // Disk
42  if(surface_info.type == 2){
43  vol->Draw();
44  surface->Draw();
45  surface_z->Draw();
46 
47  Dist->GetYaxis()->SetRangeUser(surface_info.range_min-(surface_info.range_max-surface_info.range_min)/10,
48  surface_info.range_max+(surface_info.range_max-surface_info.range_min)/10);
49 
50  // Position of the disk surface
51  line_pos = new TLine(surface_info.pos,surface_info.range_min,surface_info.pos,surface_info.range_max);
52  line_pos->SetLineColor(kRed);
53  line_pos->Draw("Same");
54  }
55 
56  // Cylinder
57  if(surface_info.type == 1){
58  vol->Draw();
59  surface->Draw();
60  surface_r->Draw();
61 
62  Dist->GetYaxis()->SetRangeUser(surface_info.range_min-(surface_info.range_max-surface_info.range_min)/20,
63  surface_info.range_max+(surface_info.range_max-surface_info.range_min)/20);
64 
65  // Position of the cylinder surface
66  line_pos = new TLine(-1*surface_info.range_max, surface_info.pos, surface_info.range_max, surface_info.pos);
67  line_pos->SetLineColor(kRed);
68  line_pos->Draw("Same");
69  }
70 
71  c->Print( (out_name+"_Dist.pdf").c_str());
72  //c->Print( (out_name+"_Dist.root").c_str());
73 
74  delete c;
75 
76  delete vol;
77  delete surface;
78  delete surface_z;
79  delete line_pos;
80 }
81 
82 
84 
85 void Initialise_hist(TGraph*& surface_hist,
86  const std::pair<std::vector<float>,std::vector<float>>& surface_pos, const sinfo& surface_info){
87 
88  if(surface_info.type != -1){
89  TGraph * Dist = new TGraph(surface_pos.first.size(), &surface_pos.second[0], &surface_pos.first[0]);
90  Dist->Draw();
91  Dist->GetXaxis()->SetTitle("Z [mm]");
92  Dist->GetYaxis()->SetTitle("R [mm]");
93  surface_hist = Dist;
94  }
95 }
96 
98 
99 void Fill(std::map<uint64_t,TGraph*>& surface_hist, std::map<uint64_t,sinfo>& surface_info,
100  const std::string& input_file, const std::string& json_surface_file, const int& nbprocess){
101 
102  nlohmann::json Det;
103  std::map<std::string,std::string> surface_name;
104 
105  if(json_surface_file != ""){
106  std::ifstream lfile(json_surface_file.c_str());
107  lfile >> Det;
108 
109  Parse_Json(Det, surface_name);
110  }
111 
112  std::map<uint64_t,std::pair<std::vector<float>,std::vector<float>>> surface_pos;
113 
114  //Get file, tree and set top branch address
115  TFile *tfile = new TFile(input_file.c_str());
116  TTree *tree = (TTree*)tfile->Get("material-tracks");
117 
118  std::vector<float> *mat_x = 0;
119  std::vector<float> *mat_y = 0;
120  std::vector<float> *mat_z = 0;
121 
122  std::vector<uint64_t> *sur_id = 0;
123  std::vector<int32_t> *sur_type = 0;
124  std::vector<float> *sur_x = 0;
125  std::vector<float> *sur_y = 0;
126  std::vector<float> *sur_z = 0;
127  std::vector<float> *sur_range_min = 0;
128  std::vector<float> *sur_range_max = 0;
129 
130  tree->SetBranchAddress("mat_x",&mat_x);
131  tree->SetBranchAddress("mat_y",&mat_y);
132  tree->SetBranchAddress("mat_z",&mat_z);
133 
134  tree->SetBranchAddress("sur_id",&sur_id);
135  tree->SetBranchAddress("sur_type",&sur_type);
136  tree->SetBranchAddress("sur_x",&sur_x);
137  tree->SetBranchAddress("sur_y",&sur_y);
138  tree->SetBranchAddress("sur_z",&sur_z);
139  tree->SetBranchAddress("sur_range_min",&sur_range_min);
140  tree->SetBranchAddress("sur_range_max",&sur_range_max);
141 
142  int nentries = tree->GetEntries();
143  if(nentries > nbprocess && nbprocess != -1) nentries = nbprocess;
144  // Limit the number of event processed event to 10000
145  // more could lead to errors with the TGraphs
146  if(nentries > 10000){
147  nentries = 10000;
148  std::cout << "Number of event reduced to 10000" << std::endl;
149  }
150  // Loop over all the material tracks.
151  for (Long64_t i=0;i<nentries; i++) {
152  if(i%1000==0) std::cout << "processed " << i << " events out of " << nentries << std::endl;
153  tree->GetEntry(i);
154 
155  // loop over all the material hits.
156  for(int j=0; j<mat_x->size(); j++ ){
157 
158  // Ignore surface of incorrect type
159  if(sur_type->at(j) == -1) continue;
160 
161  // If a surface was never encountered initialise the position info
162  if(surface_hist.find(sur_id->at(j))==surface_hist.end()){
163 
164  float pos;
165  float range;
166  if(sur_type->at(j) == 1){
167  pos = sqrt(sur_x->at(j)*sur_x->at(j)+sur_y->at(j)*sur_y->at(j));
168  }
169  if(sur_type->at(j) == 2){
170  pos = sur_z->at(j);
171  }
172  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));
173  }
174  // Fill the vector of positions for each layer.
175  surface_pos[sur_id->at(j)].first.push_back(sqrt(mat_y->at(j)*mat_y->at(j)+mat_x->at(j)*mat_x->at(j)));
176  surface_pos[sur_id->at(j)].second.push_back(mat_z->at(j));
177 
178  }
179  }
180  // Use the vector of positions to create the TGraphs
181  for (auto pos_it = surface_pos.begin(); pos_it != surface_pos.end(); pos_it++){
182  Initialise_hist(surface_hist[pos_it->first], pos_it->second, surface_info[pos_it->first]);
183  }
184 }
185 
186 
193 
194 void Mat_map_surface_plot_dist(std::string input_file = "", std::string json_surface_file = "", int nbprocess = -1, std::string name = ""){
195 
196  gStyle->SetOptStat(0);
197  gStyle->SetOptTitle(0);
198 
199  std::map<uint64_t,TGraph*> surface_hist;
200  std::map<uint64_t,sinfo> surface_info;
201  Fill(surface_hist, surface_info, input_file, json_surface_file, nbprocess);
202  for (auto hist_it = surface_hist.begin(); hist_it != surface_hist.end(); hist_it++){
203  if(hist_it->second)
204  plot(hist_it->second, surface_info[hist_it->first], name);
205  }
206 }