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 = AliMagF::kConvLHC;
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 -> positive Bz
41 positive Dip current -> negative Bx
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.
47 -----------------------------------------------
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}
55 The same applies to integral in cylindrical coordinates:
57 GetTPCIntRatCyl(rphiz,b)
58 They accept the R,Phi,Z coordinate (-pi<phi<pi) and return the field
59 integrals in cyl. coordinates.
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)
64 Note: the integrals are defined for the range -300<Z<300 and 0<R<300
66 //_______________________________________________________________________
72 fBeamType(kNoBeamField),
89 // Default constructor
93 //_______________________________________________________________________
94 AliMagF::AliMagF(const char *name, const char* title, Double_t factorSol, Double_t factorDip,
95 BMap_t maptype, BeamType_t bt, Double_t be,Int_t integ, Double_t fmax, const char* path):
96 TVirtualMagField(name),
108 fDipoleOFF(factorDip==0.),
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
122 if(integ<0 || integ > 2) {
123 AliWarning(Form("Invalid magnetic field flag: %5d; Helix tracking chosen instead",integ));
126 if (fInteg == 0) fPrecInteg = 0;
128 if (fBeamEnergy<=0 && fBeamType!=kNoBeamField) {
129 if (fBeamType == kBeamTypepp) fBeamEnergy = 7000.; // max proton energy
130 else if (fBeamType == kBeamTypeAA) fBeamEnergy = 2760; // max PbPb energy
131 else if (fBeamType == kBeamTypepA || fBeamType == kBeamTypeAp) fBeamEnergy = 2760; // same rigitiy max PbPb energy
132 AliInfo("Maximim possible beam energy for requested beam is assumed");
134 const char* parname = 0;
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));
141 SetDataFileName(path);
142 SetParamName(parname);
144 LoadParameterization();
145 InitMachineField(fBeamType,fBeamEnergy);
146 double xyz[3]={0.,0.,0.};
147 fSolenoid = GetBz(xyz);
148 SetFactorSol(factorSol);
149 SetFactorDip(factorDip);
153 //_______________________________________________________________________
154 AliMagF::AliMagF(const AliMagF &src):
155 TVirtualMagField(src),
157 fMapType(src.fMapType),
158 fSolenoid(src.fSolenoid),
159 fBeamType(src.fBeamType),
160 fBeamEnergy(src.fBeamEnergy),
162 fPrecInteg(src.fPrecInteg),
163 fFactorSol(src.fFactorSol),
164 fFactorDip(src.fFactorDip),
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)
174 if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
177 //_______________________________________________________________________
183 //_______________________________________________________________________
184 Bool_t AliMagF::LoadParameterization()
187 AliFatal(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
190 char* fname = gSystem->ExpandPathName(GetDataFileName());
191 TFile* file = TFile::Open(fname);
193 AliFatal(Form("Failed to open magnetic field data file %s\n",fname));
196 fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
198 AliFatal(Form("Did not find field %s in %s\n",GetParamName(),fname));
206 //_______________________________________________________________________
207 void AliMagF::Field(const Double_t *xyz, Double_t *b)
209 // Method to calculate the field at point xyz
211 // b[0]=b[1]=b[2]=0.0;
212 if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
213 fMeasuredMap->Field(xyz,b);
214 if (xyz[2]>fgkSol2DipZ || fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
215 else for (int i=3;i--;) b[i] *= fFactorDip;
217 else MachineField(xyz, b);
221 //_______________________________________________________________________
222 Double_t AliMagF::GetBz(const Double_t *xyz) const
224 // Method to calculate the field at point xyz
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;
233 //_______________________________________________________________________
234 AliMagF& AliMagF::operator=(const AliMagF& src)
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;
244 fPrecInteg = src.fPrecInteg;
245 fFactorSol = src.fFactorSol;
246 fFactorDip = src.fFactorDip;
248 fDipoleOFF = src.fDipoleOFF;
249 fParNames = src.fParNames;
254 //_______________________________________________________________________
255 void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
257 if (btype==kNoBeamField) {
258 fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
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
264 if (btype == kBeamTypeAA || kBeamTypepA || kBeamTypeAp) rigScale *= 208./82.;
266 fQuadGradient = 22.0002*rigScale;
267 fDipoleField = 37.8781*rigScale;
270 fCCorrField = -9.6980;
272 fACorr1Field = -13.2247;
273 fACorr2Field = 11.7905;
277 //_______________________________________________________________________
278 void AliMagF::MachineField(const Double_t *x, Double_t *b) const
280 // ---- This is the ZDC part
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;
285 const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
286 const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
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;
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;
294 double rad2 = x[0] * x[0] + x[1] * x[1];
296 b[0] = b[1] = b[2] = 0;
298 // SIDE C **************************************************
300 if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
301 b[0] = fCCorrField*fFactorDip;
303 else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
304 b[0] = fQuadGradient*x[1];
305 b[1] = fQuadGradient*x[0];
307 else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
308 b[0] = -fQuadGradient*x[1];
309 b[1] = -fQuadGradient*x[0];
311 else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
312 b[0] = -fQuadGradient*x[1];
313 b[1] = -fQuadGradient*x[0];
315 else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
316 b[0] = fQuadGradient*x[1];
317 b[1] = fQuadGradient*x[0];
319 else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
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) {
325 b[1] = -fDipoleField;
330 // SIDE A **************************************************
332 if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
333 // Compensator magnet at z = 1075 m
334 b[0] = fACorr1Field*fFactorDip;
337 if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
338 b[0] = fACorr2Field*fFactorDip;
340 else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
341 b[0] = -fQuadGradient*x[1];
342 b[1] = -fQuadGradient*x[0];
344 else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
345 b[0] = fQuadGradient*x[1];
346 b[1] = fQuadGradient*x[0];
348 else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
349 b[0] = fQuadGradient*x[1];
350 b[1] = fQuadGradient*x[0];
352 else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
353 b[0] = -fQuadGradient*x[1];
354 b[1] = -fQuadGradient*x[0];
356 else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
357 b[1] = -fDipoleField;
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) {
369 //_______________________________________________________________________
370 void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
372 // Method to calculate the integral_0^z of br,bt,bz
375 fMeasuredMap->GetTPCInt(xyz,b);
376 for (int i=3;i--;) b[i] *= fFactorSol;
380 //_______________________________________________________________________
381 void AliMagF::GetTPCRatInt(const Double_t *xyz, Double_t *b) const
383 // Method to calculate the integral_0^z of bx/bz,by/bz and (bx/bz)^2+(by/bz)^2
386 fMeasuredMap->GetTPCRatInt(xyz,b);
391 //_______________________________________________________________________
392 void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
394 // Method to calculate the integral_0^z of br,bt,bz
395 // in cylindrical coordiates ( -pi<phi<pi convention )
398 fMeasuredMap->GetTPCIntCyl(rphiz,b);
399 for (int i=3;i--;) b[i] *= fFactorSol;
403 //_______________________________________________________________________
404 void AliMagF::GetTPCRatIntCyl(const Double_t *rphiz, Double_t *b) const
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 )
410 fMeasuredMap->GetTPCRatIntCyl(rphiz,b);
415 //_______________________________________________________________________
416 void AliMagF::SetFactorSol(Float_t fc)
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;
426 //_______________________________________________________________________
427 void AliMagF::SetFactorDip(Float_t fc)
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;
437 //_______________________________________________________________________
438 Double_t AliMagF::GetFactorSol() const
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;
448 //_______________________________________________________________________
449 Double_t AliMagF::GetFactorDip() const
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;
459 //_____________________________________________________________________________
460 AliMagF* AliMagF::CreateFieldMap(Float_t l3Cur, Float_t diCur, Int_t convention, Bool_t uniform,
461 Float_t beamenergy, const Char_t *beamtype, const Char_t *path)
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)
475 const Float_t tolerance=0.03; // relative current tolerance
476 const Float_t zero=77.; // "zero" current (A)
481 Float_t l3Pol = l3Cur > 0 ? 1:-1;
482 Float_t diPol = diCur > 0 ? 1:-1;
484 l3Cur = TMath::Abs(l3Cur);
485 diCur = TMath::Abs(diCur);
487 if (TMath::Abs((sclDip=diCur/diNominalCurrent)-1.) > tolerance && !uniform) {
488 if (diCur <= zero) sclDip = 0.; // some small current.. -> Dipole OFF
490 AliFatalGeneral("AliMagF",Form("Wrong dipole current (%f A)!",diCur));
495 // special treatment of special MC with uniform mag field (normalized to 0.5 T)
496 // no check for scaling/polarities are done
498 sclL3 = l3Cur/l3NominalCurrent1;
501 if (TMath::Abs((sclL3=l3Cur/l3NominalCurrent1)-1.) < tolerance) map = k5kG;
502 else if (TMath::Abs((sclL3=l3Cur/l3NominalCurrent2)-1.) < tolerance) map = k2kG;
503 else if (l3Cur <= zero && diCur<=zero) { sclL3=0; sclDip=0; map = k5kGUniform;}
505 AliFatalGeneral("AliMagF",Form("Wrong L3 current (%f A)!",l3Cur));
509 if (sclDip!=0 && map!=k5kGUniform) {
510 if ( (l3Cur<=zero) || ((convention==kConvLHC && l3Pol!=diPol) || (convention==kConvDCS2008 && l3Pol==diPol)) ) {
511 AliFatalGeneral("AliMagF",Form("Wrong combination for L3/Dipole polarities (%c/%c) for convention %d",
512 l3Pol>0?'+':'-',diPol>0?'+':'-',GetPolarityConvention()));
516 if (l3Pol<0) sclL3 = -sclL3;
517 if (diPol<0) sclDip = -sclDip;
519 BeamType_t btype = kNoBeamField;
520 TString btypestr = beamtype;
522 TPRegexp protonBeam("(proton|p)\\s*-?\\s*\\1");
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)");
526 if (btypestr.Contains(ionBeam)) btype = kBeamTypeAA;
527 else if (btypestr.Contains(protonBeam)) btype = kBeamTypepp;
528 else if (btypestr.Contains(protonionBeam)) btype = kBeamTypepA;
529 else if (btypestr.Contains(ionprotonBeam)) btype = kBeamTypeAp;
530 else AliInfoGeneral("AliMagF",Form("Assume no LHC magnet field for the beam type %s, ",beamtype));
532 snprintf(ttl,79,"L3: %+5d Dip: %+4d kA; %s | Polarities in %s convention",(int)TMath::Sign(l3Cur,float(sclL3)),
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;
538 return new AliMagF("MagneticFieldMap", ttl,sclL3,sclDip,map,btype,beamenergy,2,10.,path);
542 //_____________________________________________________________________________
543 const char* AliMagF::GetBeamTypeText() const
545 const char *beamNA = "No Beam";
546 const char *beamPP = "p-p";
547 const char *beamPbPb= "A-A";
548 const char *beamPPb = "p-A";
549 const char *beamPbP = "A-p";
550 switch ( fBeamType ) {
551 case kBeamTypepp : return beamPP;
552 case kBeamTypeAA : return beamPbPb;
553 case kBeamTypepA : return beamPPb;
554 case kBeamTypeAp : return beamPbP;
556 default: return beamNA;
560 //_____________________________________________________________________________
561 void AliMagF::Print(Option_t *opt) const
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",
572 fBeamEnergy,fQuadGradient,fDipoleField));
573 AliInfo(Form("Uses %s of %s",GetParamName(),GetDataFileName()));