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