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