Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Field.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Field.cc
1 /*
2  Authors: Haiwang Yu
3 */
4 #include "Field.h"
5 
6 #include <phfield/PHField.h>
7 
8 #include <TVector3.h> // for TVector3
9 
10 #include <CLHEP/Units/SystemOfUnits.h>
11 
12 #include <cassert>
13 #include <limits>
14 
15 namespace genfit
16 {
18  : field_(field)
19 {
20  assert(field_);
21 }
22 
23 //
24 //bool Field::initialize(std::string inname)
25 //{
26 // field_map_r_ = new TH2D("field_map_r_", "B_{r} [kGauss]; z [cm]; r [cm]", 401, -401, 401, 151, -1, 301);
27 // field_map_r_->SetDirectory(0);
28 // field_map_z_ = new TH2D("field_map_z_", "B_{z} [kGauss]; z [cm]; r [cm]", 401, -401, 401, 151, -1, 301);
29 // field_map_z_->SetDirectory(0);
30 // TFile* fin = TFile::Open(inname.data(), "READ");
31 // if (!fin)
32 // {
33 // LogERROR("Input TFile is invalid!");
34 // return false;
35 // }
36 // TNtuple* T = (TNtuple*) fin->Get("fieldmap");
37 // if (!T)
38 // {
39 // LogERROR("Input filed map NTuple not found!");
40 // ;
41 // return false;
42 // }
43 //
44 // Float_t r;
45 // Float_t z;
46 // Float_t br;
47 // Float_t bz;
48 // T->SetBranchAddress("r", &r);
49 // T->SetBranchAddress("z", &z);
50 // T->SetBranchAddress("br", &br);
51 // T->SetBranchAddress("bz", &bz);
52 //
53 // for (long ientry = 0; ientry < T->GetEntries(); ientry++)
54 // {
55 // //if (T->GetEntry(ientry) < 0) continue;
56 // T->GetEntry(ientry);
57 // //std::cout<<r <<" ,"<<z <<" ,"<<br<<" ,"<<bz<<"\n";
58 // field_map_r_->SetBinContent(
59 // field_map_r_->GetXaxis()->FindBin(z),
60 // field_map_r_->GetYaxis()->FindBin(r),
61 // br * 10);
62 // field_map_z_->SetBinContent(
63 // field_map_z_->GetXaxis()->FindBin(z),
64 // field_map_z_->GetYaxis()->FindBin(r),
65 // bz * 10);
66 // }
67 // fin->Close();
68 //
69 // //field_map_z_->Print("all");
70 //
71 // return true;
72 //}
73 
74 //void Field::plot(std::string option){
75 //
78 // TH2D *hbr = new TH2D("hbr","|B_{r}| [kGauss]; z [cm]; r [cm]",300, -350, 350, 200, 0, 290);
79 // TH2D *hbz = new TH2D("hbz","|B_{z}| [kGauss]; z [cm]; r [cm]",300, -350, 350, 200, 0, 290);
80 // for (double r = 0; r < 300; r+=1)
81 // for (double z = 400; z > -400; z-=1) {
82 // double bx;
83 // double by;
84 // double bz;
85 // get(r, 0, z, bx, by, bz);
86 //
88 //
89 // //std::cout <<"DEBUG: "<<__LINE__<<":"<< r << " ," << z << " ," << br << " ," << bz << "\n";
90 //
91 // hbr->SetBinContent(
92 // hbr->GetXaxis()->FindBin(z),
93 // hbr->GetYaxis()->FindBin(r),
94 // bz*10);
95 // hbz->SetBinContent(
96 // hbz->GetXaxis()->FindBin(z),
97 // hbz->GetYaxis()->FindBin(r),
98 // bz*10);
99 // }
100 //
101 // TCanvas *c0 = new TCanvas("c0","c0");
102 // c0->Divide(1,2);
103 // c0->cd(1);
104 // hbr->SetStats(0);
105 // hbr->Draw("colz");
106 // c0->cd(2);
107 // hbz->SetStats(0);
108 // hbz->Draw("colz");
109 // c0->Update();
110 //
111 // c0->SaveAs("Field_plot.root");
112 // c0->SaveAs("Field_plot.pdf");
113 //
114 // delete c0;
115 // delete hbr;
116 // delete hbz;
117 //
118 // return;
119 //}
120 
121 TVector3 Field::get(const TVector3& v) const
122 {
123  double x = v.x();
124  double y = v.y();
125  double z = v.z();
126  double Bx;
127  double By;
128  double Bz;
129  get(x, y, z, Bx, By, Bz);
130  return TVector3(Bx, By, Bz);
131 }
132 
133 void Field::get(const double& x, const double& y, const double& z, double& Bx, double& By, double& Bz) const
134 {
135  assert(field_);
136 
137  const double Point[] = {x * CLHEP::cm, y * CLHEP::cm, z * CLHEP::cm, 0};
138  double Bfield[] = {std::numeric_limits<double>::signaling_NaN(),
139  std::numeric_limits<double>::signaling_NaN(),
140  std::numeric_limits<double>::signaling_NaN(),
141  std::numeric_limits<double>::signaling_NaN(),
142  std::numeric_limits<double>::signaling_NaN(),
143  std::numeric_limits<double>::signaling_NaN()};
144 
145  field_->GetFieldValue(Point, Bfield);
146 
147  Bx = Bfield[0] / CLHEP::kilogauss;
148  By = Bfield[1] / CLHEP::kilogauss;
149  Bz = Bfield[2] / CLHEP::kilogauss;
150 
151  // double r = sqrt(x * x + y * y);
152  //
153  // if (fabs(r) > 300 || fabs(z) > 400)
154  // {
155  // Bx = 0;
156  // By = 0;
157  // Bz = 0;
158  // return;
159  // }
160  //
161  // int bin_z = field_map_r_->GetXaxis()->FindBin(z);
162  // int bin_r = field_map_r_->GetYaxis()->FindBin(r);
163  //
164  // Bz = field_map_z_->GetBinContent(bin_z, bin_r);
165  // double Br = field_map_r_->GetBinContent(bin_z, bin_r);
166  //
167  // Bx = x / r * Br;
168  // By = y / r * Br;
169  //
170  // if (r == 0)
171  // {
172  // Bx = 0;
173  // By = 0;
174  // }
175 
176  //std::cout<<"DEBUG: "<<__LINE__<<": "<<z<<","<<r<<","<<bin_z<<","<<bin_r<<": "<<Br<<","<<Bz<<"\n";
177 }
178 
179 } /* End of namespace genfit */