ATLAS Offline Software
Loading...
Searching...
No Matches
BFieldSolenoid.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
3*/
4
5//
6// BFieldSolenoid.cxx
7//
8#include "MagFieldUtils/BFieldSolenoid.h"
9#include <string>
10#include <cmath>
11#include <algorithm>
12#include "TTree.h"
13using namespace std;
14
15//
16// read an ASCII field map from istream
17// return 0 if successful
18//
19int
20BFieldSolenoid::readMap( istream& input )
21{
22 const double meter(1000.0); // meter in mm
23 const double gauss(1.0e-7); // gauss in kT
24 const string myname("BFieldSolenoid::readMap()");
25 if ( m_orig == m_tilt ) delete m_orig;
26 else { delete m_orig; delete m_tilt; }
27 m_orig = m_tilt = new BFieldMesh<double>;
28 // first line contains version
29 int version;
30 input >> version;
31 if ( version != 4 ) {
32 cerr << myname << ": version number is " << version << " instead of 4" << endl;
33 return 1;
34 }
35 // second line contains number of bins
36 int nz;
37
38 int nr;
39
40 int nphi;
41 double unused;
42 input >> nz >> nr >> nphi >> unused >> unused >> unused;
43 m_orig->reserve( nz, nr, nphi+1 ); // extra phi at 2pi
44 // read mesh definitions
45 double x;
46 for ( int i = 0; i < nr; i++ ) { input >> x; m_orig->appendMesh(1,x*meter); }
47 for ( int i = 0; i < nz; i++ ) { input >> x; m_orig->appendMesh(0,x*meter); }
48 for ( int i = 0; i < nphi; i++ ) { input >> x; m_orig->appendMesh(2,x); }
49 m_orig->appendMesh(2,2*M_PI); // extra phi at 2pi
50 double zmin = m_orig->mesh(0,0);
51 double zmax = m_orig->mesh(0,nz-1);
52 double rmin = m_orig->mesh(1,0);
53 double rmax = m_orig->mesh(1,nr-1);
54 double phimin = m_orig->mesh(2,0);
55 double phimax = m_orig->mesh(2,nphi);
56 m_orig->setRange( zmin, zmax, rmin, rmax, phimin, phimax );
57 // skip one record
58 input >> unused >> unused >> unused >> unused >> unused;
59 // read field values
60 // the index order is opposite to the toroid - have to reorder
61 // also convert from gauss to Tesla here.
62 vector<double> Bz;
63
64 vector<double> Br;
65
66 vector<double> Bphi;
67 double b;
68 for ( int k = 0; k < nphi; k++ ) {
69 for ( int j = 0; j < nr; j++ ) {
70 for ( int i = 0; i < nz; i++ ) { input >> b; Bz.push_back(b*gauss); }
71 }
72 for ( int j = 0; j < nr; j++ ) {
73 for ( int i = 0; i < nz; i++ ) { input >> b; Br.push_back(b*gauss); }
74 }
75 for ( int j = 0; j < nr; j++ ) {
76 for ( int i = 0; i < nz; i++ ) { input >> b; Bphi.push_back(b*gauss); }
77 }
78 }
79 for ( int i = 0; i < nz; i++ ) {
80 for ( int j = 0; j < nr; j++ ) {
81 for ( int k = 0; k < nphi; k++ ) {
82 int index = i + nz*(j + nr*k);
83 BFieldVector<double> field( Bz[index], Br[index], Bphi[index] );
84 m_orig->appendField( field );
85 }
86 // close phi at 2pi
87 int index = i + nz*j;
88 BFieldVector<double> field( Bz[index], Br[index], Bphi[index] );
89 m_orig->appendField( field );
90 }
91 }
92 // build the LUT
93 m_orig->buildLUT();
94
95 return 0;
96}
97
98//
99// wrire the map to a ROOT file
100// if tilted = true, write the moved-and-tilted map.
101// ohterwise, write the original map.
102//
103void
104BFieldSolenoid::writeMap( TFile* rootfile, bool tilted )
105{
106 BFieldMesh<double> *map = tilted ? m_tilt : m_orig;
107 if ( map == nullptr ) return; // no map to write
108 if ( rootfile == nullptr ) return; // no file
109 if ( !rootfile->cd() ) return; // could not make it current directory
110 // define the tree
111 TTree* tree = new TTree( "BFieldSolenoid", "BFieldSolenoid version 4" );
112 double zmin = map->zmin();
113 double zmax = map->zmax();
114 double rmin = map->rmin();
115 double rmax = map->rmax();
116 double phimin = map->phimin();
117 double phimax = map->phimax();
118 int nmeshz = map->nmesh(0);
119 int nmeshr = map->nmesh(1);
120 int nmeshphi = map->nmesh(2);
121 double *meshz;
122
123 double *meshr;
124
125 double *meshphi;
126 int nfield = nmeshz*nmeshr*nmeshphi;
127 double *fieldz;
128
129 double *fieldr;
130
131 double *fieldphi;
132 meshz = new double[nmeshz];
133 meshr = new double[nmeshr];
134 meshphi = new double[nmeshphi];
135 fieldz = new double[nfield];
136 fieldr = new double[nfield];
137 fieldphi = new double[nfield];
138 // define the tree branches
139 tree->Branch( "zmin", &zmin, "zmin/D" );
140 tree->Branch( "zmax", &zmax, "zmax/D" );
141 tree->Branch( "rmin", &rmin, "rmin/D" );
142 tree->Branch( "rmax", &rmax, "rmax/D" );
143 tree->Branch( "phimin", &phimin, "phimin/D" );
144 tree->Branch( "phimax", &phimax, "phimax/D" );
145 tree->Branch( "nmeshz", &nmeshz, "nmeshz/I" );
146 tree->Branch( "meshz", meshz, "meshz[nmeshz]/D" );
147 tree->Branch( "nmeshr", &nmeshr, "nmeshr/I" );
148 tree->Branch( "meshr", meshr, "meshr[nmeshr]/D" );
149 tree->Branch( "nmeshphi", &nmeshphi, "nmeshphi/I" );
150 tree->Branch( "meshphi", meshphi, "meshphi[nmeshphi]/D" );
151 tree->Branch( "nfield", &nfield, "nfield/I" );
152 tree->Branch( "fieldz", fieldz, "fieldz[nfield]/D" );
153 tree->Branch( "fieldr", fieldr, "fieldr[nfield]/D" );
154 tree->Branch( "fieldphi", fieldphi, "fieldphi[nfield]/D" );
155 // fill the mesh and field arrays
156 for ( int j = 0; j < nmeshz; j++ ) {
157 meshz[j] = map->mesh(0,j);
158 }
159 for ( int j = 0; j < nmeshr; j++ ) {
160 meshr[j] = map->mesh(1,j);
161 }
162 for ( int j = 0; j < nmeshphi; j++ ) {
163 meshphi[j] = map->mesh(2,j);
164 }
165 for ( int j = 0; j < nfield; j++ ) {
166 const BFieldVector<double> f = map->field(j);
167 fieldz[j] = f.z();
168 fieldr[j] = f.r();
169 fieldphi[j] = f.phi();
170 }
171 // write
172 tree->Fill();
173 rootfile->Write();
174 // clean up
175 delete[] meshz;
176 delete[] meshr;
177 delete[] meshphi;
178 delete[] fieldz;
179 delete[] fieldr;
180 delete[] fieldphi;
181}
182
183//
184// read the map from a ROOT file.
185// returns 0 if successful.
186//
187int
188BFieldSolenoid::readMap( TFile* rootfile )
189{
190 if ( rootfile == nullptr ) return 1; // no file
191 if ( !rootfile->cd() ) return 2; // could not make it current directory
192 if ( m_orig == m_tilt ) delete m_orig;
193 else { delete m_orig; delete m_tilt; }
194 m_orig = m_tilt = new BFieldMesh<double>;
195 // open the tree
196 TTree* tree = (TTree*)rootfile->Get("BFieldSolenoid");
197 if ( tree == nullptr ) return 3; // no tree
198 double zmin;
199
200 double zmax;
201
202 double rmin;
203
204 double rmax;
205
206 double phimin;
207
208 double phimax;
209 int nmeshz;
210
211 int nmeshr;
212
213 int nmeshphi;
214 double *meshz;
215
216 double *meshr;
217
218 double *meshphi;
219 int nfield;
220 double *fieldz;
221
222 double *fieldr;
223
224 double *fieldphi;
225 //unsigned char *fbyte;
226 // define the fixed-sized branches first
227 tree->SetBranchAddress( "zmin", &zmin );
228 tree->SetBranchAddress( "zmax", &zmax );
229 tree->SetBranchAddress( "rmin", &rmin );
230 tree->SetBranchAddress( "rmax", &rmax );
231 tree->SetBranchAddress( "phimin", &phimin );
232 tree->SetBranchAddress( "phimax", &phimax );
233 tree->SetBranchAddress( "nmeshz", &nmeshz );
234 tree->SetBranchAddress( "nmeshr", &nmeshr );
235 tree->SetBranchAddress( "nmeshphi", &nmeshphi );
236 tree->SetBranchAddress( "nfield", &nfield );
237 // prepare arrays - need to know the maximum sizes
238 tree->GetEntry(0);
239 meshz = new double[nmeshz];
240 meshr = new double[nmeshr];
241 meshphi = new double[nmeshphi];
242 fieldz = new double[nfield];
243 fieldr = new double[nfield];
244 fieldphi = new double[nfield];
245 // define the variable length branches
246 tree->SetBranchAddress( "meshz", meshz );
247 tree->SetBranchAddress( "meshr", meshr );
248 tree->SetBranchAddress( "meshphi", meshphi );
249 tree->SetBranchAddress( "fieldz", fieldz );
250 tree->SetBranchAddress( "fieldr", fieldr );
251 tree->SetBranchAddress( "fieldphi", fieldphi );
252 // read again, and copy data
253 tree->GetEntry(0);
254 m_orig->setRange( zmin, zmax, rmin, rmax, phimin, phimax );
255 m_orig->reserve( nmeshz, nmeshr, nmeshphi );
256 for ( int j = 0; j < nmeshz; j++ ) {
257 m_orig->appendMesh( 0, meshz[j] );
258 }
259 for ( int j = 0; j < nmeshr; j++ ) {
260 m_orig->appendMesh( 1, meshr[j] );
261 }
262 for ( int j = 0; j < nmeshphi; j++ ) {
263 m_orig->appendMesh( 2, meshphi[j] );
264 }
265 for ( int j = 0; j < nfield; j++ ) {
266 BFieldVector<double> field( fieldz[j], fieldr[j], fieldphi[j] );
267 m_orig->appendField( field );
268 }
269 // clean up
270 tree->Delete();
271 delete[] meshz;
272 delete[] meshr;
273 delete[] meshphi;
274 delete[] fieldz;
275 delete[] fieldr;
276 delete[] fieldphi;
277 // build the LUTs
278 m_orig->buildLUT();
279
280 return 0;
281}
282
283//
284// Returns the magnetic field at any position.
285//
286void
287BFieldSolenoid::getB( const double *xyz, double *B, double *deriv ) const
288{
289 double z = xyz[2];
290 double r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
291 double phi = atan2( xyz[1], xyz[0] );
292 if ( m_tilt && m_tilt->inside( z, r, phi ) ) {
293 m_tilt->getB( xyz, B, deriv );
294 } else {
295 B[0] = B[1] = B[2] = 0.0;
296 if ( deriv ) for ( int i = 0; i < 9; i++ ) deriv[i] = 0.0;
297 }
298}
299
300//
301// Move and tilt the solenoid.
302// Modify the m_tilt copy and keep the m_orig copy.
303//
304void
305BFieldSolenoid::moveMap( double dx, double dy, double dz, double ax, double ay )
306{
307 if ( m_orig==nullptr ) {
308 cerr << "BFieldSolenoid::moveMap() : original map has not been read" << endl;
309 return;
310 }
311 if ( m_tilt != m_orig ) delete m_tilt;
312 //
313 // copy the mesh and make it a bit smaller
314 //
315 const double zlim = 2820.; // mm
316 const double rlim = 1075.; // mm
317 m_tilt = new BFieldMesh<double>( -zlim, zlim, 0.0, rlim, 0.0, 2*M_PI, 1.0 );
318 // z
319 m_tilt->appendMesh( 0, -zlim );
320 for ( unsigned i = 0; i < m_orig->nmesh(0); i++ ) {
321 if ( abs(m_orig->mesh(0,i)) < zlim ) m_tilt->appendMesh( 0, m_orig->mesh(0,i) );
322 }
323 m_tilt->appendMesh( 0, zlim );
324 // r
325 for ( unsigned i = 0; i < m_orig->nmesh(1); i++ ) {
326 if ( m_orig->mesh(1,i) < rlim ) m_tilt->appendMesh( 1, m_orig->mesh(1,i) );
327 }
328 m_tilt->appendMesh( 1, rlim );
329 // phi (no change)
330 for ( unsigned i = 0; i < m_orig->nmesh(2); i++ ) {
331 m_tilt->appendMesh( 2, m_orig->mesh(2,i) );
332 }
333 //
334 // loop over the new mesh, and compute the field at the
335 // corresponding location in the original map
336 //
337 double sinax( sin(ax) );
338 double cosax( cos(ax) );
339 double sinay( sin(ay) );
340 double cosay( cos(ay) );
341 for ( unsigned i = 0; i < m_tilt->nmesh(0); i++ ) {
342 double z0( m_tilt->mesh(0,i) );
343 for ( unsigned j = 0; j < m_tilt->nmesh(1); j++ ) {
344 double r0( m_tilt->mesh(1,j) );
345 for ( unsigned k = 0; k < m_tilt->nmesh(2); k++ ) {
346 double phi0( m_tilt->mesh(2,k) );
347 double x0( r0*cos(phi0) );
348 double y0( r0*sin(phi0) );
349 // shift
350 double x1( x0 - dx );
351 double y1( y0 - dy );
352 double z1( z0 - dz );
353 // rotate around x by -ax
354 double x2( x1 );
355 double y2( y1*cosax + z1*sinax );
356 double z2( z1*cosax - y1*sinax );
357 // rotate around y by -ay
358 double xyz3[3] = { x2*cosay-z2*sinay, y2, z2*cosay + x2*sinay };
359 // get (Bx,By,Bz) in the original frame
360 double B[3];
361 m_orig->getB( xyz3, B );
362 // rotate around y by +ay
363 double Bx1( B[0]*cosay + B[2]*sinay );
364 double By1( B[1] );
365 double Bz1( B[2]*cosay - B[0]*sinay );
366 // rotate around x by +ax
367 double Bx2( Bx1 );
368 double By2( By1*cosax - Bz1*sinax );
369 double Bz2( Bz1*cosax + By1*sinax );
370 // convert to cylindrical
371 double cosphi0( cos(phi0) );
372 double sinphi0( sin(phi0) );
373 double Br( Bx2*cosphi0 + By2*sinphi0 );
374 double Bphi( -Bx2*sinphi0 + By2*cosphi0 );
375 BFieldVector<double> field( Bz2, Br, Bphi );
376 m_tilt->appendField( field );
377 }
378 }
379 }
380 m_tilt->buildLUT();
381}
382
M_PI
#define M_PI
Definition ActiveFraction.h:11
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
xyz
#define xyz
x
#define x
z
#define z
unused
void unused(Args &&...)
Definition VP1ExpertSettings.cxx:39
BFieldMesh
A 3-d mesh (axes z,r,phi) representing a simple field map.
Definition BFieldMesh.h:43
BFieldSolenoid::m_orig
BFieldMesh< double > * m_orig
Definition BFieldSolenoid.h:39
BFieldSolenoid::getB
void getB(const double *xyz, double *B, double *deriv=nullptr) const
Definition BFieldSolenoid.cxx:287
BFieldSolenoid::readMap
int readMap(std::istream &input)
Definition BFieldSolenoid.cxx:20
BFieldSolenoid::m_tilt
BFieldMesh< double > * m_tilt
Definition BFieldSolenoid.h:40
BFieldSolenoid::writeMap
void writeMap(TFile *rootfile, bool tilted=false)
Definition BFieldSolenoid.cxx:104
BFieldSolenoid::moveMap
void moveMap(double dx, double dy, double dz, double ax, double ay)
Definition BFieldSolenoid.cxx:305
r
int r
Definition globals.cxx:22
rootfile
static std::vector< std::string > rootfile
Definition iLumiCalc.h:30
index
Definition index.py:1
std
STL namespace.
tree
TChain * tree
Definition tile_monitor.h:30
Generated on for ATLAS Offline Software by doxygen 1.14.0