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