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