No ISO/ANSI compliance for the old gcc compilers (2.95.2,2.96)
[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//
84737f5e 20// Author: Andreas Morsch <andreas.morsch@cern.ch>
90dd3681 21//
116cbefd 22//-----------------------------------------------------------------------
84737f5e 23
116cbefd 24#include <TSystem.h>
84737f5e 25#include <TVector.h>
116cbefd 26
84737f5e 27#include "AliFieldMap.h"
84737f5e 28
29ClassImp(AliFieldMap)
30
e2afb3b6 31//_______________________________________________________________________
32AliFieldMap::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)
84737f5e 50{
51 //
52 // Standard constructor
53 //
7b6cddfa 54 SetWriteEnable();
84737f5e 55}
56
e2afb3b6 57//_______________________________________________________________________
58AliFieldMap::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)
84737f5e 77{
78 //
79 // Standard constructor
80 //
84737f5e 81 ReadField();
7b6cddfa 82 SetWriteEnable();
84737f5e 83}
84
e2afb3b6 85//_______________________________________________________________________
84737f5e 86AliFieldMap::~AliFieldMap()
87{
e2afb3b6 88 //
89 // Destructor
90 //
84737f5e 91 delete fB;
92}
93
e2afb3b6 94//_______________________________________________________________________
95AliFieldMap::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)
84737f5e 114{
115 //
116 // Copy constructor
117 //
118 map.Copy(*this);
119}
120
e2afb3b6 121//_______________________________________________________________________
84737f5e 122void AliFieldMap::ReadField()
123{
124 //
125 // Method to read the magnetic field map from file
126 //
127 FILE* magfile;
e2afb3b6 128 // FILE* endf = fopen("end.table", "r");
129 // FILE* out = fopen("out", "w");
84737f5e 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
e2afb3b6 153 // Float_t x,y,z,b;
84737f5e 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
7b6cddfa 164 if (iz == -1) {
165// fscanf(endf,"%f %f %f", &bx,&by,&bz);
166 } else if (iz > -1) {
84737f5e 167 fscanf(magfile," %f %f %f", &bx, &by, &bz);
7b6cddfa 168 } else {
84737f5e 169 continue;
7b6cddfa 170 }
171
84737f5e 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
e2afb3b6 207//_______________________________________________________________________
5d8718b8 208void AliFieldMap::Field(Float_t *x, Float_t *b)
209{
210 //
211 // Use simple interpolation to obtain field at point x
90dd3681 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;
84737f5e 263}
264
e2afb3b6 265//_______________________________________________________________________
6c4904c2 266void AliFieldMap::Copy(TObject & /* magf */) const
84737f5e 267{
268 //
269 // Copy *this onto magf -- Not implemented
270 //
271 Fatal("Copy","Not implemented!\n");
272}
273
e2afb3b6 274//_______________________________________________________________________
84737f5e 275AliFieldMap & AliFieldMap::operator =(const AliFieldMap &magf)
276{
277 magf.Copy(*this);
278 return *this;
279}
7b6cddfa 280
e2afb3b6 281//_______________________________________________________________________
7b6cddfa 282void AliFieldMap::Streamer(TBuffer &R__b)
283{
284 // Stream an object of class AliFieldMap.
f9aab227 285 TVector* save = 0;
b4cfe01e 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 }
7b6cddfa 297}