9 #include <TTreeReader.h>
10 #include <TTreeReaderValue.h>
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);
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());
72 TTreeReader reader(treeName.c_str(), &inputFile);
74 double x = 0.,
y = 0.,
z = 0.,
r = 0.;
75 double Bx = 0., By = 0.,
Bz = 0., Br = 0.;
78 bool cylinderCoordinates =
false;
79 if (tree->FindBranch(
"r")) {
80 cylinderCoordinates =
true;
81 tree->SetBranchAddress(
"r", &
r);
82 tree->SetBranchAddress(
"Br", &Br);
84 tree->SetBranchAddress(
"x", &x);
85 tree->SetBranchAddress(
"y", &
y);
86 tree->SetBranchAddress(
"Bx", &Bx);
87 tree->SetBranchAddress(
"By", &By);
90 tree->SetBranchAddress(
"z", &
z);
91 tree->SetBranchAddress(
"Bz", &
Bz);
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");
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]");
115 for (
int i = 0; i < entries; i++) {
117 if (cylinderCoordinates) {
118 float bFieldValue = sqrt(Br * Br +
Bz *
Bz);
119 bField_rz->Fill(
z / 1000.,
r / 1000., bFieldValue);
121 float bFieldValue = sqrt(Bx * Bx + By * By +
Bz *
Bz);
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);
130 if (!cylinderCoordinates) {