Fixes for bug #52499: Field polarities inconsistiency
[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>
fe4da5cc 20
21#include "AliMagF.h"
db83d72f 22#include "AliMagWrapCheb.h"
23#include "AliLog.h"
972ca52f 24
fe4da5cc 25ClassImp(AliMagF)
26
9251fceb 27const Double_t AliMagF::fgkSol2DipZ = -700.;
1dd3d90e 28const UShort_t AliMagF::fgkPolarityConvention = kConvLHC;
db83d72f 29
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
40 positive L3 current -> negative Bz
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.
46*/
e2afb3b6 47//_______________________________________________________________________
48AliMagF::AliMagF():
db83d72f 49 TVirtualMagField(),
50 fMeasuredMap(0),
51 fMapType(k5kG),
52 fSolenoid(0),
53 fBeamType(kNoBeamField),
54 fBeamEnergy(0),
db83d72f 55 //
e2afb3b6 56 fInteg(0),
db83d72f 57 fPrecInteg(0),
58 fFactorSol(1.),
59 fFactorDip(1.),
60 fMax(15),
61 fDipoleOFF(kFALSE),
e2afb3b6 62 //
db83d72f 63 fQuadGradient(0),
64 fDipoleField(0),
65 fCCorrField(0),
66 fACorr1Field(0),
67 fACorr2Field(0),
68 fParNames("","")
69{
e2afb3b6 70 // Default constructor
71 //
72}
73
74//_______________________________________________________________________
db83d72f 75AliMagF::AliMagF(const char *name, const char* title, Int_t integ,
76 Double_t factorSol, Double_t factorDip,
77 Double_t fmax, BMap_t maptype, const char* path,
9251fceb 78 BeamType_t bt, Double_t be):
db83d72f 79 TVirtualMagField(name),
80 fMeasuredMap(0),
81 fMapType(maptype),
82 fSolenoid(0),
83 fBeamType(bt),
84 fBeamEnergy(be),
db83d72f 85 //
86 fInteg(integ),
604e0531 87 fPrecInteg(1),
db83d72f 88 fFactorSol(1.),
89 fFactorDip(1.),
972ca52f 90 fMax(fmax),
db83d72f 91 fDipoleOFF(factorDip==0.),
92 //
93 fQuadGradient(0),
94 fDipoleField(0),
95 fCCorrField(0),
96 fACorr1Field(0),
97 fACorr2Field(0),
98 fParNames("","")
fe4da5cc 99{
9251fceb 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
aee8290b 103 //
db83d72f 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;
aee8290b 110 //
db83d72f 111 const char* parname = 0;
112 //
f04e7f5f 113 if (fMapType == k2kG) parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
114 else if (fMapType == k5kG) parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
115 else if (fMapType == k5kGUniform) parname = "Sol30_Dip6_Uniform";
116 else AliFatal(Form("Unknown field identifier %d is requested\n",fMapType));
db83d72f 117 //
118 SetDataFileName(path);
119 SetParamName(parname);
120 //
db83d72f 121 LoadParameterization();
122 InitMachineField(fBeamType,fBeamEnergy);
f04e7f5f 123 double xyz[3]={0.,0.,0.};
124 fSolenoid = GetBz(xyz);
125 SetFactorSol(factorSol);
126 SetFactorDip(factorDip);
fe4da5cc 127}
128
e2afb3b6 129//_______________________________________________________________________
eeda4611 130AliMagF::AliMagF(const AliMagF &src):
db83d72f 131 TVirtualMagField(src),
132 fMeasuredMap(0),
133 fMapType(src.fMapType),
134 fSolenoid(src.fSolenoid),
135 fBeamType(src.fBeamType),
136 fBeamEnergy(src.fBeamEnergy),
eeda4611 137 fInteg(src.fInteg),
138 fPrecInteg(src.fPrecInteg),
db83d72f 139 fFactorSol(src.fFactorSol),
140 fFactorDip(src.fFactorDip),
eeda4611 141 fMax(src.fMax),
db83d72f 142 fDipoleOFF(src.fDipoleOFF),
143 fQuadGradient(src.fQuadGradient),
144 fDipoleField(src.fDipoleField),
145 fCCorrField(src.fCCorrField),
146 fACorr1Field(src.fACorr1Field),
147 fACorr2Field(src.fACorr2Field),
148 fParNames(src.fParNames)
eeda4611 149{
db83d72f 150 if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
eeda4611 151}
152
153//_______________________________________________________________________
db83d72f 154AliMagF::~AliMagF()
ff66b122 155{
db83d72f 156 delete fMeasuredMap;
157}
158
159//_______________________________________________________________________
160Bool_t AliMagF::LoadParameterization()
161{
162 if (fMeasuredMap) {
163 AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
164 return kTRUE;
165 }
ff66b122 166 //
db83d72f 167 char* fname = gSystem->ExpandPathName(GetDataFileName());
168 TFile* file = TFile::Open(fname);
169 if (!file) {
170 AliError(Form("Failed to open magnetic field data file %s\n",fname));
171 return kFALSE;
172 }
ff66b122 173 //
db83d72f 174 fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
175 if (!fMeasuredMap) {
176 AliError(Form("Did not find field %s in %s\n",GetParamName(),fname));
177 return kFALSE;
178 }
179 file->Close();
180 delete file;
181 return kTRUE;
ff66b122 182}
183
db83d72f 184
ff66b122 185//_______________________________________________________________________
db83d72f 186void AliMagF::Field(const Double_t *xyz, Double_t *b)
fe4da5cc 187{
db83d72f 188 // Method to calculate the field at point xyz
aee8290b 189 //
9251fceb 190 // b[0]=b[1]=b[2]=0.0;
191 if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
db83d72f 192 fMeasuredMap->Field(xyz,b);
193 if (xyz[2]>fgkSol2DipZ || fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
9251fceb 194 else for (int i=3;i--;) b[i] *= fFactorDip;
db83d72f 195 }
9251fceb 196 else MachineField(xyz, b);
aee8290b 197 //
fe4da5cc 198}
eeda4611 199
200//_______________________________________________________________________
db83d72f 201Double_t AliMagF::GetBz(const Double_t *xyz) const
eeda4611 202{
db83d72f 203 // Method to calculate the field at point xyz
204 //
9251fceb 205 if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
206 double bz = fMeasuredMap->GetBz(xyz);
207 return (xyz[2]>fgkSol2DipZ || fDipoleOFF) ? bz*fFactorSol : bz*fFactorDip;
db83d72f 208 }
9251fceb 209 else return 0.;
eeda4611 210}
211
212//_______________________________________________________________________
db83d72f 213AliMagF& AliMagF::operator=(const AliMagF& src)
eeda4611 214{
db83d72f 215 if (this != &src && src.fMeasuredMap) {
216 if (fMeasuredMap) delete fMeasuredMap;
217 fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
218 SetName(src.GetName());
219 fSolenoid = src.fSolenoid;
220 fBeamType = src.fBeamType;
221 fBeamEnergy = src.fBeamEnergy;
db83d72f 222 fInteg = src.fInteg;
223 fPrecInteg = src.fPrecInteg;
224 fFactorSol = src.fFactorSol;
225 fFactorDip = src.fFactorDip;
226 fMax = src.fMax;
227 fDipoleOFF = src.fDipoleOFF;
228 fParNames = src.fParNames;
229 }
230 return *this;
eeda4611 231}
232
233//_______________________________________________________________________
db83d72f 234void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
eeda4611 235{
9251fceb 236 if (btype==kNoBeamField || benergy<1.) {
db83d72f 237 fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
9251fceb 238 return;
db83d72f 239 }
240 //
9251fceb 241 double rigScale = benergy/7000.; // scale according to ratio of E/Enominal
242 // for ions assume PbPb (with energy provided per nucleon) and account for A/Z
243 if (btype == kBeamTypeAA) rigScale *= 208./82.;
244 //
245 fQuadGradient = 22.0002*rigScale;
246 fDipoleField = 37.8781*rigScale;
247 //
248 // SIDE C
249 fCCorrField = -9.6980;
250 // SIDE A
251 fACorr1Field = -13.2247;
252 fACorr2Field = 11.7905;
db83d72f 253 //
eeda4611 254}
eed8a1a2 255
db83d72f 256//_______________________________________________________________________
257void AliMagF::MachineField(const Double_t *x, Double_t *b) const
eed8a1a2 258{
db83d72f 259 // ---- This is the ZDC part
9251fceb 260 // Compansators for Alice Muon Arm Dipole
261 const Double_t kBComp1CZ = 1075., kBComp1hDZ = 260./2., kBComp1SqR = 4.0*4.0;
262 const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5;
263 //
264 const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
90ae20c9 265 const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
9251fceb 266 const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5;
267 const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5;
db83d72f 268 //
9251fceb 269 const Double_t kDip1CZ = 6310.8, kDip1hDZ = 945./2., kDip1SqRC = 4.5*4.5, kDip1SqRA = 3.375*3.375;
270 const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.5*4.5, kDip2SqRA = 3.75*3.75;
271 const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4;
db83d72f 272 //
273 double rad2 = x[0] * x[0] + x[1] * x[1];
274 //
9251fceb 275 b[0] = b[1] = b[2] = 0;
276 //
db83d72f 277 // SIDE C **************************************************
278 if(x[2]<0.){
9251fceb 279 if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
280 b[0] = fCCorrField*fFactorDip;
db83d72f 281 }
9251fceb 282 else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
db83d72f 283 b[0] = fQuadGradient*x[1];
284 b[1] = fQuadGradient*x[0];
db83d72f 285 }
9251fceb 286 else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
db83d72f 287 b[0] = -fQuadGradient*x[1];
288 b[1] = -fQuadGradient*x[0];
db83d72f 289 }
9251fceb 290 else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
db83d72f 291 b[0] = -fQuadGradient*x[1];
292 b[1] = -fQuadGradient*x[0];
db83d72f 293 }
9251fceb 294 else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
db83d72f 295 b[0] = fQuadGradient*x[1];
296 b[1] = fQuadGradient*x[0];
db83d72f 297 }
9251fceb 298 else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
db83d72f 299 b[1] = fDipoleField;
db83d72f 300 }
9251fceb 301 else if(TMath::Abs(x[2]+kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRC) {
302 double dxabs = TMath::Abs(x[0])-kDip2DXC;
303 if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRC) {
db83d72f 304 b[1] = -fDipoleField;
db83d72f 305 }
306 }
307 }
308 //
309 // SIDE A **************************************************
310 else{
9251fceb 311 if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
db83d72f 312 // Compensator magnet at z = 1075 m
9251fceb 313 b[0] = fACorr1Field*fFactorDip;
db83d72f 314 }
315 //
9251fceb 316 if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
317 b[0] = fACorr2Field*fFactorDip;
318 }
319 else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
db83d72f 320 b[0] = -fQuadGradient*x[1];
321 b[1] = -fQuadGradient*x[0];
eed8a1a2 322 }
9251fceb 323 else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
324 b[0] = fQuadGradient*x[1];
325 b[1] = fQuadGradient*x[0];
eed8a1a2 326 }
9251fceb 327 else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
328 b[0] = fQuadGradient*x[1];
329 b[1] = fQuadGradient*x[0];
db83d72f 330 }
9251fceb 331 else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
db83d72f 332 b[0] = -fQuadGradient*x[1];
333 b[1] = -fQuadGradient*x[0];
db83d72f 334 }
9251fceb 335 else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
db83d72f 336 b[1] = -fDipoleField;
db83d72f 337 }
9251fceb 338 else if(TMath::Abs(x[2]-kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRA) {
339 double dxabs = TMath::Abs(x[0])-kDip2DXA;
340 if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRA) {
db83d72f 341 b[1] = fDipoleField;
342 }
343 }
344 }
9251fceb 345 //
db83d72f 346}
347
348//_______________________________________________________________________
349void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
350{
351 // Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
352 b[0]=b[1]=b[2]=0.0;
353 if (fMeasuredMap) {
354 fMeasuredMap->GetTPCInt(xyz,b);
355 for (int i=3;i--;) b[i] *= fFactorSol;
356 }
357}
358
359//_______________________________________________________________________
360void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
361{
362 // Method to calculate the integral of magnetic integral from point to nearest cathode plane
363 // in cylindrical coordiates ( -pi<phi<pi convention )
364 b[0]=b[1]=b[2]=0.0;
365 if (fMeasuredMap) {
366 fMeasuredMap->GetTPCIntCyl(rphiz,b);
367 for (int i=3;i--;) b[i] *= fFactorSol;
368 }
eed8a1a2 369}
1dd3d90e 370
371//_______________________________________________________________________
372void AliMagF::SetFactorSol(Float_t fc)
373{
374 // set the sign/scale of the current in the L3 according to fgkPolarityConvention
375 switch (fgkPolarityConvention) {
376 case kConvDCS2008: fFactorSol = -fc; break;
377 case kConvLHC : fFactorSol = -fc; break;
378 default : fFactorSol = fc; break; // case kConvMap2005: fFactorSol = fc; break;
379 }
380}
381
382//_______________________________________________________________________
383void AliMagF::SetFactorDip(Float_t fc)
384{
385 // set the sign*scale of the current in the Dipole according to fgkPolarityConvention
386 switch (fgkPolarityConvention) {
387 case kConvDCS2008: fFactorDip = fc; break;
388 case kConvLHC : fFactorDip = -fc; break;
389 default : fFactorDip = fc; break; // case kConvMap2005: fFactorDip = fc; break;
390 }
391}
392
393//_______________________________________________________________________
394Double_t AliMagF::GetFactorSol() const
395{
396 // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe
397 switch (fgkPolarityConvention) {
398 case kConvDCS2008: return -fFactorSol;
399 case kConvLHC : return -fFactorSol;
400 default : return fFactorSol; // case kConvMap2005: return fFactorSol;
401 }
402}
403
404//_______________________________________________________________________
405Double_t AliMagF::GetFactorDip() const
406{
407 // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe
408 switch (fgkPolarityConvention) {
409 case kConvDCS2008: return fFactorDip;
410 case kConvLHC : return -fFactorDip;
411 default : return fFactorDip; // case kConvMap2005: return fFactorDip;
412 }
413}