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