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