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