[u/mrichter/AliRoot.git] / STEER / AliMagF.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


#include <TClass.h>
#include <TFile.h>
#include <TSystem.h>
#include <TPRegexp.h>

#include "AliMagF.h"
#include "AliMagWrapCheb.h"
#include "AliLog.h"

ClassImp(AliMagF)

const Double_t AliMagF::fgkSol2DipZ    =  -700.;  
const UShort_t AliMagF::fgkPolarityConvention = kConvLHC;

/*
 Explanation for polarity conventions: these are the mapping between the
 current signs and main field components in L3 (Bz) and Dipole (Bx) (in Alice frame)
 1) kConvMap2005: used for the field mapping in 2005
 positive L3  current -> negative Bz
 positive Dip current -> positive Bx 
 2) kConvMapDCS2008: defined by the microswitches/cabling of power converters as of 2008 - 1st half 2009
 positive L3  current -> positive Bz
 positive Dip current -> positive Bx
 3) kConvLHC : defined by LHC
 positive L3  current -> negative Bz
 positive Dip current -> negative Bx
 
 Note: only "negative Bz(L3) with postive Bx(Dipole)" and its inverse was mapped in 2005. Hence 
 the GRP Manager will reject the runs with the current combinations (in the convention defined by the
 static Int_t AliMagF::GetPolarityConvention()) which do not lead to such field polarities.
*/
//_______________________________________________________________________
AliMagF::AliMagF():
  TVirtualMagField(),
  fMeasuredMap(0),
  fMapType(k5kG),
  fSolenoid(0),
  fBeamType(kNoBeamField),
  fBeamEnergy(0),
  //
  fInteg(0),
  fPrecInteg(0),
  fFactorSol(1.),
  fFactorDip(1.),
  fMax(15),
  fDipoleOFF(kFALSE),
  //
  fQuadGradient(0),
  fDipoleField(0),
  fCCorrField(0), 
  fACorr1Field(0),
  fACorr2Field(0),
  fParNames("","")
{
  // Default constructor
  //
}

//_______________________________________________________________________
AliMagF::AliMagF(const char *name, const char* title, Double_t factorSol, Double_t factorDip, 
		 BMap_t maptype, BeamType_t btype, Double_t benergy, Int_t integ, Double_t fmax, 
		 const char* path):
  TVirtualMagField(name),
  fMeasuredMap(0),
  fMapType(maptype),
  fSolenoid(0),
  fBeamType(btype),
  fBeamEnergy(benergy),
  //
  fInteg(integ),
  fPrecInteg(1),
  fFactorSol(1.),
  fFactorDip(1.),
  fMax(fmax),
  fDipoleOFF(factorDip==0.),
  //
  fQuadGradient(0),
  fDipoleField(0),
  fCCorrField(0), 
  fACorr1Field(0),
  fACorr2Field(0),
  fParNames("","")
{
  // Initialize the field with Geant integration option "integ" and max field "fmax,
  // Impose scaling of parameterized L3 field by factorSol and of dipole by factorDip.
  // The "be" is the energy of the beam in GeV/nucleon
  //
  SetTitle(title);
  if(integ<0 || integ > 2) {
    AliWarning(Form("Invalid magnetic field flag: %5d; Helix tracking chosen instead",integ));
    fInteg = 2;
  }
  if (fInteg == 0) fPrecInteg = 0;
  //
  if (fBeamEnergy<=0 && fBeamType!=kNoBeamField) {
    if      (fBeamType == kBeamTypepp) fBeamEnergy = 7000.; // max proton energy
    else if (fBeamType == kBeamTypeAA) fBeamEnergy = 5500;  // max PbPb energy
    AliInfo("Maximim possible beam energy for requested beam is assumed");
  } 
  const char* parname = 0;
  //  
  if      (fMapType == k2kG) parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
  else if (fMapType == k5kG) parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
  else if (fMapType == k5kGUniform) parname = "Sol30_Dip6_Uniform";
  else AliFatal(Form("Unknown field identifier %d is requested\n",fMapType));
  //
  SetDataFileName(path);
  SetParamName(parname);
  //
  LoadParameterization();
  InitMachineField(fBeamType,fBeamEnergy);
  double xyz[3]={0.,0.,0.};
  fSolenoid = GetBz(xyz);
  SetFactorSol(factorSol);
  SetFactorDip(factorDip);
  AliInfo(Form("Alice   B fields: Solenoid (%+.2f*)%.0f kG, Dipole %s (%+.2f) %s",
	       factorSol,(fMapType==k5kG||fMapType==k5kGUniform)?5.:2.,
	       fDipoleOFF ? "OFF":"ON",factorDip,fMapType==k5kGUniform?" |Constant Field!":""));
  AliInfo(Form("Machine B fields for %s beam (%.0f GeV): QGrad: %.4f Dipole: %.4f",
	       fBeamType==kBeamTypeAA ? "A-A":(fBeamType==kBeamTypepp ? "p-p":"OFF"),
	       fBeamEnergy,fQuadGradient,fDipoleField));
}

