Examples/Scripts/momentumDistributions.C

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Utilities/AngleHelpers.hpp"
0010
0011 #include <numbers>
0012
0013 #include "TFile.h"
0014 #include "TH1F.h"
0015 #include "TH2F.h"
0016 #include "TProfile.h"
0017 #include "TROOT.h"
0018 #include "TTree.h"
0019
0020 // This root script creates different momentum distributions with the inFile
0021 // created from the ExtrapolationTest.
0022 // To plot two momentum distributions of this kind  in one canvas the root
0023 // script "compareDistributions.C" can be used.
0024
0025 void momentumDistributions(std::string inFile, std::string treeName,
0026                            std::string outFile, int nBins, float r, float zMin,
0027                            float zMax, float etaMin, float etaMax,
0028                            float thetaMin = 0., float thetaMax = std::numbers::pi_v<float>) {
0029   std::cout << "Opening file: " << inFile << std::endl;
0030   TFile inputFile(inFile.c_str());
0031   std::cout << "Reading tree: " << treeName << std::endl;
0032   TTree* tree = (TTree*)inputFile.Get(treeName.c_str());
0033
0034   int nHits;
0035   float eta;
0036
0037   std::vector<float>* x = new std::vector<float>;
0038   std::vector<float>* y = new std::vector<float>;
0039   std::vector<float>* z = new std::vector<float>;
0040
0041   tree->SetBranchAddress("nHits", &nHits);
0042   tree->SetBranchAddress("Eta", &eta);
0043   tree->SetBranchAddress("StepX", &x);
0044   tree->SetBranchAddress("StepY", &y);
0045   tree->SetBranchAddress("StepZ", &z);
0046
0047   Int_t entries = tree->GetEntries();
0048   std::cout << "Creating new output file: " << outFile
0049             << " and writing "
0050                "material maps"
0051             << std::endl;
0052   TFile outputFile(outFile.c_str(), "recreate");
0053
0054   // distributions of the number of hits versus momentum coordinates
0055   TProfile* nHits_eta = new TProfile("nHits_eta", "Hits in sensitive Material",
0056                                      nBins, etaMin, etaMax);
0057   nHits_eta->GetXaxis()->SetTitle("#eta");
0058   nHits_eta->GetYaxis()->SetTitle("#hits");
0059   TProfile* nHits_theta = new TProfile(
0060       "nHits_theta", "Hits in sensitive Material", nBins, thetaMin, thetaMax);
0061   nHits_theta->GetXaxis()->SetTitle("#theta [rad]");
0062   nHits_theta->GetYaxis()->SetTitle("#hits");
0063   TProfile* nHits_z =
0064       new TProfile("nHits_z", "Hits in sensitive Material", nBins, zMin, zMax);
0065   nHits_z->GetXaxis()->SetTitle("z coordinate of momentum [mm]");
0066   nHits_z->GetYaxis()->SetTitle("#hits");
0067
0068   // distributions of the momentum coordinates
0069   TH1F* Eta = new TH1F("eta", "Distribution of #eta", nBins, etaMin, etaMax);
0070   Eta->GetXaxis()->SetTitle("#eta");
0071   Eta->GetYaxis()->SetTitle("#events");
0072   // distributions of the momentum coordinates calculated from eta - since in
0073   // the extrapolation test eta is flat randomly generated and theta and z are
0074   // calculated from eta.
0075   TH1F* Theta =
0076       new TH1F("theta", "Distribution of #theta", nBins, thetaMin, thetaMax);
0077   Theta->GetXaxis()->SetTitle("#theta [rad]");
0078   Theta->GetYaxis()->SetTitle("#events");
0079   TH1F* Z = new TH1F("z", "Distribution of z coordinate of the momentum", nBins,
0080                      zMin, zMax);
0081   Z->GetXaxis()->SetTitle("z coordinate of momentum [mm]");
0082   Z->GetYaxis()->SetTitle("#events");
0083
0084   // hit distributions
0085   TH1F* hitsEta =
0086       new TH1F("hitsEta", "Sensitive Hit Distribution", nBins, etaMin, etaMax);
0087   hitsEta->GetXaxis()->SetTitle("#eta");
0088   hitsEta->GetYaxis()->SetTitle("#hits");
0089   TH1F* hitsTheta = new TH1F("hitsTheta", "Sensitive Hit Distribution", nBins,
0090                              thetaMin, thetaMax);
0091   hitsTheta->GetXaxis()->SetTitle("#theta");
0092   hitsTheta->GetYaxis()->SetTitle("#hits");
0093   TH1F* hitsZ =
0094       new TH1F("hitsZ", "Sensitive Hit Distribution", nBins, zMin, zMax);
0095   hitsZ->GetXaxis()->SetTitle("z [mm]");
0096   hitsZ->GetYaxis()->SetTitle("#hits");
0097
0098   for (int i = 0; i < entries; i++) {
0099     tree->GetEvent(i);
0100     double theta = 2. * atan(exp(-eta));
0101     double zDir = r / tan(theta);
0102
0103     nHits_eta->Fill(eta, nHits);
0104     nHits_theta->Fill(theta, nHits);
0105     nHits_z->Fill(zDir, nHits);
0106
0107     Eta->Fill(eta, 1);
0108     Theta->Fill(theta);
0109     Z->Fill(zDir, 1);
0110
0111     for (int j = 0; j < x->size(); j++) {
0112       float hitTheta = std::atan2(std::hypot(x->at(j), y->at(j)), z->at(j));
0113       hitsEta->Fill(Acts::AngleHelpers::etaFromTheta(hitTheta));
0114       hitsTheta->Fill(hitTheta);
0115       hitsZ->Fill(z->at(j));
0116     }
0117   }
0118   inputFile.Close();
0119
0120   nHits_eta->Write();
0121   nHits_theta->Write();
0122   nHits_z->Write();
0123
0124   Eta->Write();
0125   Theta->Write();
0126   Z->Write();
0127
0128   hitsEta->Write();
0129   hitsTheta->Write();
0130   hitsZ->Write();
0131
0132   delete nHits_eta;
0133   delete nHits_theta;
0134   delete nHits_z;
0135
0136   delete Eta;
0137   delete Theta;
0138   delete Z;
0139
0140   delete hitsEta;
0141   delete hitsTheta;
0142   delete hitsZ;
0143
0144   delete x;
0145   delete y;
0146   delete z;
0147
0148   outputFile.Close();
0149 }