Replacing Header with Id
[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 // Author: Andreas Morsch <andreas.morsch@cern.ch>
21 //
22 //-----------------------------------------------------------------------
23
24 #include <TSystem.h>
25 #include <TVector.h>
26
27 #include "AliFieldMap.h"
28
29 ClassImp(AliFieldMap)
30
31 //_______________________________________________________________________
32 AliFieldMap::AliFieldMap():
33   fXbeg(0),
34   fYbeg(0),
35   fZbeg(0),
36   fXend(0),
37   fYend(0),
38   fZend(0),
39   fXdel(0),
40   fYdel(0),
41   fZdel(0),
42   fXdeli(0),
43   fYdeli(0),
44   fZdeli(0),
45   fXn(0),
46   fYn(0),
47   fZn(0),
48   fWriteEnable(0),
49   fB(0)
50 {
51   //
52   // Standard constructor
53   //
54   SetWriteEnable();
55 }
56
57 //_______________________________________________________________________
58 AliFieldMap::AliFieldMap(const char *name, const char *title):
59   TNamed(name,title),
60   fXbeg(0),
61   fYbeg(0),
62   fZbeg(0),
63   fXend(0),
64   fYend(0),
65   fZend(0),
66   fXdel(0),
67   fYdel(0),
68   fZdel(0),
69   fXdeli(0),
70   fYdeli(0),
71   fZdeli(0),
72   fXn(0),
73   fYn(0),
74   fZn(0),
75   fWriteEnable(0),
76   fB(0)
77 {
78   //
79   // Standard constructor
80   //
81   ReadField();
82   SetWriteEnable();
83 }
84
85 //_______________________________________________________________________
86 AliFieldMap::~AliFieldMap()
87 {
88   //
89   // Destructor
90   //  
91   delete fB;
92 }
93
94 //_______________________________________________________________________
95 AliFieldMap::AliFieldMap(const AliFieldMap &map):
96   TNamed(map),
97   fXbeg(0),
98   fYbeg(0),
99   fZbeg(0),
100   fXend(0),
101   fYend(0),
102   fZend(0),
103   fXdel(0),
104   fYdel(0),
105   fZdel(0),
106   fXdeli(0),
107   fYdeli(0),
108   fZdeli(0),
109   fXn(0),
110   fYn(0),
111   fZn(0),
112   fWriteEnable(0),
113   fB(0)
114 {
115   //
116   // Copy constructor
117   //
118   map.Copy(*this);
119 }
120
121 //_______________________________________________________________________
122 void AliFieldMap::ReadField()
123 {
124   // 
125   // Method to read the magnetic field map from file
126   //
127   FILE* magfile;
128   //  FILE* endf = fopen("end.table", "r");
129   //  FILE* out  = fopen("out", "w");
130   
131   Int_t   ix, iy, iz, ipx, ipy, ipz;
132   Float_t bx, by, bz;
133   char *fname = 0;
134   printf("%s: Reading Magnetic Field Map %s from file %s\n",
135          ClassName(),fName.Data(),fTitle.Data());
136
137   fname   = gSystem->ExpandPathName(fTitle.Data());
138   magfile = fopen(fname,"r");
139   delete [] fname;
140
141   fscanf(magfile,"%d %d %d %f %f %f %f %f %f", 
142          &fXn, &fYn, &fZn, &fXbeg, &fYbeg, &fZbeg, &fXdel, &fYdel, &fZdel);
143  
144   fXdeli = 1./fXdel;
145   fYdeli = 1./fYdel;    
146   fZdeli = 1./fZdel;    
147   fXend  = fXbeg + (fXn-1)*fXdel;
148   fYend  = fYbeg + (fYn-1)*fYdel;
149   fZend  = fZbeg + (fZn-1)*fZdel;
150   
151   Int_t nDim   = fXn*fYn*fZn;
152
153   //  Float_t x,y,z,b;
154
155   fB = new TVector(3*nDim);
156   if (magfile) {
157       for (ix = 0; ix < fXn; ix++) {
158           ipx=ix*3*(fZn*fYn);
159           for (iy = 0; iy < fYn; iy++) {
160               ipy=ipx+iy*3*fZn;
161               for (iz = 0; iz < fZn; iz++) {
162                   ipz=ipy+iz*3;
163
164                   if (iz == -1) {
165 //                    fscanf(endf,"%f %f %f", &bx,&by,&bz);
166                   } else if (iz > -1) {
167                       fscanf(magfile," %f %f %f", &bx, &by, &bz);
168                   } else {
169                       continue;
170                   }
171                   
172 //                fscanf(magfile,"%f %f %f %f %f %f %f ",
173 //                       &x, &y, &z, &bx,&by,&bz, &b);
174 //                fprintf(out, "%15.8e %15.8e %15.8e \n", bx, by, bz);
175
176                   (*fB)(ipz+2) = 10.*bz;
177                   (*fB)(ipz+1) = 10.*by;
178                   (*fB)(ipz  ) = 10.*bx;
179               } //iz
180           } // iy
181       } // ix
182 /*
183 //
184 // this part for interpolation between z = 700 and 720 cm to get
185 // z = 710 cm
186 //      
187       for (ix = 0; ix < fXn; ix++) {
188           ipx=ix*3*(fZn*fYn);
189           for (iy = 0; iy < fYn; iy++) {
190               ipy=ipx+iy*3*fZn;
191               Float_t bxx = (Bx(ix,iy,0) + Bx(ix,iy,2))/2.;
192               Float_t byy = (By(ix,iy,0) + By(ix,iy,2))/2.;
193               Float_t bzz = (Bz(ix,iy,0) + Bz(ix,iy,2))/2.;           
194               ipz=ipy+3;
195               (*fB)(ipz+2) = bzz;
196               (*fB)(ipz+1) = byy;
197               (*fB)(ipz  ) = bxx;
198           } // iy
199       } // ix
200 */      
201   } else { 
202       printf("%s: File %s not found !\n",ClassName(),fTitle.Data());
203       exit(1);
204   } // if mafile
205 }
206
207 //_______________________________________________________________________
208 void AliFieldMap::Field(Float_t *x, Float_t *b)
209 {
210   //
211   // Use simple interpolation to obtain field at point x
212   //
213     Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1, 
214         bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
215     const Double_t kone=1;
216     Int_t ix, iy, iz;
217     b[0]=b[1]=b[2]=0;
218     //
219     
220     xl[0]=TMath::Abs(x[0])-fXbeg;
221     xl[1]=TMath::Abs(x[1])-fYbeg;
222     xl[2]=x[2]-fZbeg;
223     
224     hix=xl[0]*fXdeli;
225     ratx=hix-int(hix);
226     ix=int(hix);
227     
228     hiy=xl[1]*fYdeli;
229     raty=hiy-int(hiy);
230     iy=int(hiy);
231     
232     hiz=xl[2]*fZdeli;
233     ratz=hiz-int(hiz);
234     iz=int(hiz);
235
236     ratx1=kone-ratx;
237     raty1=kone-raty;
238     ratz1=kone-ratz;
239
240     bhyhz = Bx(ix  ,iy+1,iz+1)*ratx1+Bx(ix+1,iy+1,iz+1)*ratx;
241     bhylz = Bx(ix  ,iy+1,iz  )*ratx1+Bx(ix+1,iy+1,iz  )*ratx;
242     blyhz = Bx(ix  ,iy  ,iz+1)*ratx1+Bx(ix+1,iy  ,iz+1)*ratx;
243     blylz = Bx(ix  ,iy  ,iz  )*ratx1+Bx(ix+1,iy  ,iz  )*ratx;
244     bhz   = blyhz             *raty1+bhyhz             *raty;
245     blz   = blylz             *raty1+bhylz             *raty;
246     b[0]  = blz               *ratz1+bhz               *ratz;
247     //
248     bhyhz = By(ix  ,iy+1,iz+1)*ratx1+By(ix+1,iy+1,iz+1)*ratx;
249     bhylz = By(ix  ,iy+1,iz  )*ratx1+By(ix+1,iy+1,iz  )*ratx;
250     blyhz = By(ix  ,iy  ,iz+1)*ratx1+By(ix+1,iy  ,iz+1)*ratx;
251     blylz = By(ix  ,iy  ,iz  )*ratx1+By(ix+1,iy  ,iz  )*ratx;
252     bhz   = blyhz             *raty1+bhyhz             *raty;
253     blz   = blylz             *raty1+bhylz             *raty;
254     b[1]  = blz               *ratz1+bhz               *ratz;
255     //
256     bhyhz = Bz(ix  ,iy+1,iz+1)*ratx1+Bz(ix+1,iy+1,iz+1)*ratx;
257     bhylz = Bz(ix  ,iy+1,iz  )*ratx1+Bz(ix+1,iy+1,iz  )*ratx;
258     blyhz = Bz(ix  ,iy  ,iz+1)*ratx1+Bz(ix+1,iy  ,iz+1)*ratx;
259     blylz = Bz(ix  ,iy  ,iz  )*ratx1+Bz(ix+1,iy  ,iz  )*ratx;
260     bhz   = blyhz             *raty1+bhyhz             *raty;
261     blz   = blylz             *raty1+bhylz             *raty;
262     b[2]  = blz               *ratz1+bhz               *ratz;
263 }
264
265 //_______________________________________________________________________
266 void AliFieldMap::Copy(AliFieldMap & /* magf */) const
267 {
268   //
269   // Copy *this onto magf -- Not implemented
270   //
271   Fatal("Copy","Not implemented!\n");
272 }
273
274 //_______________________________________________________________________
275 AliFieldMap & AliFieldMap::operator =(const AliFieldMap &magf)
276 {
277   magf.Copy(*this);
278   return *this;
279 }
280
281 //_______________________________________________________________________
282 void AliFieldMap::Streamer(TBuffer &R__b)
283 {
284    // Stream an object of class AliFieldMap.
285     TVector* save = 0;
286     
287     if (R__b.IsReading()) {
288         AliFieldMap::Class()->ReadBuffer(R__b, this);
289     } else {
290         if (!fWriteEnable) {
291             save = fB;
292             fB = 0;
293         }
294         AliFieldMap::Class()->WriteBuffer(R__b, this);
295         if (!fWriteEnable) fB = save;
296     }
297 }