//_______________________________________________________________________
AliMagF::AliMagF(const AliMagF &src):
  TVirtualMagField(src),
  fMeasuredMap(0),
  fMapType(src.fMapType),
  fSolenoid(src.fSolenoid),
  fBeamType(src.fBeamType),
  fBeamEnergy(src.fBeamEnergy),
  fInteg(src.fInteg),
  fPrecInteg(src.fPrecInteg),
  fFactorSol(src.fFactorSol),
  fFactorDip(src.fFactorDip),
  fMax(src.fMax),
  fDipoleOFF(src.fDipoleOFF),
  fQuadGradient(src.fQuadGradient),
  fDipoleField(src.fDipoleField),
  fCCorrField(src.fCCorrField), 
  fACorr1Field(src.fACorr1Field),
  fACorr2Field(src.fACorr2Field),
  fParNames(src.fParNames)
{
  if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
}

//_______________________________________________________________________
AliMagF::~AliMagF()
{
  delete fMeasuredMap;
}

//_______________________________________________________________________
Bool_t AliMagF::LoadParameterization()
{
  if (fMeasuredMap) {
    AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
    return kTRUE;
  }
  //
  char* fname = gSystem->ExpandPathName(GetDataFileName());
  TFile* file = TFile::Open(fname);
  if (!file) {
    AliError(Form("Failed to open magnetic field data file %s\n",fname)); 
    return kFALSE;
  }
  //
  fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
  if (!fMeasuredMap) {
    AliError(Form("Did not find field %s in %s\n",GetParamName(),fname)); 
    return kFALSE;
  }
  file->Close();
  delete file;
  return kTRUE;
}


//_______________________________________________________________________
void AliMagF::Field(const Double_t *xyz, Double_t *b)
{
  // Method to calculate the field at point  xyz
  //
  //  b[0]=b[1]=b[2]=0.0;
  if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
    fMeasuredMap->Field(xyz,b);
    if (xyz[2]>fgkSol2DipZ || fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
    else                                  for (int i=3;i--;) b[i] *= fFactorDip;    
  }
  else MachineField(xyz, b);
  //
}

//_______________________________________________________________________
Double_t AliMagF::GetBz(const Double_t *xyz) const
{
  // Method to calculate the field at point  xyz
  //
  if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
    double bz = fMeasuredMap->GetBz(xyz);
    return (xyz[2]>fgkSol2DipZ || fDipoleOFF) ? bz*fFactorSol : bz*fFactorDip;    
  }
  else return 0.;
}

//_______________________________________________________________________
AliMagF& AliMagF::operator=(const AliMagF& src)
{
  if (this != &src && src.fMeasuredMap) { 
    if (fMeasuredMap) delete fMeasuredMap;
    fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
    SetName(src.GetName());
    fSolenoid    = src.fSolenoid;
    fBeamType    = src.fBeamType;
    fBeamEnergy  = src.fBeamEnergy;
    fInteg       = src.fInteg;
    fPrecInteg   = src.fPrecInteg;
    fFactorSol   = src.fFactorSol;
    fFactorDip   = src.fFactorDip;
    fMax         = src.fMax;
    fDipoleOFF   = src.fDipoleOFF;
    fParNames    = src.fParNames;
  }
  return *this;
}

//_______________________________________________________________________
void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
{
  if (btype==kNoBeamField) {
    fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
    return;
  }
  //
  double rigScale = benergy/7000.;   // scale according to ratio of E/Enominal
  // for ions assume PbPb (with energy provided per nucleon) and account for A/Z
  if (btype == kBeamTypeAA) rigScale *= 208./82.;
  //
  fQuadGradient = 22.0002*rigScale;
  fDipoleField  = 37.8781*rigScale;
  //
  // SIDE C
  fCCorrField   = -9.6980;
  // SIDE A
  fACorr1Field  = -13.2247;
  fACorr2Field  =  11.7905;
  //
}

