Examples/Scripts/printHits.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 <TTreeReader.h>
0010 #include <TTreeReaderValue.h>
0011 #include "TFile.h"
0012 #include "TH1F.h"
0013 #include "TH2F.h"
0014 #include "TProfile.h"
0015 #include "TROOT.h"
0016 #include "TTree.h"
0017
0018 // This root script creates different histograms displaying the sensitive
0019 // material, the boundaries and the material of the detector in different views.
0020 // The in file needs to be in the format created form the ExtrapolationTest.
0021 // To plot two or three histograms of this kind  in one canvas the root script
0022 // "compareHistogram.C" can be used.
0023
0024 void
0025 printHits(std::string inFile,
0026           std::string treeName,
0027           std::string outFile,
0028           float       rmin,
0029           float       rmax,
0030           float       zmin,
0031           float       zmax,
0032           int         nBins)
0033 {
0034   std::cout << "Opening file: " << inFile << std::endl;
0035   TFile inputFile(inFile.c_str());
0036   std::cout << "Reading tree: " << treeName << std::endl;
0037   TTree* tree = (TTree*)inputFile.Get(treeName.c_str());
0038
0039   TTreeReader reader(treeName.c_str(), &inputFile);
0040
0041   int   nHits;
0042   float eta;
0043   float theta;
0044   float zDir;
0045
0046   std::vector<float>* x      = new std::vector<float>;
0047   std::vector<float>* y      = new std::vector<float>;
0048   std::vector<float>* z      = new std::vector<float>;
0049   std::vector<float>* sens   = new std::vector<float>;
0050   std::vector<float>* mat    = new std::vector<float>;
0051   std::vector<float>* bounds = new std::vector<float>;
0052
0053   tree->SetBranchAddress("step_x", &x);
0054   tree->SetBranchAddress("step_y", &y);
0055   tree->SetBranchAddress("step_z", &z);
0056   tree->SetBranchAddress("sensitive", &sens);
0057
0058   if (tree->FindBranch("boundary")) {
0059     std::cout << "No BoundarySteps are given." << std::endl;
0060     tree->SetBranchAddress("boundary", &bounds);
0061   }
0062   if (tree->FindBranch("material")) {
0063     std::cout << "No MaterialSteps are given." << std::endl;
0064     tree->SetBranchAddress("material", &mat);
0065   }
0066   Int_t entries = tree->GetEntries();
0067   std::cout << "Creating new output file: " << outFile
0068             << " and writing out histograms. " << std::endl;
0069   TFile outputFile(outFile.c_str(), "recreate");
0070
0071   // full
0072   TH2F* Full_xy = new TH2F(
0073       "Full_xy", "Full material", nBins, rmin, rmax, nBins, rmin, rmax);
0074   TH2F* Full_zr = new TH2F(
0075       "Full_zr", "Full material", nBins, zmin, zmax, nBins, 0., rmax);
0076
0077   // sensitive
0078   TH2F* Sens_xy = new TH2F(
0079       "Sens_xy", "Sensitive material", nBins, rmin, rmax, nBins, rmin, rmax);
0080   Sens_xy->GetXaxis()->SetTitle("x [mm]");
0081   Sens_xy->GetYaxis()->SetTitle("y [mm]");
0082   TH2F* Sens_zx = new TH2F(
0083       "Sens_zx", "Sensitive material", nBins, zmin, zmax, nBins, rmin, rmax);
0084   Sens_zx->GetXaxis()->SetTitle("z [mm]");
0085   Sens_zx->GetYaxis()->SetTitle("x [mm]");
0086   TH2F* Sens_zy = new TH2F(
0087       "Sens_zy", "Sensitive material", nBins, zmin, zmax, nBins, rmin, rmax);
0088   Sens_zy->GetXaxis()->SetTitle("z [mm]");
0089   Sens_zy->GetYaxis()->SetTitle("y [mm]");
0090   TH2F* Sens_zr = new TH2F(
0091       "Sens_zr", "Sensitive material", nBins, zmin, zmax, nBins, 0., rmax);
0092   Sens_zr->GetXaxis()->SetTitle("z [mm]");
0093   Sens_zr->GetYaxis()->SetTitle("r [mm]");
0094
0095   // boundaries
0096   TH2F* Bounds_xy = new TH2F(
0097       "Bounds_xy", "Boundaries", nBins, rmin, rmax, nBins, rmin, rmax);
0098   Bounds_xy->GetXaxis()->SetTitle("x [mm]");
0099   Bounds_xy->GetYaxis()->SetTitle("y [mm]");
0100   TH2F* Bounds_zx = new TH2F(
0101       "Bounds_zx", "Boundaries", nBins, zmin, zmax, nBins, rmin, rmax);
0102   Bounds_zx->GetXaxis()->SetTitle("z [mm]");
0103   Bounds_zx->GetYaxis()->SetTitle("x [mm]");
0104   TH2F* Bounds_zy = new TH2F(
0105       "Bounds_zy", "Boundaries", nBins, zmin, zmax, nBins, rmin, rmax);
0106   Bounds_zy->GetXaxis()->SetTitle("z [mm]");
0107   Bounds_zy->GetYaxis()->SetTitle("y [mm]");
0108   TH2F* Bounds_zr
0109       = new TH2F("Bounds_zr", "Boundaries", nBins, zmin, zmax, nBins, 0., rmax);
0110   Bounds_zr->GetXaxis()->SetTitle("z [mm]");
0111   Bounds_zr->GetYaxis()->SetTitle("r [mm]");
0112
0113   // material
0114   TH2F* Mat_xy
0115       = new TH2F("Mat_xy", "Material", nBins, rmin, rmax, nBins, rmin, rmax);
0116   Mat_xy->GetXaxis()->SetTitle("x [mm]");
0117   Mat_xy->GetYaxis()->SetTitle("y [mm]");
0118   TH2F* Mat_zx
0119       = new TH2F("Mat_zx", "Material", nBins, zmin, zmax, nBins, rmin, rmax);
0120   Mat_zx->GetXaxis()->SetTitle("z [mm]");
0121   Mat_zx->GetYaxis()->SetTitle("x [mm]");
0122   TH2F* Mat_zy
0123       = new TH2F("Mat_zy", "Material", nBins, zmin, zmax, nBins, rmin, rmax);
0124   Mat_zy->GetXaxis()->SetTitle("z [mm]");
0125   Mat_zy->GetYaxis()->SetTitle("y [mm]");
0126   TH2F* Mat_zr
0127       = new TH2F("Mat_zr", "Material", nBins, zmin, zmax, nBins, 0., rmax);
0128   Mat_zr->GetXaxis()->SetTitle("z [mm]");
0129   Mat_zr->GetYaxis()->SetTitle("r [mm]");
0130
0131   for (int i = 0; i < entries; i++) {
0132     tree->GetEvent(i);
0133
0134     for (int j = 0; j < x->size(); j++) {
0135       // sensitive
0136       if (z->at(j) >= zmin && z->at(j) <= zmax) {
0137         Sens_xy->Fill(x->at(j), y->at(j), sens->at(j));
0138         Full_xy->Fill(x->at(j), y->at(j));
0139       }
0140       Sens_zx->Fill(z->at(j), x->at(j), sens->at(j));
0141       Sens_zy->Fill(z->at(j), y->at(j), sens->at(j));
0142       Sens_zr->Fill(z->at(j), std::hypot(x->at(j), y->at(j)));
0143       Full_zr->Fill(z->at(j), std::hypot(x->at(j), y->at(j)));
0144
0145       // boundaries
0146       if (tree->FindBranch("boundary")) {
0147         if (z->at(j) >= zmin && z->at(j) <= zmax)
0148           Bounds_xy->Fill(x->at(j), y->at(j), bounds->at(j));
0149         Bounds_zx->Fill(z->at(j), x->at(j), bounds->at(j));
0150         Bounds_zy->Fill(z->at(j), y->at(j), bounds->at(j));
0151         Bounds_zr->Fill(z->at(j),
0152                         std::hypot(x->at(j), y->at(j)),
0153                         bounds->at(j));
0154       }
0155       // material
0156       if (tree->FindBranch("material")) {
0157         if (z->at(j) >= zmin && z->at(j) <= zmax)
0158           Mat_xy->Fill(x->at(j), y->at(j), mat->at(j));
0159         Mat_zx->Fill(z->at(j), x->at(j), mat->at(j));
0160         Mat_zy->Fill(z->at(j), y->at(j), mat->at(j));
0161         Mat_zr->Fill(z->at(j),
0162                      std::hypot(x->at(j), y->at(j)),
0163                      mat->at(j));
0164       }
0165     }
0166   }
0167   inputFile.Close();
0168
0169   // full
0170   Full_xy->Write();
0171   delete Full_xy;
0172   Full_zr->Write();
0173   delete Full_zr;
0174
0175   // sensitive
0176   Sens_xy->Write();
0177   delete Sens_xy;
0178   Sens_zx->Write();
0179   delete Sens_zx;
0180   Sens_zy->Write();
0181   delete Sens_zy;
0182   Sens_zr->Write();
0183   delete Sens_zr;
0184
0185   // boundaries
0186   Bounds_xy->Write();
0187   delete Bounds_xy;
0188   Bounds_zx->Write();
0189   delete Bounds_zx;
0190   Bounds_zy->Write();
0191   delete Bounds_zy;
0192   Bounds_zr->Write();
0193   delete Bounds_zr;
0194
0195   // material
0196   Mat_xy->Write();
0197   delete Mat_xy;
0198   Mat_zx->Write();
0199   delete Mat_zx;
0200   Mat_zy->Write();
0201   delete Mat_zy;
0202   Mat_zr->Write();
0203   delete Mat_zr;
0204
0205   delete x;
0206   delete y;
0207   delete z;
0208   delete sens;
0209   delete bounds;
0210   delete mat;
0211
0212   outputFile.Close();
0213 }