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