//_______________________________________________________________________
void AliMagF::MachineField(const Double_t *x, Double_t *b) const
{
  // ---- This is the ZDC part
  // Compansators for Alice Muon Arm Dipole
  const Double_t kBComp1CZ = 1075., kBComp1hDZ = 260./2., kBComp1SqR = 4.0*4.0; 
  const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5; 
  //  
  const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
  const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
  const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5;
  const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5;
  //
  const Double_t kDip1CZ = 6310.8,  kDip1hDZ = 945./2., kDip1SqRC = 4.5*4.5, kDip1SqRA = 3.375*3.375;
  const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.5*4.5, kDip2SqRA = 3.75*3.75;
  const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4;
  //
  double rad2 = x[0] * x[0] + x[1] * x[1];
  //
  b[0] = b[1] = b[2] = 0;
  //
  // SIDE C **************************************************
  if(x[2]<0.){  
    if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
      b[0] = fCCorrField*fFactorDip;
    } 
    else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
      b[0] = fQuadGradient*x[1];
      b[1] = fQuadGradient*x[0];
    }
    else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
      b[0] = -fQuadGradient*x[1];
      b[1] = -fQuadGradient*x[0];
    }
    else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
      b[0] = -fQuadGradient*x[1];
      b[1] = -fQuadGradient*x[0];
    }
    else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
      b[0] = fQuadGradient*x[1];
      b[1] = fQuadGradient*x[0];
    }
    else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
      b[1] = fDipoleField;
    }
    else if(TMath::Abs(x[2]+kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRC) {
      double dxabs = TMath::Abs(x[0])-kDip2DXC;
      if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRC) {
	b[1] = -fDipoleField;
      }
    }
  }
  //
  // SIDE A **************************************************
  else{        
    if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
      // Compensator magnet at z = 1075 m 
      b[0] = fACorr1Field*fFactorDip;
    }
    //
    if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
      b[0] = fACorr2Field*fFactorDip;
    }
    else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
      b[0] = -fQuadGradient*x[1];
      b[1] = -fQuadGradient*x[0];
    }
    else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
      b[0] =  fQuadGradient*x[1];
      b[1] =  fQuadGradient*x[0];
    }
    else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
      b[0] =  fQuadGradient*x[1];
      b[1] =  fQuadGradient*x[0];
    }
    else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
      b[0] = -fQuadGradient*x[1];
      b[1] = -fQuadGradient*x[0];
    }
    else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
      b[1] = -fDipoleField;
    }
    else if(TMath::Abs(x[2]-kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRA) {
      double dxabs = TMath::Abs(x[0])-kDip2DXA;
      if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRA) {
	b[1] = fDipoleField;
      }
    }
  }
  //
}

//_______________________________________________________________________
void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
{
  // Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
  b[0]=b[1]=b[2]=0.0;
  if (fMeasuredMap) {
    fMeasuredMap->GetTPCInt(xyz,b);
    for (int i=3;i--;) b[i] *= fFactorSol;
  }
}

//_______________________________________________________________________
void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
{
  // Method to calculate the integral of magnetic integral from point to nearest cathode plane
  // in cylindrical coordiates ( -pi<phi<pi convention )
  b[0]=b[1]=b[2]=0.0;
  if (fMeasuredMap) {
    fMeasuredMap->GetTPCIntCyl(rphiz,b);
    for (int i=3;i--;) b[i] *= fFactorSol;
  }
}

//_______________________________________________________________________
void AliMagF::SetFactorSol(Float_t fc)
{
  // set the sign/scale of the current in the L3 according to fgkPolarityConvention
  switch (fgkPolarityConvention) {
  case kConvDCS2008: fFactorSol = -fc; break;
  case kConvLHC    : fFactorSol = -fc; break;
  default          : fFactorSol =  fc; break;  // case kConvMap2005: fFactorSol =  fc; break;
  }
}

//_______________________________________________________________________
void AliMagF::SetFactorDip(Float_t fc)
{
  // set the sign*scale of the current in the Dipole according to fgkPolarityConvention
  switch (fgkPolarityConvention) {
  case kConvDCS2008: fFactorDip =  fc; break;
  case kConvLHC    : fFactorDip = -fc; break;
  default          : fFactorDip =  fc; break;  // case kConvMap2005: fFactorDip =  fc; break;
  }
}

//_______________________________________________________________________
Double_t AliMagF::GetFactorSol() const
{
  // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe 
  switch (fgkPolarityConvention) {
  case kConvDCS2008: return -fFactorSol;
  case kConvLHC    : return -fFactorSol;
  default          : return  fFactorSol;       //  case kConvMap2005: return  fFactorSol;
  }
}

//_______________________________________________________________________
Double_t AliMagF::GetFactorDip() const
{
  // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe 
  switch (fgkPolarityConvention) {
  case kConvDCS2008: return  fFactorDip;
  case kConvLHC    : return -fFactorDip;
  default          : return  fFactorDip;       //  case kConvMap2005: return  fFactorDip;
  }
}

