Examples/Scripts/fullMaterial.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 /*
0010  * fullMaterial.cxx
0011  *
0012  *  Created on: 22 Aug 2016
0013  *      Author: jhrdinka
0014  */
0015
0016 #include "Acts/Utilities/AngleHelpers.hpp"
0017
0018 #include "TFile.h"
0019 #include "TH1F.h"
0020 #include "TH2F.h"
0021 #include "TROOT.h"
0022 #include "TTree.h"
0023
0024 // This root script writes out the full material described in
0025 // MaterialTracks. It sums up all the material accumulated along one track
0026 // and averages the material over all tracks. It generates Histograms of the
0027 // different MaterialSlab along phi, theta,z and eta.
0028 // It is forseen to be used with an input file generated by the
0029 // materialTrackWriter.
0030
0031 void
0032 fullMaterial(std::string inFile,
0033              std::string treename,
0034              std::string outFile,
0035              int         binsZ,
0036              float       minZ,
0037              float       maxZ,
0038              int         binsPhi,
0039              float       minPhi,
0040              float       maxPhi,
0041              int         binsTheta,
0042              float       minTheta,
0043              float       maxTheta,
0044              int         binsEta,
0045              float       minEta,
0046              float       maxEta)
0047 {
0048   std::cout << "Opening file: " << inFile << std::endl;
0049   TFile inputFile(inFile.c_str());
0050   std::cout << "Reading tree: " << treename << std::endl;
0051   TTreeReader reader(treename.c_str(), &inputFile);
0052
0053   // get the MaterialTrack entities
0054   std::vector<Acts::MaterialTrack> mrecords;
0055   std::cout << "Accessing Branch 'MaterialTracks'" << std::endl;
0056   TTreeReaderValue<Acts::MaterialTrack> mtRecord(reader, "MaterialTracks");
0057   while (reader.Next()) { mrecords.push_back(*mtRecord); }
0058   inputFile.Close();
0059
0060   std::cout << "Creating new output file: " << outFile
0061             << " and writing "
0062                "material maps"
0063             << std::endl;
0064   TFile outputFile(outFile.c_str(), "recreate");
0065   // thickness
0066   TProfile* t_phi = new TProfile("t_phi", "t_phi", binsPhi, minPhi, maxPhi);
0067   TProfile* t_eta = new TProfile("t_eta", "t_eta", binsEta, minEta, maxEta);
0068   TProfile* t_theta
0069       = new TProfile("t_theta", "t_theta", binsTheta, minTheta, maxTheta);
0070
0071   // thickness in X0
0072   TProfile* tInX0_phi
0073       = new TProfile("tInX0_phi", "tInX0_phi", binsPhi, minPhi, maxPhi);
0074   TProfile* tInX0_eta
0075       = new TProfile("tInX0_eta", "tInX0_eta", binsEta, minEta, maxEta);
0076   TProfile* tInX0_theta = new TProfile(
0077       "tInX0_theta", "tInX0_theta", binsTheta, minTheta, maxTheta);
0078   // thickness in L0
0079   TProfile* tInL0_phi
0080       = new TProfile("tInL0_phi", "tInL0_phi", binsPhi, minPhi, maxPhi);
0081   TProfile* tInL0_eta
0082       = new TProfile("tInL0_eta", "tInL0_eta", binsEta, minEta, maxEta);
0083   TProfile* tInL0_theta = new TProfile(
0084       "tInL0_theta", "tInL0_theta", binsTheta, minTheta, maxTheta);
0085   // A
0086   TProfile* A_phi = new TProfile("A_phi", "A_phi", binsPhi, minPhi, maxPhi);
0087   TProfile* A_eta = new TProfile("A_eta", "A_eta", binsEta, minEta, maxEta);
0088   TProfile* A_theta
0089       = new TProfile("A_theta", "A_theta", binsTheta, minTheta, maxTheta);
0090   // Z
0091   TProfile* Z_phi = new TProfile("Z_phi", "Z_phi", binsPhi, minPhi, maxPhi);
0092   TProfile* Z_eta = new TProfile("Z_eta", "Z_eta", binsEta, minEta, maxEta);
0093   TProfile* Z_theta
0094       = new TProfile("Z_theta", "Z_theta", binsTheta, minTheta, maxTheta);
0095   // Rho
0096   TProfile* rho_phi
0097       = new TProfile("rho_phi", "rho_phi", binsPhi, minPhi, maxPhi);
0098   TProfile* rho_eta
0099       = new TProfile("rho_eta", "rho_eta", binsEta, minEta, maxEta);
0100   TProfile* rho_theta
0101       = new TProfile("rho_theta", "rho_theta", binsTheta, minTheta, maxTheta);
0102   // x0
0103   TProfile* x0_phi = new TProfile("x0_phi", "x0_phi", binsPhi, minPhi, maxPhi);
0104   TProfile* x0_eta = new TProfile("x0_eta", "x0_eta", binsEta, minEta, maxEta);
0105   TProfile* x0_theta
0106       = new TProfile("x0_theta", "x0_theta", binsTheta, minTheta, maxTheta);
0107   // l0
0108   TProfile* l0_phi = new TProfile("l0_phi", "l0_phi", binsPhi, minPhi, maxPhi);
0109   TProfile* l0_eta = new TProfile("l0_eta", "l0_eta", binsEta, minEta, maxEta);
0110   TProfile* l0_theta
0111       = new TProfile("l0_theta", "l0_theta", binsTheta, minTheta, maxTheta);
0112
0113   for (auto& mrecord : mrecords) {
0114     std::vector<Acts::MaterialStep> steps = mrecord.materialSteps();
0115     float                           theta = mrecord.theta();
0116     float                           eta   = Acts::AngleHelpers::etaFromTheta(theta);
0117     float                           phi   = VectorHelpers::phi(mrecord);
0118
0119     float thickness = 0.;
0120     float rho       = 0.;
0121     float A         = 0.;
0122     float Z         = 0.;
0123     float tInX0     = 0.;
0124     float tInL0     = 0.;
0125
0126     for (auto& step : steps) {
0127       //  std::cout << "t" << step.material().thickness() << std::endl;
0128       float t = step.material().thickness();
0129       float r = step.material().material().massDensity();
0130       thickness += step.material().thickness();
0131       rho += r * t;
0132       A += step.material().material().Ar() * r * t;
0133       Z += step.material().material().Z() * r * t;
0134       tInX0 += (t != 0. && step.material().x0() != 0.)
0135           ? t / step.material().x0()
0136           : 0.;
0137
0138       tInL0 += (t != 0. && step.material().x0() != 0.)
0139           ? t / step.material().l0()
0140           : 0.;
0141     }
0142     if (rho != 0.) {
0143       A /= rho;
0144       Z /= rho;
0145     }
0146     if (thickness != 0.) rho /= thickness;
0147
0148     t_phi->Fill(phi, thickness);
0149     t_theta->Fill(theta, thickness);
0150     t_eta->Fill(eta, thickness);
0151
0152     if (tInX0 != 0.) {
0153       x0_phi->Fill(phi, thickness / tInX0);
0154       x0_theta->Fill(theta, thickness / tInX0);
0155       x0_eta->Fill(eta, thickness / tInX0);
0156     }
0157
0158     if (tInL0 != 0.) {
0159       l0_phi->Fill(phi, thickness / tInL0);
0160       l0_theta->Fill(theta, thickness / tInL0);
0161       l0_eta->Fill(eta, thickness / tInL0);
0162     }
0163
0164     tInX0_phi->Fill(phi, tInX0);
0165     tInX0_theta->Fill(theta, tInX0);
0166     tInX0_eta->Fill(eta, tInX0);
0167
0168     tInL0_phi->Fill(phi, tInL0);
0169     tInL0_theta->Fill(theta, tInL0);
0170     tInL0_eta->Fill(eta, tInL0);
0171
0172     A_phi->Fill(phi, A);
0173     A_theta->Fill(theta, A);
0174     A_eta->Fill(eta, A);
0175
0176     Z_phi->Fill(phi, Z);
0177     Z_theta->Fill(theta, Z);
0178     Z_eta->Fill(eta, Z);
0179
0180     rho_phi->Fill(phi, rho);
0181     rho_theta->Fill(theta, rho);
0182     rho_eta->Fill(eta, rho);
0183   }
0184   // thickness
0185   t_phi->Write();
0186   t_theta->Write();
0187   t_eta->Write();
0188   delete t_theta;
0189   delete t_eta;
0190   delete t_phi;
0191   // thickness in X0
0192   tInX0_phi->Write();
0193   tInX0_theta->Write();
0194   tInX0_eta->Write();
0195   delete tInX0_phi;
0196   delete tInX0_eta;
0197   delete tInX0_theta;
0198   // thickness in L0
0199   tInL0_phi->Write();
0200   tInL0_theta->Write();
0201   tInL0_eta->Write();
0202   delete tInL0_phi;
0203   delete tInL0_eta;
0204   delete tInL0_theta;
0205   // A
0206   A_phi->Write();
0207   A_eta->Write();
0208   A_theta->Write();
0209   delete A_phi;
0210   delete A_eta;
0211   delete A_theta;
0212   // Z
0213   Z_phi->Write();
0214   Z_eta->Write();
0215   Z_theta->Write();
0216   delete Z_phi;
0217   delete Z_eta;
0218   delete Z_theta;
0219   // rho
0220   rho_phi->Write();
0221   rho_eta->Write();
0222   rho_theta->Write();
0223   delete rho_phi;
0224   delete rho_eta;
0225   delete rho_theta;
0226   // x0
0227   x0_phi->Write();
0228   x0_eta->Write();
0229   x0_theta->Write();
0230   delete x0_phi;
0231   delete x0_eta;
0232   delete x0_theta;
0233   // l0
0234   l0_phi->Write();
0235   l0_eta->Write();
0236   l0_theta->Write();
0237   delete l0_phi;
0238   delete l0_eta;
0239   delete l0_theta;
0240
0241   outputFile.Close();
0242 }