Added to AliMagF the definition (const) of the polarities conventions.
[u/mrichter/AliRoot.git] / STEER / AliMagF.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
17 #include <TClass.h>
18 #include <TFile.h>
19 #include <TSystem.h>
20
21 #include "AliMagF.h"
22 #include "AliMagWrapCheb.h"
23 #include "AliLog.h"
24
25 ClassImp(AliMagF)
26
27 const Double_t AliMagF::fgkSol2DipZ    =  -700.;  
28 const UShort_t AliMagF::fgkPolarityConvention = kConvLHC;
29
30 /*
31  Explanation for polarity conventions: these are the mapping between the
32  current signs and main field components in L3 (Bz) and Dipole (Bx) (in Alice frame)
33  1) kConvMap2005: used for the field mapping in 2005
34  positive L3  current -> negative Bz
35  positive Dip current -> positive Bx 
36  2) kConvMapDCS2008: defined by the microswitches/cabling of power converters as of 2008 - 1st half 2009
37  positive L3  current -> positive Bz
38  positive Dip current -> positive Bx
39  3) kConvLHC : defined by LHC
40  positive L3  current -> negative Bz
41  positive Dip current -> negative Bx
42  
43  Note: only "negative Bz(L3) with postive Bx(Dipole)" and its inverse was mapped in 2005. Hence 
44  the GRP Manager will reject the runs with the current combinations (in the convention defined by the
45  static Int_t AliMagF::GetPolarityConvention()) which do not lead to such field polarities.
46 */
47 //_______________________________________________________________________
48 AliMagF::AliMagF():
49   TVirtualMagField(),
50   fMeasuredMap(0),
51   fMapType(k5kG),
52   fSolenoid(0),
53   fBeamType(kNoBeamField),
54   fBeamEnergy(0),
55   //
56   fInteg(0),
57   fPrecInteg(0),
58   fFactorSol(1.),
59   fFactorDip(1.),
60   fMax(15),
61   fDipoleOFF(kFALSE),
62   //
63   fQuadGradient(0),
64   fDipoleField(0),
65   fCCorrField(0), 
66   fACorr1Field(0),
67   fACorr2Field(0),
68   fParNames("","")
69 {
70   // Default constructor
71   //
72 }
73
74 //_______________________________________________________________________
75 AliMagF::AliMagF(const char *name, const char* title, Int_t integ, 
76                  Double_t factorSol, Double_t factorDip, 
77                  Double_t fmax, BMap_t maptype, const char* path,
78                  BeamType_t bt, Double_t be):
79   TVirtualMagField(name),
80   fMeasuredMap(0),
81   fMapType(maptype),
82   fSolenoid(0),
83   fBeamType(bt),
84   fBeamEnergy(be),
85   //
86   fInteg(integ),
87   fPrecInteg(1),
88   fFactorSol(1.),
89   fFactorDip(1.),
90   fMax(fmax),
91   fDipoleOFF(factorDip==0.),
92   //
93   fQuadGradient(0),
94   fDipoleField(0),
95   fCCorrField(0), 
96   fACorr1Field(0),
97   fACorr2Field(0),
98   fParNames("","")
99 {
100   // Initialize the field with Geant integration option "integ" and max field "fmax,
101   // Impose scaling of parameterized L3 field by factorSol and of dipole by factorDip.
102   // The "be" is the energy of the beam in GeV/nucleon
103   //
104   SetTitle(title);
105   if(integ<0 || integ > 2) {
106     AliWarning(Form("Invalid magnetic field flag: %5d; Helix tracking chosen instead",integ));
107     fInteg = 2;
108   }
109   if (fInteg == 0) fPrecInteg = 0;
110   //
111   const char* parname = 0;
112   //  
113   if      (fMapType == k2kG) parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
114   else if (fMapType == k5kG) parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
115   else if (fMapType == k5kGUniform) parname = "Sol30_Dip6_Uniform";
116   else AliFatal(Form("Unknown field identifier %d is requested\n",fMapType));
117   //
118   SetDataFileName(path);
119   SetParamName(parname);
120   //
121   LoadParameterization();
122   InitMachineField(fBeamType,fBeamEnergy);
123   double xyz[3]={0.,0.,0.};
124   fSolenoid = GetBz(xyz);
125   SetFactorSol(factorSol);
126   SetFactorDip(factorDip);
127 }
128
129 //_______________________________________________________________________
130 AliMagF::AliMagF(const AliMagF &src):
131   TVirtualMagField(src),
132   fMeasuredMap(0),
133   fMapType(src.fMapType),
134   fSolenoid(src.fSolenoid),
135   fBeamType(src.fBeamType),
136   fBeamEnergy(src.fBeamEnergy),
137   fInteg(src.fInteg),
138   fPrecInteg(src.fPrecInteg),
139   fFactorSol(src.fFactorSol),
140   fFactorDip(src.fFactorDip),
141   fMax(src.fMax),
142   fDipoleOFF(src.fDipoleOFF),
143   fQuadGradient(src.fQuadGradient),
144   fDipoleField(src.fDipoleField),
145   fCCorrField(src.fCCorrField), 
146   fACorr1Field(src.fACorr1Field),
147   fACorr2Field(src.fACorr2Field),
148   fParNames(src.fParNames)
149 {
150   if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
151 }
152
153 //_______________________________________________________________________
154 AliMagF::~AliMagF()
155 {
156   delete fMeasuredMap;
157 }
158
159 //_______________________________________________________________________
160 Bool_t AliMagF::LoadParameterization()
161 {
162   if (fMeasuredMap) {
163     AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
164     return kTRUE;
165   }
166   //
167   char* fname = gSystem->ExpandPathName(GetDataFileName());
168   TFile* file = TFile::Open(fname);
169   if (!file) {
170     AliError(Form("Failed to open magnetic field data file %s\n",fname)); 
171     return kFALSE;
172   }
173   //
174   fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
175   if (!fMeasuredMap) {
176     AliError(Form("Did not find field %s in %s\n",GetParamName(),fname)); 
177     return kFALSE;
178   }
179   file->Close();
180   delete file;
181   return kTRUE;
182 }
183
184
185 //_______________________________________________________________________
186 void AliMagF::Field(const Double_t *xyz, Double_t *b)
187 {
188   // Method to calculate the field at point  xyz
189   //
190   //  b[0]=b[1]=b[2]=0.0;
191   if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
192     fMeasuredMap->Field(xyz,b);
193     if (xyz[2]>fgkSol2DipZ || fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
194     else                                  for (int i=3;i--;) b[i] *= fFactorDip;    
195   }
196   else MachineField(xyz, b);
197   //
198 }
199
200 //_______________________________________________________________________
201 Double_t AliMagF::GetBz(const Double_t *xyz) const
202 {
203   // Method to calculate the field at point  xyz
204   //
205   if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
206     double bz = fMeasuredMap->GetBz(xyz);
207     return (xyz[2]>fgkSol2DipZ || fDipoleOFF) ? bz*fFactorSol : bz*fFactorDip;    
208   }
209   else return 0.;
210 }
211
212 //_______________________________________________________________________
213 AliMagF& AliMagF::operator=(const AliMagF& src)
214 {
215   if (this != &src && src.fMeasuredMap) { 
216     if (fMeasuredMap) delete fMeasuredMap;
217     fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
218     SetName(src.GetName());
219     fSolenoid    = src.fSolenoid;
220     fBeamType    = src.fBeamType;
221     fBeamEnergy  = src.fBeamEnergy;
222     fInteg       = src.fInteg;
223     fPrecInteg   = src.fPrecInteg;
224     fFactorSol   = src.fFactorSol;
225     fFactorDip   = src.fFactorDip;
226     fMax         = src.fMax;
227     fDipoleOFF   = src.fDipoleOFF;
228     fParNames    = src.fParNames;
229   }
230   return *this;
231 }
232
233 //_______________________________________________________________________
234 void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
235 {
236   if (btype==kNoBeamField || benergy<1.) {
237     fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
238     return;
239   }
240   //
241   double rigScale = benergy/7000.;   // scale according to ratio of E/Enominal
242   // for ions assume PbPb (with energy provided per nucleon) and account for A/Z
243   if (btype == kBeamTypeAA) rigScale *= 208./82.;
244   //
245   fQuadGradient = 22.0002*rigScale;
246   fDipoleField  = 37.8781*rigScale;
247   //
248   // SIDE C
249   fCCorrField   = -9.6980;
250   // SIDE A
251   fACorr1Field  = -13.2247;
252   fACorr2Field  =  11.7905;
253   //
254 }
255
256 //_______________________________________________________________________
257 void AliMagF::MachineField(const Double_t *x, Double_t *b) const
258 {
259   // ---- This is the ZDC part
260   // Compansators for Alice Muon Arm Dipole
261   const Double_t kBComp1CZ = 1075., kBComp1hDZ = 260./2., kBComp1SqR = 4.0*4.0; 
262   const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5; 
263   //  
264   const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
265   const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
266   const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5;
267   const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5;
268   //
269   const Double_t kDip1CZ = 6310.8,  kDip1hDZ = 945./2., kDip1SqRC = 4.5*4.5, kDip1SqRA = 3.375*3.375;
270   const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.5*4.5, kDip2SqRA = 3.75*3.75;
271   const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4;
272   //
273   double rad2 = x[0] * x[0] + x[1] * x[1];
274   //
275   b[0] = b[1] = b[2] = 0;
276   //
277   // SIDE C **************************************************
278   if(x[2]<0.){  
279     if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
280       b[0] = fCCorrField*fFactorDip;
281     } 
282     else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
283       b[0] = fQuadGradient*x[1];
284       b[1] = fQuadGradient*x[0];
285     }
286     else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
287       b[0] = -fQuadGradient*x[1];
288       b[1] = -fQuadGradient*x[0];
289     }
290     else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
291       b[0] = -fQuadGradient*x[1];
292       b[1] = -fQuadGradient*x[0];
293     }
294     else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
295       b[0] = fQuadGradient*x[1];
296       b[1] = fQuadGradient*x[0];
297     }
298     else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
299       b[1] = fDipoleField;
300     }
301     else if(TMath::Abs(x[2]+kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRC) {
302       double dxabs = TMath::Abs(x[0])-kDip2DXC;
303       if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRC) {
304         b[1] = -fDipoleField;
305       }
306     }
307   }
308   //
309   // SIDE A **************************************************
310   else{        
311     if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
312       // Compensator magnet at z = 1075 m 
313       b[0] = fACorr1Field*fFactorDip;
314     }
315     //
316     if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
317       b[0] = fACorr2Field*fFactorDip;
318     }
319     else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
320       b[0] = -fQuadGradient*x[1];
321       b[1] = -fQuadGradient*x[0];
322     }
323     else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
324       b[0] =  fQuadGradient*x[1];
325       b[1] =  fQuadGradient*x[0];
326     }
327     else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
328       b[0] =  fQuadGradient*x[1];
329       b[1] =  fQuadGradient*x[0];
330     }
331     else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
332       b[0] = -fQuadGradient*x[1];
333       b[1] = -fQuadGradient*x[0];
334     }
335     else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
336       b[1] = -fDipoleField;
337     }
338     else if(TMath::Abs(x[2]-kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRA) {
339       double dxabs = TMath::Abs(x[0])-kDip2DXA;
340       if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRA) {
341         b[1] = fDipoleField;
342       }
343     }
344   }
345   //
346 }
347
348 //_______________________________________________________________________
349 void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
350 {
351   // Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
352   b[0]=b[1]=b[2]=0.0;
353   if (fMeasuredMap) {
354     fMeasuredMap->GetTPCInt(xyz,b);
355     for (int i=3;i--;) b[i] *= fFactorSol;
356   }
357 }
358
359 //_______________________________________________________________________
360 void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
361 {
362   // Method to calculate the integral of magnetic integral from point to nearest cathode plane
363   // in cylindrical coordiates ( -pi<phi<pi convention )
364   b[0]=b[1]=b[2]=0.0;
365   if (fMeasuredMap) {
366     fMeasuredMap->GetTPCIntCyl(rphiz,b);
367     for (int i=3;i--;) b[i] *= fFactorSol;
368   }
369 }
370
371 //_______________________________________________________________________
372 void AliMagF::SetFactorSol(Float_t fc)
373 {
374   // set the sign/scale of the current in the L3 according to fgkPolarityConvention
375   switch (fgkPolarityConvention) {
376   case kConvDCS2008: fFactorSol = -fc; break;
377   case kConvLHC    : fFactorSol = -fc; break;
378   default          : fFactorSol =  fc; break;  // case kConvMap2005: fFactorSol =  fc; break;
379   }
380 }
381
382 //_______________________________________________________________________
383 void AliMagF::SetFactorDip(Float_t fc)
384 {
385   // set the sign*scale of the current in the Dipole according to fgkPolarityConvention
386   switch (fgkPolarityConvention) {
387   case kConvDCS2008: fFactorDip =  fc; break;
388   case kConvLHC    : fFactorDip = -fc; break;
389   default          : fFactorDip =  fc; break;  // case kConvMap2005: fFactorDip =  fc; break;
390   }
391 }
392
393 //_______________________________________________________________________
394 Double_t AliMagF::GetFactorSol() const
395 {
396   // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe 
397   switch (fgkPolarityConvention) {
398   case kConvDCS2008: return -fFactorSol;
399   case kConvLHC    : return -fFactorSol;
400   default          : return  fFactorSol;       //  case kConvMap2005: return  fFactorSol;
401   }
402 }
403
404 //_______________________________________________________________________
405 Double_t AliMagF::GetFactorDip() const
406 {
407   // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe 
408   switch (fgkPolarityConvention) {
409   case kConvDCS2008: return  fFactorDip;
410   case kConvLHC    : return -fFactorDip;
411   default          : return  fFactorDip;       //  case kConvMap2005: return  fFactorDip;
412   }
413 }