Analysis Software
Documentation for sPHENIX simulation 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 || surface_info.type == 4){
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 int& nbprocess){
101  std::map<std::string,std::string> surface_name;
102  std::map<uint64_t,std::pair<std::vector<float>,std::vector<float>>> surface_pos;
103 
104  //Get file, tree and set top branch address
105  TFile *tfile = new TFile(input_file.c_str());
106  TTree *tree = (TTree*)tfile->Get("material-tracks");
107 
108  std::vector<float> *mat_x = 0;
109  std::vector<float> *mat_y = 0;
110  std::vector<float> *mat_z = 0;
111 
112  std::vector<uint64_t> *sur_id = 0;
113  std::vector<int32_t> *sur_type = 0;
114  std::vector<float> *sur_x = 0;
115  std::vector<float> *sur_y = 0;
116  std::vector<float> *sur_z = 0;
117  std::vector<float> *sur_range_min = 0;
118  std::vector<float> *sur_range_max = 0;
119 
120  tree->SetBranchAddress("mat_x",&mat_x);
121  tree->SetBranchAddress("mat_y",&mat_y);
122  tree->SetBranchAddress("mat_z",&mat_z);
123 
124  tree->SetBranchAddress("sur_id",&sur_id);
125  tree->SetBranchAddress("sur_type",&sur_type);
126  tree->SetBranchAddress("sur_x",&sur_x);
127  tree->SetBranchAddress("sur_y",&sur_y);
128  tree->SetBranchAddress("sur_z",&sur_z);
129  tree->SetBranchAddress("sur_range_min",&sur_range_min);
130  tree->SetBranchAddress("sur_range_max",&sur_range_max);
131 
132  int nentries = tree->GetEntries();
133  if(nentries > nbprocess && nbprocess != -1) nentries = nbprocess;
134  // Limit the number of event processed event to 10000
135  // more could lead to errors with the TGraphs
136  if(nentries > 10000){
137  nentries = 10000;
138  std::cout << "Number of event reduced to 10000" << std::endl;
139  }
140  // Loop over all the material tracks.
141  for (Long64_t i=0;i<nentries; i++) {
142  if(i%1000==0) std::cout << "processed " << i << " events out of " << nentries << std::endl;
143  tree->GetEntry(i);
144 
145  // loop over all the material hits.
146  for(int j=0; j<mat_x->size(); j++ ){
147 
148  // Ignore surface of incorrect type
149  if(sur_type->at(j) == -1) continue;
150 
151  // If a surface was never encountered initialise the position info
152  if(surface_hist.find(sur_id->at(j))==surface_hist.end()){
153 
154  float pos;
155  float range;
156  if(sur_type->at(j) == 1){
157  pos = sqrt(sur_x->at(j)*sur_x->at(j)+sur_y->at(j)*sur_y->at(j));
158  }
159  if(sur_type->at(j) == 2 || sur_type->at(j) == 4){
160  pos = sur_z->at(j);
161  }
162  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));
163  }
164  // Fill the vector of positions for each layer.
165  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)));
166  surface_pos[sur_id->at(j)].second.push_back(mat_z->at(j));
167 
168  }
169  }
170  // Use the vector of positions to create the TGraphs
171  for (auto pos_it = surface_pos.begin(); pos_it != surface_pos.end(); pos_it++){
172  Initialise_hist(surface_hist[pos_it->first], pos_it->second, surface_info[pos_it->first]);
173  }
174 }
175 
176 
180 
181 
182 void Mat_map_surface_plot_dist(std::string input_file = "", int nbprocess = -1, std::string name = ""){
183 
184  gStyle->SetOptStat(0);
185  gStyle->SetOptTitle(0);
186 
187  std::map<uint64_t,TGraph*> surface_hist;
188  std::map<uint64_t,sinfo> surface_info;
189  Fill(surface_hist, surface_info, input_file, nbprocess);
190  for (auto hist_it = surface_hist.begin(); hist_it != surface_hist.end(); hist_it++){
191  if(hist_it->second)
192  plot(hist_it->second, surface_info[hist_it->first], name);
193  }
194 }