ATLAS Offline Software
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 "BFieldSolenoid.h"
9 #include <string>
10 #include <cmath>
11 #include <algorithm>
12 #include "TTree.h"
13 using namespace std;
14 
15 //
16 // read an ASCII field map from istream
17 // return 0 if successful
18 //
19 int
20 BFieldSolenoid::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 //
103 void
104 BFieldSolenoid::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 //
187 int
188 BFieldSolenoid::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 //
286 void
287 BFieldSolenoid::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 //
304 void
305 BFieldSolenoid::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 
beamspotman.r
def r
Definition: beamspotman.py:676
BFieldSolenoid.h
plotBeamSpotCompare.x1
x1
Definition: plotBeamSpotCompare.py:216
BFieldSolenoid::getB
void getB(const double *xyz, double *B, double *deriv=nullptr) const
Definition: BFieldSolenoid.cxx:287
python.CaloRecoConfig.f
f
Definition: CaloRecoConfig.py:127
fitman.ax
ax
Definition: fitman.py:522
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:64
BFieldMesh::phimin
double phimin() const
minimun in phi
BFieldSolenoid::moveMap
void moveMap(double dx, double dy, double dz, double ax, double ay)
Definition: BFieldSolenoid.cxx:305
PixelAthClusterMonAlgCfg.zmin
zmin
Definition: PixelAthClusterMonAlgCfg.py:176
index
Definition: index.py:1
InDetAccessor::phi0
@ phi0
Definition: InDetAccessor.h:33
plotBeamSpotCompare.x2
x2
Definition: plotBeamSpotCompare.py:218
tree
TChain * tree
Definition: tile_monitor.h:30
M_PI
#define M_PI
Definition: ActiveFraction.h:11
xyz
#define xyz
ReadOfcFromCool.field
field
Definition: ReadOfcFromCool.py:48
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
BFieldSolenoid::readMap
int readMap(std::istream &input)
Definition: BFieldSolenoid.cxx:20
x
#define x
LArCalib_HVScale2NtupleConfig.rootfile
string rootfile
Definition: LArCalib_HVScale2NtupleConfig.py:74
BFieldSolenoid::writeMap
void writeMap(TFile *rootfile, bool tilted=false)
Definition: BFieldSolenoid.cxx:104
makeTRTBarrelCans.y1
tuple y1
Definition: makeTRTBarrelCans.py:15
Br
Definition: VP1BPhysConvertor.h:14
BFieldMesh::phimax
double phimax() const
maximum in phi
python.SystemOfUnits.meter
int meter
Definition: SystemOfUnits.py:61
lumiFormat.i
int i
Definition: lumiFormat.py:92
z
#define z
BFieldMesh< double >
makeTRTBarrelCans.y2
tuple y2
Definition: makeTRTBarrelCans.py:18
BFieldMesh::zmax
double zmax() const
maximum in z
PlotPulseshapeFromCool.input
input
Definition: PlotPulseshapeFromCool.py:106
TRT_PAI_physicsConstants::r0
const double r0
electron radius{cm}
Definition: TRT_PAI_physicsConstants.h:20
BFieldMesh::mesh
double mesh(size_t axis, size_t index) const
coordinate along axis (0 (z), 1 (r), 2 (phi)) of the cell at index (0 to nmesh-1)
PixelAthClusterMonAlgCfg.zmax
zmax
Definition: PixelAthClusterMonAlgCfg.py:176
BFieldMesh::nmesh
unsigned nmesh(size_t axis) const
number of cells along each axis = 0 (z), 1 (r), 2 (phi)
BFieldMesh::field
const BFieldVector< T > & field(size_t index) const
field vector at cell corner at index
TRT::Track::z0
@ z0
Definition: InnerDetector/InDetCalibEvent/TRT_CalibData/TRT_CalibData/TrackInfo.h:63
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
dqt_zlumi_alleff_HIST.B
B
Definition: dqt_zlumi_alleff_HIST.py:110
get_generator_info.version
version
Definition: get_generator_info.py:33
makeTRTBarrelCans.dy
tuple dy
Definition: makeTRTBarrelCans.py:21
unused
void unused(Args &&...)
Definition: VP1ExpertSettings.cxx:39
BFieldMesh::zmin
double zmin() const
minimum in z
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
python.SystemOfUnits.gauss
int gauss
Definition: SystemOfUnits.py:230
BFieldVector< double >
makeTRTBarrelCans.dx
tuple dx
Definition: makeTRTBarrelCans.py:20
BFieldMesh::rmax
double rmax() const
maximum in r
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
fitman.k
k
Definition: fitman.py:528
BFieldMesh::rmin
double rmin() const
minimun in r
fitman.ay
ay
Definition: fitman.py:525
Generated on Fri Jul 5 2024 21:11:10 for ATLAS Offline Software by   doxygen 1.8.18