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)
237 if (src.fMeasuredMap) {
238 if (fMeasuredMap) delete fMeasuredMap;
239 fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
241 SetName(src.GetName());
242 fSolenoid = src.fSolenoid;
243 fBeamType = src.fBeamType;
244 fBeamEnergy = src.fBeamEnergy;
246 fPrecInteg = src.fPrecInteg;
247 fFactorSol = src.fFactorSol;
248 fFactorDip = src.fFactorDip;
250 fDipoleOFF = src.fDipoleOFF;
251 fParNames = src.fParNames;
256 //_______________________________________________________________________
257 void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
259 if (btype==kNoBeamField) {
260 fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
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
266 if (btype==kBeamTypeAA || btype==kBeamTypepA || btype==kBeamTypeAp) rigScale *= 208./82.;
268 fQuadGradient = 22.0002*rigScale;
269 fDipoleField = 37.8781*rigScale;
272 fCCorrField = -9.6980;
274 fACorr1Field = -13.2247;
275 fACorr2Field = 11.7905;
279 //_______________________________________________________________________
280 void AliMagF::MachineField(const Double_t *x, Double_t *b) const
282 // ---- This is the ZDC part
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;
287 const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
288 const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
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;
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;
296 double rad2 = x[0] * x[0] + x[1] * x[1];
298 b[0] = b[1] = b[2] = 0;
300 // SIDE C **************************************************
302 if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
303 b[0] = fCCorrField*fFactorDip;
305 else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
306 b[0] = fQuadGradient*x[1];
307 b[1] = fQuadGradient*x[0];
309 else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
310 b[0] = -fQuadGradient*x[1];
311 b[1] = -fQuadGradient*x[0];
313 else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
314 b[0] = -fQuadGradient*x[1];
315 b[1] = -fQuadGradient*x[0];
317 else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
318 b[0] = fQuadGradient*x[1];
319 b[1] = fQuadGradient*x[0];
321 else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
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) {
327 b[1] = -fDipoleField;
332 // SIDE A **************************************************
334 if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
335 // Compensator magnet at z = 1075 m
336 b[0] = fACorr1Field*fFactorDip;
339 if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
340 b[0] = fACorr2Field*fFactorDip;
342 else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
343 b[0] = -fQuadGradient*x[1];
344 b[1] = -fQuadGradient*x[0];
346 else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
347 b[0] = fQuadGradient*x[1];
348 b[1] = fQuadGradient*x[0];
350 else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
351 b[0] = fQuadGradient*x[1];
352 b[1] = fQuadGradient*x[0];
354 else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
355 b[0] = -fQuadGradient*x[1];
356 b[1] = -fQuadGradient*x[0];
358 else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
359 b[1] = -fDipoleField;
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) {
371 //_______________________________________________________________________
372 void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
374 // Method to calculate the integral_0^z of br,bt,bz
377 fMeasuredMap->GetTPCInt(xyz,b);
378 for (int i=3;i--;) b[i] *= fFactorSol;
382 //_______________________________________________________________________
383 void AliMagF::GetTPCRatInt(const Double_t *xyz, Double_t *b) const
385 // Method to calculate the integral_0^z of bx/bz,by/bz and (bx/bz)^2+(by/bz)^2
388 fMeasuredMap->GetTPCRatInt(xyz,b);
393 //_______________________________________________________________________
394 void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
396 // Method to calculate the integral_0^z of br,bt,bz
397 // in cylindrical coordiates ( -pi<phi<pi convention )
400 fMeasuredMap->GetTPCIntCyl(rphiz,b);
401 for (int i=3;i--;) b[i] *= fFactorSol;
405 //_______________________________________________________________________
406 void AliMagF::GetTPCRatIntCyl(const Double_t *rphiz, Double_t *b) const
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 )
412 fMeasuredMap->GetTPCRatIntCyl(rphiz,b);
417 //_______________________________________________________________________
418 void AliMagF::SetFactorSol(Float_t fc)
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;
428 //_______________________________________________________________________
429 void AliMagF::SetFactorDip(Float_t fc)
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;
439 //_______________________________________________________________________
440 Double_t AliMagF::GetFactorSol() const
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;
450 //_______________________________________________________________________
451 Double_t AliMagF::GetFactorDip() const
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;
461 //_____________________________________________________________________________
462 AliMagF* AliMagF::CreateFieldMap(Float_t l3Cur, Float_t diCur, Int_t convention, Bool_t uniform,
463 Float_t beamenergy, const Char_t *beamtype, const Char_t *path)
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)
477 const Float_t tolerance=0.03; // relative current tolerance
478 const Float_t zero=77.; // "zero" current (A)
483 Float_t l3Pol = l3Cur > 0 ? 1:-1;
484 Float_t diPol = diCur > 0 ? 1:-1;
486 l3Cur = TMath::Abs(l3Cur);
487 diCur = TMath::Abs(diCur);
489 if (TMath::Abs((sclDip=diCur/diNominalCurrent)-1.) > tolerance && !uniform) {
490 if (diCur <= zero) sclDip = 0.; // some small current.. -> Dipole OFF
492 AliFatalGeneral("AliMagF",Form("Wrong dipole current (%f A)!",diCur));
497 // special treatment of special MC with uniform mag field (normalized to 0.5 T)
498 // no check for scaling/polarities are done
500 sclL3 = l3Cur/l3NominalCurrent1;
503 if (TMath::Abs((sclL3=l3Cur/l3NominalCurrent1)-1.) < tolerance) map = k5kG;
504 else if (TMath::Abs((sclL3=l3Cur/l3NominalCurrent2)-1.) < tolerance) map = k2kG;
505 else if (l3Cur <= zero && diCur<=zero) { sclL3=0; sclDip=0; map = k5kGUniform;}
507 AliFatalGeneral("AliMagF",Form("Wrong L3 current (%f A)!",l3Cur));
511 if (sclDip!=0 && map!=k5kGUniform) {
512 if ( (l3Cur<=zero) || ((convention==kConvLHC && l3Pol!=diPol) || (convention==kConvDCS2008 && l3Pol==diPol)) ) {
513 AliFatalGeneral("AliMagF",Form("Wrong combination for L3/Dipole polarities (%c/%c) for convention %d",
514 l3Pol>0?'+':'-',diPol>0?'+':'-',GetPolarityConvention()));
518 if (l3Pol<0) sclL3 = -sclL3;
519 if (diPol<0) sclDip = -sclDip;
521 BeamType_t btype = kNoBeamField;
522 TString btypestr = beamtype;
524 TPRegexp protonBeam("(proton|p)\\s*-?\\s*\\1");
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)");
528 if (btypestr.Contains(ionBeam)) btype = kBeamTypeAA;
529 else if (btypestr.Contains(protonBeam)) btype = kBeamTypepp;
530 else if (btypestr.Contains(protonionBeam)) btype = kBeamTypepA;
531 else if (btypestr.Contains(ionprotonBeam)) btype = kBeamTypeAp;
532 else AliInfoGeneral("AliMagF",Form("Assume no LHC magnet field for the beam type %s, ",beamtype));
534 snprintf(ttl,79,"L3: %+5d Dip: %+4d kA; %s | Polarities in %s convention",(int)TMath::Sign(l3Cur,float(sclL3)),
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;
540 return new AliMagF("MagneticFieldMap", ttl,sclL3,sclDip,map,btype,beamenergy,2,10.,path);
544 //_____________________________________________________________________________
545 const char* AliMagF::GetBeamTypeText() const
547 const char *beamNA = "No Beam";
548 const char *beamPP = "p-p";
549 const char *beamPbPb= "A-A";
550 const char *beamPPb = "p-A";
551 const char *beamPbP = "A-p";
552 switch ( fBeamType ) {
553 case kBeamTypepp : return beamPP;
554 case kBeamTypeAA : return beamPbPb;
555 case kBeamTypepA : return beamPPb;
556 case kBeamTypeAp : return beamPbP;
558 default: return beamNA;
562 //_____________________________________________________________________________
563 void AliMagF::Print(Option_t *opt) const
565 // print short or long info
566 TString opts = opt; opts.ToLower();
567 AliInfo(Form("%s:%s",GetName(),GetTitle()));
568 AliInfo(Form("Solenoid (%+.2f*)%.0f kG, Dipole %s (%+.2f) %s",
569 GetFactorSol(),(fMapType==k5kG||fMapType==k5kGUniform)?5.:2.,
570 fDipoleOFF ? "OFF":"ON",GetFactorDip(),fMapType==k5kGUniform?" |Constant Field!":""));
571 if (opts.Contains("a")) {
572 AliInfo(Form("Machine B fields for %s beam (%.0f GeV): QGrad: %.4f Dipole: %.4f",
574 fBeamEnergy,fQuadGradient,fDipoleField));
575 AliInfo(Form("Uses %s of %s",GetParamName(),GetDataFileName()));