//_____________________________________________________________________________
AliMagF* AliMagF::CreateFieldMap(Float_t l3Cur, Float_t diCur, Int_t convention, Bool_t uniform,
				 Float_t beamenergy, const Char_t *beamtype, const Char_t *path) 
{
  //------------------------------------------------
  // The magnetic field map, defined externally...
  // L3 current 30000 A  -> 0.5 T
  // L3 current 12000 A  -> 0.2 T
  // dipole current 6000 A
  // The polarities must match the convention (LHC or DCS2008) 
  // unless the special uniform map was used for MC
  //------------------------------------------------
  const Float_t l3NominalCurrent1=30000.; // (A)
  const Float_t l3NominalCurrent2=12000.; // (A)
  const Float_t diNominalCurrent =6000. ; // (A)

  const Float_t tolerance=0.03; // relative current tolerance
  const Float_t zero=77.;       // "zero" current (A)
  //
  BMap_t map;
  double sclL3,sclDip;
  //
  Float_t l3Pol = l3Cur > 0 ? 1:-1;
  Float_t diPol = diCur > 0 ? 1:-1;
 
  l3Cur = TMath::Abs(l3Cur);
  diCur = TMath::Abs(diCur);
  //
  if (TMath::Abs((sclDip=diCur/diNominalCurrent)-1.) > tolerance && !uniform) {
    if (diCur <= zero) sclDip = 0.; // some small current.. -> Dipole OFF
    else {
      AliErrorGeneral("AliMagF",Form("Wrong dipole current (%f A)!",diCur));
      return 0;
    }
  }
  //
  if (uniform) { 
    // special treatment of special MC with uniform mag field (normalized to 0.5 T)
    // no check for scaling/polarities are done
    map   = k5kGUniform;
    sclL3 = l3Cur/l3NominalCurrent1; 
  }
  else {
    if      (TMath::Abs((sclL3=l3Cur/l3NominalCurrent1)-1.) < tolerance) map  = k5kG;
    else if (TMath::Abs((sclL3=l3Cur/l3NominalCurrent2)-1.) < tolerance) map  = k2kG;
    else if (l3Cur <= zero)                                { sclL3 = 0;  map  = k5kGUniform;}
    else {
      AliErrorGeneral("AliMagF",Form("Wrong L3 current (%f A)!",l3Cur));
      return 0;
    }
  }
  //
  if (sclDip!=0 && (map==k5kG || map==k2kG) &&
      ((convention==kConvLHC     && l3Pol!=diPol) ||
       (convention==kConvDCS2008 && l3Pol==diPol)) ) { 
    AliErrorGeneral("AliMagF",Form("Wrong combination for L3/Dipole polarities (%c/%c) for convention %d",
				   l3Pol>0?'+':'-',diPol>0?'+':'-',GetPolarityConvention()));
    return 0;
  }
  //
  if (l3Pol<0) sclL3  = -sclL3;
  if (diPol<0) sclDip = -sclDip;
  //
  BeamType_t btype = kNoBeamField;
  TString btypestr = beamtype;
  btypestr.ToLower();
  TPRegexp protonBeam("(proton|p)\\s*-?\\s*\\1");
  TPRegexp ionBeam("(lead|pb|ion|a)\\s*-?\\s*\\1");
  if (btypestr.Contains(ionBeam)) btype = kBeamTypeAA;
  else if (btypestr.Contains(protonBeam)) btype = kBeamTypepp;
  else AliInfoGeneral("AliMagF",Form("Assume no LHC magnet field for the beam type %s, ",beamtype));
  char ttl[80];
  sprintf(ttl,"L3: %+5d Dip: %+4d kA; %s | Polarities in %s convention",(int)TMath::Sign(l3Cur,float(sclL3)),
	  (int)TMath::Sign(diCur,float(sclDip)),uniform ? " Constant":"",
	  convention==kConvLHC ? "LHC":"DCS2008");
  // LHC and DCS08 conventions have opposite dipole polarities
  if ( GetPolarityConvention() != convention) sclDip = -sclDip;
  //
  return new AliMagF("MagneticFieldMap", ttl,sclL3,sclDip,map,btype,beamenergy,2,10.,path);
  //
}

//_____________________________________________________________________________
const char*  AliMagF::GetBeamTypeText() const
{
  const char *beamNA  = "No Beam";
  const char *beamPP  = "p-p";
  const char *beamPbPb= "Pb-Pb";
  switch ( fBeamType ) {
  case kBeamTypepp : return beamPP;
  case kBeamTypeAA : return beamPbPb;
  case kNoBeamField: 
  default:           return beamNA;
  }
}
Commit	Line	Data
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
fe4da5cc	22	#include "AliMagF.h"
db83d72f	23	#include "AliMagWrapCheb.h"
db83d72f	24	#include "AliLog.h"
972ca52f	25
fe4da5cc	26	ClassImp(AliMagF)
fe4da5cc	27
9251fceb	28	const Double_t AliMagF::fgkSol2DipZ = -700.;
1dd3d90e	29	const UShort_t AliMagF::fgkPolarityConvention = kConvLHC;
db83d72f	30
1dd3d90e	31	/*
	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 -> negative Bz
	42	positive Dip current -> negative Bx
	43
	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.
	47	*/
e2afb3b6	48	//_______________________________________________________________________
e2afb3b6	49	AliMagF::AliMagF():
db83d72f	50	TVirtualMagField(),
	51	fMeasuredMap(0),
	52	fMapType(k5kG),
	53	fSolenoid(0),
	54	fBeamType(kNoBeamField),
	55	fBeamEnergy(0),
db83d72f	56	//
e2afb3b6	57	fInteg(0),
db83d72f	58	fPrecInteg(0),
	59	fFactorSol(1.),
	60	fFactorDip(1.),
	61	fMax(15),
	62	fDipoleOFF(kFALSE),
e2afb3b6	63	//
db83d72f	64	fQuadGradient(0),
	65	fDipoleField(0),
	66	fCCorrField(0),
	67	fACorr1Field(0),
	68	fACorr2Field(0),
	69	fParNames("","")
	70	{
e2afb3b6	71	// Default constructor
	72	//
	73	}
	74
	75	//_______________________________________________________________________
