[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 -> positive 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 bt, Double_t be,Int_t integ, Double_t fmax, const char* path):
  TVirtualMagField(name),
  fMeasuredMap(0),
  fMapType(maptype),
  fSolenoid(0),
  fBeamType(bt),
  fBeamEnergy(be),
  //
  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 = 2750;  // 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",
	       bt==kBeamTypeAA ? "A-A":(bt==kBeamTypepp ? "p-p":"OFF"),be,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 sqrts, 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,sqrts/2,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
e3eadfac	41	positive L3 current -> positive Bz
1dd3d90e	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,
02233f2b	77	BMap_t maptype, BeamType_t bt, Double_t be,Int_t integ, Double_t fmax, const char* path):
db83d72f	78	TVirtualMagField(name),
	79	fMeasuredMap(0),
	80	fMapType(maptype),
	81	fSolenoid(0),
02233f2b	82	fBeamType(bt),
02233f2b	83	fBeamEnergy(be),
db83d72f	84	//
db83d72f	85	fInteg(integ),
604e0531	86	fPrecInteg(1),
db83d72f	87	fFactorSol(1.),
db83d72f	88	fFactorDip(1.),
972ca52f	89	fMax(fmax),
db83d72f	90	fDipoleOFF(factorDip==0.),
	91	//
	92	fQuadGradient(0),
	93	fDipoleField(0),
	94	fCCorrField(0),
	95	fACorr1Field(0),
	96	fACorr2Field(0),
	97	fParNames("","")
fe4da5cc	98	{
9251fceb	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
aee8290b	102	//
db83d72f	103	SetTitle(title);
	104	if(integ<0 \|\| integ > 2) {
	105	AliWarning(Form("Invalid magnetic field flag: %5d; Helix tracking chosen instead",integ));
	106	fInteg = 2;
	107	}
	108	if (fInteg == 0) fPrecInteg = 0;
aee8290b	109	//
4642ac4b	110	if (fBeamEnergy<=0 && fBeamType!=kNoBeamField) {
4642ac4b	111	if (fBeamType == kBeamTypepp) fBeamEnergy = 7000.; // max proton energy
02233f2b	112	else if (fBeamType == kBeamTypeAA) fBeamEnergy = 2750; // max PbPb energy
4642ac4b	113	AliInfo("Maximim possible beam energy for requested beam is assumed");
4642ac4b	114	}
db83d72f	115	const char* parname = 0;
db83d72f	116	//
f04e7f5f	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));
db83d72f	121	//
	122	SetDataFileName(path);
	123	SetParamName(parname);
	124	//
db83d72f	125	LoadParameterization();
db83d72f	126	InitMachineField(fBeamType,fBeamEnergy);
f04e7f5f	127	double xyz[3]={0.,0.,0.};
	128	fSolenoid = GetBz(xyz);
	129	SetFactorSol(factorSol);
	130	SetFactorDip(factorDip);
e86708b3	131	AliInfo(Form("Alice B fields: Solenoid (%+.2f*)%.0f kG, Dipole %s (%+.2f) %s",
e86708b3	132	factorSol,(fMapType==k5kG\|\|fMapType==k5kGUniform)?5.:2.,
439b5096	133	fDipoleOFF ? "OFF":"ON",factorDip,fMapType==k5kGUniform?" \|Constant Field!":""));
439b5096	134	AliInfo(Form("Machine B fields for %s beam (%.0f GeV): QGrad: %.4f Dipole: %.4f",
02233f2b	135	bt==kBeamTypeAA ? "A-A":(bt==kBeamTypepp ? "p-p":"OFF"),be,fQuadGradient,fDipoleField));
fe4da5cc	136	}
fe4da5cc	137
eeda4611	138	//_______________________________________________________________________
eeda4611	139	AliMagF::AliMagF(const AliMagF &src):
db83d72f	140	TVirtualMagField(src),
	141	fMeasuredMap(0),
	142	fMapType(src.fMapType),
	143	fSolenoid(src.fSolenoid),
	144	fBeamType(src.fBeamType),
	145	fBeamEnergy(src.fBeamEnergy),
