Examples/Scripts/momentumDistributions.C

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