4642ac4b	76	AliMagF::AliMagF(const char name, const char title, Double_t factorSol, Double_t factorDip,
	77	BMap_t maptype, BeamType_t btype, Double_t benergy, Int_t integ, Double_t fmax,
	78	const char* path):
db83d72f	79	TVirtualMagField(name),
	80	fMeasuredMap(0),
	81	fMapType(maptype),
	82	fSolenoid(0),
4642ac4b	83	fBeamType(btype),
4642ac4b	84	fBeamEnergy(benergy),
db83d72f	85	//
db83d72f	86	fInteg(integ),
604e0531	87	fPrecInteg(1),
db83d72f	88	fFactorSol(1.),
db83d72f	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	//
4642ac4b	111	if (fBeamEnergy<=0 && fBeamType!=kNoBeamField) {
	112	if (fBeamType == kBeamTypepp) fBeamEnergy = 7000.; // max proton energy
	113	else if (fBeamType == kBeamTypeAA) fBeamEnergy = 5500; // max PbPb energy
	114	AliInfo("Maximim possible beam energy for requested beam is assumed");
	115	}
db83d72f	116	const char* parname = 0;
db83d72f	117	//
f04e7f5f	118	if (fMapType == k2kG) parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
	119	else if (fMapType == k5kG) parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
	120	else if (fMapType == k5kGUniform) parname = "Sol30_Dip6_Uniform";
	121	else AliFatal(Form("Unknown field identifier %d is requested\n",fMapType));
db83d72f	122	//
	123	SetDataFileName(path);
	124	SetParamName(parname);
	125	//
db83d72f	126	LoadParameterization();
db83d72f	127	InitMachineField(fBeamType,fBeamEnergy);
f04e7f5f	128	double xyz[3]={0.,0.,0.};
	129	fSolenoid = GetBz(xyz);
	130	SetFactorSol(factorSol);
	131	SetFactorDip(factorDip);
e86708b3	132	AliInfo(Form("Alice B fields: Solenoid (%+.2f*)%.0f kG, Dipole %s (%+.2f) %s",
e86708b3	133	factorSol,(fMapType==k5kG\|\|fMapType==k5kGUniform)?5.:2.,
439b5096	134	fDipoleOFF ? "OFF":"ON",factorDip,fMapType==k5kGUniform?" \|Constant Field!":""));
439b5096	135	AliInfo(Form("Machine B fields for %s beam (%.0f GeV): QGrad: %.4f Dipole: %.4f",
4642ac4b	136	fBeamType==kBeamTypeAA ? "A-A":(fBeamType==kBeamTypepp ? "p-p":"OFF"),
4642ac4b	137	fBeamEnergy,fQuadGradient,fDipoleField));
fe4da5cc	138	}
fe4da5cc	139
eeda4611	140	//_______________________________________________________________________
eeda4611	141	AliMagF::AliMagF(const AliMagF &src):
db83d72f	142	TVirtualMagField(src),
	143	fMeasuredMap(0),
	144	fMapType(src.fMapType),
	145	fSolenoid(src.fSolenoid),
	146	fBeamType(src.fBeamType),
	147	fBeamEnergy(src.fBeamEnergy),
eeda4611	148	fInteg(src.fInteg),
eeda4611	149	fPrecInteg(src.fPrecInteg),
db83d72f	150	fFactorSol(src.fFactorSol),
db83d72f	151	fFactorDip(src.fFactorDip),
eeda4611	152	fMax(src.fMax),
db83d72f	153	fDipoleOFF(src.fDipoleOFF),
	154	fQuadGradient(src.fQuadGradient),
	155	fDipoleField(src.fDipoleField),
	156	fCCorrField(src.fCCorrField),
	157	fACorr1Field(src.fACorr1Field),
	158	fACorr2Field(src.fACorr2Field),
	159	fParNames(src.fParNames)
