Add directory structure (needed if new in cvmfs)
[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//-----------------------------------------------------------------------
90dd3681 19//
af7ba10c 20// Class to handle the field
21// I/O and interpolation
22// of the field map to return the value
23// of the magnetic field in an arbitrary position
84737f5e 24// Author: Andreas Morsch <andreas.morsch@cern.ch>
90dd3681 25//
116cbefd 26//-----------------------------------------------------------------------
84737f5e 27
b09247a2 28#include <cstdlib>
a1e17193 29#include <TClass.h>
116cbefd 30#include <TSystem.h>
116cbefd 31
21bf7095 32#include "AliLog.h"
84737f5e 33#include "AliFieldMap.h"
84737f5e 34
35ClassImp(AliFieldMap)
36
e2afb3b6 37//_______________________________________________________________________
38AliFieldMap::AliFieldMap():
39 fXbeg(0),
40 fYbeg(0),
41 fZbeg(0),
42 fXend(0),
43 fYend(0),
44 fZend(0),
45 fXdel(0),
46 fYdel(0),
47 fZdel(0),
48 fXdeli(0),
49 fYdeli(0),
50 fZdeli(0),
51 fXn(0),
52 fYn(0),
53 fZn(0),
54 fWriteEnable(0),
55 fB(0)
84737f5e 56{
57 //
58 // Standard constructor
59 //
7b6cddfa 60 SetWriteEnable();
84737f5e 61}
62
e2afb3b6 63//_______________________________________________________________________
64AliFieldMap::AliFieldMap(const char *name, const char *title):
65 TNamed(name,title),
66 fXbeg(0),
67 fYbeg(0),
68 fZbeg(0),
69 fXend(0),
70 fYend(0),
71 fZend(0),
72 fXdel(0),
73 fYdel(0),
74 fZdel(0),
75 fXdeli(0),
76 fYdeli(0),
77 fZdeli(0),
78 fXn(0),
79 fYn(0),
80 fZn(0),
81 fWriteEnable(0),
82 fB(0)
84737f5e 83{
84 //
85 // Standard constructor
86 //
84737f5e 87 ReadField();
7b6cddfa 88 SetWriteEnable();
84737f5e 89}
90
e2afb3b6 91//_______________________________________________________________________
84737f5e 92AliFieldMap::~AliFieldMap()
93{
e2afb3b6 94 //
95 // Destructor
96 //
84737f5e 97 delete fB;
98}
99
e2afb3b6 100//_______________________________________________________________________
101AliFieldMap::AliFieldMap(const AliFieldMap &map):
102 TNamed(map),
103 fXbeg(0),
104 fYbeg(0),
105 fZbeg(0),
106 fXend(0),
107 fYend(0),
108 fZend(0),
109 fXdel(0),
110 fYdel(0),
111 fZdel(0),
112 fXdeli(0),
113 fYdeli(0),
114 fZdeli(0),
115 fXn(0),
116 fYn(0),
117 fZn(0),
118 fWriteEnable(0),
119 fB(0)
84737f5e 120{
121 //
122 // Copy constructor
123 //
124 map.Copy(*this);
125}
126
e2afb3b6 127//_______________________________________________________________________
84737f5e 128void AliFieldMap::ReadField()
129{
130 //
131 // Method to read the magnetic field map from file
132 //
133 FILE* magfile;
e2afb3b6 134 // FILE* endf = fopen("end.table", "r");
135 // FILE* out = fopen("out", "w");
84737f5e 136
137 Int_t ix, iy, iz, ipx, ipy, ipz;
138 Float_t bx, by, bz;
139 char *fname = 0;
21bf7095 140 AliInfo(Form("Reading Magnetic Field Map %s from file %s",
141 fName.Data(),fTitle.Data()));
84737f5e 142
143 fname = gSystem->ExpandPathName(fTitle.Data());
144 magfile = fopen(fname,"r");
145 delete [] fname;
146
147 fscanf(magfile,"%d %d %d %f %f %f %f %f %f",
148 &fXn, &fYn, &fZn, &fXbeg, &fYbeg, &fZbeg, &fXdel, &fYdel, &fZdel);
149
150 fXdeli = 1./fXdel;
151 fYdeli = 1./fYdel;
152 fZdeli = 1./fZdel;
153 fXend = fXbeg + (fXn-1)*fXdel;
154 fYend = fYbeg + (fYn-1)*fYdel;
155 fZend = fZbeg + (fZn-1)*fZdel;
156
157 Int_t nDim = fXn*fYn*fZn;
158
e2afb3b6 159 // Float_t x,y,z,b;
84737f5e 160
161 fB = new TVector(3*nDim);
162 if (magfile) {
163 for (ix = 0; ix < fXn; ix++) {
164 ipx=ix*3*(fZn*fYn);
165 for (iy = 0; iy < fYn; iy++) {
166 ipy=ipx+iy*3*fZn;
167 for (iz = 0; iz < fZn; iz++) {
168 ipz=ipy+iz*3;
169
7b6cddfa 170 if (iz == -1) {
171// fscanf(endf,"%f %f %f", &bx,&by,&bz);
172 } else if (iz > -1) {
84737f5e 173 fscanf(magfile," %f %f %f", &bx, &by, &bz);
7b6cddfa 174 } else {
84737f5e 175 continue;
7b6cddfa 176 }
177
84737f5e 178// fscanf(magfile,"%f %f %f %f %f %f %f ",
179// &x, &y, &z, &bx,&by,&bz, &b);
180// fprintf(out, "%15.8e %15.8e %15.8e \n", bx, by, bz);
181
182 (*fB)(ipz+2) = 10.*bz;
183 (*fB)(ipz+1) = 10.*by;
184 (*fB)(ipz ) = 10.*bx;
185 } //iz
186 } // iy
187 } // ix
188/*
189//
190// this part for interpolation between z = 700 and 720 cm to get
191// z = 710 cm
192//
193 for (ix = 0; ix < fXn; ix++) {
194 ipx=ix*3*(fZn*fYn);
195 for (iy = 0; iy < fYn; iy++) {
196 ipy=ipx+iy*3*fZn;
197 Float_t bxx = (Bx(ix,iy,0) + Bx(ix,iy,2))/2.;
198 Float_t byy = (By(ix,iy,0) + By(ix,iy,2))/2.;
199 Float_t bzz = (Bz(ix,iy,0) + Bz(ix,iy,2))/2.;
200 ipz=ipy+3;
201 (*fB)(ipz+2) = bzz;
202 (*fB)(ipz+1) = byy;
203 (*fB)(ipz ) = bxx;
204 } // iy
205 } // ix
206*/
207 } else {
208 printf("%s: File %s not found !\n",ClassName(),fTitle.Data());
209 exit(1);
210 } // if mafile
211}
212
e2afb3b6 213//_______________________________________________________________________
ff66b122 214void AliFieldMap::Field(float *x, float *b) const
215{
216 //
217 // Use simple interpolation to obtain field at point x
218 //
219 float ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
220 bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
221 const float kone=1;
222 Int_t ix, iy, iz;
223 b[0]=b[1]=b[2]=0;
224 //
225
226 xl[0] = x[0] - fXbeg;
227 xl[1] = x[1] - fYbeg;
228 xl[2] = x[2] - fZbeg;
229
230 hix=TMath::Max(0.,TMath::Min(xl[0]*fXdeli,fXn-1.0001));
231 ratx=hix-int(hix);
232 ix=int(hix);
233
234 hiy=TMath::Max(0.,TMath::Min(xl[1]*fYdeli,fYn-1.0001));
235 raty=hiy-int(hiy);
236 iy=int(hiy);
237
238 hiz=TMath::Max(0.,TMath::Min(xl[2]*fZdeli,fZn-1.0001));
239 ratz=hiz-int(hiz);
240 iz=int(hiz);
241
242 ratx1=kone-ratx;
243 raty1=kone-raty;
244 ratz1=kone-ratz;
245
246 if (!fB) return;
247
248 bhyhz = Bx(ix ,iy+1,iz+1)*ratx1+Bx(ix+1,iy+1,iz+1)*ratx;
249 bhylz = Bx(ix ,iy+1,iz )*ratx1+Bx(ix+1,iy+1,iz )*ratx;
250 blyhz = Bx(ix ,iy ,iz+1)*ratx1+Bx(ix+1,iy ,iz+1)*ratx;
251 blylz = Bx(ix ,iy ,iz )*ratx1+Bx(ix+1,iy ,iz )*ratx;
252 bhz = blyhz *raty1+bhyhz *raty;
253 blz = blylz *raty1+bhylz *raty;
254 b[0] = blz *ratz1+bhz *ratz;
255 //
256 bhyhz = By(ix ,iy+1,iz+1)*ratx1+By(ix+1,iy+1,iz+1)*ratx;
257 bhylz = By(ix ,iy+1,iz )*ratx1+By(ix+1,iy+1,iz )*ratx;
258 blyhz = By(ix ,iy ,iz+1)*ratx1+By(ix+1,iy ,iz+1)*ratx;
259 blylz = By(ix ,iy ,iz )*ratx1+By(ix+1,iy ,iz )*ratx;
260 bhz = blyhz *raty1+bhyhz *raty;
261 blz = blylz *raty1+bhylz *raty;
262 b[1] = blz *ratz1+bhz *ratz;
263 //
264 bhyhz = Bz(ix ,iy+1,iz+1)*ratx1+Bz(ix+1,iy+1,iz+1)*ratx;
265 bhylz = Bz(ix ,iy+1,iz )*ratx1+Bz(ix+1,iy+1,iz )*ratx;
266 blyhz = Bz(ix ,iy ,iz+1)*ratx1+Bz(ix+1,iy ,iz+1)*ratx;
267 blylz = Bz(ix ,iy ,iz )*ratx1+Bz(ix+1,iy ,iz )*ratx;
268 bhz = blyhz *raty1+bhyhz *raty;
269 blz = blylz *raty1+bhylz *raty;
270 b[2] = blz *ratz1+bhz *ratz;
271}
272
273//_______________________________________________________________________
274void AliFieldMap::Field(double *x, double *b) const
5d8718b8 275{
276 //
277 // Use simple interpolation to obtain field at point x
90dd3681 278 //
279 Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
280 bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
281 const Double_t kone=1;
282 Int_t ix, iy, iz;
283 b[0]=b[1]=b[2]=0;
284 //
285
415bda4f 286 xl[0] = x[0] - fXbeg;
287 xl[1] = x[1] - fYbeg;
288 xl[2] = x[2] - fZbeg;
90dd3681 289
cb3dd940 290 hix=TMath::Max(0.,TMath::Min(xl[0]*fXdeli,fXn-1.0001));
90dd3681 291 ratx=hix-int(hix);
292 ix=int(hix);
293
cb3dd940 294 hiy=TMath::Max(0.,TMath::Min(xl[1]*fYdeli,fYn-1.0001));
90dd3681 295 raty=hiy-int(hiy);
296 iy=int(hiy);
297
cb3dd940 298 hiz=TMath::Max(0.,TMath::Min(xl[2]*fZdeli,fZn-1.0001));
90dd3681 299 ratz=hiz-int(hiz);
300 iz=int(hiz);
301
302 ratx1=kone-ratx;
303 raty1=kone-raty;
304 ratz1=kone-ratz;
305
39b41d10 306 if (!fB) return;
307
90dd3681 308 bhyhz = Bx(ix ,iy+1,iz+1)*ratx1+Bx(ix+1,iy+1,iz+1)*ratx;
309 bhylz = Bx(ix ,iy+1,iz )*ratx1+Bx(ix+1,iy+1,iz )*ratx;
310 blyhz = Bx(ix ,iy ,iz+1)*ratx1+Bx(ix+1,iy ,iz+1)*ratx;
311 blylz = Bx(ix ,iy ,iz )*ratx1+Bx(ix+1,iy ,iz )*ratx;
312 bhz = blyhz *raty1+bhyhz *raty;
313 blz = blylz *raty1+bhylz *raty;
314 b[0] = blz *ratz1+bhz *ratz;
315 //
316 bhyhz = By(ix ,iy+1,iz+1)*ratx1+By(ix+1,iy+1,iz+1)*ratx;
317 bhylz = By(ix ,iy+1,iz )*ratx1+By(ix+1,iy+1,iz )*ratx;
318 blyhz = By(ix ,iy ,iz+1)*ratx1+By(ix+1,iy ,iz+1)*ratx;
319 blylz = By(ix ,iy ,iz )*ratx1+By(ix+1,iy ,iz )*ratx;
320 bhz = blyhz *raty1+bhyhz *raty;
321 blz = blylz *raty1+bhylz *raty;
322 b[1] = blz *ratz1+bhz *ratz;
323 //
324 bhyhz = Bz(ix ,iy+1,iz+1)*ratx1+Bz(ix+1,iy+1,iz+1)*ratx;
325 bhylz = Bz(ix ,iy+1,iz )*ratx1+Bz(ix+1,iy+1,iz )*ratx;
326 blyhz = Bz(ix ,iy ,iz+1)*ratx1+Bz(ix+1,iy ,iz+1)*ratx;
327 blylz = Bz(ix ,iy ,iz )*ratx1+Bz(ix+1,iy ,iz )*ratx;
328 bhz = blyhz *raty1+bhyhz *raty;
329 blz = blylz *raty1+bhylz *raty;
330 b[2] = blz *ratz1+bhz *ratz;
84737f5e 331}
332
e2afb3b6 333//_______________________________________________________________________
6c4904c2 334void AliFieldMap::Copy(TObject & /* magf */) const
84737f5e 335{
336 //
337 // Copy *this onto magf -- Not implemented
338 //
21bf7095 339 AliFatal("Not implemented!");
84737f5e 340}
341
e2afb3b6 342//_______________________________________________________________________
84737f5e 343AliFieldMap & AliFieldMap::operator =(const AliFieldMap &magf)
344{
345 magf.Copy(*this);
346 return *this;
347}
7b6cddfa 348
e2afb3b6 349//_______________________________________________________________________
7b6cddfa 350void AliFieldMap::Streamer(TBuffer &R__b)
351{
352 // Stream an object of class AliFieldMap.
f9aab227 353 TVector* save = 0;
b4cfe01e 354
355 if (R__b.IsReading()) {
356 AliFieldMap::Class()->ReadBuffer(R__b, this);
357 } else {
358 if (!fWriteEnable) {
359 save = fB;
360 fB = 0;
361 }
362 AliFieldMap::Class()->WriteBuffer(R__b, this);
363 if (!fWriteEnable) fB = save;
364 }
7b6cddfa 365}