Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
layerMaterial.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file layerMaterial.C
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2017 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 //
10 // layerMaterial.C
11 //
12 //
13 // Created by Julia Hrdinka on 18/08/16.
14 //
15 //
16 
17 #include <algorithm>
18 #include <iostream>
19 #include <vector>
20 #include "TFile.h"
21 #include "TH1F.h"
22 #include "TObject.h"
23 #include "TROOT.h"
24 #include "TTree.h"
25 
26 // This root script generates material histograms and maps of a layer.
27 // It is foreseen to be used for output generated by the RootMaterialWriter or
28 // the RootMaterialStepWriter. For the indicated Layer it produces the
29 // histograms and maps of the material properties along the local coordinates of
30 // the layer.
31 // It also writes the assigned and real global positions into histograms, if
32 // given.
33 
34 void
35 layerMaterial(std::string inFile,
36  std::string outFile,
37  int binsLoc1,
38  int binsLoc2,
39  int binsZ,
40  float minGlobZ,
41  float maxGlobZ,
42  int binsR,
43  float minGlobR,
44  float maxGlobR,
45  std::string layerName)
46 {
47  std::cout << "Opening file: " << inFile << std::endl;
48  TFile inputFile(inFile.c_str());
49  TTree* layer = (TTree*)inputFile.Get(layerName.c_str());
50  std::cout << "Reading tree: " << layerName << std::endl;
51 
52  std::vector<float>* loc0 = new std::vector<float>;
53  std::vector<float>* loc1 = new std::vector<float>;
54  std::vector<float>* A = new std::vector<float>;
55  std::vector<float>* Z = new std::vector<float>;
56  std::vector<float>* x0 = new std::vector<float>;
57  std::vector<float>* l0 = new std::vector<float>;
58  std::vector<float>* d = new std::vector<float>;
59  std::vector<float>* rho = new std::vector<float>;
60  std::vector<float>* dInX0 = new std::vector<float>;
61  std::vector<float>* dInL0 = new std::vector<float>;
62  std::vector<float>* t = new std::vector<float>;
63 
64  std::vector<float>* globX = new std::vector<float>;
65  std::vector<float>* globY = new std::vector<float>;
66  std::vector<float>* globZ = new std::vector<float>;
67  std::vector<float>* globR = new std::vector<float>;
68 
69  std::vector<float>* assignedGlobX = new std::vector<float>;
70  std::vector<float>* assignedGlobY = new std::vector<float>;
71  std::vector<float>* assignedGlobZ = new std::vector<float>;
72  std::vector<float>* assignedGlobR = new std::vector<float>;
73 
74  layer->SetBranchAddress("loc0", &loc0);
75  layer->SetBranchAddress("loc1", &loc1);
76  layer->SetBranchAddress("A", &A);
77  layer->SetBranchAddress("Z", &Z);
78  layer->SetBranchAddress("x0", &x0);
79  layer->SetBranchAddress("l0", &l0);
80  layer->SetBranchAddress("thickness", &d);
81  layer->SetBranchAddress("rho", &rho);
82  layer->SetBranchAddress("tInX0", &dInX0);
83  layer->SetBranchAddress("tInL0", &dInL0);
84  layer->SetBranchAddress("thickness", &t);
85 
86  if (layer->FindBranch("globX") && layer->FindBranch("globY")
87  && layer->FindBranch("globZ") && layer->FindBranch("globR")) {
88  layer->SetBranchAddress("globX", &globX);
89  layer->SetBranchAddress("globY", &globY);
90  layer->SetBranchAddress("globZ", &globZ);
91  layer->SetBranchAddress("globR", &globR);
92  }
93 
94  if (layer->FindBranch("assignedGlobX") && layer->FindBranch("assignedGlobY")
95  && layer->FindBranch("assignedGlobZ")
96  && layer->FindBranch("assignedGlobR")) {
97  layer->SetBranchAddress("assignedGlobX", &assignedGlobX);
98  layer->SetBranchAddress("assignedGlobY", &assignedGlobY);
99  layer->SetBranchAddress("assignedGlobZ", &assignedGlobZ);
100  layer->SetBranchAddress("assignedGlobR", &assignedGlobR);
101  }
102  layer->GetEntry(0);
103 
104  // get minima and maxima for different layers
105  auto minmax0 = std::minmax_element(loc0->begin(), loc0->end());
106  float min0 = *minmax0.first;
107  float max0 = *minmax0.second;
108 
109  auto minmax1 = std::minmax_element(loc1->begin(), loc1->end());
110  float min1 = *minmax1.first;
111  float max1 = *minmax1.second;
112 
113  inputFile.Close();
114  std::cout << "Creating new output file: " << outFile
115  << " and writing "
116  "material maps"
117  << std::endl;
118  TFile outputFile(outFile.c_str(), "update");
119  TDirectory* dir = outputFile.mkdir(layerName.c_str());
120  dir->cd();
121  // thickness in X0
122  TProfile* dInX0_loc1
123  = new TProfile("dInX0_loc1", "dInX0_loc1", binsLoc1, min1, max1);
124  TProfile* dInX0_loc0
125  = new TProfile("dInX0_loc0", "dInX0_loc0", binsLoc1, min0, max0);
126  TProfile2D* dInX0_map = new TProfile2D(
127  "dInX0", "dInX0", binsLoc1, min1, max1, binsLoc1, min0, max0);
128  // thickness in L0
129  TProfile* dInL0_loc1
130  = new TProfile("dInL0_loc1", "dInL0_loc1", binsLoc2, min1, max1);
131  TProfile* dInL0_loc0
132  = new TProfile("dInL0_loc0", "dInL0_loc0", binsLoc1, min0, max0);
133  TProfile2D* dInL0_map = new TProfile2D(
134  "dInL0", "dInL0", binsLoc2, min1, max1, binsLoc1, min0, max0);
135  // A
136  TProfile* A_loc1 = new TProfile("A_loc1", "A_loc1", binsLoc2, min1, max1);
137  TProfile* A_loc0 = new TProfile("A_loc0", "A_loc0", binsLoc1, min0, max0);
138  TProfile2D* A_map
139  = new TProfile2D("A", "A", binsLoc2, min1, max1, binsLoc1, min0, max0);
140  // Z
141  TProfile* Z_loc1 = new TProfile("Z_loc1", "Z_loc1", binsLoc2, min1, max1);
142  TProfile* Z_loc0 = new TProfile("Z_loc0", "Z_loc0", binsLoc1, min0, max0);
143  TProfile2D* Z_map
144  = new TProfile2D("Z", "Z", binsLoc2, min1, max1, binsLoc1, min0, max0);
145  // Rho
146  TProfile* rho_loc1
147  = new TProfile("rho_loc1", "rho_loc1", binsLoc2, min1, max1);
148  TProfile* rho_loc0
149  = new TProfile("rho_loc0", "rho_loc0", binsLoc1, min0, max0);
150  TProfile2D* rho_map = new TProfile2D(
151  "rho", "rho", binsLoc2, min1, max1, binsLoc1, min0, max0);
152  // x0
153  TProfile* x0_loc1 = new TProfile("x0_loc1", "x0_loc1", binsLoc2, min1, max1);
154  TProfile* x0_loc0 = new TProfile("x0_loc0", "x0_loc0", binsLoc1, min0, max0);
155  TProfile2D* x0_map
156  = new TProfile2D("x0", "x0", binsLoc2, min1, max1, binsLoc1, min0, max0);
157  // l0
158  TProfile* l0_loc1 = new TProfile("l0_loc1", "l0_loc1", binsLoc2, min1, max1);
159  TProfile* l0_loc0 = new TProfile("l0_loc0", "l0_loc0", binsLoc1, min0, max0);
160  TProfile2D* l0_map
161  = new TProfile2D("l0", "l0", binsLoc2, min1, max1, binsLoc1, min0, max0);
162 
163  // thickness
164  TProfile* t_loc1 = new TProfile("t_loc1", "t_loc1", binsLoc2, min1, max1);
165  TProfile* t_loc0 = new TProfile("t_loc0", "t_loc0", binsLoc1, min0, max0);
166  TProfile2D* t_map
167  = new TProfile2D("t", "t", binsLoc2, min1, max1, binsLoc1, min0, max0);
168 
169  // global
170  TH2F* glob_r_z = new TH2F(
171  "r_z", "r_z", binsZ, minGlobZ, maxGlobZ, binsR, minGlobR, maxGlobR);
172  TH2F* assigned_r_z = new TH2F("r_z_assigned",
173  "r_z_assigned",
174  binsZ,
175  minGlobZ,
176  maxGlobZ,
177  binsR,
178  minGlobR,
179  maxGlobR);
180 
181  TH2F* glob_x_y = new TH2F(
182  "x_y", "x_y", binsR, -maxGlobR, maxGlobR, binsR, -maxGlobR, maxGlobR);
183  TH2F* assigned_x_y = new TH2F("x_y_assigned",
184  "x_y_assigned",
185  binsR,
186  -maxGlobR,
187  maxGlobR,
188  binsR,
189  -maxGlobR,
190  maxGlobR);
191 
192  size_t nEntries = loc1->size();
193  for (int i = 0; i < nEntries; i++) {
194  // A
195  A_loc1->Fill(loc1->at(i), A->at(i));
196  A_loc0->Fill(loc0->at(i), A->at(i));
197  A_map->Fill(loc1->at(i), loc0->at(i), A->at(i));
198  // Z
199  Z_loc1->Fill(loc1->at(i), Z->at(i));
200  Z_loc0->Fill(loc0->at(i), Z->at(i));
201  Z_map->Fill(loc1->at(i), loc0->at(i), Z->at(i));
202  // x0
203  x0_loc1->Fill(loc1->at(i), x0->at(i));
204  x0_loc0->Fill(loc0->at(i), x0->at(i));
205  x0_map->Fill(loc1->at(i), loc0->at(i), x0->at(i));
206  // l0
207  l0_loc1->Fill(loc1->at(i), l0->at(i));
208  l0_loc0->Fill(loc0->at(i), l0->at(i));
209  l0_map->Fill(loc1->at(i), loc0->at(i), l0->at(i));
210  // rho
211  rho_loc1->Fill(loc1->at(i), rho->at(i));
212  rho_loc0->Fill(loc0->at(i), rho->at(i));
213  rho_map->Fill(loc1->at(i), loc0->at(i), rho->at(i));
214  // thickness in X0
215  dInX0_loc1->Fill(loc1->at(i), dInX0->at(i));
216  dInX0_loc0->Fill(loc0->at(i), dInX0->at(i));
217  dInX0_map->Fill(loc1->at(i), loc0->at(i), dInX0->at(i));
218  // thickness in L0
219  dInL0_loc1->Fill(loc1->at(i), dInL0->at(i));
220  dInL0_loc0->Fill(loc0->at(i), dInL0->at(i));
221  dInL0_map->Fill(loc1->at(i), loc0->at(i), dInL0->at(i));
222  // thickness
223  t_loc1->Fill(loc1->at(i), t->at(i));
224  t_loc0->Fill(loc0->at(i), t->at(i));
225  t_map->Fill(loc1->at(i), loc0->at(i), t->at(i));
226 
227  // fill global r/z
228  if (globZ->size() && globR->size())
229  glob_r_z->Fill(globZ->at(i), globR->at(i));
230 
231  if (assignedGlobZ->size() && assignedGlobR->size())
232  assigned_r_z->Fill(assignedGlobZ->at(i), assignedGlobR->at(i));
233 
234  // fill global x/y
235  if (globX->size() && globY->size())
236  glob_x_y->Fill(globX->at(i), globY->at(i));
237 
238  if (assignedGlobX->size() && assignedGlobY->size())
239  assigned_x_y->Fill(assignedGlobX->at(i), assignedGlobY->at(i));
240  }
241 
242  gStyle->SetOptStat(0);
243  // A
244  A_loc1->Write();
245  A_loc0->Write();
246  A_map->Write();
247  delete A_loc1;
248  delete A_loc0;
249  delete A_map;
250  // Z
251  Z_loc1->Write();
252  Z_loc0->Write();
253  Z_map->Write();
254  delete Z_loc1;
255  delete Z_loc0;
256  delete Z_map;
257  // x0
258  x0_loc1->Write();
259  x0_loc0->Write();
260  x0_map->Write();
261  delete x0_loc1;
262  delete x0_loc0;
263  delete x0_map;
264  // l0
265  l0_loc1->Write();
266  l0_loc0->Write();
267  l0_map->Write();
268  delete l0_loc1;
269  delete l0_loc0;
270  delete l0_map;
271  // rho
272  rho_loc1->Write();
273  rho_loc0->Write();
274  rho_map->Write();
275  delete rho_loc1;
276  delete rho_loc0;
277  delete rho_map;
278  // thickness in X0
279  dInX0_loc1->Write();
280  dInX0_loc0->Write();
281  dInX0_map->Write();
282  delete dInX0_loc1;
283  delete dInX0_loc0;
284  delete dInX0_map;
285  // thickness in L0
286  dInL0_loc1->Write();
287  dInL0_loc0->Write();
288  dInL0_map->Write();
289  delete dInL0_loc1;
290  delete dInL0_loc0;
291  delete dInL0_map;
292  // thickness
293  t_loc1->Write();
294  t_loc0->Write();
295  t_map->Write();
296 
297  delete loc0;
298  delete loc1;
299  delete A;
300  delete Z;
301  delete x0;
302  delete l0;
303  delete d;
304  delete rho;
305  delete dInX0;
306  delete dInL0;
307  delete t;
308 
309  glob_r_z->Write();
310  assigned_r_z->Write();
311  glob_x_y->Write();
312  assigned_x_y->Write();
313  delete glob_r_z;
314  delete assigned_r_z;
315  delete glob_x_y;
316  delete assigned_x_y;
317  outputFile.Close();
318 }