eeda4611	160	{
db83d72f	161	if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
eeda4611	162	}
eeda4611	163
e2afb3b6	164	//_______________________________________________________________________
db83d72f	165	AliMagF::~AliMagF()
ff66b122	166	{
db83d72f	167	delete fMeasuredMap;
	168	}
	169
	170	//_______________________________________________________________________
	171	Bool_t AliMagF::LoadParameterization()
	172	{
	173	if (fMeasuredMap) {
	174	AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
	175	return kTRUE;
	176	}
ff66b122	177	//
db83d72f	178	char* fname = gSystem->ExpandPathName(GetDataFileName());
	179	TFile* file = TFile::Open(fname);
	180	if (!file) {
	181	AliError(Form("Failed to open magnetic field data file %s\n",fname));
	182	return kFALSE;
	183	}
ff66b122	184	//
db83d72f	185	fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
	186	if (!fMeasuredMap) {
	187	AliError(Form("Did not find field %s in %s\n",GetParamName(),fname));
	188	return kFALSE;
	189	}
	190	file->Close();
	191	delete file;
	192	return kTRUE;
ff66b122	193	}
ff66b122	194
db83d72f	195
ff66b122	196	//_______________________________________________________________________
db83d72f	197	void AliMagF::Field(const Double_t xyz, Double_t b)
fe4da5cc	198	{
db83d72f	199	// Method to calculate the field at point xyz
aee8290b	200	//
9251fceb	201	// b[0]=b[1]=b[2]=0.0;
9251fceb	202	if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
db83d72f	203	fMeasuredMap->Field(xyz,b);
db83d72f	204	if (xyz[2]>fgkSol2DipZ \|\| fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
9251fceb	205	else for (int i=3;i--;) b[i] *= fFactorDip;
db83d72f	206	}
9251fceb	207	else MachineField(xyz, b);
aee8290b	208	//
fe4da5cc	209	}
eeda4611	210
eeda4611	211	//_______________________________________________________________________
db83d72f	212	Double_t AliMagF::GetBz(const Double_t *xyz) const
eeda4611	213	{
db83d72f	214	// Method to calculate the field at point xyz
db83d72f	215	//
9251fceb	216	if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
	217	double bz = fMeasuredMap->GetBz(xyz);
	218	return (xyz[2]>fgkSol2DipZ \|\| fDipoleOFF) ? bzfFactorSol : bzfFactorDip;
db83d72f	219	}
9251fceb	220	else return 0.;
eeda4611	221	}
	222
	223	//_______________________________________________________________________
db83d72f	224	AliMagF& AliMagF::operator=(const AliMagF& src)
eeda4611	225	{
db83d72f	226	if (this != &src && src.fMeasuredMap) {
	227	if (fMeasuredMap) delete fMeasuredMap;
	228	fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
	229	SetName(src.GetName());
	230	fSolenoid = src.fSolenoid;
	231	fBeamType = src.fBeamType;
	232	fBeamEnergy = src.fBeamEnergy;
db83d72f	233	fInteg = src.fInteg;
	234	fPrecInteg = src.fPrecInteg;
	235	fFactorSol = src.fFactorSol;
	236	fFactorDip = src.fFactorDip;
	237	fMax = src.fMax;
	238	fDipoleOFF = src.fDipoleOFF;
	239	fParNames = src.fParNames;
	240	}
	241	return *this;
eeda4611	242	}
	243
	244	//_______________________________________________________________________
db83d72f	245	void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
eeda4611	246	{
4642ac4b	247	if (btype==kNoBeamField) {
db83d72f	248	fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
9251fceb	249	return;
db83d72f	250	}
db83d72f	251	//
9251fceb	252	double rigScale = benergy/7000.; // scale according to ratio of E/Enominal
	253	// for ions assume PbPb (with energy provided per nucleon) and account for A/Z
	254	if (btype == kBeamTypeAA) rigScale *= 208./82.;
	255	//
	256	fQuadGradient = 22.0002*rigScale;
	257	fDipoleField = 37.8781*rigScale;
	258	//
	259	// SIDE C
	260	fCCorrField = -9.6980;
	261	// SIDE A
	262	fACorr1Field = -13.2247;
	263	fACorr2Field = 11.7905;
db83d72f	264	//
eeda4611	265	}
eed8a1a2	266
db83d72f	267	//_______________________________________________________________________
db83d72f	268	void AliMagF::MachineField(const Double_t x, Double_t b) const
eed8a1a2	269	{
db83d72f	270	// ---- This is the ZDC part
9251fceb	271	// Compansators for Alice Muon Arm Dipole
	272	const Double_t kBComp1CZ = 1075., kBComp1hDZ = 260./2., kBComp1SqR = 4.0*4.0;
	273	const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5;
	274	//
	275	const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
90ae20c9	276	const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
9251fceb	277	const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5;
9251fceb	278	const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5;
db83d72f	279	//
9251fceb	280	const Double_t kDip1CZ = 6310.8, kDip1hDZ = 945./2., kDip1SqRC = 4.54.5, kDip1SqRA = 3.3753.375;
	281	const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.54.5, kDip2SqRA = 3.753.75;
	282	const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4;
db83d72f	283	//
	284	double rad2 = x[0] * x[0] + x[1] * x[1];
	285	//
9251fceb	286	b[0] = b[1] = b[2] = 0;
9251fceb	287	//
db83d72f	288	// SIDE C **************************************************
db83d72f	289	if(x[2]<0.){
9251fceb	290	if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
9251fceb	291	b[0] = fCCorrField*fFactorDip;
db83d72f	292	}
9251fceb	293	else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
db83d72f	294	b[0] = fQuadGradient*x[1];
db83d72f	295	b[1] = fQuadGradient*x[0];
db83d72f	296	}
9251fceb	297	else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
db83d72f	298	b[0] = -fQuadGradient*x[1];
db83d72f	299	b[1] = -fQuadGradient*x[0];
db83d72f	300	}
9251fceb	301	else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
db83d72f	302	b[0] = -fQuadGradient*x[1];
db83d72f	303	b[1] = -fQuadGradient*x[0];
db83d72f	304	}
9251fceb	305	else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
db83d72f	306	b[0] = fQuadGradient*x[1];
db83d72f	307	b[1] = fQuadGradient*x[0];
db83d72f	308	}
9251fceb	309	else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
db83d72f	310	b[1] = fDipoleField;
db83d72f	311	}
9251fceb	312	else if(TMath::Abs(x[2]+kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRC) {
	313	double dxabs = TMath::Abs(x[0])-kDip2DXC;
	314	if ( (dxabsdxabs + x[1]x[1])<kDip2SqRC) {
db83d72f	315	b[1] = -fDipoleField;
db83d72f	316	}
	317	}
	318	}
	319	//
	320	// SIDE A **************************************************
	321	else{
9251fceb	322	if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
db83d72f	323	// Compensator magnet at z = 1075 m
9251fceb	324	b[0] = fACorr1Field*fFactorDip;
db83d72f	325	}
db83d72f	326	//
9251fceb	327	if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
	328	b[0] = fACorr2Field*fFactorDip;
	329	}
	330	else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