eeda4611	146	fInteg(src.fInteg),
eeda4611	147	fPrecInteg(src.fPrecInteg),
db83d72f	148	fFactorSol(src.fFactorSol),
db83d72f	149	fFactorDip(src.fFactorDip),
eeda4611	150	fMax(src.fMax),
db83d72f	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)
eeda4611	158	{
db83d72f	159	if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
eeda4611	160	}
eeda4611	161
e2afb3b6	162	//_______________________________________________________________________
db83d72f	163	AliMagF::~AliMagF()
ff66b122	164	{
db83d72f	165	delete fMeasuredMap;
	166	}
	167
	168	//_______________________________________________________________________
	169	Bool_t AliMagF::LoadParameterization()
	170	{
	171	if (fMeasuredMap) {
	172	AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
	173	return kTRUE;
	174	}
ff66b122	175	//
db83d72f	176	char* fname = gSystem->ExpandPathName(GetDataFileName());
	177	TFile* file = TFile::Open(fname);
	178	if (!file) {
	179	AliError(Form("Failed to open magnetic field data file %s\n",fname));
	180	return kFALSE;
	181	}
ff66b122	182	//
db83d72f	183	fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
	184	if (!fMeasuredMap) {
	185	AliError(Form("Did not find field %s in %s\n",GetParamName(),fname));
	186	return kFALSE;
	187	}
	188	file->Close();
	189	delete file;
	190	return kTRUE;
ff66b122	191	}
ff66b122	192
db83d72f	193
ff66b122	194	//_______________________________________________________________________
db83d72f	195	void AliMagF::Field(const Double_t xyz, Double_t b)
fe4da5cc	196	{
db83d72f	197	// Method to calculate the field at point xyz
aee8290b	198	//
9251fceb	199	// b[0]=b[1]=b[2]=0.0;
9251fceb	200	if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
db83d72f	201	fMeasuredMap->Field(xyz,b);
db83d72f	202	if (xyz[2]>fgkSol2DipZ \|\| fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
9251fceb	203	else for (int i=3;i--;) b[i] *= fFactorDip;
db83d72f	204	}
9251fceb	205	else MachineField(xyz, b);
aee8290b	206	//
fe4da5cc	207	}
eeda4611	208
eeda4611	209	//_______________________________________________________________________
db83d72f	210	Double_t AliMagF::GetBz(const Double_t *xyz) const
eeda4611	211	{
db83d72f	212	// Method to calculate the field at point xyz
db83d72f	213	//
9251fceb	214	if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
	215	double bz = fMeasuredMap->GetBz(xyz);
	216	return (xyz[2]>fgkSol2DipZ \|\| fDipoleOFF) ? bzfFactorSol : bzfFactorDip;
db83d72f	217	}
9251fceb	218	else return 0.;
eeda4611	219	}
	220
	221	//_______________________________________________________________________
db83d72f	222	AliMagF& AliMagF::operator=(const AliMagF& src)
eeda4611	223	{
db83d72f	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;
db83d72f	231	fInteg = src.fInteg;
	232	fPrecInteg = src.fPrecInteg;
	233	fFactorSol = src.fFactorSol;
	234	fFactorDip = src.fFactorDip;
	235	fMax = src.fMax;
	236	fDipoleOFF = src.fDipoleOFF;
	237	fParNames = src.fParNames;
	238	}
	239	return *this;
eeda4611	240	}
	241
	242	//_______________________________________________________________________
