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