db83d72f	331	b[0] = -fQuadGradient*x[1];
db83d72f	332	b[1] = -fQuadGradient*x[0];
eed8a1a2	333	}
9251fceb	334	else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
	335	b[0] = fQuadGradient*x[1];
	336	b[1] = fQuadGradient*x[0];
eed8a1a2	337	}
9251fceb	338	else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
	339	b[0] = fQuadGradient*x[1];
	340	b[1] = fQuadGradient*x[0];
db83d72f	341	}
9251fceb	342	else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
db83d72f	343	b[0] = -fQuadGradient*x[1];
db83d72f	344	b[1] = -fQuadGradient*x[0];
db83d72f	345	}
9251fceb	346	else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
db83d72f	347	b[1] = -fDipoleField;
db83d72f	348	}
9251fceb	349	else if(TMath::Abs(x[2]-kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRA) {
	350	double dxabs = TMath::Abs(x[0])-kDip2DXA;
	351	if ( (dxabsdxabs + x[1]x[1])<kDip2SqRA) {
db83d72f	352	b[1] = fDipoleField;
	353	}
	354	}
	355	}
9251fceb	356	//
db83d72f	357	}
	358
	359	//_______________________________________________________________________
	360	void AliMagF::GetTPCInt(const Double_t xyz, Double_t b) const
	361	{
	362	// Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
	363	b[0]=b[1]=b[2]=0.0;
	364	if (fMeasuredMap) {
	365	fMeasuredMap->GetTPCInt(xyz,b);
	366	for (int i=3;i--;) b[i] *= fFactorSol;
	367	}
	368	}
	369
	370	//_______________________________________________________________________
	371	void AliMagF::GetTPCIntCyl(const Double_t rphiz, Double_t b) const
	372	{
	373	// Method to calculate the integral of magnetic integral from point to nearest cathode plane
	374	// in cylindrical coordiates ( -pi<phi<pi convention )
	375	b[0]=b[1]=b[2]=0.0;
	376	if (fMeasuredMap) {
	377	fMeasuredMap->GetTPCIntCyl(rphiz,b);
	378	for (int i=3;i--;) b[i] *= fFactorSol;
	379	}
eed8a1a2	380	}
1dd3d90e	381
	382	//_______________________________________________________________________
	383	void AliMagF::SetFactorSol(Float_t fc)
	384	{
	385	// set the sign/scale of the current in the L3 according to fgkPolarityConvention
	386	switch (fgkPolarityConvention) {
	387	case kConvDCS2008: fFactorSol = -fc; break;
	388	case kConvLHC : fFactorSol = -fc; break;
	389	default : fFactorSol = fc; break; // case kConvMap2005: fFactorSol = fc; break;
	390	}
	391	}
	392
	393	//_______________________________________________________________________
	394	void AliMagF::SetFactorDip(Float_t fc)
	395	{
	396	// set the sign*scale of the current in the Dipole according to fgkPolarityConvention
	397	switch (fgkPolarityConvention) {
	398	case kConvDCS2008: fFactorDip = fc; break;
	399	case kConvLHC : fFactorDip = -fc; break;
	400	default : fFactorDip = fc; break; // case kConvMap2005: fFactorDip = fc; break;
	401	}
	402	}
	403
	404	//_______________________________________________________________________
	405	Double_t AliMagF::GetFactorSol() const
	406	{
	407	// return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe
	408	switch (fgkPolarityConvention) {
	409	case kConvDCS2008: return -fFactorSol;
	410	case kConvLHC : return -fFactorSol;
	411	default : return fFactorSol; // case kConvMap2005: return fFactorSol;
	412	}
	413	}
	414
	415	//_______________________________________________________________________
	416	Double_t AliMagF::GetFactorDip() const
	417	{
	418	// return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe
	419	switch (fgkPolarityConvention) {
	420	case kConvDCS2008: return fFactorDip;
	421	case kConvLHC : return -fFactorDip;
	422	default : return fFactorDip; // case kConvMap2005: return fFactorDip;
	423	}
	424	}
