1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
23 #include "AliMagWrapCheb.h"
28 const Double_t AliMagF::fgkSol2DipZ = -700.;
29 const UShort_t AliMagF::fgkPolarityConvention = kConvLHC;
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
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.
48 //_______________________________________________________________________
54 fBeamType(kNoBeamField),
71 // Default constructor
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),
90 fDipoleOFF(factorDip==0.),
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
104 if(integ<0 || integ > 2) {
105 AliWarning(Form("Invalid magnetic field flag: %5d; Helix tracking chosen instead",integ));
108 if (fInteg == 0) fPrecInteg = 0;
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");
115 const char* parname = 0;
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));
122 SetDataFileName(path);
123 SetParamName(parname);
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));
138 //_______________________________________________________________________
139 AliMagF::AliMagF(const AliMagF &src):
140 TVirtualMagField(src),
142 fMapType(src.fMapType),
143 fSolenoid(src.fSolenoid),
144 fBeamType(src.fBeamType),
145 fBeamEnergy(src.fBeamEnergy),
147 fPrecInteg(src.fPrecInteg),
148 fFactorSol(src.fFactorSol),
149 fFactorDip(src.fFactorDip),
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)
159 if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
162 //_______________________________________________________________________
168 //_______________________________________________________________________
169 Bool_t AliMagF::LoadParameterization()
172 AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
176 char* fname = gSystem->ExpandPathName(GetDataFileName());
177 TFile* file = TFile::Open(fname);
179 AliError(Form("Failed to open magnetic field data file %s\n",fname));
183 fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
185 AliError(Form("Did not find field %s in %s\n",GetParamName(),fname));
194 //_______________________________________________________________________
195 void AliMagF::Field(const Double_t *xyz, Double_t *b)
197 // Method to calculate the field at point xyz
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;
205 else MachineField(xyz, b);
209 //_______________________________________________________________________
210 Double_t AliMagF::GetBz(const Double_t *xyz) const
212 // Method to calculate the field at point xyz
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;
221 //_______________________________________________________________________
222 AliMagF& AliMagF::operator=(const AliMagF& src)
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;
232 fPrecInteg = src.fPrecInteg;
233 fFactorSol = src.fFactorSol;
234 fFactorDip = src.fFactorDip;
236 fDipoleOFF = src.fDipoleOFF;
237 fParNames = src.fParNames;
242 //_______________________________________________________________________
243 void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
245 if (btype==kNoBeamField) {
246 fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
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.;
254 fQuadGradient = 22.0002*rigScale;
255 fDipoleField = 37.8781*rigScale;
258 fCCorrField = -9.6980;
260 fACorr1Field = -13.2247;
261 fACorr2Field = 11.7905;
265 //_______________________________________________________________________
266 void AliMagF::MachineField(const Double_t *x, Double_t *b) const
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;
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;
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;
282 double rad2 = x[0] * x[0] + x[1] * x[1];
284 b[0] = b[1] = b[2] = 0;
286 // SIDE C **************************************************
288 if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
289 b[0] = fCCorrField*fFactorDip;
291 else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
292 b[0] = fQuadGradient*x[1];
293 b[1] = fQuadGradient*x[0];
295 else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
296 b[0] = -fQuadGradient*x[1];
297 b[1] = -fQuadGradient*x[0];
299 else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
300 b[0] = -fQuadGradient*x[1];
301 b[1] = -fQuadGradient*x[0];
303 else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
304 b[0] = fQuadGradient*x[1];
305 b[1] = fQuadGradient*x[0];
307 else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
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;
318 // SIDE A **************************************************
320 if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
321 // Compensator magnet at z = 1075 m
322 b[0] = fACorr1Field*fFactorDip;
325 if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
326 b[0] = fACorr2Field*fFactorDip;
328 else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
329 b[0] = -fQuadGradient*x[1];
330 b[1] = -fQuadGradient*x[0];
332 else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
333 b[0] = fQuadGradient*x[1];
334 b[1] = fQuadGradient*x[0];
336 else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
337 b[0] = fQuadGradient*x[1];
338 b[1] = fQuadGradient*x[0];
340 else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
341 b[0] = -fQuadGradient*x[1];
342 b[1] = -fQuadGradient*x[0];
344 else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
345 b[1] = -fDipoleField;
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) {
357 //_______________________________________________________________________
358 void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
360 // Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
363 fMeasuredMap->GetTPCInt(xyz,b);
364 for (int i=3;i--;) b[i] *= fFactorSol;
368 //_______________________________________________________________________
369 void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
371 // Method to calculate the integral of magnetic integral from point to nearest cathode plane
372 // in cylindrical coordiates ( -pi<phi<pi convention )
375 fMeasuredMap->GetTPCIntCyl(rphiz,b);
376 for (int i=3;i--;) b[i] *= fFactorSol;
380 //_______________________________________________________________________
381 void AliMagF::SetFactorSol(Float_t fc)
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;
391 //_______________________________________________________________________
392 void AliMagF::SetFactorDip(Float_t fc)
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;
402 //_______________________________________________________________________
403 Double_t AliMagF::GetFactorSol() const
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;
413 //_______________________________________________________________________
414 Double_t AliMagF::GetFactorDip() const
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;
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)
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)
440 const Float_t tolerance=0.03; // relative current tolerance
441 const Float_t zero=77.; // "zero" current (A)
446 Float_t l3Pol = l3Cur > 0 ? 1:-1;
447 Float_t diPol = diCur > 0 ? 1:-1;
449 l3Cur = TMath::Abs(l3Cur);
450 diCur = TMath::Abs(diCur);
452 if (TMath::Abs((sclDip=diCur/diNominalCurrent)-1.) > tolerance && !uniform) {
453 if (diCur <= zero) sclDip = 0.; // some small current.. -> Dipole OFF
455 AliErrorGeneral("AliMagF",Form("Wrong dipole current (%f A)!",diCur));
461 // special treatment of special MC with uniform mag field (normalized to 0.5 T)
462 // no check for scaling/polarities are done
464 sclL3 = l3Cur/l3NominalCurrent1;
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 && diCur<=zero) { sclL3=0; sclDip=0; map = k5kGUniform;}
471 AliErrorGeneral("AliMagF",Form("Wrong L3 current (%f A)!",l3Cur));
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()));
484 if (l3Pol<0) sclL3 = -sclL3;
485 if (diPol<0) sclDip = -sclDip;
487 BeamType_t btype = kNoBeamField;
488 TString btypestr = beamtype;
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));
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;
502 return new AliMagF("MagneticFieldMap", ttl,sclL3,sclDip,map,btype,beamenergy,2,10.,path);
506 //_____________________________________________________________________________
507 const char* AliMagF::GetBeamTypeText() const
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;
516 default: return beamNA;