ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
momentumDistributions.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file momentumDistributions.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 "TFile.h"
10 #include "TH1F.h"
11 #include "TH2F.h"
12 #include "TProfile.h"
13 #include "TROOT.h"
14 #include "TTree.h"
15 
16 // This root script creates different momentum distributions with the inFile
17 // created from the ExtrapolationTest.
18 // To plot two momentum distributions of this kind in one canvas the root
19 // script "compareDistributions.C" can be used.
20 
21 void
22 momentumDistributions(std::string inFile,
23  std::string treeName,
24  std::string outFile,
25  int nBins,
26  float r,
27  float zMin,
28  float zMax,
29  float etaMin,
30  float etaMax,
31  float thetaMin = 0.,
32  float thetaMax = M_PI)
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  int nHits;
40  float eta;
41 
42  std::vector<float>* x = new std::vector<float>;
43  std::vector<float>* y = new std::vector<float>;
44  std::vector<float>* z = new std::vector<float>;
45 
46  tree->SetBranchAddress("nHits", &nHits);
47  tree->SetBranchAddress("Eta", &eta);
48  tree->SetBranchAddress("StepX", &x);
49  tree->SetBranchAddress("StepY", &y);
50  tree->SetBranchAddress("StepZ", &z);
51 
52  Int_t entries = tree->GetEntries();
53  std::cout << "Creating new output file: " << outFile
54  << " and writing "
55  "material maps"
56  << std::endl;
57  TFile outputFile(outFile.c_str(), "recreate");
58 
59  // distributions of the number of hits versus momentum coordinates
60  TProfile* nHits_eta = new TProfile(
61  "nHits_eta", "Hits in sensitive Material", nBins, etaMin, etaMax);
62  nHits_eta->GetXaxis()->SetTitle("#eta");
63  nHits_eta->GetYaxis()->SetTitle("#hits");
64  TProfile* nHits_theta = new TProfile(
65  "nHits_theta", "Hits in sensitive Material", nBins, thetaMin, thetaMax);
66  nHits_theta->GetXaxis()->SetTitle("#theta [rad]");
67  nHits_theta->GetYaxis()->SetTitle("#hits");
68  TProfile* nHits_z = new TProfile(
69  "nHits_z", "Hits in sensitive Material", nBins, zMin, zMax);
70  nHits_z->GetXaxis()->SetTitle("z coordinate of momentum [mm]");
71  nHits_z->GetYaxis()->SetTitle("#hits");
72 
73  // distributions of the momentum coordinates
74  TH1F* Eta = new TH1F("eta", "Distribution of #eta", nBins, etaMin, etaMax);
75  Eta->GetXaxis()->SetTitle("#eta");
76  Eta->GetYaxis()->SetTitle("#events");
77  // distributions of the momentum coordinates calculated from eta - since in
78  // the extrapolation test eta is flat randomly generated and theta and z are
79  // calculated from eta.
80  TH1F* Theta
81  = new TH1F("theta", "Distribution of #theta", nBins, thetaMin, thetaMax);
82  Theta->GetXaxis()->SetTitle("#theta [rad]");
83  Theta->GetYaxis()->SetTitle("#events");
84  TH1F* Z = new TH1F(
85  "z", "Distribution of z coordinate of the momentum", nBins, zMin, zMax);
86  Z->GetXaxis()->SetTitle("z coordinate of momentum [mm]");
87  Z->GetYaxis()->SetTitle("#events");
88 
89  // hit distributions
90  TH1F* hitsEta = new TH1F(
91  "hitsEta", "Sensitive Hit Distribution", nBins, etaMin, etaMax);
92  hitsEta->GetXaxis()->SetTitle("#eta");
93  hitsEta->GetYaxis()->SetTitle("#hits");
94  TH1F* hitsTheta = new TH1F(
95  "hitsTheta", "Sensitive Hit Distribution", nBins, thetaMin, thetaMax);
96  hitsTheta->GetXaxis()->SetTitle("#theta");
97  hitsTheta->GetYaxis()->SetTitle("#hits");
98  TH1F* hitsZ
99  = new TH1F("hitsZ", "Sensitive Hit Distribution", nBins, zMin, zMax);
100  hitsZ->GetXaxis()->SetTitle("z [mm]");
101  hitsZ->GetYaxis()->SetTitle("#hits");
102 
103  for (int i = 0; i < entries; i++) {
104  tree->GetEvent(i);
105  double theta = 2. * atan(exp(-eta));
106  double zDir = r / tan(theta);
107 
108  nHits_eta->Fill(eta, nHits);
109  nHits_theta->Fill(theta, nHits);
110  nHits_z->Fill(zDir, nHits);
111 
112  Eta->Fill(eta, 1);
113  Theta->Fill(theta);
114  Z->Fill(zDir, 1);
115 
116  for (int j = 0; j < x->size(); j++) {
117  float hitTheta
118  = atan2(sqrt(x->at(j) * x->at(j) + y->at(j) * y->at(j)), z->at(j));
119  hitsEta->Fill(-log(tan(hitTheta * 0.5)));
120  hitsTheta->Fill(hitTheta);
121  hitsZ->Fill(z->at(j));
122  }
123  }
124  inputFile.Close();
125 
126  nHits_eta->Write();
127  nHits_theta->Write();
128  nHits_z->Write();
129 
130  Eta->Write();
131  Theta->Write();
132  Z->Write();
133 
134  hitsEta->Write();
135  hitsTheta->Write();
136  hitsZ->Write();
137 
138  delete nHits_eta;
139  delete nHits_theta;
140  delete nHits_z;
141 
142  delete Eta;
143  delete Theta;
144  delete Z;
145 
146  delete hitsEta;
147  delete hitsTheta;
148  delete hitsZ;
149 
150  delete x;
151  delete y;
152  delete z;
153 
154  outputFile.Close();
155 }