33fe5eb1	425
	426	//_____________________________________________________________________________
	427	AliMagF* AliMagF::CreateFieldMap(Float_t l3Cur, Float_t diCur, Int_t convention, Bool_t uniform,
	428	Float_t beamenergy, const Char_t beamtype, const Char_t path)
	429	{
	430	//------------------------------------------------
	431	// The magnetic field map, defined externally...
	432	// L3 current 30000 A -> 0.5 T
	433	// L3 current 12000 A -> 0.2 T
	434	// dipole current 6000 A
	435	// The polarities must match the convention (LHC or DCS2008)
	436	// unless the special uniform map was used for MC
	437	//------------------------------------------------
	438	const Float_t l3NominalCurrent1=30000.; // (A)
	439	const Float_t l3NominalCurrent2=12000.; // (A)
	440	const Float_t diNominalCurrent =6000. ; // (A)
	441
	442	const Float_t tolerance=0.03; // relative current tolerance
	443	const Float_t zero=77.; // "zero" current (A)
	444	//
	445	BMap_t map;
	446	double sclL3,sclDip;
	447	//
	448	Float_t l3Pol = l3Cur > 0 ? 1:-1;
	449	Float_t diPol = diCur > 0 ? 1:-1;
	450
	451	l3Cur = TMath::Abs(l3Cur);
	452	diCur = TMath::Abs(diCur);
	453	//
	454	if (TMath::Abs((sclDip=diCur/diNominalCurrent)-1.) > tolerance && !uniform) {
	455	if (diCur <= zero) sclDip = 0.; // some small current.. -> Dipole OFF
	456	else {
	457	AliErrorGeneral("AliMagF",Form("Wrong dipole current (%f A)!",diCur));
	458	return 0;
	459	}
	460	}
	461	//
	462	if (uniform) {
	463	// special treatment of special MC with uniform mag field (normalized to 0.5 T)
	464	// no check for scaling/polarities are done
	465	map = k5kGUniform;
	466	sclL3 = l3Cur/l3NominalCurrent1;
	467	}
	468	else {
	469	if (TMath::Abs((sclL3=l3Cur/l3NominalCurrent1)-1.) < tolerance) map = k5kG;
	470	else if (TMath::Abs((sclL3=l3Cur/l3NominalCurrent2)-1.) < tolerance) map = k2kG;
	471	else if (l3Cur <= zero) { sclL3 = 0; map = k5kGUniform;}
	472	else {
	473	AliErrorGeneral("AliMagF",Form("Wrong L3 current (%f A)!",l3Cur));
	474	return 0;
	475	}
	476	}
	477	//
	478	if (sclDip!=0 && (map==k5kG \|\| map==k2kG) &&
	479	((convention==kConvLHC && l3Pol!=diPol) \|\|
	480	(convention==kConvDCS2008 && l3Pol==diPol)) ) {
	481	AliErrorGeneral("AliMagF",Form("Wrong combination for L3/Dipole polarities (%c/%c) for convention %d",
	482	l3Pol>0?'+':'-',diPol>0?'+':'-',GetPolarityConvention()));
	483	return 0;
	484	}
	485	//
	486	if (l3Pol<0) sclL3 = -sclL3;
	487	if (diPol<0) sclDip = -sclDip;
	488	//
489	BeamType_t btype = kNoBeamField;
490	TString btypestr = beamtype;
491	btypestr.ToLower();
492	TPRegexp protonBeam("(proton\|p)\\s-?\\s\\1");
493	TPRegexp ionBeam("(lead\|pb\|ion\|a)\\s-?\\s\\1");
494	if (btypestr.Contains(ionBeam)) btype = kBeamTypeAA;
495	else if (btypestr.Contains(protonBeam)) btype = kBeamTypepp;
496	else AliInfoGeneral("AliMagF",Form("Assume no LHC magnet field for the beam type %s, ",beamtype));
497	char ttl[80];
498	sprintf(ttl,"L3: %+5d Dip: %+4d kA; %s \| Polarities in %s convention",(int)TMath::Sign(l3Cur,float(sclL3)),
499	(int)TMath::Sign(diCur,float(sclDip)),uniform ? " Constant":"",
500	convention==kConvLHC ? "LHC":"DCS2008");
501	// LHC and DCS08 conventions have opposite dipole polarities
502	if ( GetPolarityConvention() != convention) sclDip = -sclDip;
503	//
4642ac4b	504	return new AliMagF("MagneticFieldMap", ttl,sclL3,sclDip,map,btype,beamenergy,2,10.,path);
33fe5eb1	505	//
	506	}
	507
	508	//_____________________________________________________________________________
	509	const char* AliMagF::GetBeamTypeText() const
	510	{
	511	const char *beamNA = "No Beam";
	512	const char *beamPP = "p-p";
	513	const char *beamPbPb= "Pb-Pb";
	514	switch ( fBeamType ) {
	515	case kBeamTypepp : return beamPP;
	516	case kBeamTypeAA : return beamPbPb;
	517	case kNoBeamField:
	518	default: return beamNA;
	519	}
	520	}
	521