4 #include <TDirectory.h>
8 #include <Geant4/G4SystemOfUnits.hh>
12 #include <boost/tuple/tuple_comparison.hpp>
33 cout <<
" ------------- PHField2D::PHField2D() ------------------" << endl;
36 TFile *rootinput = TFile::Open(filename.c_str());
39 cout <<
" could not open " << filename <<
" exiting now" << endl;
43 if (
Verbosity() > 0) cout <<
" Field grid file: " << filename << endl;
49 TNtuple *field_map = (TNtuple *) gDirectory->Get(
"fieldmap");
56 field_map = (TNtuple *) gDirectory->Get(
"map");
59 cout <<
"PHField2D: could not locate ntuple of name map or fieldmap, exiting now" << endl;
64 field_map->SetBranchAddress(
"z", &ROOT_Z);
65 field_map->SetBranchAddress(
"r", &ROOT_R);
66 field_map->SetBranchAddress(
"bz", &ROOT_BZ);
67 field_map->SetBranchAddress(
"br", &ROOT_BR);
71 nz = field_map->GetEntries(
"z>-1e6");
72 nr = field_map->GetEntries(
"r>-1e6");
73 static const int NENTRIES = field_map->GetEntries();
76 if (
Verbosity() > 0) cout <<
" The field grid contained " << NENTRIES <<
" entries" << endl;
79 cout <<
"\n NENTRIES should be the same as the following values:"
80 <<
"\n [ Number of values r,z: "
81 << nr <<
" " << nz <<
" ]! " << endl;
86 cout <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
87 <<
"\n The file you entered is not a \"table\" of values"
88 <<
"\n Something very likely went oh so wrong"
89 <<
"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
93 std::set<float> z_set, r_set;
102 cout <<
" --> Sorting Entries..." << endl;
104 std::map<trio, trio> sorted_map;
105 for (
int i = 0; i < field_map->GetEntries(); i++)
107 field_map->GetEntry(i);
108 trio coord_key(ROOT_Z, ROOT_R);
109 trio field_val(ROOT_BZ, ROOT_BR);
110 sorted_map[coord_key] = field_val;
112 z_set.insert(ROOT_Z *
cm);
113 r_set.insert(ROOT_R * cm);
119 map<trio, trio>::iterator
it = sorted_map.begin();
121 float last_z = it->first.get<0>();
122 for (it = sorted_map.begin(); it != sorted_map.end(); ++
it)
124 if (it->first.get<0>() != last_z)
126 last_z = it->first.get<0>();
134 cout <<
" --> Putting entries into containers... " << endl;
138 minz_ = *(z_set.begin());
139 std::set<float>::iterator ziter = z_set.end();
153 BFieldR_.resize(nz, vector<float>(nr, 0));
154 BFieldZ_.resize(nz, vector<float>(nr, 0));
157 unsigned int ir = 0, iz = 0;
158 map<trio, trio>::iterator iter = sorted_map.begin();
159 for (; iter != sorted_map.end(); ++iter)
162 float z = iter->first.get<0>() *
cm;
163 float r = iter->first.get<1>() *
cm;
182 if (iz > 0 && z <
z_map_[iz - 1])
184 cout <<
"!!!!!!!!! Your map isn't ordered.... z: " << z <<
" zprev: " <<
z_map_[iz - 1] << endl;
193 if (fabs(z) < 10 && ir < 10 &&
Verbosity() > 3)
208 if (
Verbosity() > 0) cout <<
" Mag field z boundaries (min,max): (" <<
minz_ /
cm <<
", " <<
maxz_ /
cm <<
") cm" << endl;
209 if (
Verbosity() > 0) cout <<
" Mag field r max boundary: " <<
r_map_.back() /
cm <<
" cm" << endl;
212 cout <<
" -----------------------------------------------------------" << endl;
218 cout <<
"\nPHField2D::GetFieldValue" << endl;
222 double r = sqrt(x * x + y * y);
225 if (phi < 0) phi += 2 *
M_PI;
231 double cylpoint[4] = {
z,
r,
phi, 0};
237 Bfield[0] = cos(phi) * BFieldCyl[1] - sin(phi) * BFieldCyl[2];
240 Bfield[1] = sin(phi) * BFieldCyl[1] + cos(phi) * BFieldCyl[2];
243 Bfield[2] = BFieldCyl[0];
251 cout <<
"!!!!!!!!!! Field point not in defined region (outside of z bounds)" << endl;
256 cout <<
"END PHField2D::GetFieldValue\n"
257 <<
" ---> {Bx, By, Bz} : "
258 <<
"< " << Bfield[0] <<
", " << Bfield[1] <<
", " << Bfield[2] <<
" >" << endl;
266 float z = CylPoint[0];
267 float r = CylPoint[1];
274 cout <<
"GetFieldCyl@ <z,r>: {" << z <<
"," << r <<
"}" << endl;
279 cout <<
"!!!! Point not in defined region (radius too large in specific z-plane)" << endl;
290 if (!((r >
r_map_[r_index0]) && (r <
r_map_[r_index1])))
293 vector<float>::const_iterator riter = upper_bound(
r_map_.begin(),
r_map_.end(),
r);
294 r_index0 = distance(
r_map_.begin(), riter) - 1;
295 if (r_index0 >=
r_map_.size())
298 cout <<
"!!!! Point not in defined region (radius too large in specific z-plane)" << endl;
302 r_index1 = r_index0 + 1;
303 if (r_index1 >=
r_map_.size())
306 cout <<
"!!!! Point not in defined region (radius too large in specific z-plane)" << endl;
318 if (!((z >
z_map_[z_index0]) && (z <
z_map_[z_index1])))
321 vector<float>::const_iterator ziter = upper_bound(
z_map_.begin(),
z_map_.end(),
z);
322 z_index0 = distance(
z_map_.begin(), ziter) - 1;
323 z_index1 = z_index0 + 1;
324 if (z_index1 >=
z_map_.size())
327 cout <<
"!!!! Point not in defined region (z too large in specific r-plane)" << endl;
336 double Br000 =
BFieldR_[z_index0][r_index0];
337 double Br010 =
BFieldR_[z_index0][r_index1];
338 double Br100 =
BFieldR_[z_index1][r_index0];
339 double Br110 =
BFieldR_[z_index1][r_index1];
341 double Bz000 =
BFieldZ_[z_index0][r_index0];
342 double Bz100 =
BFieldZ_[z_index1][r_index0];
343 double Bz010 =
BFieldZ_[z_index0][r_index1];
344 double Bz110 =
BFieldZ_[z_index1][r_index1];
346 double zweight = z -
z_map_[z_index0];
347 double zspacing = z_map_[z_index1] - z_map_[z_index0];
350 double rweight = r -
r_map_[r_index0];
351 double rspacing = r_map_[r_index1] - r_map_[r_index0];
356 (1 - zweight) * ((1 - rweight) * Bz000 +
358 zweight * ((1 - rweight) * Bz100 +
363 (1 - zweight) * ((1 - rweight) * Br000 +
365 zweight * ((1 - rweight) * Br100 +
373 cout <<
"End GFCyl Call: <bz,br,bphi> : {"
374 << BfieldCyl[0] /
gauss <<
"," << BfieldCyl[1] /
gauss <<
"," << BfieldCyl[2] /
gauss <<
"}"
385 << it->first.get<0>() /
cm <<
","
386 << it->first.get<1>() /
cm <<
">"