Corrected protection.
[u/mrichter/AliRoot.git] / STEER / AliMagF.cxx
CommitLineData
4c039060 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
4c039060 16
db83d72f 17#include <TClass.h>
18#include <TFile.h>
19#include <TSystem.h>
33fe5eb1 20#include <TPRegexp.h>
fe4da5cc 21
22#include "AliMagF.h"
db83d72f 23#include "AliMagWrapCheb.h"
24#include "AliLog.h"
972ca52f 25
fe4da5cc 26ClassImp(AliMagF)
27
9251fceb 28const Double_t AliMagF::fgkSol2DipZ = -700.;
47c3d315 29const UShort_t AliMagF::fgkPolarityConvention = AliMagF::kConvLHC;
1dd3d90e 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
e3eadfac 40 positive L3 current -> positive Bz
1dd3d90e 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.
704fefba 46
47 -----------------------------------------------
48
49 Explanation on integrals in the TPC region
50 GetTPCInt(xyz,b) and GetTPCRatInt(xyz,b) give integrals from point (x,y,z) to point (x,y,0)
51 (irrespectively of the z sign) of the following:
52 TPCInt: b contains int{bx}, int{by}, int{bz}
53 TPCRatInt: b contains int{bx/bz}, int{by/bz}, int{(bx/bz)^2+(by/bz)^2}
54
55 The same applies to integral in cylindrical coordinates:
56 GetTPCIntCyl(rphiz,b)
57 GetTPCIntRatCyl(rphiz,b)
58 They accept the R,Phi,Z coordinate (-pi<phi<pi) and return the field
59 integrals in cyl. coordinates.
60
61 Thus, to compute the integral from arbitrary xy_z1 to xy_z2, one should take
62 b = b1-b2 with b1 and b2 coming from GetTPCInt(xy_z1,b1) and GetTPCInt(xy_z2,b2)
63
64 Note: the integrals are defined for the range -300<Z<300 and 0<R<300
1dd3d90e 65*/
e2afb3b6 66//_______________________________________________________________________
67AliMagF::AliMagF():
db83d72f 68 TVirtualMagField(),
69 fMeasuredMap(0),
70 fMapType(k5kG),
71 fSolenoid(0),
72 fBeamType(kNoBeamField),
73 fBeamEnergy(0),
db83d72f 74 //
e2afb3b6 75 fInteg(0),
db83d72f 76 fPrecInteg(0),
77 fFactorSol(1.),
78 fFactorDip(1.),
79 fMax(15),
80 fDipoleOFF(kFALSE),
e2afb3b6 81 //
db83d72f 82 fQuadGradient(0),
83 fDipoleField(0),
84 fCCorrField(0),
85 fACorr1Field(0),
86 fACorr2Field(0),
87 fParNames("","")
88{
e2afb3b6 89 // Default constructor
90 //
91}
92
93//_______________________________________________________________________
4642ac4b 94AliMagF::AliMagF(const char *name, const char* title, Double_t factorSol, Double_t factorDip,
02233f2b 95 BMap_t maptype, BeamType_t bt, Double_t be,Int_t integ, Double_t fmax, const char* path):
db83d72f 96 TVirtualMagField(name),
97 fMeasuredMap(0),
98 fMapType(maptype),
99 fSolenoid(0),
02233f2b 100 fBeamType(bt),
101 fBeamEnergy(be),
db83d72f 102 //
103 fInteg(integ),
604e0531 104 fPrecInteg(1),
db83d72f 105 fFactorSol(1.),
106 fFactorDip(1.),
972ca52f 107 fMax(fmax),
db83d72f 108 fDipoleOFF(factorDip==0.),
109 //
110 fQuadGradient(0),
111 fDipoleField(0),
112 fCCorrField(0),
113 fACorr1Field(0),
114 fACorr2Field(0),
115 fParNames("","")
fe4da5cc 116{
9251fceb 117 // Initialize the field with Geant integration option "integ" and max field "fmax,
118 // Impose scaling of parameterized L3 field by factorSol and of dipole by factorDip.
119 // The "be" is the energy of the beam in GeV/nucleon
aee8290b 120 //
db83d72f 121 SetTitle(title);
122 if(integ<0 || integ > 2) {
123 AliWarning(Form("Invalid magnetic field flag: %5d; Helix tracking chosen instead",integ));
124 fInteg = 2;
125 }
126 if (fInteg == 0) fPrecInteg = 0;
aee8290b 127 //
4642ac4b 128 if (fBeamEnergy<=0 && fBeamType!=kNoBeamField) {
129 if (fBeamType == kBeamTypepp) fBeamEnergy = 7000.; // max proton energy
02233f2b 130 else if (fBeamType == kBeamTypeAA) fBeamEnergy = 2750; // max PbPb energy
4642ac4b 131 AliInfo("Maximim possible beam energy for requested beam is assumed");
132 }
db83d72f 133 const char* parname = 0;
134 //
f04e7f5f 135 if (fMapType == k2kG) parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
136 else if (fMapType == k5kG) parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
137 else if (fMapType == k5kGUniform) parname = "Sol30_Dip6_Uniform";
138 else AliFatal(Form("Unknown field identifier %d is requested\n",fMapType));
db83d72f 139 //
140 SetDataFileName(path);
141 SetParamName(parname);
142 //
db83d72f 143 LoadParameterization();
144 InitMachineField(fBeamType,fBeamEnergy);
f04e7f5f 145 double xyz[3]={0.,0.,0.};
146 fSolenoid = GetBz(xyz);
147 SetFactorSol(factorSol);
148 SetFactorDip(factorDip);
77c9a262 149 Print("a");
fe4da5cc 150}
151
e2afb3b6 152//_______________________________________________________________________
eeda4611 153AliMagF::AliMagF(const AliMagF &src):
db83d72f 154 TVirtualMagField(src),
155 fMeasuredMap(0),
156 fMapType(src.fMapType),
157 fSolenoid(src.fSolenoid),
158 fBeamType(src.fBeamType),
159 fBeamEnergy(src.fBeamEnergy),
eeda4611 160 fInteg(src.fInteg),
161 fPrecInteg(src.fPrecInteg),
db83d72f 162 fFactorSol(src.fFactorSol),
163 fFactorDip(src.fFactorDip),
eeda4611 164 fMax(src.fMax),
db83d72f 165 fDipoleOFF(src.fDipoleOFF),
166 fQuadGradient(src.fQuadGradient),
167 fDipoleField(src.fDipoleField),
168 fCCorrField(src.fCCorrField),
169 fACorr1Field(src.fACorr1Field),
170 fACorr2Field(src.fACorr2Field),
171 fParNames(src.fParNames)
eeda4611 172{
db83d72f 173 if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
eeda4611 174}
175
176//_______________________________________________________________________
db83d72f 177AliMagF::~AliMagF()
ff66b122 178{
db83d72f 179 delete fMeasuredMap;
180}
181
182//_______________________________________________________________________
183Bool_t AliMagF::LoadParameterization()
184{
185 if (fMeasuredMap) {
706eaf0b 186 AliFatal(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
db83d72f 187 }
ff66b122 188 //
db83d72f 189 char* fname = gSystem->ExpandPathName(GetDataFileName());
190 TFile* file = TFile::Open(fname);
191 if (!file) {
af19d8b9 192 AliFatal(Form("Failed to open magnetic field data file %s\n",fname));
db83d72f 193 }
ff66b122 194 //
db83d72f 195 fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
196 if (!fMeasuredMap) {
af19d8b9 197 AliFatal(Form("Did not find field %s in %s\n",GetParamName(),fname));
db83d72f 198 }
199 file->Close();
200 delete file;
201 return kTRUE;
ff66b122 202}
203
db83d72f 204
ff66b122 205//_______________________________________________________________________
db83d72f 206void AliMagF::Field(const Double_t *xyz, Double_t *b)
fe4da5cc 207{
db83d72f 208 // Method to calculate the field at point xyz
aee8290b 209 //
9251fceb 210 // b[0]=b[1]=b[2]=0.0;
211 if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
db83d72f 212 fMeasuredMap->Field(xyz,b);
213 if (xyz[2]>fgkSol2DipZ || fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
9251fceb 214 else for (int i=3;i--;) b[i] *= fFactorDip;
db83d72f 215 }
9251fceb 216 else MachineField(xyz, b);
aee8290b 217 //
fe4da5cc 218}
eeda4611 219
220//_______________________________________________________________________
db83d72f 221Double_t AliMagF::GetBz(const Double_t *xyz) const
eeda4611 222{
db83d72f 223 // Method to calculate the field at point xyz
224 //
9251fceb 225 if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
226 double bz = fMeasuredMap->GetBz(xyz);
227 return (xyz[2]>fgkSol2DipZ || fDipoleOFF) ? bz*fFactorSol : bz*fFactorDip;
db83d72f 228 }
9251fceb 229 else return 0.;
eeda4611 230}
231
232//_______________________________________________________________________
db83d72f 233AliMagF& AliMagF::operator=(const AliMagF& src)
eeda4611 234{
db83d72f 235 if (this != &src && src.fMeasuredMap) {
236 if (fMeasuredMap) delete fMeasuredMap;
237 fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
238 SetName(src.GetName());
239 fSolenoid = src.fSolenoid;
240 fBeamType = src.fBeamType;
241 fBeamEnergy = src.fBeamEnergy;
db83d72f 242 fInteg = src.fInteg;
243 fPrecInteg = src.fPrecInteg;
244 fFactorSol = src.fFactorSol;
245 fFactorDip = src.fFactorDip;
246 fMax = src.fMax;
247 fDipoleOFF = src.fDipoleOFF;
248 fParNames = src.fParNames;
249 }
250 return *this;
eeda4611 251}
252
253//_______________________________________________________________________
db83d72f 254void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
eeda4611 255{
4642ac4b 256 if (btype==kNoBeamField) {
db83d72f 257 fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
9251fceb 258 return;
db83d72f 259 }
260 //
9251fceb 261 double rigScale = benergy/7000.; // scale according to ratio of E/Enominal
262 // for ions assume PbPb (with energy provided per nucleon) and account for A/Z
263 if (btype == kBeamTypeAA) rigScale *= 208./82.;
264 //
265 fQuadGradient = 22.0002*rigScale;
266 fDipoleField = 37.8781*rigScale;
267 //
268 // SIDE C
269 fCCorrField = -9.6980;
270 // SIDE A
271 fACorr1Field = -13.2247;
272 fACorr2Field = 11.7905;
db83d72f 273 //
eeda4611 274}
eed8a1a2 275
db83d72f 276//_______________________________________________________________________
277void AliMagF::MachineField(const Double_t *x, Double_t *b) const
eed8a1a2 278{
db83d72f 279 // ---- This is the ZDC part
9251fceb 280 // Compansators for Alice Muon Arm Dipole
281 const Double_t kBComp1CZ = 1075., kBComp1hDZ = 260./2., kBComp1SqR = 4.0*4.0;
282 const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5;
283 //
284 const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
90ae20c9 285 const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
9251fceb 286 const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5;
287 const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5;
db83d72f 288 //
9251fceb 289 const Double_t kDip1CZ = 6310.8, kDip1hDZ = 945./2., kDip1SqRC = 4.5*4.5, kDip1SqRA = 3.375*3.375;
290 const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.5*4.5, kDip2SqRA = 3.75*3.75;
291 const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4;
db83d72f 292 //
293 double rad2 = x[0] * x[0] + x[1] * x[1];
294 //
9251fceb 295 b[0] = b[1] = b[2] = 0;
296 //
db83d72f 297 // SIDE C **************************************************
298 if(x[2]<0.){
9251fceb 299 if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
300 b[0] = fCCorrField*fFactorDip;
db83d72f 301 }
9251fceb 302 else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
db83d72f 303 b[0] = fQuadGradient*x[1];
304 b[1] = fQuadGradient*x[0];
db83d72f 305 }
9251fceb 306 else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
db83d72f 307 b[0] = -fQuadGradient*x[1];
308 b[1] = -fQuadGradient*x[0];
db83d72f 309 }
9251fceb 310 else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
db83d72f 311 b[0] = -fQuadGradient*x[1];
312 b[1] = -fQuadGradient*x[0];
db83d72f 313 }
9251fceb 314 else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
db83d72f 315 b[0] = fQuadGradient*x[1];
316 b[1] = fQuadGradient*x[0];
db83d72f 317 }
9251fceb 318 else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
db83d72f 319 b[1] = fDipoleField;
db83d72f 320 }
9251fceb 321 else if(TMath::Abs(x[2]+kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRC) {
322 double dxabs = TMath::Abs(x[0])-kDip2DXC;
323 if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRC) {
db83d72f 324 b[1] = -fDipoleField;
db83d72f 325 }
326 }
327 }
328 //
329 // SIDE A **************************************************
330 else{
9251fceb 331 if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
db83d72f 332 // Compensator magnet at z = 1075 m
9251fceb 333 b[0] = fACorr1Field*fFactorDip;
db83d72f 334 }
335 //
9251fceb 336 if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
337 b[0] = fACorr2Field*fFactorDip;
338 }
339 else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
db83d72f 340 b[0] = -fQuadGradient*x[1];
341 b[1] = -fQuadGradient*x[0];
eed8a1a2 342 }
9251fceb 343 else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
344 b[0] = fQuadGradient*x[1];
345 b[1] = fQuadGradient*x[0];
eed8a1a2 346 }
9251fceb 347 else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
348 b[0] = fQuadGradient*x[1];
349 b[1] = fQuadGradient*x[0];
db83d72f 350 }
9251fceb 351 else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
db83d72f 352 b[0] = -fQuadGradient*x[1];
353 b[1] = -fQuadGradient*x[0];
db83d72f 354 }
9251fceb 355 else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
db83d72f 356 b[1] = -fDipoleField;
db83d72f 357 }
9251fceb 358 else if(TMath::Abs(x[2]-kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRA) {
359 double dxabs = TMath::Abs(x[0])-kDip2DXA;
360 if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRA) {
db83d72f 361 b[1] = fDipoleField;
362 }
363 }
364 }
9251fceb 365 //
db83d72f 366}
367
368//_______________________________________________________________________
369void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
370{
47c3d315 371 // Method to calculate the integral_0^z of br,bt,bz
db83d72f 372 b[0]=b[1]=b[2]=0.0;
373 if (fMeasuredMap) {
374 fMeasuredMap->GetTPCInt(xyz,b);
375 for (int i=3;i--;) b[i] *= fFactorSol;
376 }
377}
378
379//_______________________________________________________________________
47c3d315 380void AliMagF::GetTPCRatInt(const Double_t *xyz, Double_t *b) const
381{
382 // Method to calculate the integral_0^z of bx/bz,by/bz and (bx/bz)^2+(by/bz)^2
383 b[0]=b[1]=b[2]=0.0;
384 if (fMeasuredMap) {
385 fMeasuredMap->GetTPCRatInt(xyz,b);
386 b[2] /= 100;
387 }
388}
389
390//_______________________________________________________________________
db83d72f 391void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
392{
47c3d315 393 // Method to calculate the integral_0^z of br,bt,bz
db83d72f 394 // in cylindrical coordiates ( -pi<phi<pi convention )
395 b[0]=b[1]=b[2]=0.0;
396 if (fMeasuredMap) {
397 fMeasuredMap->GetTPCIntCyl(rphiz,b);
398 for (int i=3;i--;) b[i] *= fFactorSol;
399 }
eed8a1a2 400}
1dd3d90e 401
402//_______________________________________________________________________
47c3d315 403void AliMagF::GetTPCRatIntCyl(const Double_t *rphiz, Double_t *b) const
404{
405 // Method to calculate the integral_0^z of bx/bz,by/bz and (bx/bz)^2+(by/bz)^2
406 // in cylindrical coordiates ( -pi<phi<pi convention )
407 b[0]=b[1]=b[2]=0.0;
408 if (fMeasuredMap) {
409 fMeasuredMap->GetTPCRatIntCyl(rphiz,b);
410 b[2] /= 100;
411 }
412}
413
414//_______________________________________________________________________
1dd3d90e 415void AliMagF::SetFactorSol(Float_t fc)
416{
417 // set the sign/scale of the current in the L3 according to fgkPolarityConvention
418 switch (fgkPolarityConvention) {
419 case kConvDCS2008: fFactorSol = -fc; break;
420 case kConvLHC : fFactorSol = -fc; break;
421 default : fFactorSol = fc; break; // case kConvMap2005: fFactorSol = fc; break;
422 }
423}
424
425//_______________________________________________________________________
426void AliMagF::SetFactorDip(Float_t fc)
427{
428 // set the sign*scale of the current in the Dipole according to fgkPolarityConvention
429 switch (fgkPolarityConvention) {
430 case kConvDCS2008: fFactorDip = fc; break;
431 case kConvLHC : fFactorDip = -fc; break;
432 default : fFactorDip = fc; break; // case kConvMap2005: fFactorDip = fc; break;
433 }
434}
435
436//_______________________________________________________________________
437Double_t AliMagF::GetFactorSol() const
438{
439 // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe
440 switch (fgkPolarityConvention) {
441 case kConvDCS2008: return -fFactorSol;
442 case kConvLHC : return -fFactorSol;
443 default : return fFactorSol; // case kConvMap2005: return fFactorSol;
444 }
445}
446
447//_______________________________________________________________________
448Double_t AliMagF::GetFactorDip() const
449{
450 // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe
451 switch (fgkPolarityConvention) {
452 case kConvDCS2008: return fFactorDip;
453 case kConvLHC : return -fFactorDip;
454 default : return fFactorDip; // case kConvMap2005: return fFactorDip;
455 }
456}
33fe5eb1 457
458//_____________________________________________________________________________
459AliMagF* AliMagF::CreateFieldMap(Float_t l3Cur, Float_t diCur, Int_t convention, Bool_t uniform,
5cf76849 460 Float_t beamenergy, const Char_t *beamtype, const Char_t *path)
33fe5eb1 461{
462 //------------------------------------------------
463 // The magnetic field map, defined externally...
464 // L3 current 30000 A -> 0.5 T
465 // L3 current 12000 A -> 0.2 T
466 // dipole current 6000 A
467 // The polarities must match the convention (LHC or DCS2008)
468 // unless the special uniform map was used for MC
469 //------------------------------------------------
470 const Float_t l3NominalCurrent1=30000.; // (A)
471 const Float_t l3NominalCurrent2=12000.; // (A)
472 const Float_t diNominalCurrent =6000. ; // (A)
473
474 const Float_t tolerance=0.03; // relative current tolerance
475 const Float_t zero=77.; // "zero" current (A)
476 //
3b94b44f 477 BMap_t map = k5kG;
33fe5eb1 478 double sclL3,sclDip;
479 //
480 Float_t l3Pol = l3Cur > 0 ? 1:-1;
481 Float_t diPol = diCur > 0 ? 1:-1;
482
483 l3Cur = TMath::Abs(l3Cur);
484 diCur = TMath::Abs(diCur);
485 //
486 if (TMath::Abs((sclDip=diCur/diNominalCurrent)-1.) > tolerance && !uniform) {
487 if (diCur <= zero) sclDip = 0.; // some small current.. -> Dipole OFF
488 else {
706eaf0b 489 AliFatalGeneral("AliMagF",Form("Wrong dipole current (%f A)!",diCur));
33fe5eb1 490 }
491 }
492 //
493 if (uniform) {
494 // special treatment of special MC with uniform mag field (normalized to 0.5 T)
495 // no check for scaling/polarities are done
496 map = k5kGUniform;
497 sclL3 = l3Cur/l3NominalCurrent1;
498 }
499 else {
500 if (TMath::Abs((sclL3=l3Cur/l3NominalCurrent1)-1.) < tolerance) map = k5kG;
501 else if (TMath::Abs((sclL3=l3Cur/l3NominalCurrent2)-1.) < tolerance) map = k2kG;
e0418f7d 502 else if (l3Cur <= zero && diCur<=zero) { sclL3=0; sclDip=0; map = k5kGUniform;}
33fe5eb1 503 else {
706eaf0b 504 AliFatalGeneral("AliMagF",Form("Wrong L3 current (%f A)!",l3Cur));
33fe5eb1 505 }
506 }
507 //
17c30c5b 508 if (sclDip!=0 && map!=k5kGUniform) {
509 if ( (l3Cur<=zero) || ((convention==kConvLHC && l3Pol!=diPol) || (convention==kConvDCS2008 && l3Pol==diPol)) ) {
706eaf0b 510 AliFatalGeneral("AliMagF",Form("Wrong combination for L3/Dipole polarities (%c/%c) for convention %d",
17c30c5b 511 l3Pol>0?'+':'-',diPol>0?'+':'-',GetPolarityConvention()));
17c30c5b 512 }
33fe5eb1 513 }
514 //
515 if (l3Pol<0) sclL3 = -sclL3;
516 if (diPol<0) sclDip = -sclDip;
517 //
518 BeamType_t btype = kNoBeamField;
519 TString btypestr = beamtype;
520 btypestr.ToLower();
521 TPRegexp protonBeam("(proton|p)\\s*-?\\s*\\1");
522 TPRegexp ionBeam("(lead|pb|ion|a)\\s*-?\\s*\\1");
523 if (btypestr.Contains(ionBeam)) btype = kBeamTypeAA;
524 else if (btypestr.Contains(protonBeam)) btype = kBeamTypepp;
525 else AliInfoGeneral("AliMagF",Form("Assume no LHC magnet field for the beam type %s, ",beamtype));
526 char ttl[80];
527 sprintf(ttl,"L3: %+5d Dip: %+4d kA; %s | Polarities in %s convention",(int)TMath::Sign(l3Cur,float(sclL3)),
528 (int)TMath::Sign(diCur,float(sclDip)),uniform ? " Constant":"",
529 convention==kConvLHC ? "LHC":"DCS2008");
530 // LHC and DCS08 conventions have opposite dipole polarities
531 if ( GetPolarityConvention() != convention) sclDip = -sclDip;
532 //
5cf76849 533 return new AliMagF("MagneticFieldMap", ttl,sclL3,sclDip,map,btype,beamenergy,2,10.,path);
33fe5eb1 534 //
535}
536
537//_____________________________________________________________________________
538const char* AliMagF::GetBeamTypeText() const
539{
540 const char *beamNA = "No Beam";
541 const char *beamPP = "p-p";
542 const char *beamPbPb= "Pb-Pb";
543 switch ( fBeamType ) {
544 case kBeamTypepp : return beamPP;
545 case kBeamTypeAA : return beamPbPb;
546 case kNoBeamField:
547 default: return beamNA;
548 }
549}
550
77c9a262 551//_____________________________________________________________________________
552void AliMagF::Print(Option_t *opt) const
553{
554 // print short or long info
555 TString opts = opt; opts.ToLower();
556 AliInfo(Form("%s:%s",GetName(),GetTitle()));
557 AliInfo(Form("Solenoid (%+.2f*)%.0f kG, Dipole %s (%+.2f) %s",
558 GetFactorSol(),(fMapType==k5kG||fMapType==k5kGUniform)?5.:2.,
559 fDipoleOFF ? "OFF":"ON",GetFactorDip(),fMapType==k5kGUniform?" |Constant Field!":""));
560 if (opts.Contains("a")) {
561 AliInfo(Form("Machine B fields for %s beam (%.0f GeV): QGrad: %.4f Dipole: %.4f",
562 fBeamType==kBeamTypeAA ? "A-A":(fBeamType==kBeamTypepp ? "p-p":"OFF"),
563 fBeamEnergy,fQuadGradient,fDipoleField));
564 AliInfo(Form("Uses %s of %s",GetParamName(),GetDataFileName()));
565 }
566}