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