update for AD from M.Broz, see ALIROOT-5766
[u/mrichter/AliRoot.git] / STEER / AliFieldMap.cxx
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
16 /* $Id$ */
17
18 //-----------------------------------------------------------------------
19 //
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
24 // Author: Andreas Morsch <andreas.morsch@cern.ch>
25 //
26 //-----------------------------------------------------------------------
27
28 #include <cstdlib>
29 #include <TClass.h>
30 #include <TSystem.h>
31
32 #include "AliLog.h"
33 #include "AliFieldMap.h"
34
35 ClassImp(AliFieldMap)
36
37 //_______________________________________________________________________
38 AliFieldMap::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)
56 {
57   //
58   // Standard constructor
59   //
60   SetWriteEnable();
61 }
62
63 //_______________________________________________________________________
64 AliFieldMap::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)
83 {
84   //
85   // Standard constructor
86   //
87   ReadField();
88   SetWriteEnable();
89 }
90
91 //_______________________________________________________________________
92 AliFieldMap::~AliFieldMap()
93 {
94   //
95   // Destructor
96   //  
97   delete fB;
98 }
99
100 //_______________________________________________________________________
101 AliFieldMap::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)
120 {
121   //
122   // Copy constructor
123   //
124   map.Copy(*this);
125 }
126
127 //_______________________________________________________________________
128 void AliFieldMap::ReadField()
129 {
130   // 
131   // Method to read the magnetic field map from file
132   //
133   FILE* magfile;
134   //  FILE* endf = fopen("end.table", "r");
135   //  FILE* out  = fopen("out", "w");
136   
137   Int_t   ix, iy, iz, ipx, ipy, ipz;
138   Float_t bx, by, bz;
139   char *fname = 0;
140   AliInfo(Form("Reading Magnetic Field Map %s from file %s",
141                fName.Data(),fTitle.Data()));
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
159   //  Float_t x,y,z,b;
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
170                   if (iz == -1) {
171 //                    fscanf(endf,"%f %f %f", &bx,&by,&bz);
172                   } else if (iz > -1) {
173                       fscanf(magfile," %f %f %f", &bx, &by, &bz);
174                   } else {
175                       continue;
176                   }
177                   
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
213 //_______________________________________________________________________
214 void 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 //_______________________________________________________________________
274 void AliFieldMap::Field(double *x, double *b) const
275 {
276   //
277   // Use simple interpolation to obtain field at point x
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     
286     xl[0] = x[0] - fXbeg;
287     xl[1] = x[1] - fYbeg;
288     xl[2] = x[2] - fZbeg;
289     
290     hix=TMath::Max(0.,TMath::Min(xl[0]*fXdeli,fXn-1.0001));
291     ratx=hix-int(hix);
292     ix=int(hix);
293     
294     hiy=TMath::Max(0.,TMath::Min(xl[1]*fYdeli,fYn-1.0001));
295     raty=hiy-int(hiy);
296     iy=int(hiy);
297     
298     hiz=TMath::Max(0.,TMath::Min(xl[2]*fZdeli,fZn-1.0001));
299     ratz=hiz-int(hiz);
300     iz=int(hiz);
301
302     ratx1=kone-ratx;
303     raty1=kone-raty;
304     ratz1=kone-ratz;
305
306     if (!fB) return;
307
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;
331 }
332
333 //_______________________________________________________________________
334 void AliFieldMap::Copy(TObject & /* magf */) const
335 {
336   //
337   // Copy *this onto magf -- Not implemented
338   //
339   AliFatal("Not implemented!");
340 }
341
342 //_______________________________________________________________________
343 AliFieldMap & AliFieldMap::operator =(const AliFieldMap &magf)
344 {
345   magf.Copy(*this);
346   return *this;
347 }
348
349 //_______________________________________________________________________
350 void AliFieldMap::Streamer(TBuffer &R__b)
351 {
352    // Stream an object of class AliFieldMap.
353     TVector* save = 0;
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     }
365 }