changes in the MagF constructor
[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 #include <TPRegexp.h>
21
22 #include "AliMagF.h"
23 #include "AliMagWrapCheb.h"
24 #include "AliLog.h"
25
26 ClassImp(AliMagF)
27
28 const Double_t AliMagF::fgkSol2DipZ    =  -700.;  
29 const UShort_t AliMagF::fgkPolarityConvention = kConvLHC;
30
31 /*
32  Explanation for polarity conventions: these are the mapping between the
33  current signs and main field components in L3 (Bz) and Dipole (Bx) (in Alice frame)
34  1) kConvMap2005: used for the field mapping in 2005
35  positive L3  current -> negative Bz
36  positive Dip current -> positive Bx 
37  2) kConvMapDCS2008: defined by the microswitches/cabling of power converters as of 2008 - 1st half 2009
38  positive L3  current -> positive Bz
39  positive Dip current -> positive Bx
40  3) kConvLHC : defined by LHC
41  positive L3  current -> negative Bz
42  positive Dip current -> negative Bx
43  
44  Note: only "negative Bz(L3) with postive Bx(Dipole)" and its inverse was mapped in 2005. Hence 
45  the GRP Manager will reject the runs with the current combinations (in the convention defined by the
46  static Int_t AliMagF::GetPolarityConvention()) which do not lead to such field polarities.
47 */
48 //_______________________________________________________________________
49 AliMagF::AliMagF():
50   TVirtualMagField(),
51   fMeasuredMap(0),
52   fMapType(k5kG),
53   fSolenoid(0),
54   fBeamType(kNoBeamField),
55   fBeamEnergy(0),
56   //
57   fInteg(0),
58   fPrecInteg(0),
59   fFactorSol(1.),
60   fFactorDip(1.),
61   fMax(15),
62   fDipoleOFF(kFALSE),
63   //
64   fQuadGradient(0),
65   fDipoleField(0),
66   fCCorrField(0), 
67   fACorr1Field(0),
68   fACorr2Field(0),
69   fParNames("","")
70 {
71   // Default constructor
72   //
73 }
74
75 //_______________________________________________________________________
76 AliMagF::AliMagF(const char *name, const char* title, Double_t factorSol, Double_t factorDip, 
77                  BMap_t maptype, BeamType_t btype, Double_t benergy, Int_t integ, Double_t fmax, 
78                  const char* path):
79   TVirtualMagField(name),
80   fMeasuredMap(0),
81   fMapType(maptype),
82   fSolenoid(0),
83   fBeamType(btype),
84   fBeamEnergy(benergy),
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   if (fBeamEnergy<=0 && fBeamType!=kNoBeamField) {
112     if      (fBeamType == kBeamTypepp) fBeamEnergy = 7000.; // max proton energy
113     else if (fBeamType == kBeamTypeAA) fBeamEnergy = 5500;  // max PbPb energy
114     AliInfo("Maximim possible beam energy for requested beam is assumed");
115   } 
116   const char* parname = 0;
117   //  
118   if      (fMapType == k2kG) parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
119   else if (fMapType == k5kG) parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
120   else if (fMapType == k5kGUniform) parname = "Sol30_Dip6_Uniform";
121   else AliFatal(Form("Unknown field identifier %d is requested\n",fMapType));
122   //
123   SetDataFileName(path);
124   SetParamName(parname);
125   //
126   LoadParameterization();
127   InitMachineField(fBeamType,fBeamEnergy);
128   double xyz[3]={0.,0.,0.};
129   fSolenoid = GetBz(xyz);
130   SetFactorSol(factorSol);
131   SetFactorDip(factorDip);
132   AliInfo(Form("Alice   B fields: Solenoid (%+.2f*)%.0f kG, Dipole %s (%+.2f) %s",
133                factorSol,(fMapType==k5kG||fMapType==k5kGUniform)?5.:2.,
134                fDipoleOFF ? "OFF":"ON",factorDip,fMapType==k5kGUniform?" |Constant Field!":""));
135   AliInfo(Form("Machine B fields for %s beam (%.0f GeV): QGrad: %.4f Dipole: %.4f",
136                fBeamType==kBeamTypeAA ? "A-A":(fBeamType==kBeamTypepp ? "p-p":"OFF"),
137                fBeamEnergy,fQuadGradient,fDipoleField));
138 }
139
140 //_______________________________________________________________________
141 AliMagF::AliMagF(const AliMagF &src):
142   TVirtualMagField(src),
143   fMeasuredMap(0),
144   fMapType(src.fMapType),
145   fSolenoid(src.fSolenoid),
146   fBeamType(src.fBeamType),
147   fBeamEnergy(src.fBeamEnergy),
148   fInteg(src.fInteg),
149   fPrecInteg(src.fPrecInteg),
150   fFactorSol(src.fFactorSol),
151   fFactorDip(src.fFactorDip),
152   fMax(src.fMax),
153   fDipoleOFF(src.fDipoleOFF),
154   fQuadGradient(src.fQuadGradient),
155   fDipoleField(src.fDipoleField),
156   fCCorrField(src.fCCorrField), 
157   fACorr1Field(src.fACorr1Field),
158   fACorr2Field(src.fACorr2Field),
159   fParNames(src.fParNames)
160 {
161   if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
162 }
163
164 //_______________________________________________________________________
165 AliMagF::~AliMagF()
166 {
167   delete fMeasuredMap;
168 }
169
170 //_______________________________________________________________________
171 Bool_t AliMagF::LoadParameterization()
172 {
173   if (fMeasuredMap) {
174     AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
175     return kTRUE;
176   }
177   //
178   char* fname = gSystem->ExpandPathName(GetDataFileName());
179   TFile* file = TFile::Open(fname);
180   if (!file) {
181     AliError(Form("Failed to open magnetic field data file %s\n",fname)); 
182     return kFALSE;
183   }
184   //
185   fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
186   if (!fMeasuredMap) {
187     AliError(Form("Did not find field %s in %s\n",GetParamName(),fname)); 
188     return kFALSE;
189   }
190   file->Close();
191   delete file;
192   return kTRUE;
193 }
194
195
196 //_______________________________________________________________________
197 void AliMagF::Field(const Double_t *xyz, Double_t *b)
198 {
199   // Method to calculate the field at point  xyz
200   //
201   //  b[0]=b[1]=b[2]=0.0;
202   if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
203     fMeasuredMap->Field(xyz,b);
204     if (xyz[2]>fgkSol2DipZ || fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
205     else                                  for (int i=3;i--;) b[i] *= fFactorDip;    
206   }
207   else MachineField(xyz, b);
208   //
209 }
210
211 //_______________________________________________________________________
212 Double_t AliMagF::GetBz(const Double_t *xyz) const
213 {
214   // Method to calculate the field at point  xyz
215   //
216   if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
217     double bz = fMeasuredMap->GetBz(xyz);
218     return (xyz[2]>fgkSol2DipZ || fDipoleOFF) ? bz*fFactorSol : bz*fFactorDip;    
219   }
220   else return 0.;
221 }
222
223 //_______________________________________________________________________
224 AliMagF& AliMagF::operator=(const AliMagF& src)
225 {
226   if (this != &src && src.fMeasuredMap) { 
227     if (fMeasuredMap) delete fMeasuredMap;
228     fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
229     SetName(src.GetName());
230     fSolenoid    = src.fSolenoid;
231     fBeamType    = src.fBeamType;
232     fBeamEnergy  = src.fBeamEnergy;
233     fInteg       = src.fInteg;
234     fPrecInteg   = src.fPrecInteg;
235     fFactorSol   = src.fFactorSol;
236     fFactorDip   = src.fFactorDip;
237     fMax         = src.fMax;
238     fDipoleOFF   = src.fDipoleOFF;
239     fParNames    = src.fParNames;
240   }
241   return *this;
242 }
243
244 //_______________________________________________________________________
245 void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
246 {
247   if (btype==kNoBeamField) {
248     fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
249     return;
250   }
251   //
252   double rigScale = benergy/7000.;   // scale according to ratio of E/Enominal
253   // for ions assume PbPb (with energy provided per nucleon) and account for A/Z
254   if (btype == kBeamTypeAA) rigScale *= 208./82.;
255   //
256   fQuadGradient = 22.0002*rigScale;
257   fDipoleField  = 37.8781*rigScale;
258   //
259   // SIDE C
260   fCCorrField   = -9.6980;
261   // SIDE A
262   fACorr1Field  = -13.2247;
263   fACorr2Field  =  11.7905;
264   //
265 }
266
267 //_______________________________________________________________________
268 void AliMagF::MachineField(const Double_t *x, Double_t *b) const
269 {
270   // ---- This is the ZDC part
271   // Compansators for Alice Muon Arm Dipole
272   const Double_t kBComp1CZ = 1075., kBComp1hDZ = 260./2., kBComp1SqR = 4.0*4.0; 
273   const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5; 
274   //  
275   const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
276   const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
277   const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5;
278   const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5;
279   //
280   const Double_t kDip1CZ = 6310.8,  kDip1hDZ = 945./2., kDip1SqRC = 4.5*4.5, kDip1SqRA = 3.375*3.375;
281   const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.5*4.5, kDip2SqRA = 3.75*3.75;
282   const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4;
283   //
284   double rad2 = x[0] * x[0] + x[1] * x[1];
285   //
286   b[0] = b[1] = b[2] = 0;
287   //
288   // SIDE C **************************************************
289   if(x[2]<0.){  
290     if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
291       b[0] = fCCorrField*fFactorDip;
292     } 
293     else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
294       b[0] = fQuadGradient*x[1];
295       b[1] = fQuadGradient*x[0];
296     }
297     else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
298       b[0] = -fQuadGradient*x[1];
299       b[1] = -fQuadGradient*x[0];
300     }
301     else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
302       b[0] = -fQuadGradient*x[1];
303       b[1] = -fQuadGradient*x[0];
304     }
305     else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
306       b[0] = fQuadGradient*x[1];
307       b[1] = fQuadGradient*x[0];
308     }
309     else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
310       b[1] = fDipoleField;
311     }
312     else if(TMath::Abs(x[2]+kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRC) {
313       double dxabs = TMath::Abs(x[0])-kDip2DXC;
314       if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRC) {
315         b[1] = -fDipoleField;
316       }
317     }
318   }
319   //
320   // SIDE A **************************************************
321   else{        
322     if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
323       // Compensator magnet at z = 1075 m 
324       b[0] = fACorr1Field*fFactorDip;
325     }
326     //
327     if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
328       b[0] = fACorr2Field*fFactorDip;
329     }
330     else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
331       b[0] = -fQuadGradient*x[1];
332       b[1] = -fQuadGradient*x[0];
333     }
334     else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
335       b[0] =  fQuadGradient*x[1];
336       b[1] =  fQuadGradient*x[0];
337     }
338     else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
339       b[0] =  fQuadGradient*x[1];
340       b[1] =  fQuadGradient*x[0];
341     }
342     else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
343       b[0] = -fQuadGradient*x[1];
344       b[1] = -fQuadGradient*x[0];
345     }
346     else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
347       b[1] = -fDipoleField;
348     }
349     else if(TMath::Abs(x[2]-kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRA) {
350       double dxabs = TMath::Abs(x[0])-kDip2DXA;
351       if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRA) {
352         b[1] = fDipoleField;
353       }
354     }
355   }
356   //
357 }
358
359 //_______________________________________________________________________
360 void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
361 {
362   // Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
363   b[0]=b[1]=b[2]=0.0;
364   if (fMeasuredMap) {
365     fMeasuredMap->GetTPCInt(xyz,b);
366     for (int i=3;i--;) b[i] *= fFactorSol;
367   }
368 }
369
370 //_______________________________________________________________________
371 void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
372 {
373   // Method to calculate the integral of magnetic integral from point to nearest cathode plane
374   // in cylindrical coordiates ( -pi<phi<pi convention )
375   b[0]=b[1]=b[2]=0.0;
376   if (fMeasuredMap) {
377     fMeasuredMap->GetTPCIntCyl(rphiz,b);
378     for (int i=3;i--;) b[i] *= fFactorSol;
379   }
380 }
381
382 //_______________________________________________________________________
383 void AliMagF::SetFactorSol(Float_t fc)
384 {
385   // set the sign/scale of the current in the L3 according to fgkPolarityConvention
386   switch (fgkPolarityConvention) {
387   case kConvDCS2008: fFactorSol = -fc; break;
388   case kConvLHC    : fFactorSol = -fc; break;
389   default          : fFactorSol =  fc; break;  // case kConvMap2005: fFactorSol =  fc; break;
390   }
391 }
392
393 //_______________________________________________________________________
394 void AliMagF::SetFactorDip(Float_t fc)
395 {
396   // set the sign*scale of the current in the Dipole according to fgkPolarityConvention
397   switch (fgkPolarityConvention) {
398   case kConvDCS2008: fFactorDip =  fc; break;
399   case kConvLHC    : fFactorDip = -fc; break;
400   default          : fFactorDip =  fc; break;  // case kConvMap2005: fFactorDip =  fc; break;
401   }
402 }
403
404 //_______________________________________________________________________
405 Double_t AliMagF::GetFactorSol() const
406 {
407   // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe 
408   switch (fgkPolarityConvention) {
409   case kConvDCS2008: return -fFactorSol;
410   case kConvLHC    : return -fFactorSol;
411   default          : return  fFactorSol;       //  case kConvMap2005: return  fFactorSol;
412   }
413 }
414
415 //_______________________________________________________________________
416 Double_t AliMagF::GetFactorDip() const
417 {
418   // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe 
419   switch (fgkPolarityConvention) {
420   case kConvDCS2008: return  fFactorDip;
421   case kConvLHC    : return -fFactorDip;
422   default          : return  fFactorDip;       //  case kConvMap2005: return  fFactorDip;
423   }
424 }
425
426 //_____________________________________________________________________________
427 AliMagF* AliMagF::CreateFieldMap(Float_t l3Cur, Float_t diCur, Int_t convention, Bool_t uniform,
428                                  Float_t beamenergy, const Char_t *beamtype, const Char_t *path) 
429 {
430   //------------------------------------------------
431   // The magnetic field map, defined externally...
432   // L3 current 30000 A  -> 0.5 T
433   // L3 current 12000 A  -> 0.2 T
434   // dipole current 6000 A
435   // The polarities must match the convention (LHC or DCS2008) 
436   // unless the special uniform map was used for MC
437   //------------------------------------------------
438   const Float_t l3NominalCurrent1=30000.; // (A)
439   const Float_t l3NominalCurrent2=12000.; // (A)
440   const Float_t diNominalCurrent =6000. ; // (A)
441
442   const Float_t tolerance=0.03; // relative current tolerance
443   const Float_t zero=77.;       // "zero" current (A)
444   //
445   BMap_t map;
446   double sclL3,sclDip;
447   //
448   Float_t l3Pol = l3Cur > 0 ? 1:-1;
449   Float_t diPol = diCur > 0 ? 1:-1;
450  
451   l3Cur = TMath::Abs(l3Cur);
452   diCur = TMath::Abs(diCur);
453   //
454   if (TMath::Abs((sclDip=diCur/diNominalCurrent)-1.) > tolerance && !uniform) {
455     if (diCur <= zero) sclDip = 0.; // some small current.. -> Dipole OFF
456     else {
457       AliErrorGeneral("AliMagF",Form("Wrong dipole current (%f A)!",diCur));
458       return 0;
459     }
460   }
461   //
462   if (uniform) { 
463     // special treatment of special MC with uniform mag field (normalized to 0.5 T)
464     // no check for scaling/polarities are done
465     map   = k5kGUniform;
466     sclL3 = l3Cur/l3NominalCurrent1; 
467   }
468   else {
469     if      (TMath::Abs((sclL3=l3Cur/l3NominalCurrent1)-1.) < tolerance) map  = k5kG;
470     else if (TMath::Abs((sclL3=l3Cur/l3NominalCurrent2)-1.) < tolerance) map  = k2kG;
471     else if (l3Cur <= zero)                                { sclL3 = 0;  map  = k5kGUniform;}
472     else {
473       AliErrorGeneral("AliMagF",Form("Wrong L3 current (%f A)!",l3Cur));
474       return 0;
475     }
476   }
477   //
478   if (sclDip!=0 && (map==k5kG || map==k2kG) &&
479       ((convention==kConvLHC     && l3Pol!=diPol) ||
480        (convention==kConvDCS2008 && l3Pol==diPol)) ) { 
481     AliErrorGeneral("AliMagF",Form("Wrong combination for L3/Dipole polarities (%c/%c) for convention %d",
482                                    l3Pol>0?'+':'-',diPol>0?'+':'-',GetPolarityConvention()));
483     return 0;
484   }
485   //
486   if (l3Pol<0) sclL3  = -sclL3;
487   if (diPol<0) sclDip = -sclDip;
488   //
489   BeamType_t btype = kNoBeamField;
490   TString btypestr = beamtype;
491   btypestr.ToLower();
492   TPRegexp protonBeam("(proton|p)\\s*-?\\s*\\1");
493   TPRegexp ionBeam("(lead|pb|ion|a)\\s*-?\\s*\\1");
494   if (btypestr.Contains(ionBeam)) btype = kBeamTypeAA;
495   else if (btypestr.Contains(protonBeam)) btype = kBeamTypepp;
496   else AliInfoGeneral("AliMagF",Form("Assume no LHC magnet field for the beam type %s, ",beamtype));
497   char ttl[80];
498   sprintf(ttl,"L3: %+5d Dip: %+4d kA; %s | Polarities in %s convention",(int)TMath::Sign(l3Cur,float(sclL3)),
499           (int)TMath::Sign(diCur,float(sclDip)),uniform ? " Constant":"",
500           convention==kConvLHC ? "LHC":"DCS2008");
501   // LHC and DCS08 conventions have opposite dipole polarities
502   if ( GetPolarityConvention() != convention) sclDip = -sclDip;
503   //
504   return new AliMagF("MagneticFieldMap", ttl,sclL3,sclDip,map,btype,beamenergy,2,10.,path);
505   //
506 }
507
508 //_____________________________________________________________________________
509 const char*  AliMagF::GetBeamTypeText() const
510 {
511   const char *beamNA  = "No Beam";
512   const char *beamPP  = "p-p";
513   const char *beamPbPb= "Pb-Pb";
514   switch ( fBeamType ) {
515   case kBeamTypepp : return beamPP;
516   case kBeamTypeAA : return beamPbPb;
517   case kNoBeamField: 
518   default:           return beamNA;
519   }
520 }
521