Introducing Header instead of Log
[u/mrichter/AliRoot.git] / STEER / AliMagFCM.cxx
CommitLineData
aee8290b 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
fb17acd4 16/* $Header$ */
e2afb3b6 17
fb17acd4 18#include "TSystem.h"
65fb704d 19#include "TVector.h"
aee8290b 20
21#include "AliMagFCM.h"
aee8290b 22
23ClassImp(AliMagFCM)
24
e2afb3b6 25//_______________________________________________________________________
26AliMagFCM::AliMagFCM():
27 fXbeg(0),
28 fYbeg(0),
29 fZbeg(0),
30 fXdel(0),
31 fYdel(0),
32 fZdel(0),
33 fSolenoid(0),
34 fXdeli(0),
35 fYdeli(0),
36 fZdeli(0),
37 fXn(0),
38 fYn(0),
39 fZn(0),
40 fB(0)
41{
42 //
43 // Standard constructor
44 //
45 fType = kConMesh;
46 fMap = 2;
47 SetSolenoidField();
48}
49
50//_______________________________________________________________________
aee8290b 51AliMagFCM::AliMagFCM(const char *name, const char *title, const Int_t integ,
e2afb3b6 52 const Float_t factor, const Float_t fmax):
53 AliMagF(name,title,integ,factor,fmax),
54 fXbeg(0),
55 fYbeg(0),
56 fZbeg(0),
57 fXdel(0),
58 fYdel(0),
59 fZdel(0),
60 fSolenoid(0),
61 fXdeli(0),
62 fYdeli(0),
63 fZdeli(0),
64 fXn(0),
65 fYn(0),
66 fZn(0),
67 fB(0)
aee8290b 68{
69 //
70 // Standard constructor
71 //
72 fType = kConMesh;
d8408e76 73 fMap = 2;
4cc8933f 74 SetSolenoidField();
75
e2afb3b6 76 if(fDebug>-1) Info("ctor",
77 "%s: Constant Mesh Field %s created: map= %d, factor= %f, file= %s\n",
9e1a0ddb 78 ClassName(),fName.Data(), fMap, factor,fTitle.Data());
aee8290b 79}
80
e2afb3b6 81//_______________________________________________________________________
82AliMagFCM::AliMagFCM(const AliMagFCM &magf):
83 AliMagF(magf),
84 fXbeg(0),
85 fYbeg(0),
86 fZbeg(0),
87 fXdel(0),
88 fYdel(0),
89 fZdel(0),
90 fSolenoid(0),
91 fXdeli(0),
92 fYdeli(0),
93 fZdeli(0),
94 fXn(0),
95 fYn(0),
96 fZn(0),
97 fB(0)
aee8290b 98{
99 //
100 // Copy constructor
101 //
102 magf.Copy(*this);
103}
104
e2afb3b6 105//_______________________________________________________________________
aee8290b 106void AliMagFCM::Field(Float_t *x, Float_t *b)
107{
108 //
109 // Method to calculate the magnetic field
110 //
111 Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
112 bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
113 const Double_t kone=1;
114 Int_t ix, iy, iz;
115
116 // --- find the position in the grid ---
117
118 b[0]=b[1]=b[2]=0;
119 if(-700<x[2] && x[2]<fZbeg && x[0]*x[0]+(x[1]+30)*(x[1]+30) < 560*560) {
4cc8933f 120 b[2]= fSolenoid;
aee8290b 121 } else {
122 Bool_t infield=(fZbeg<=x[2] && x[2]<fZbeg+fZdel*(fZn-1)
123 && ( fXbeg <= TMath::Abs(x[0]) && TMath::Abs(x[0]) < fXbeg+fXdel*(fXn-1) )
124 && ( fYbeg <= TMath::Abs(x[1]) && TMath::Abs(x[1]) < fYbeg+fYdel*(fYn-1) ));
125 if(infield) {
126 xl[0]=TMath::Abs(x[0])-fXbeg;
127 xl[1]=TMath::Abs(x[1])-fYbeg;
128 xl[2]=x[2]-fZbeg;
129
130 // --- start with x
131
132 hix=xl[0]*fXdeli;
133 ratx=hix-int(hix);
134 ix=int(hix);
135
136 hiy=xl[1]*fYdeli;
137 raty=hiy-int(hiy);
138 iy=int(hiy);
139
140 hiz=xl[2]*fZdeli;
141 ratz=hiz-int(hiz);
142 iz=int(hiz);
143
144 if(fMap==2) {
145 // ... simple interpolation
146 ratx1=kone-ratx;
147 raty1=kone-raty;
148 ratz1=kone-ratz;
149 bhyhz = Bx(ix ,iy+1,iz+1)*ratx1+Bx(ix+1,iy+1,iz+1)*ratx;
150 bhylz = Bx(ix ,iy+1,iz )*ratx1+Bx(ix+1,iy+1,iz )*ratx;
151 blyhz = Bx(ix ,iy ,iz+1)*ratx1+Bx(ix+1,iy ,iz+1)*ratx;
152 blylz = Bx(ix ,iy ,iz )*ratx1+Bx(ix+1,iy ,iz )*ratx;
153 bhz = blyhz *raty1+bhyhz *raty;
154 blz = blylz *raty1+bhylz *raty;
155 b[0] = blz *ratz1+bhz *ratz;
156 //
157 bhyhz = By(ix ,iy+1,iz+1)*ratx1+By(ix+1,iy+1,iz+1)*ratx;
158 bhylz = By(ix ,iy+1,iz )*ratx1+By(ix+1,iy+1,iz )*ratx;
159 blyhz = By(ix ,iy ,iz+1)*ratx1+By(ix+1,iy ,iz+1)*ratx;
160 blylz = By(ix ,iy ,iz )*ratx1+By(ix+1,iy ,iz )*ratx;
161 bhz = blyhz *raty1+bhyhz *raty;
162 blz = blylz *raty1+bhylz *raty;
163 b[1] = blz *ratz1+bhz *ratz;
164 //
165 bhyhz = Bz(ix ,iy+1,iz+1)*ratx1+Bz(ix+1,iy+1,iz+1)*ratx;
166 bhylz = Bz(ix ,iy+1,iz )*ratx1+Bz(ix+1,iy+1,iz )*ratx;
167 blyhz = Bz(ix ,iy ,iz+1)*ratx1+Bz(ix+1,iy ,iz+1)*ratx;
168 blylz = Bz(ix ,iy ,iz )*ratx1+Bz(ix+1,iy ,iz )*ratx;
169 bhz = blyhz *raty1+bhyhz *raty;
170 blz = blylz *raty1+bhylz *raty;
171 b[2] = blz *ratz1+bhz *ratz;
172 //printf("ratx,raty,ratz,b[0],b[1],b[2] %f %f %f %f %f %f\n",
173 //ratx,raty,ratz,b[0],b[1],b[2]);
174 //
175 // ... use the dipole symmetry
176 if (x[0]*x[1] < 0) b[1]=-b[1];
177 if (x[0]<0) b[2]=-b[2];
178 } else {
179 printf("Invalid field map for constant mesh %d\n",fMap);
180 }
181 } else {
d95ea23f 182//This is the ZDC part
183 Float_t rad2=x[0]*x[0]+x[1]*x[1];
006cb30c 184 if(x[2]>kCORBEG2 && x[2]<kCOREND2){
d95ea23f 185 if(rad2<kCOR2RA2){
186 b[0] = kFCORN2;
187 }
188 }
189 else if(x[2]>kZ1BEG && x[2]<kZ1END){
190 if(rad2<kZ1RA2){
191 b[0] = -kG1*x[1];
192 b[1] = -kG1*x[0];
193 }
194 }
195 else if(x[2]>kZ2BEG && x[2]<kZ2END){
196 if(rad2<kZ2RA2){
197 b[0] = kG1*x[1];
198 b[1] = kG1*x[0];
199 }
200 }
201 else if(x[2]>kZ3BEG && x[2]<kZ3END){
202 if(rad2<kZ3RA2){
203 b[0] = kG1*x[1];
204 b[1] = kG1*x[0];
205 }
206 }
207 else if(x[2]>kZ4BEG && x[2]<kZ4END){
208 if(rad2<kZ4RA2){
209 b[0] = -kG1*x[1];
210 b[1] = -kG1*x[0];
211 }
212 }
213 else if(x[2]>kD1BEG && x[2]<kD1END){
214 if(rad2<kD1RA2){
215 b[1] = -kFDIP;
216 }
217 }
218 else if(x[2]>kD2BEG && x[2]<kD2END){
219 if(((x[0]-kXCEN1D2)*(x[0]-kXCEN1D2)+(x[1]-kYCEN1D2)*(x[1]-kYCEN1D2))<kD2RA2
220 || ((x[0]-kXCEN2D2)*(x[0]-kXCEN2D2)+(x[1]-kYCEN2D2)*(x[1]-kYCEN2D2))<kD2RA2){
221 b[1] = kFDIP;
aee8290b 222 }
223 }
d95ea23f 224
225 }
aee8290b 226 }
227 if(fFactor!=1) {
228 b[0]*=fFactor;
229 b[1]*=fFactor;
230 b[2]*=fFactor;
231 }
232}
233
e2afb3b6 234//_______________________________________________________________________
aee8290b 235void AliMagFCM::ReadField()
236{
237 //
238 // Method to read the magnetic field map from file
239 //
240 FILE *magfile;
241 Int_t ix, iy, iz, ipx, ipy, ipz;
242 Float_t bx, by, bz;
243 char *fname;
9e1a0ddb 244 if(fDebug) printf("%s: Reading Magnetic Field %s from file %s\n",ClassName(),fName.Data(),fTitle.Data());
aee8290b 245 fname = gSystem->ExpandPathName(fTitle.Data());
246 magfile=fopen(fname,"r");
247 delete [] fname;
248 if (magfile) {
249 fscanf(magfile,"%d %d %d %f %f %f %f %f %f",
250 &fXn, &fYn, &fZn, &fXdel, &fYdel, &fZdel, &fXbeg, &fYbeg, &fZbeg);
9e1a0ddb 251 if(fDebug>1) printf("%s: fXn %d, fYn %d, fZn %d, fXdel %f, fYdel %f, fZdel %f, fXbeg %f, fYbeg %f, fZbeg %f\n",
252 ClassName(),fXn, fYn, fZn, fXdel, fYdel, fZdel, fXbeg, fYbeg, fZbeg);
aee8290b 253 fXdeli=1./fXdel;
254 fYdeli=1./fYdel;
255 fZdeli=1./fZdel;
256 fB = new TVector(3*fXn*fYn*fZn);
257 for (iz=0; iz<fZn; iz++) {
258 ipz=iz*3*(fXn*fYn);
259 for (iy=0; iy<fYn; iy++) {
260 ipy=ipz+iy*3*fXn;
261 for (ix=0; ix<fXn; ix++) {
262 ipx=ipy+ix*3;
263 fscanf(magfile,"%f %f %f",&bz,&by,&bx);
264 (*fB)(ipx+2)=bz;
265 (*fB)(ipx+1)=by;
266 (*fB)(ipx )=bx;
267 }
268 }
269 }
270 } else {
9e1a0ddb 271 printf("%s: File %s not found !\n",ClassName(),fTitle.Data());
aee8290b 272 exit(1);
273 }
274}
275
e2afb3b6 276//_______________________________________________________________________
8918e700 277void AliMagFCM::Copy(AliMagFCM & /* magf */) const
aee8290b 278{
279 //
8918e700 280 // Copy *this onto magf -- Not implemented
aee8290b 281 //
282 Fatal("Copy","Not implemented!\n");
283}
284