db83d72f	243	void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
eeda4611	244	{
4642ac4b	245	if (btype==kNoBeamField) {
db83d72f	246	fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
9251fceb	247	return;
db83d72f	248	}
db83d72f	249	//
9251fceb	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.;
	253	//
	254	fQuadGradient = 22.0002*rigScale;
	255	fDipoleField = 37.8781*rigScale;
	256	//
	257	// SIDE C
	258	fCCorrField = -9.6980;
	259	// SIDE A
	260	fACorr1Field = -13.2247;
	261	fACorr2Field = 11.7905;
db83d72f	262	//
eeda4611	263	}
eed8a1a2	264
db83d72f	265	//_______________________________________________________________________
db83d72f	266	void AliMagF::MachineField(const Double_t x, Double_t b) const
eed8a1a2	267	{
db83d72f	268	// ---- This is the ZDC part
9251fceb	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;
	272	//
	273	const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
90ae20c9	274	const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
9251fceb	275	const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5;
9251fceb	276	const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5;
db83d72f	277	//
9251fceb	278	const Double_t kDip1CZ = 6310.8, kDip1hDZ = 945./2., kDip1SqRC = 4.54.5, kDip1SqRA = 3.3753.375;
	279	const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.54.5, kDip2SqRA = 3.753.75;
	280	const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4;
db83d72f	281	//
	282	double rad2 = x[0] * x[0] + x[1] * x[1];
	283	//
9251fceb	284	b[0] = b[1] = b[2] = 0;
9251fceb	285	//
db83d72f	286	// SIDE C **************************************************
db83d72f	287	if(x[2]<0.){
9251fceb	288	if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
9251fceb	289	b[0] = fCCorrField*fFactorDip;
db83d72f	290	}
9251fceb	291	else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
db83d72f	292	b[0] = fQuadGradient*x[1];
db83d72f	293	b[1] = fQuadGradient*x[0];
db83d72f	294	}
9251fceb	295	else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
db83d72f	296	b[0] = -fQuadGradient*x[1];
db83d72f	297	b[1] = -fQuadGradient*x[0];
db83d72f	298	}
9251fceb	299	else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
db83d72f	300	b[0] = -fQuadGradient*x[1];
db83d72f	301	b[1] = -fQuadGradient*x[0];
db83d72f	302	}
9251fceb	303	else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
db83d72f	304	b[0] = fQuadGradient*x[1];
db83d72f	305	b[1] = fQuadGradient*x[0];
db83d72f	306	}
9251fceb	307	else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
db83d72f	308	b[1] = fDipoleField;
db83d72f	309	}
9251fceb	310	else if(TMath::Abs(x[2]+kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRC) {
	311	double dxabs = TMath::Abs(x[0])-kDip2DXC;
	312	if ( (dxabsdxabs + x[1]x[1])<kDip2SqRC) {
db83d72f	313	b[1] = -fDipoleField;
db83d72f	314	}
	315	}
	316	}
	317	//
	318	// SIDE A **************************************************
	319	else{
9251fceb	320	if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
db83d72f	321	// Compensator magnet at z = 1075 m
9251fceb	322	b[0] = fACorr1Field*fFactorDip;
db83d72f	323	}
db83d72f	324	//
9251fceb	325	if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
	326	b[0] = fACorr2Field*fFactorDip;
	327	}
	328	else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
db83d72f	329	b[0] = -fQuadGradient*x[1];
db83d72f	330	b[1] = -fQuadGradient*x[0];
eed8a1a2	331	}
9251fceb	332	else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
	333	b[0] = fQuadGradient*x[1];
	334	b[1] = fQuadGradient*x[0];
eed8a1a2	335	}
9251fceb	336	else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
	337	b[0] = fQuadGradient*x[1];
	338	b[1] = fQuadGradient*x[0];
db83d72f	339	}
9251fceb	340	else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
db83d72f	341	b[0] = -fQuadGradient*x[1];
db83d72f	342	b[1] = -fQuadGradient*x[0];
db83d72f	343	}
9251fceb	344	else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
db83d72f	345	b[1] = -fDipoleField;
db83d72f	346	}
9251fceb	347	else if(TMath::Abs(x[2]-kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRA) {
	348	double dxabs = TMath::Abs(x[0])-kDip2DXA;
	349	if ( (dxabsdxabs + x[1]x[1])<kDip2SqRA) {
db83d72f	350	b[1] = fDipoleField;
	351	}
	352	}
	353	}
