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