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