9251fceb	354	//
db83d72f	355	}
	356
	357	//_______________________________________________________________________
	358	void AliMagF::GetTPCInt(const Double_t xyz, Double_t b) const
	359	{
	360	// Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
	361	b[0]=b[1]=b[2]=0.0;
	362	if (fMeasuredMap) {
	363	fMeasuredMap->GetTPCInt(xyz,b);
	364	for (int i=3;i--;) b[i] *= fFactorSol;
	365	}
	366	}
	367
	368	//_______________________________________________________________________
	369	void AliMagF::GetTPCIntCyl(const Double_t rphiz, Double_t b) const
	370	{
	371	// Method to calculate the integral of magnetic integral from point to nearest cathode plane
	372	// in cylindrical coordiates ( -pi<phi<pi convention )
	373	b[0]=b[1]=b[2]=0.0;
	374	if (fMeasuredMap) {
	375	fMeasuredMap->GetTPCIntCyl(rphiz,b);
	376	for (int i=3;i--;) b[i] *= fFactorSol;
	377	}
eed8a1a2	378	}
1dd3d90e	379
	380	//_______________________________________________________________________
	381	void AliMagF::SetFactorSol(Float_t fc)
	382	{
	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;
	388	}
	389	}
	390
	391	//_______________________________________________________________________
	392	void AliMagF::SetFactorDip(Float_t fc)
	393	{
	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;
	399	}
	400	}
	401
	402	//_______________________________________________________________________
	403	Double_t AliMagF::GetFactorSol() const
	404	{
	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;
	410	}
	411	}
	412
	413	//_______________________________________________________________________
	414	Double_t AliMagF::GetFactorDip() const
	415	{
	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;
	421	}
	422	}
33fe5eb1	423
	424	//_____________________________________________________________________________
	425	AliMagF* AliMagF::CreateFieldMap(Float_t l3Cur, Float_t diCur, Int_t convention, Bool_t uniform,
02233f2b	426	Float_t sqrts, const Char_t beamtype, const Char_t path)
33fe5eb1	427	{
	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)
	439
	440	const Float_t tolerance=0.03; // relative current tolerance
	441	const Float_t zero=77.; // "zero" current (A)
	442	//
	443	BMap_t map;
	444	double sclL3,sclDip;
	445	//
	446	Float_t l3Pol = l3Cur > 0 ? 1:-1;
	447	Float_t diPol = diCur > 0 ? 1:-1;
	448
	449	l3Cur = TMath::Abs(l3Cur);
	450	diCur = TMath::Abs(diCur);
	451	//
	452	if (TMath::Abs((sclDip=diCur/diNominalCurrent)-1.) > tolerance && !uniform) {
	453	if (diCur <= zero) sclDip = 0.; // some small current.. -> Dipole OFF
	454	else {
	455	AliErrorGeneral("AliMagF",Form("Wrong dipole current (%f A)!",diCur));
	456	return 0;
	457	}
	458	}
	459	//
	460	if (uniform) {
	461	// special treatment of special MC with uniform mag field (normalized to 0.5 T)
	462	// no check for scaling/polarities are done
	463	map = k5kGUniform;
	464	sclL3 = l3Cur/l3NominalCurrent1;
	465	}
	466	else {
	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) { sclL3 = 0; map = k5kGUniform;}
	470	else {
	471	AliErrorGeneral("AliMagF",Form("Wrong L3 current (%f A)!",l3Cur));
	472	return 0;
	473	}
	474	}
	475	//
	476	if (sclDip!=0 && (map==k5kG \|\| map==k2kG) &&
	477	((convention==kConvLHC && l3Pol!=diPol) \|\|
	478	(convention==kConvDCS2008 && l3Pol==diPol)) ) {
	479	AliErrorGeneral("AliMagF",Form("Wrong combination for L3/Dipole polarities (%c/%c) for convention %d",
	480	l3Pol>0?'+':'-',diPol>0?'+':'-',GetPolarityConvention()));
	481	return 0;
	482	}
	483	//
	484	if (l3Pol<0) sclL3 = -sclL3;
	485	if (diPol<0) sclDip = -sclDip;
	486	//
	487	BeamType_t btype = kNoBeamField;
	488	TString btypestr = beamtype;
	489	btypestr.ToLower();
	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));
495	char ttl[80];
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;
501	//
02233f2b	502	return new AliMagF("MagneticFieldMap", ttl,sclL3,sclDip,map,btype,sqrts/2,2,10.,path);
33fe5eb1	503	//
	504	}
	505
	506	//_____________________________________________________________________________
	507	const char* AliMagF::GetBeamTypeText() const
	508	{
	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;
	515	case kNoBeamField:
	516	default: return beamNA;
	517	}
	518	}
	519