Examples/Scripts/printBField.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 script prints the histogram of a magnetic field map.
0019 // To be used with the Output of the RootInterpolatedBFieldWriter.
0020
0021 /// to print out the FCC field map please use
0022 /// @code
0023 /// printBField("FCChhBField.root","bField","printBField_FCC.root",-20.,20.,-30.,30.,400.)
0024 /// @endcode
0025 /// ro print out the ATLAS BField map please use
0026 /// @code
0027 /// printBField("ATLASBField.root","bField","printBField_ATLAS.root",-10.,10.,-15.,15.,200.)
0028 /// @endcode
0029
0030 /// @param inFile The root input file containing the magnetic field values and
0031 /// positions either in cylinder (Branch names: 'r','z','Br','Bz') or cartesian
0032 /// coordinates (Branch names: 'x','y','z','Bx','By','Bz')
0033 /// @param the name of the tree containing the branches
0034 /// @param rMin The minimum value of the position in either r (for cylinder
0035 /// coordinates) or x/y (for cartesian coordinates) to be printed in [m]
0036 /// @param rMin The minimum value of the position in either r (for cylinder
0037 /// coordinates) or x/y (for cartesian coordinates) to be printed in [m]
0038 /// @param rMin The maximum value of the position in either r (for cylinder
0039 /// coordinates) or x/y (for cartesian coordinates) to be printed in [m]
0040 /// @param rMin The minimum value of the position in z in [m]
0041 /// @param rMin The maximum value of the position in z in [m]
0042 /// @param nBins Number of bins which should be used for the histogram (on all
0043 /// axes)
0044 /// @note This script just writes out the values which are read in from the
0045 /// given input file. It does no interpolation in between the values. This means,
0046 /// in case the binning is chosen too high, empty bins will appear.
0047 void
0048 printBField(std::string inFile,
0049             std::string treeName,
0050             std::string outFile,
0051             float       rmin,
0052             float       rmax,
0053             float       zmin,
0054             float       zmax,
0055             int         nBins)
0056 {
0057   const Int_t NRGBs        = 5;
0058   const Int_t NCont        = 255;
0059   Double_t    stops[NRGBs] = {0.00, 0.34, 0.61, 0.84, 1.00};
0060   Double_t    red[NRGBs]   = {0.00, 0.00, 0.87, 1.00, 0.51};
0061   Double_t    green[NRGBs] = {0.00, 0.81, 1.00, 0.20, 0.00};
0062   Double_t    blue[NRGBs]  = {0.51, 1.00, 0.12, 0.00, 0.00};
0063   TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
0064   gStyle->SetNumberContours(NCont);
0065   gStyle->SetOptStat(0);
0066
0067   std::cout << "Opening file: " << inFile << std::endl;
0068   TFile inputFile(inFile.c_str());
0069   std::cout << "Reading tree: " << treeName << std::endl;
0070   TTree* tree = (TTree*)inputFile.Get(treeName.c_str());
0071
0072   TTreeReader reader(treeName.c_str(), &inputFile);
0073
0074   double x = 0., y = 0., z = 0., r = 0.;
0075   double Bx = 0., By = 0., Bz = 0., Br = 0.;
0076
0077   // find out if file is given in cylinder coordinates or cartesian coordinates
0078   bool cylinderCoordinates = false;
0079   if (tree->FindBranch("r")) {
0080     cylinderCoordinates = true;
0081     tree->SetBranchAddress("r", &r);
0082     tree->SetBranchAddress("Br", &Br);
0083   } else {
0084     tree->SetBranchAddress("x", &x);
0085     tree->SetBranchAddress("y", &y);
0086     tree->SetBranchAddress("Bx", &Bx);
0087     tree->SetBranchAddress("By", &By);
0088   }
0089   // should be given for sure
0090   tree->SetBranchAddress("z", &z);
0091   tree->SetBranchAddress("Bz", &Bz);
0092
0093   Int_t entries = tree->GetEntries();
0094   std::cout << "Creating new output file: " << outFile
0095             << " and writing out histograms. " << std::endl;
0096   TFile outputFile(outFile.c_str(), "recreate");
0097
0098   TProfile2D* bField_rz = new TProfile2D(
0099       "BField_rz", "Magnetic Field", nBins, zmin, zmax, nBins * 0.5, 0., rmax);
0100   bField_rz->GetXaxis()->SetTitle("z [m]");
0101   bField_rz->GetYaxis()->SetTitle("r [m]");
0102   TProfile2D* bField_xy = new TProfile2D(
0103       "BField_xy", "Magnetic Field", nBins, rmin, rmax, nBins, rmin, rmax);
0104   bField_xy->GetXaxis()->SetTitle("x [m]");
0105   bField_xy->GetYaxis()->SetTitle("y [m]");
0106   TProfile2D* bField_yz = new TProfile2D(
0107       "BField_yz", "Magnetic Field", nBins, zmin, zmax, nBins, rmin, rmax);
0108   bField_yz->GetXaxis()->SetTitle("z [m]");
0109   bField_yz->GetYaxis()->SetTitle("y [m]");
0110   TProfile2D* bField_xz = new TProfile2D(
0111       "BField_xz", "Magnetic Field", nBins, zmin, zmax, nBins, rmin, rmax);
0112   bField_xz->GetXaxis()->SetTitle("z [m]");
0113   bField_xz->GetYaxis()->SetTitle("x [m]");
0114
0115   for (int i = 0; i < entries; i++) {
0116     tree->GetEvent(i);
0117     if (cylinderCoordinates) {
0118       float bFieldValue = std::hypot(Br, Bz);
0119       bField_rz->Fill(z / 1000., r / 1000., bFieldValue);
0120     } else {
0121       float bFieldValue = std::hypot(Bx, By, Bz);
0122
0123       bField_xy->Fill(x / 1000., y / 1000., bFieldValue);
0124       bField_yz->Fill(z / 1000., y / 1000., bFieldValue);
0125       bField_xz->Fill(z / 1000., x / 1000., bFieldValue);
0126     }
0127   }
0128   inputFile.Close();
0129
0130   if (!cylinderCoordinates) {
0131     bField_xy->Write();
0132     bField_yz->Write();
0133     bField_xz->Write();
0134   } else
0135     bField_rz->Write();
0136
0137   delete bField_rz;
0138   delete bField_xy;
0139   delete bField_yz;
0140   delete bField_xz;
0141
0142   outputFile.Close();
0143 }