alieve_main/alieve_main.cxx
[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
b09247a2 28#include <cstdlib>
a1e17193 29#include <TClass.h>
116cbefd 30#include <TSystem.h>
116cbefd 31
21bf7095 32#include "AliLog.h"
84737f5e 33#include "AliFieldMap.h"
84737f5e 34
35ClassImp(AliFieldMap)
36
e2afb3b6 37//_______________________________________________________________________
38AliFieldMap::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)
84737f5e 56{
57 //
58 // Standard constructor
59 //
7b6cddfa 60 SetWriteEnable();
84737f5e 61}
62
e2afb3b6 63//_______________________________________________________________________
64AliFieldMap::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)
84737f5e 83{
84 //
85 // Standard constructor
86 //
84737f5e 87 ReadField();
7b6cddfa 88 SetWriteEnable();
84737f5e 89}
90
e2afb3b6 91//_______________________________________________________________________
84737f5e 92AliFieldMap::~AliFieldMap()
93{
e2afb3b6 94 //
95 // Destructor
96 //
84737f5e 97 delete fB;
98}
99
e2afb3b6 100//_______________________________________________________________________
101AliFieldMap::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)
84737f5e 120{
121 //
122 // Copy constructor
123 //
124 map.Copy(*this);
125}
126
e2afb3b6 127//_______________________________________________________________________
84737f5e 128void AliFieldMap::ReadField()
129{
130 //
131 // Method to read the magnetic field map from file
132 //
133 FILE* magfile;
e2afb3b6 134 // FILE* endf = fopen("end.table", "r");
135 // FILE* out = fopen("out", "w");
84737f5e 136
137 Int_t ix, iy, iz, ipx, ipy, ipz;
138 Float_t bx, by, bz;
139 char *fname = 0;
21bf7095 140 AliInfo(Form("Reading Magnetic Field Map %s from file %s",
141 fName.Data(),fTitle.Data()));
84737f5e 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
e2afb3b6 159 // Float_t x,y,z,b;
84737f5e 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
7b6cddfa 170 if (iz == -1) {
171// fscanf(endf,"%f %f %f", &bx,&by,&bz);
172 } else if (iz > -1) {
84737f5e 173 fscanf(magfile," %f %f %f", &bx, &by, &bz);
7b6cddfa 174 } else {
84737f5e 175 continue;
7b6cddfa 176 }
177
84737f5e 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
e2afb3b6 213//_______________________________________________________________________
af7ba10c 214void AliFieldMap::Field(Float_t *x, Float_t *b) const
5d8718b8 215{
216 //
217 // Use simple interpolation to obtain field at point x
90dd3681 218 //
219 Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
220 bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
221 const Double_t kone=1;
222 Int_t ix, iy, iz;
223 b[0]=b[1]=b[2]=0;
224 //
225
415bda4f 226 xl[0] = x[0] - fXbeg;
227 xl[1] = x[1] - fYbeg;
228 xl[2] = x[2] - fZbeg;
90dd3681 229
cb3dd940 230 hix=TMath::Max(0.,TMath::Min(xl[0]*fXdeli,fXn-1.0001));
90dd3681 231 ratx=hix-int(hix);
232 ix=int(hix);
233
cb3dd940 234 hiy=TMath::Max(0.,TMath::Min(xl[1]*fYdeli,fYn-1.0001));
90dd3681 235 raty=hiy-int(hiy);
236 iy=int(hiy);
237
cb3dd940 238 hiz=TMath::Max(0.,TMath::Min(xl[2]*fZdeli,fZn-1.0001));
90dd3681 239 ratz=hiz-int(hiz);
240 iz=int(hiz);
241
242 ratx1=kone-ratx;
243 raty1=kone-raty;
244 ratz1=kone-ratz;
245
39b41d10 246 if (!fB) return;
247
90dd3681 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;
84737f5e 271}
272
e2afb3b6 273//_______________________________________________________________________
6c4904c2 274void AliFieldMap::Copy(TObject & /* magf */) const
84737f5e 275{
276 //
277 // Copy *this onto magf -- Not implemented
278 //
21bf7095 279 AliFatal("Not implemented!");
84737f5e 280}
281
e2afb3b6 282//_______________________________________________________________________
84737f5e 283AliFieldMap & AliFieldMap::operator =(const AliFieldMap &magf)
284{
285 magf.Copy(*this);
286 return *this;
287}
7b6cddfa 288
e2afb3b6 289//_______________________________________________________________________
7b6cddfa 290void AliFieldMap::Streamer(TBuffer &R__b)
291{
292 // Stream an object of class AliFieldMap.
f9aab227 293 TVector* save = 0;
b4cfe01e 294
295 if (R__b.IsReading()) {
296 AliFieldMap::Class()->ReadBuffer(R__b, this);
297 } else {
298 if (!fWriteEnable) {
299 save = fB;
300 fB = 0;
301 }
302 AliFieldMap::Class()->WriteBuffer(R__b, this);
303 if (!fWriteEnable) fB = save;
304 }
7b6cddfa 305}