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