ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
printHits.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file printHits.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 #include <TTreeReader.h>
10 #include <TTreeReaderValue.h>
11 #include "TFile.h"
12 #include "TH1F.h"
13 #include "TH2F.h"
14 #include "TProfile.h"
15 #include "TROOT.h"
16 #include "TTree.h"
17 
18 // This root script creates different histograms displaying the sensitive
19 // material, the boundaries and the material of the detector in different views.
20 // The in file needs to be in the format created form the ExtrapolationTest.
21 // To plot two or three histograms of this kind in one canvas the root script
22 // "compareHistogram.C" can be used.
23 
24 void
25 printHits(std::string inFile,
26  std::string treeName,
27  std::string outFile,
28  float rmin,
29  float rmax,
30  float zmin,
31  float zmax,
32  int nBins)
33 {
34  std::cout << "Opening file: " << inFile << std::endl;
35  TFile inputFile(inFile.c_str());
36  std::cout << "Reading tree: " << treeName << std::endl;
37  TTree* tree = (TTree*)inputFile.Get(treeName.c_str());
38 
39  TTreeReader reader(treeName.c_str(), &inputFile);
40 
41  int nHits;
42  float eta;
43  float theta;
44  float zDir;
45 
46  std::vector<float>* x = new std::vector<float>;
47  std::vector<float>* y = new std::vector<float>;
48  std::vector<float>* z = new std::vector<float>;
49  std::vector<float>* sens = new std::vector<float>;
50  std::vector<float>* mat = new std::vector<float>;
51  std::vector<float>* bounds = new std::vector<float>;
52 
53  tree->SetBranchAddress("step_x", &x);
54  tree->SetBranchAddress("step_y", &y);
55  tree->SetBranchAddress("step_z", &z);
56  tree->SetBranchAddress("sensitive", &sens);
57 
58  if (tree->FindBranch("boundary")) {
59  std::cout << "No BoundarySteps are given." << std::endl;
60  tree->SetBranchAddress("boundary", &bounds);
61  }
62  if (tree->FindBranch("material")) {
63  std::cout << "No MaterialSteps are given." << std::endl;
64  tree->SetBranchAddress("material", &mat);
65  }
66  Int_t entries = tree->GetEntries();
67  std::cout << "Creating new output file: " << outFile
68  << " and writing out histograms. " << std::endl;
69  TFile outputFile(outFile.c_str(), "recreate");
70 
71  // full
72  TH2F* Full_xy = new TH2F(
73  "Full_xy", "Full material", nBins, rmin, rmax, nBins, rmin, rmax);
74  TH2F* Full_zr = new TH2F(
75  "Full_zr", "Full material", nBins, zmin, zmax, nBins, 0., rmax);
76 
77  // sensitive
78  TH2F* Sens_xy = new TH2F(
79  "Sens_xy", "Sensitive material", nBins, rmin, rmax, nBins, rmin, rmax);
80  Sens_xy->GetXaxis()->SetTitle("x [mm]");
81  Sens_xy->GetYaxis()->SetTitle("y [mm]");
82  TH2F* Sens_zx = new TH2F(
83  "Sens_zx", "Sensitive material", nBins, zmin, zmax, nBins, rmin, rmax);
84  Sens_zx->GetXaxis()->SetTitle("z [mm]");
85  Sens_zx->GetYaxis()->SetTitle("x [mm]");
86  TH2F* Sens_zy = new TH2F(
87  "Sens_zy", "Sensitive material", nBins, zmin, zmax, nBins, rmin, rmax);
88  Sens_zy->GetXaxis()->SetTitle("z [mm]");
89  Sens_zy->GetYaxis()->SetTitle("y [mm]");
90  TH2F* Sens_zr = new TH2F(
91  "Sens_zr", "Sensitive material", nBins, zmin, zmax, nBins, 0., rmax);
92  Sens_zr->GetXaxis()->SetTitle("z [mm]");
93  Sens_zr->GetYaxis()->SetTitle("r [mm]");
94 
95  // boundaries
96  TH2F* Bounds_xy = new TH2F(
97  "Bounds_xy", "Boundaries", nBins, rmin, rmax, nBins, rmin, rmax);
98  Bounds_xy->GetXaxis()->SetTitle("x [mm]");
99  Bounds_xy->GetYaxis()->SetTitle("y [mm]");
100  TH2F* Bounds_zx = new TH2F(
101  "Bounds_zx", "Boundaries", nBins, zmin, zmax, nBins, rmin, rmax);
102  Bounds_zx->GetXaxis()->SetTitle("z [mm]");
103  Bounds_zx->GetYaxis()->SetTitle("x [mm]");
104  TH2F* Bounds_zy = new TH2F(
105  "Bounds_zy", "Boundaries", nBins, zmin, zmax, nBins, rmin, rmax);
106  Bounds_zy->GetXaxis()->SetTitle("z [mm]");
107  Bounds_zy->GetYaxis()->SetTitle("y [mm]");
108  TH2F* Bounds_zr
109  = new TH2F("Bounds_zr", "Boundaries", nBins, zmin, zmax, nBins, 0., rmax);
110  Bounds_zr->GetXaxis()->SetTitle("z [mm]");
111  Bounds_zr->GetYaxis()->SetTitle("r [mm]");
112 
113  // material
114  TH2F* Mat_xy
115  = new TH2F("Mat_xy", "Material", nBins, rmin, rmax, nBins, rmin, rmax);
116  Mat_xy->GetXaxis()->SetTitle("x [mm]");
117  Mat_xy->GetYaxis()->SetTitle("y [mm]");
118  TH2F* Mat_zx
119  = new TH2F("Mat_zx", "Material", nBins, zmin, zmax, nBins, rmin, rmax);
120  Mat_zx->GetXaxis()->SetTitle("z [mm]");
121  Mat_zx->GetYaxis()->SetTitle("x [mm]");
122  TH2F* Mat_zy
123  = new TH2F("Mat_zy", "Material", nBins, zmin, zmax, nBins, rmin, rmax);
124  Mat_zy->GetXaxis()->SetTitle("z [mm]");
125  Mat_zy->GetYaxis()->SetTitle("y [mm]");
126  TH2F* Mat_zr
127  = new TH2F("Mat_zr", "Material", nBins, zmin, zmax, nBins, 0., rmax);
128  Mat_zr->GetXaxis()->SetTitle("z [mm]");
129  Mat_zr->GetYaxis()->SetTitle("r [mm]");
130 
131  for (int i = 0; i < entries; i++) {
132  tree->GetEvent(i);
133 
134  for (int j = 0; j < x->size(); j++) {
135  // sensitive
136  if (z->at(j) >= zmin && z->at(j) <= zmax) {
137  Sens_xy->Fill(x->at(j), y->at(j), sens->at(j));
138  Full_xy->Fill(x->at(j), y->at(j));
139  }
140  Sens_zx->Fill(z->at(j), x->at(j), sens->at(j));
141  Sens_zy->Fill(z->at(j), y->at(j), sens->at(j));
142  Sens_zr->Fill(z->at(j), sqrt(x->at(j) * x->at(j) + y->at(j) * y->at(j)));
143  Full_zr->Fill(z->at(j), sqrt(x->at(j) * x->at(j) + y->at(j) * y->at(j)));
144 
145  // boundaries
146  if (tree->FindBranch("boundary")) {
147  if (z->at(j) >= zmin && z->at(j) <= zmax)
148  Bounds_xy->Fill(x->at(j), y->at(j), bounds->at(j));
149  Bounds_zx->Fill(z->at(j), x->at(j), bounds->at(j));
150  Bounds_zy->Fill(z->at(j), y->at(j), bounds->at(j));
151  Bounds_zr->Fill(z->at(j),
152  sqrt(x->at(j) * x->at(j) + y->at(j) * y->at(j)),
153  bounds->at(j));
154  }
155  // material
156  if (tree->FindBranch("material")) {
157  if (z->at(j) >= zmin && z->at(j) <= zmax)
158  Mat_xy->Fill(x->at(j), y->at(j), mat->at(j));
159  Mat_zx->Fill(z->at(j), x->at(j), mat->at(j));
160  Mat_zy->Fill(z->at(j), y->at(j), mat->at(j));
161  Mat_zr->Fill(z->at(j),
162  sqrt(x->at(j) * x->at(j) + y->at(j) * y->at(j)),
163  mat->at(j));
164  }
165  }
166  }
167  inputFile.Close();
168 
169  // full
170  Full_xy->Write();
171  delete Full_xy;
172  Full_zr->Write();
173  delete Full_zr;
174 
175  // sensitive
176  Sens_xy->Write();
177  delete Sens_xy;
178  Sens_zx->Write();
179  delete Sens_zx;
180  Sens_zy->Write();
181  delete Sens_zy;
182  Sens_zr->Write();
183  delete Sens_zr;
184 
185  // boundaries
186  Bounds_xy->Write();
187  delete Bounds_xy;
188  Bounds_zx->Write();
189  delete Bounds_zx;
190  Bounds_zy->Write();
191  delete Bounds_zy;
192  Bounds_zr->Write();
193  delete Bounds_zr;
194 
195  // material
196  Mat_xy->Write();
197  delete Mat_xy;
198  Mat_zx->Write();
199  delete Mat_zx;
200  Mat_zy->Write();
201  delete Mat_zy;
202  Mat_zr->Write();
203  delete Mat_zr;
204 
205  delete x;
206  delete y;
207  delete z;
208  delete sens;
209  delete bounds;
210  delete mat;
211 
212  outputFile.Close();
213 }