EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
printBField.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file printBField.C
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2017 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include <TTreeReader.h>
10 #include <TTreeReaderValue.h>
11 #include "TFile.h"
12 #include "TH1F.h"
13 #include "TH2F.h"
14 #include "TProfile.h"
15 #include "TROOT.h"
16 #include "TTree.h"
17 
18 // This script prints the histogram of a magnetic field map.
19 // To be used with the Output of the RootInterpolatedBFieldWriter.
20 
29 
47 void
48 printBField(std::string inFile,
49  std::string treeName,
50  std::string outFile,
51  float rmin,
52  float rmax,
53  float zmin,
54  float zmax,
55  int nBins)
56 {
57  const Int_t NRGBs = 5;
58  const Int_t NCont = 255;
59  Double_t stops[NRGBs] = {0.00, 0.34, 0.61, 0.84, 1.00};
60  Double_t red[NRGBs] = {0.00, 0.00, 0.87, 1.00, 0.51};
61  Double_t green[NRGBs] = {0.00, 0.81, 1.00, 0.20, 0.00};
62  Double_t blue[NRGBs] = {0.51, 1.00, 0.12, 0.00, 0.00};
63  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
64  gStyle->SetNumberContours(NCont);
65  gStyle->SetOptStat(0);
66 
67  std::cout << "Opening file: " << inFile << std::endl;
68  TFile inputFile(inFile.c_str());
69  std::cout << "Reading tree: " << treeName << std::endl;
70  TTree* tree = (TTree*)inputFile.Get(treeName.c_str());
71 
72  TTreeReader reader(treeName.c_str(), &inputFile);
73 
74  double x = 0., y = 0., z = 0., r = 0.;
75  double Bx = 0., By = 0., Bz = 0., Br = 0.;
76 
77  // find out if file is given in cylinder coordinates or cartesian corrdinates
78  bool cylinderCoordinates = false;
79  if (tree->FindBranch("r")) {
80  cylinderCoordinates = true;
81  tree->SetBranchAddress("r", &r);
82  tree->SetBranchAddress("Br", &Br);
83  } else {
84  tree->SetBranchAddress("x", &x);
85  tree->SetBranchAddress("y", &y);
86  tree->SetBranchAddress("Bx", &Bx);
87  tree->SetBranchAddress("By", &By);
88  }
89  // should be given for sure
90  tree->SetBranchAddress("z", &z);
91  tree->SetBranchAddress("Bz", &Bz);
92 
93  Int_t entries = tree->GetEntries();
94  std::cout << "Creating new output file: " << outFile
95  << " and writing out histograms. " << std::endl;
96  TFile outputFile(outFile.c_str(), "recreate");
97 
98  TProfile2D* bField_rz = new TProfile2D(
99  "BField_rz", "Magnetic Field", nBins, zmin, zmax, nBins * 0.5, 0., rmax);
100  bField_rz->GetXaxis()->SetTitle("z [m]");
101  bField_rz->GetYaxis()->SetTitle("r [m]");
102  TProfile2D* bField_xy = new TProfile2D(
103  "BField_xy", "Magnetic Field", nBins, rmin, rmax, nBins, rmin, rmax);
104  bField_xy->GetXaxis()->SetTitle("x [m]");
105  bField_xy->GetYaxis()->SetTitle("y [m]");
106  TProfile2D* bField_yz = new TProfile2D(
107  "BField_yz", "Magnetic Field", nBins, zmin, zmax, nBins, rmin, rmax);
108  bField_yz->GetXaxis()->SetTitle("z [m]");
109  bField_yz->GetYaxis()->SetTitle("y [m]");
110  TProfile2D* bField_xz = new TProfile2D(
111  "BField_xz", "Magnetic Field", nBins, zmin, zmax, nBins, rmin, rmax);
112  bField_xz->GetXaxis()->SetTitle("z [m]");
113  bField_xz->GetYaxis()->SetTitle("x [m]");
114 
115  for (int i = 0; i < entries; i++) {
116  tree->GetEvent(i);
117  if (cylinderCoordinates) {
118  float bFieldValue = sqrt(Br * Br + Bz * Bz);
119  bField_rz->Fill(z / 1000., r / 1000., bFieldValue);
120  } else {
121  float bFieldValue = sqrt(Bx * Bx + By * By + Bz * Bz);
122 
123  bField_xy->Fill(x / 1000., y / 1000., bFieldValue);
124  bField_yz->Fill(z / 1000., y / 1000., bFieldValue);
125  bField_xz->Fill(z / 1000., x / 1000., bFieldValue);
126  }
127  }
128  inputFile.Close();
129 
130  if (!cylinderCoordinates) {
131  bField_xy->Write();
132  bField_yz->Write();
133  bField_xz->Write();
134  } else
135  bField_rz->Write();
136 
137  delete bField_rz;
138  delete bField_xy;
139  delete bField_yz;
140  delete bField_xz;
141 
142  outputFile.Close();
143 }