coding conventions corrections
[u/mrichter/AliRoot.git] / STEER / AliFieldMap.cxx
CommitLineData
84737f5e 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
acd84897 16/* $Id$ */
84737f5e 17
116cbefd 18//-----------------------------------------------------------------------
5d8718b8 19// Class to read the magnetic field map and to calculate
20// the magnetic field in a given point. For the moment uses a simple
21// linear interpolation between the nodes of a rectangular grid
84737f5e 22// Author: Andreas Morsch <andreas.morsch@cern.ch>
116cbefd 23//-----------------------------------------------------------------------
84737f5e 24
116cbefd 25#include <TSystem.h>
84737f5e 26#include <TVector.h>
116cbefd 27
84737f5e 28#include "AliFieldMap.h"
84737f5e 29
30ClassImp(AliFieldMap)
31
e2afb3b6 32//_______________________________________________________________________
33AliFieldMap::AliFieldMap():
34 fXbeg(0),
35 fYbeg(0),
36 fZbeg(0),
37 fXend(0),
38 fYend(0),
39 fZend(0),
40 fXdel(0),
41 fYdel(0),
42 fZdel(0),
43 fXdeli(0),
44 fYdeli(0),
45 fZdeli(0),
46 fXn(0),
47 fYn(0),
48 fZn(0),
49 fWriteEnable(0),
50 fB(0)
84737f5e 51{
52 //
53 // Standard constructor
54 //
7b6cddfa 55 SetWriteEnable();
84737f5e 56}
57
e2afb3b6 58//_______________________________________________________________________
59AliFieldMap::AliFieldMap(const char *name, const char *title):
60 TNamed(name,title),
61 fXbeg(0),
62 fYbeg(0),
63 fZbeg(0),
64 fXend(0),
65 fYend(0),
66 fZend(0),
67 fXdel(0),
68 fYdel(0),
69 fZdel(0),
70 fXdeli(0),
71 fYdeli(0),
72 fZdeli(0),
73 fXn(0),
74 fYn(0),
75 fZn(0),
76 fWriteEnable(0),
77 fB(0)
84737f5e 78{
79 //
80 // Standard constructor
81 //
84737f5e 82 ReadField();
7b6cddfa 83 SetWriteEnable();
84737f5e 84}
85
e2afb3b6 86//_______________________________________________________________________
84737f5e 87AliFieldMap::~AliFieldMap()
88{
e2afb3b6 89 //
90 // Destructor
91 //
84737f5e 92 delete fB;
93}
94
e2afb3b6 95//_______________________________________________________________________
96AliFieldMap::AliFieldMap(const AliFieldMap &map):
97 TNamed(map),
98 fXbeg(0),
99 fYbeg(0),
100 fZbeg(0),
101 fXend(0),
102 fYend(0),
103 fZend(0),
104 fXdel(0),
105 fYdel(0),
106 fZdel(0),
107 fXdeli(0),
108 fYdeli(0),
109 fZdeli(0),
110 fXn(0),
111 fYn(0),
112 fZn(0),
113 fWriteEnable(0),
114 fB(0)
84737f5e 115{
116 //
117 // Copy constructor
118 //
119 map.Copy(*this);
120}
121
e2afb3b6 122//_______________________________________________________________________
84737f5e 123void AliFieldMap::ReadField()
124{
125 //
126 // Method to read the magnetic field map from file
127 //
128 FILE* magfile;
e2afb3b6 129 // FILE* endf = fopen("end.table", "r");
130 // FILE* out = fopen("out", "w");
84737f5e 131
132 Int_t ix, iy, iz, ipx, ipy, ipz;
133 Float_t bx, by, bz;
134 char *fname = 0;
135 printf("%s: Reading Magnetic Field Map %s from file %s\n",
136 ClassName(),fName.Data(),fTitle.Data());
137
138 fname = gSystem->ExpandPathName(fTitle.Data());
139 magfile = fopen(fname,"r");
140 delete [] fname;
141
142 fscanf(magfile,"%d %d %d %f %f %f %f %f %f",
143 &fXn, &fYn, &fZn, &fXbeg, &fYbeg, &fZbeg, &fXdel, &fYdel, &fZdel);
144
145 fXdeli = 1./fXdel;
146 fYdeli = 1./fYdel;
147 fZdeli = 1./fZdel;
148 fXend = fXbeg + (fXn-1)*fXdel;
149 fYend = fYbeg + (fYn-1)*fYdel;
150 fZend = fZbeg + (fZn-1)*fZdel;
151
152 Int_t nDim = fXn*fYn*fZn;
153
e2afb3b6 154 // Float_t x,y,z,b;
84737f5e 155
156 fB = new TVector(3*nDim);
157 if (magfile) {
158 for (ix = 0; ix < fXn; ix++) {
159 ipx=ix*3*(fZn*fYn);
160 for (iy = 0; iy < fYn; iy++) {
161 ipy=ipx+iy*3*fZn;
162 for (iz = 0; iz < fZn; iz++) {
163 ipz=ipy+iz*3;
164
7b6cddfa 165 if (iz == -1) {
166// fscanf(endf,"%f %f %f", &bx,&by,&bz);
167 } else if (iz > -1) {
84737f5e 168 fscanf(magfile," %f %f %f", &bx, &by, &bz);
7b6cddfa 169 } else {
84737f5e 170 continue;
7b6cddfa 171 }
172
84737f5e 173// fscanf(magfile,"%f %f %f %f %f %f %f ",
174// &x, &y, &z, &bx,&by,&bz, &b);
175// fprintf(out, "%15.8e %15.8e %15.8e \n", bx, by, bz);
176
177 (*fB)(ipz+2) = 10.*bz;
178 (*fB)(ipz+1) = 10.*by;
179 (*fB)(ipz ) = 10.*bx;
180 } //iz
181 } // iy
182 } // ix
183/*
184//
185// this part for interpolation between z = 700 and 720 cm to get
186// z = 710 cm
187//
188 for (ix = 0; ix < fXn; ix++) {
189 ipx=ix*3*(fZn*fYn);
190 for (iy = 0; iy < fYn; iy++) {
191 ipy=ipx+iy*3*fZn;
192 Float_t bxx = (Bx(ix,iy,0) + Bx(ix,iy,2))/2.;
193 Float_t byy = (By(ix,iy,0) + By(ix,iy,2))/2.;
194 Float_t bzz = (Bz(ix,iy,0) + Bz(ix,iy,2))/2.;
195 ipz=ipy+3;
196 (*fB)(ipz+2) = bzz;
197 (*fB)(ipz+1) = byy;
198 (*fB)(ipz ) = bxx;
199 } // iy
200 } // ix
201*/
202 } else {
203 printf("%s: File %s not found !\n",ClassName(),fTitle.Data());
204 exit(1);
205 } // if mafile
206}
207
e2afb3b6 208//_______________________________________________________________________
5d8718b8 209// void AliFieldMap::Field(Float_t *x, Float_t *b)
210// {
211// //
212// // Use simple interpolation to obtain field at point x
213// //
214// Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
215// bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
216// const Double_t kone=1;
217// Int_t ix, iy, iz;
218// b[0]=b[1]=b[2]=0;
219// //
84737f5e 220
5d8718b8 221// xl[0]=TMath::Abs(x[0])-fXbeg;
222// xl[1]=TMath::Abs(x[1])-fYbeg;
223// xl[2]=x[2]-fZbeg;
84737f5e 224
5d8718b8 225// hix=xl[0]*fXdeli;
226// ratx=hix-int(hix);
227// ix=int(hix);
84737f5e 228
5d8718b8 229// hiy=xl[1]*fYdeli;
230// raty=hiy-int(hiy);
231// iy=int(hiy);
84737f5e 232
5d8718b8 233// hiz=xl[2]*fZdeli;
234// ratz=hiz-int(hiz);
235// iz=int(hiz);
236
237// ratx1=kone-ratx;
238// raty1=kone-raty;
239// ratz1=kone-ratz;
240
241// bhyhz = Bx(ix ,iy+1,iz+1)*ratx1+Bx(ix+1,iy+1,iz+1)*ratx;
242// bhylz = Bx(ix ,iy+1,iz )*ratx1+Bx(ix+1,iy+1,iz )*ratx;
243// blyhz = Bx(ix ,iy ,iz+1)*ratx1+Bx(ix+1,iy ,iz+1)*ratx;
244// blylz = Bx(ix ,iy ,iz )*ratx1+Bx(ix+1,iy ,iz )*ratx;
245// bhz = blyhz *raty1+bhyhz *raty;
246// blz = blylz *raty1+bhylz *raty;
247// b[0] = blz *ratz1+bhz *ratz;
248// //
249// bhyhz = By(ix ,iy+1,iz+1)*ratx1+By(ix+1,iy+1,iz+1)*ratx;
250// bhylz = By(ix ,iy+1,iz )*ratx1+By(ix+1,iy+1,iz )*ratx;
251// blyhz = By(ix ,iy ,iz+1)*ratx1+By(ix+1,iy ,iz+1)*ratx;
252// blylz = By(ix ,iy ,iz )*ratx1+By(ix+1,iy ,iz )*ratx;
253// bhz = blyhz *raty1+bhyhz *raty;
254// blz = blylz *raty1+bhylz *raty;
255// b[1] = blz *ratz1+bhz *ratz;
256// //
257// bhyhz = Bz(ix ,iy+1,iz+1)*ratx1+Bz(ix+1,iy+1,iz+1)*ratx;
258// bhylz = Bz(ix ,iy+1,iz )*ratx1+Bz(ix+1,iy+1,iz )*ratx;
259// blyhz = Bz(ix ,iy ,iz+1)*ratx1+Bz(ix+1,iy ,iz+1)*ratx;
260// blylz = Bz(ix ,iy ,iz )*ratx1+Bz(ix+1,iy ,iz )*ratx;
261// bhz = blyhz *raty1+bhyhz *raty;
262// blz = blylz *raty1+bhylz *raty;
263// b[2] = blz *ratz1+bhz *ratz;
264// }
265//_______________________________________________________________________
266void AliFieldMap::Field(Float_t *x, Float_t *b)
267{
268 //
269 // Use simple interpolation to obtain field at point x
270 // Faster version which uses direct access to the TVector
271 // instead of Bx(...), By(...), Bz(...) interfaces.
272 // The interpolation is done first in Z (4 points),
273 // then in Y (2 points), and finally in X direction (one point).
274
275 Double_t h=(x[2]-fZbeg)*fZdeli;
276 Int_t i=int(h); // ix
277 register Int_t index=i;
278 Double_t ratz=h-i;
279 Double_t ratz1=1-ratz;
280
281 h=(TMath::Abs(x[1])-fYbeg)*fYdeli;
282 i=int(h); // iy
283 index+=(i*fZn);
284 Double_t raty=h-i;
285 Double_t raty1=1-raty;
286
287 h=(TMath::Abs(x[0])-fXbeg)*fXdeli;
288 i=int(h); // iz
289 index+=(i*fZn*fYn);
290 index*=3; //index now points to the beginning of the field data (ix,iy,iz)
291 Double_t ratx=h-i;
292 Double_t ratx1=1-ratx;
293
294 // The way we retrive data from fB should be optimized
295 // with respect to the CPU cache
296
297 Float_t bx00=(*fB)[index++]*ratz1;
298 Float_t by00=(*fB)[index++]*ratz1;
299 Float_t bz00=(*fB)[index++]*ratz1;
300 bx00+=(*fB)[index++]*ratz;
301 by00+=(*fB)[index++]*ratz;
302 bz00+=(*fB)[index]*ratz;
303
304 index+=(3*fZn-5);
305
306 Float_t bx01=(*fB)[index++]*ratz1;
307 Float_t by01=(*fB)[index++]*ratz1;
308 Float_t bz01=(*fB)[index++]*ratz1;
309 bx01+=(*fB)[index++]*ratz;
310 by01+=(*fB)[index++]*ratz;
311 bz01+=(*fB)[index]*ratz;
312
313 index+=(3*fYn*fZn-5);
314
315 Float_t bx11=(*fB)[index++]*ratz1;
316 Float_t by11=(*fB)[index++]*ratz1;
317 Float_t bz11=(*fB)[index++]*ratz1;
318 bx11+=(*fB)[index++]*ratz;
319 by11+=(*fB)[index++]*ratz;
320 bz11+=(*fB)[index]*ratz;
321
322 index-=(3*fZn+5);
323
324 Float_t bx10=(*fB)[index++]*ratz1;
325 Float_t by10=(*fB)[index++]*ratz1;
326 Float_t bz10=(*fB)[index++]*ratz1;
327 bx10+=(*fB)[index++]*ratz;
328 by10+=(*fB)[index++]*ratz;
329 bz10+=(*fB)[index]*ratz;
330
331 b[0]= (bx00*raty1+bx01*raty)*ratx1 +(bx10*raty1+bx11*raty)*ratx;
332 b[1]= (by00*raty1+by01*raty)*ratx1 +(by10*raty1+by11*raty)*ratx;
333 b[2]= (bz00*raty1+bz01*raty)*ratx1 +(bz10*raty1+bz11*raty)*ratx;
84737f5e 334}
335
e2afb3b6 336//_______________________________________________________________________
84737f5e 337void AliFieldMap::Copy(AliFieldMap & /* magf */) const
338{
339 //
340 // Copy *this onto magf -- Not implemented
341 //
342 Fatal("Copy","Not implemented!\n");
343}
344
e2afb3b6 345//_______________________________________________________________________
84737f5e 346AliFieldMap & AliFieldMap::operator =(const AliFieldMap &magf)
347{
348 magf.Copy(*this);
349 return *this;
350}
7b6cddfa 351
e2afb3b6 352//_______________________________________________________________________
7b6cddfa 353void AliFieldMap::Streamer(TBuffer &R__b)
354{
355 // Stream an object of class AliFieldMap.
f9aab227 356 TVector* save = 0;
b4cfe01e 357
358 if (R__b.IsReading()) {
359 AliFieldMap::Class()->ReadBuffer(R__b, this);
360 } else {
361 if (!fWriteEnable) {
362 save = fB;
363 fB = 0;
364 }
365 AliFieldMap::Class()->WriteBuffer(R__b, this);
366 if (!fWriteEnable) fB = save;
367 }
7b6cddfa 368}