ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fullMaterial.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file fullMaterial.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  * fullMaterial.cxx
11  *
12  * Created on: 22 Aug 2016
13  * Author: jhrdinka
14  */
15 
16 #include "TFile.h"
17 #include "TH1F.h"
18 #include "TH2F.h"
19 #include "TROOT.h"
20 #include "TTree.h"
21 
22 // This root script writes out the full material described in
23 // MaterialTracks. It sums up all the material accumulated along one track
24 // and averages the material over all tracks. It generates Histograms of the
25 // different MaterialProperties along phi, theta,z and eta.
26 // It is forseen to be used with an input file generated by the
27 // materialTrackWriter.
28 
29 void
30 fullMaterial(std::string inFile,
31  std::string treename,
32  std::string outFile,
33  int binsZ,
34  float minZ,
35  float maxZ,
36  int binsPhi,
37  float minPhi,
38  float maxPhi,
39  int binsTheta,
40  float minTheta,
41  float maxTheta,
42  int binsEta,
43  float minEta,
44  float maxEta)
45 {
46  std::cout << "Opening file: " << inFile << std::endl;
47  TFile inputFile(inFile.c_str());
48  std::cout << "Reading tree: " << treename << std::endl;
49  TTreeReader reader(treename.c_str(), &inputFile);
50 
51  // get the MaterialTrack entities
52  std::vector<Acts::MaterialTrack> mrecords;
53  std::cout << "Accessing Branch 'MaterialTracks'" << std::endl;
54  TTreeReaderValue<Acts::MaterialTrack> mtRecord(reader, "MaterialTracks");
55  while (reader.Next()) { mrecords.push_back(*mtRecord); }
56  inputFile.Close();
57 
58  std::cout << "Creating new output file: " << outFile
59  << " and writing "
60  "material maps"
61  << std::endl;
62  TFile outputFile(outFile.c_str(), "recreate");
63  // thickness
64  TProfile* t_phi = new TProfile("t_phi", "t_phi", binsPhi, minPhi, maxPhi);
65  TProfile* t_eta = new TProfile("t_eta", "t_eta", binsEta, minEta, maxEta);
66  TProfile* t_theta
67  = new TProfile("t_theta", "t_theta", binsTheta, minTheta, maxTheta);
68 
69  // thickness in X0
70  TProfile* tInX0_phi
71  = new TProfile("tInX0_phi", "tInX0_phi", binsPhi, minPhi, maxPhi);
72  TProfile* tInX0_eta
73  = new TProfile("tInX0_eta", "tInX0_eta", binsEta, minEta, maxEta);
74  TProfile* tInX0_theta = new TProfile(
75  "tInX0_theta", "tInX0_theta", binsTheta, minTheta, maxTheta);
76  // thickness in L0
77  TProfile* tInL0_phi
78  = new TProfile("tInL0_phi", "tInL0_phi", binsPhi, minPhi, maxPhi);
79  TProfile* tInL0_eta
80  = new TProfile("tInL0_eta", "tInL0_eta", binsEta, minEta, maxEta);
81  TProfile* tInL0_theta = new TProfile(
82  "tInL0_theta", "tInL0_theta", binsTheta, minTheta, maxTheta);
83  // A
84  TProfile* A_phi = new TProfile("A_phi", "A_phi", binsPhi, minPhi, maxPhi);
85  TProfile* A_eta = new TProfile("A_eta", "A_eta", binsEta, minEta, maxEta);
86  TProfile* A_theta
87  = new TProfile("A_theta", "A_theta", binsTheta, minTheta, maxTheta);
88  // Z
89  TProfile* Z_phi = new TProfile("Z_phi", "Z_phi", binsPhi, minPhi, maxPhi);
90  TProfile* Z_eta = new TProfile("Z_eta", "Z_eta", binsEta, minEta, maxEta);
91  TProfile* Z_theta
92  = new TProfile("Z_theta", "Z_theta", binsTheta, minTheta, maxTheta);
93  // Rho
94  TProfile* rho_phi
95  = new TProfile("rho_phi", "rho_phi", binsPhi, minPhi, maxPhi);
96  TProfile* rho_eta
97  = new TProfile("rho_eta", "rho_eta", binsEta, minEta, maxEta);
98  TProfile* rho_theta
99  = new TProfile("rho_theta", "rho_theta", binsTheta, minTheta, maxTheta);
100  // x0
101  TProfile* x0_phi = new TProfile("x0_phi", "x0_phi", binsPhi, minPhi, maxPhi);
102  TProfile* x0_eta = new TProfile("x0_eta", "x0_eta", binsEta, minEta, maxEta);
103  TProfile* x0_theta
104  = new TProfile("x0_theta", "x0_theta", binsTheta, minTheta, maxTheta);
105  // l0
106  TProfile* l0_phi = new TProfile("l0_phi", "l0_phi", binsPhi, minPhi, maxPhi);
107  TProfile* l0_eta = new TProfile("l0_eta", "l0_eta", binsEta, minEta, maxEta);
108  TProfile* l0_theta
109  = new TProfile("l0_theta", "l0_theta", binsTheta, minTheta, maxTheta);
110 
111  for (auto& mrecord : mrecords) {
112  std::vector<Acts::MaterialStep> steps = mrecord.materialSteps();
113  float theta = mrecord.theta();
114  float eta = -log(tan(theta * 0.5));
115  float phi = VectorHelpers::phi(mrecord);
116 
117  float thickness = 0.;
118  float rho = 0.;
119  float A = 0.;
120  float Z = 0.;
121  float tInX0 = 0.;
122  float tInL0 = 0.;
123 
124  for (auto& step : steps) {
125  // std::cout << "t" << step.material().thickness() << std::endl;
126  float t = step.material().thickness();
127  float r = step.material().material().massDensity();
128  thickness += step.material().thickness();
129  rho += r * t;
130  A += step.material().material().Ar() * r * t;
131  Z += step.material().material().Z() * r * t;
132  tInX0 += (t != 0. && step.material().x0() != 0.)
133  ? t / step.material().x0()
134  : 0.;
135 
136  tInL0 += (t != 0. && step.material().x0() != 0.)
137  ? t / step.material().l0()
138  : 0.;
139  }
140  if (rho != 0.) {
141  A /= rho;
142  Z /= rho;
143  }
144  if (thickness != 0.) rho /= thickness;
145 
146  t_phi->Fill(phi, thickness);
147  t_theta->Fill(theta, thickness);
148  t_eta->Fill(eta, thickness);
149 
150  if (tInX0 != 0.) {
151  x0_phi->Fill(phi, thickness / tInX0);
152  x0_theta->Fill(theta, thickness / tInX0);
153  x0_eta->Fill(eta, thickness / tInX0);
154  }
155 
156  if (tInL0 != 0.) {
157  l0_phi->Fill(phi, thickness / tInL0);
158  l0_theta->Fill(theta, thickness / tInL0);
159  l0_eta->Fill(eta, thickness / tInL0);
160  }
161 
162  tInX0_phi->Fill(phi, tInX0);
163  tInX0_theta->Fill(theta, tInX0);
164  tInX0_eta->Fill(eta, tInX0);
165 
166  tInL0_phi->Fill(phi, tInL0);
167  tInL0_theta->Fill(theta, tInL0);
168  tInL0_eta->Fill(eta, tInL0);
169 
170  A_phi->Fill(phi, A);
171  A_theta->Fill(theta, A);
172  A_eta->Fill(eta, A);
173 
174  Z_phi->Fill(phi, Z);
175  Z_theta->Fill(theta, Z);
176  Z_eta->Fill(eta, Z);
177 
178  rho_phi->Fill(phi, rho);
179  rho_theta->Fill(theta, rho);
180  rho_eta->Fill(eta, rho);
181  }
182  // thickness
183  t_phi->Write();
184  t_theta->Write();
185  t_eta->Write();
186  delete t_theta;
187  delete t_eta;
188  delete t_phi;
189  // thickness in X0
190  tInX0_phi->Write();
191  tInX0_theta->Write();
192  tInX0_eta->Write();
193  delete tInX0_phi;
194  delete tInX0_eta;
195  delete tInX0_theta;
196  // thickness in L0
197  tInL0_phi->Write();
198  tInL0_theta->Write();
199  tInL0_eta->Write();
200  delete tInL0_phi;
201  delete tInL0_eta;
202  delete tInL0_theta;
203  // A
204  A_phi->Write();
205  A_eta->Write();
206  A_theta->Write();
207  delete A_phi;
208  delete A_eta;
209  delete A_theta;
210  // Z
211  Z_phi->Write();
212  Z_eta->Write();
213  Z_theta->Write();
214  delete Z_phi;
215  delete Z_eta;
216  delete Z_theta;
217  // rho
218  rho_phi->Write();
219  rho_eta->Write();
220  rho_theta->Write();
221  delete rho_phi;
222  delete rho_eta;
223  delete rho_theta;
224  // x0
225  x0_phi->Write();
226  x0_eta->Write();
227  x0_theta->Write();
228  delete x0_phi;
229  delete x0_eta;
230  delete x0_theta;
231  // l0
232  l0_phi->Write();
233  l0_eta->Write();
234  l0_theta->Write();
235  delete l0_phi;
236  delete l0_eta;
237  delete l0_theta;
238 
239  outputFile.Close();
240 }