MagneticField.cxx

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
 
#include "Field/MagneticField.h"
#include <TFile.h>             // for TFile
#include <TPRegexp.h>          // for TPRegexp
#include <TString.h>           // for TString
#include <TSystem.h>           // for TSystem, gSystem
#include <fairlogger/Logger.h> // for FairLogger
#include "FairParamList.h"
#include "FairRun.h"
#include "FairRuntimeDb.h"
 
using namespace o2::field;
 
ClassImp(MagneticField);
 
const Double_t MagneticField::sSolenoidToDipoleZ = -700.;
 
const UShort_t MagneticField::sPolarityConvention = MagneticField::kConvLHC;
 
MagneticField::MagneticField()
  : FairField(),
    mMeasuredMap(nullptr),
    mFastField(nullptr),
    mMapType(MagFieldParam::k5kG),
    mSolenoid(0),
    mBeamType(MagFieldParam::kNoBeamField),
    mBeamEnergy(0),
    mDefaultIntegration(0),
    mPrecisionInteg(0),
    mMultipicativeFactorSolenoid(1.),
    mMultipicativeFactorDipole(1.),
    mMaxField(15),
    mDipoleOnOffFlag(kFALSE),
    mQuadrupoleGradient(0),
    mDipoleField(0),
    mCompensatorField2C(0),
    mCompensatorField1A(0),
    mCompensatorField2A(0),
    mParameterNames("", "")
{
  /*
   * Default constructor
   */
  fType = 2; // flag non-constant field
}
 
MagneticField::MagneticField(const char* name, const char* title, Double_t factorSol, Double_t factorDip,
                             MagFieldParam::BMap_t maptype, MagFieldParam::BeamType_t bt, Double_t be, Int_t integ,
                             Double_t fmax, const std::string path)
  : FairField(name, title),
    mMeasuredMap(nullptr),
    mFastField(nullptr),
    mMapType(maptype),
    mSolenoid(0),
    mBeamType(bt),
    mBeamEnergy(be),
    mDefaultIntegration(integ),
    mPrecisionInteg(1),
    mMultipicativeFactorSolenoid(factorSol),
    mMultipicativeFactorDipole(factorDip),
    mMaxField(fmax),
    mDipoleOnOffFlag(factorDip == 0.),
    mQuadrupoleGradient(0),
    mDipoleField(0),
    mCompensatorField2C(0),
    mCompensatorField1A(0),
    mCompensatorField2A(0),
    mParameterNames("", "")
{
  /*
   * Constructor for human readable params
   */
 
  setDataFileName(path.c_str());
  CreateField();
}
 
MagneticField::MagneticField(const MagFieldParam& param)
  : FairField(param.GetName(), param.GetTitle()),
    mMeasuredMap(nullptr),
    mFastField(nullptr),
    mMapType(param.GetMapType()),
    mSolenoid(0),
    mBeamType(param.GetBeamType()),
    mBeamEnergy(param.GetBeamEnergy()),
    mDefaultIntegration(param.GetDefInt()),
    mPrecisionInteg(1),
    mMultipicativeFactorSolenoid(param.GetFactorSol()), // temporary
    mMultipicativeFactorDipole(param.GetFactorDip()),   // temporary
    mMaxField(param.GetMaxField()),
    mDipoleOnOffFlag(param.GetFactorDip() == 0.),
    mQuadrupoleGradient(0),
    mDipoleField(0),
    mCompensatorField2C(0),
    mCompensatorField1A(0),
    mCompensatorField2A(0),
    mParameterNames("", "")
{
  /*
   * Constructor for FairParam derived params
   */
 
  setDataFileName(param.GetMapPath());
  CreateField();
}
 
MagneticField* MagneticField::createNominalField(int fld, bool uniform)
{
  float fldCoeffL3, fldCoeffDip;
  o2::field::MagFieldParam::BMap_t fldType;
  if (uniform) {
    fldCoeffL3 = float(fld) / 5.;
    fldCoeffDip = fld > 0 ? 1. : -1;
    fldType = o2::field::MagFieldParam::k5kGUniform;
  } else {
    switch (std::abs(fld)) {
      case 5:
        fldType = o2::field::MagFieldParam::k5kG;
        fldCoeffL3 = fldCoeffDip = fld > 0 ? 1. : -1;
        break;
      case 0:
        fldType = o2::field::MagFieldParam::k5kG;
        fldCoeffL3 = fldCoeffDip = 0;
        break;
      case 2:
        fldType = o2::field::MagFieldParam::k2kG;
        fldCoeffL3 = fldCoeffDip = fld > 0 ? 1. : -1;
        break;
      default:
        LOG(fatal) << "Field option " << fld << " is not supported, use +-2, +-5 or 0 or <int_kilogauss>U";
    };
  }
  return new o2::field::MagneticField("Maps", "Maps", fldCoeffL3, fldCoeffDip, fldType);
}
 
void MagneticField::CreateField()
{
  /*
   * field initialization
   */
 
  fType = 2; // flag non-constant field
 
  // does real creation of the field
  if (mDefaultIntegration < 0 || mDefaultIntegration > 2) {
    LOG(warning) << "MagneticField::CreateField: Invalid magnetic field flag: " << mDefaultIntegration
                 << "; Helix tracking chosen instead";
    mDefaultIntegration = 2;
  }
  if (mDefaultIntegration == 0) {
    mPrecisionInteg = 0;
  }
 
  if (mBeamEnergy <= 0 && mBeamType != MagFieldParam::kNoBeamField) {
    if (mBeamType == MagFieldParam::kBeamTypepp) {
      mBeamEnergy = 7000.; // max proton energy
    } else if (mBeamType == MagFieldParam::kBeamTypeAA) {
      mBeamEnergy = 2760; // max PbPb energy
    } else if (mBeamType == MagFieldParam::kBeamTypepA || mBeamType == MagFieldParam::kBeamTypeAp) {
      mBeamEnergy = 2760; // same rigitiy max PbPb energy
    }
    //
    LOG(info) << "MagneticField::CreateField: Maximim possible beam energy for requested beam is assumed";
  }
 
  const char* parname = nullptr;
 
  if (mMapType == MagFieldParam::k2kG) {
    parname = mDipoleOnOffFlag ? "Sol12_Dip0_Hole" : "Sol12_Dip6_Hole";
  } else if (mMapType == MagFieldParam::k5kG) {
    parname = mDipoleOnOffFlag ? "Sol30_Dip0_Hole" : "Sol30_Dip6_Hole";
  } else if (mMapType == MagFieldParam::k5kGUniform) {
    parname = "Sol30_Dip6_Uniform";
  } else {
    LOG(fatal) << "MagneticField::CreateField: Unknown field identifier " << mMapType << " is requested\n";
  }
 
  setParameterName(parname);
 
  loadParameterization();
  initializeMachineField(mBeamType, mBeamEnergy);
  setFactorSolenoid(mMultipicativeFactorSolenoid);
  setFactorDipole(mMultipicativeFactorDipole);
  double xyz[3] = {0., 0., 0.};
  mSolenoid = getBz(xyz);
  Print("a");
  //
}
 
Bool_t MagneticField::loadParameterization()
{
  /*
   * load parametrization for measured field
   */
 
  if (mMeasuredMap) {
    LOG(fatal) << "MagneticField::loadParameterization: Field data " << getParameterName()
               << " are already loaded from " << getDataFileName() << "\n";
  }
  TString fname = getDataFileName();
  gSystem->ExpandPathName(fname);
  TFile* file = TFile::Open(fname);
  if (!file) {
    LOG(fatal) << "MagneticField::loadParameterization: Failed to open magnetic field data file " << fname << "\n";
  }
 
  mMeasuredMap =
    std::unique_ptr<MagneticWrapperChebyshev>(dynamic_cast<MagneticWrapperChebyshev*>(file->Get(getParameterName())));
  if (!mMeasuredMap) {
    LOG(fatal) << "MagneticField::loadParameterization: Did not find field " << getParameterName() << " in " << fname
               << "%s\n";
  }
  file->Close();
  delete file;
  return kTRUE;
}
 
void MagneticField::Field(const Double_t* __restrict__ xyz, Double_t* __restrict__ b)
{
  /*
   * query field value at point
   */
 
  //  b[0]=b[1]=b[2]=0.0;
  if (mFastField && mFastField->Field(xyz, b)) {
    return;
  }
 
  if (mMeasuredMap && xyz[2] > mMeasuredMap->getMinZ() && xyz[2] < mMeasuredMap->getMaxZ()) {
    mMeasuredMap->Field(xyz, b);
    if (xyz[2] > sSolenoidToDipoleZ || mDipoleOnOffFlag) {
      for (int i = 3; i--;) {
        b[i] *= mMultipicativeFactorSolenoid;
      }
    } else {
      for (int i = 3; i--;) {
        b[i] *= mMultipicativeFactorDipole;
      }
    }
  } else {
    MachineField(xyz, b);
  }
}
 
Double_t MagneticField::getBz(const Double_t* xyz) const
{
  /*
   * query field Bz component at point
   */
 
  if (mFastField) {
    double bz = 0;
    if (mFastField->GetBz(xyz, bz)) {
      return bz;
    }
  }
  if (mMeasuredMap && xyz[2] > mMeasuredMap->getMinZ() && xyz[2] < mMeasuredMap->getMaxZ()) {
    double bz = mMeasuredMap->getBz(xyz);
    return (xyz[2] > sSolenoidToDipoleZ || mDipoleOnOffFlag) ? bz * mMultipicativeFactorSolenoid
                                                             : bz * mMultipicativeFactorDipole;
  } else {
    return 0.;
  }
}
 
MagneticField& MagneticField::operator=(const MagneticField& src)
{
  /*
   * assignment operator
   */
 
  if (this != &src) {
    if (src.mMeasuredMap) {
      mMeasuredMap.reset(new MagneticWrapperChebyshev(*src.getMeasuredMap()));
    }
    SetName(src.GetName());
    mSolenoid = src.mSolenoid;
    mBeamType = src.mBeamType;
    mBeamEnergy = src.mBeamEnergy;
    mDefaultIntegration = src.mDefaultIntegration;
    mPrecisionInteg = src.mPrecisionInteg;
    mMultipicativeFactorSolenoid = src.mMultipicativeFactorSolenoid;
    mMultipicativeFactorDipole = src.mMultipicativeFactorDipole;
    mMaxField = src.mMaxField;
    mDipoleOnOffFlag = src.mDipoleOnOffFlag;
    mParameterNames = src.mParameterNames;
    mFastField.reset(src.mFastField ? new MagFieldFast(*src.getFastField()) : nullptr);
  }
  return *this;
}
 
void MagneticField::initializeMachineField(MagFieldParam::BeamType_t btype, Double_t benergy)
{
  if (btype == MagFieldParam::kNoBeamField) {
    mQuadrupoleGradient = mDipoleField = mCompensatorField2C = mCompensatorField1A = mCompensatorField2A = 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 == MagFieldParam::kBeamTypeAA /* || btype==kBeamTypepA || btype==kBeamTypeAp */) {
    rigScale *= 208. / 82.;
  }
  // Attention: in p-Pb the energy recorded in the GRP is the PROTON energy, no rigidity
  // rescaling is needed
 
  mQuadrupoleGradient = 22.0002 * rigScale;
  mDipoleField = 37.8781 * rigScale;
 
  // SIDE C
  mCompensatorField2C = -9.6980;
  // SIDE A
  mCompensatorField1A = -13.2247;
  mCompensatorField2A = 11.7905;
}
 
void MagneticField::MachineField(const Double_t* __restrict__ x, Double_t* __restrict__ 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] = mCompensatorField2C * mMultipicativeFactorDipole;
    } else if (TMath::Abs(x[2] + kTripQ1CZ) < kTripQ1hDZ && rad2 < kTripQ1SqR) {
      b[0] = mQuadrupoleGradient * x[1];
      b[1] = mQuadrupoleGradient * x[0];
    } else if (TMath::Abs(x[2] + kTripQ2CZ) < kTripQ2hDZ && rad2 < kTripQ2SqR) {
      b[0] = -mQuadrupoleGradient * x[1];
      b[1] = -mQuadrupoleGradient * x[0];
    } else if (TMath::Abs(x[2] + kTripQ3CZ) < kTripQ3hDZ && rad2 < kTripQ3SqR) {
      b[0] = -mQuadrupoleGradient * x[1];
      b[1] = -mQuadrupoleGradient * x[0];
    } else if (TMath::Abs(x[2] + kTripQ4CZ) < kTripQ4hDZ && rad2 < kTripQ4SqR) {
      b[0] = mQuadrupoleGradient * x[1];
      b[1] = mQuadrupoleGradient * x[0];
    } else if (TMath::Abs(x[2] + kDip1CZ) < kDip1hDZ && rad2 < kDip1SqRC) {
      b[1] = mDipoleField;
    } 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] = -mDipoleField;
      }
    }
  }
 
  // SIDE A
  else {
    if (TMath::Abs(x[2] - kBComp1CZ) < kBComp1hDZ && rad2 < kBComp1SqR) {
      // Compensator magnet at z = 1075 m
      b[0] = mCompensatorField1A * mMultipicativeFactorDipole;
    }
 
    if (TMath::Abs(x[2] - kBComp2CZ) < kBComp2hDZ && rad2 < kBComp2SqR) {
      b[0] = mCompensatorField2A * mMultipicativeFactorDipole;
    } else if (TMath::Abs(x[2] - kTripQ1CZ) < kTripQ1hDZ && rad2 < kTripQ1SqR) {
      b[0] = -mQuadrupoleGradient * x[1];
      b[1] = -mQuadrupoleGradient * x[0];
    } else if (TMath::Abs(x[2] - kTripQ2CZ) < kTripQ2hDZ && rad2 < kTripQ2SqR) {
      b[0] = mQuadrupoleGradient * x[1];
      b[1] = mQuadrupoleGradient * x[0];
    } else if (TMath::Abs(x[2] - kTripQ3CZ) < kTripQ3hDZ && rad2 < kTripQ3SqR) {
      b[0] = mQuadrupoleGradient * x[1];
      b[1] = mQuadrupoleGradient * x[0];
    } else if (TMath::Abs(x[2] - kTripQ4CZ) < kTripQ4hDZ && rad2 < kTripQ4SqR) {
      b[0] = -mQuadrupoleGradient * x[1];
      b[1] = -mQuadrupoleGradient * x[0];
    } else if (TMath::Abs(x[2] - kDip1CZ) < kDip1hDZ && rad2 < kDip1SqRA) {
      b[1] = -mDipoleField;
    } 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] = mDipoleField;
      }
    }
  }
}
 
void MagneticField::getTPCIntegral(const Double_t* xyz, Double_t* b) const
{
  b[0] = b[1] = b[2] = 0.0;
  if (mMeasuredMap) {
    mMeasuredMap->getTPCIntegral(xyz, b);
    for (int i = 3; i--;) {
      b[i] *= mMultipicativeFactorSolenoid;
    }
  }
}
 
void MagneticField::getTPCRatIntegral(const Double_t* xyz, Double_t* b) const
{
  b[0] = b[1] = b[2] = 0.0;
  if (mMeasuredMap) {
    mMeasuredMap->getTPCRatIntegral(xyz, b);
    b[2] /= 100;
  }
}
 
void MagneticField::getTPCIntegralCylindrical(const Double_t* rphiz, Double_t* b) const
{
  b[0] = b[1] = b[2] = 0.0;
  if (mMeasuredMap) {
    mMeasuredMap->getTPCIntegralCylindrical(rphiz, b);
    for (int i = 3; i--;) {
      b[i] *= mMultipicativeFactorSolenoid;
    }
  }
}
 
void MagneticField::getTPCRatIntegralCylindrical(const Double_t* rphiz, Double_t* b) const
{
  b[0] = b[1] = b[2] = 0.0;
  if (mMeasuredMap) {
    mMeasuredMap->getTPCRatIntegralCylindrical(rphiz, b);
    b[2] /= 100;
  }
}
 
void MagneticField::setFactorSolenoid(Float_t fc)
{
  switch (sPolarityConvention) {
    case kConvDCS2008:
      mMultipicativeFactorSolenoid = -fc;
      break;
    case kConvLHC:
      mMultipicativeFactorSolenoid = -fc;
      break;
    default:
      mMultipicativeFactorSolenoid = fc;
      break; // case kConvMap2005: mMultipicativeFactorSolenoid =  fc; break;
  }
  if (mFastField) {
    mFastField->setFactorSol(getFactorSolenoid());
  }
}
 
void MagneticField::setFactorDipole(Float_t fc)
{
  switch (sPolarityConvention) {
    case kConvDCS2008:
      mMultipicativeFactorDipole = fc;
      break;
    case kConvLHC:
      mMultipicativeFactorDipole = -fc;
      break;
    default:
      mMultipicativeFactorDipole = fc;
      break; // case kConvMap2005: mMultipicativeFactorDipole =  fc; break;
  }
}
 
Double_t MagneticField::getFactorSolenoid() const
{
  switch (sPolarityConvention) {
    case kConvDCS2008:
      return -mMultipicativeFactorSolenoid;
    case kConvLHC:
      return -mMultipicativeFactorSolenoid;
    default:
      return mMultipicativeFactorSolenoid; //  case kConvMap2005: return  mMultipicativeFactorSolenoid;
  }
}
 
Double_t MagneticField::getFactorDipole() const
{
  switch (sPolarityConvention) {
    case kConvDCS2008:
      return mMultipicativeFactorDipole;
    case kConvLHC:
      return -mMultipicativeFactorDipole;
    default:
      return mMultipicativeFactorDipole; //  case kConvMap2005: return  mMultipicativeFactorDipole;
  }
}
 
void MagneticField::rescaleField(float l3Cur, float diCur, bool uniform, int convention)
{
  // this function taks as input magnet currents and rescales existing field if the map is compatible
  float sclL3 = l3Cur, sclDip = diCur;
  MagFieldParam::BMap_t map = getFieldMapScale(sclL3, sclDip, uniform);
  if (map != mMapType) {
    LOGP(fatal, "Provided L3current={} DipCurrent={} uniform={} leads to map type {}, incompatible with loaded {}", l3Cur, diCur, uniform, (int)map, (int)mMapType);
  }
  setFactorSolenoid(sclL3);
  setFactorDipole(sclDip);
  LOGP(info, "Updating magnetic field: L3current={} DipCurrent={} uniform={}", l3Cur, diCur, uniform);
}
 
MagFieldParam::BMap_t MagneticField::getFieldMapScale(float& l3, float& dip, bool uniform, int convention)
{
  // this function taks as input magnet currents and returns the field type and scalings for L3 and dipole
  const float l3NominalCurrent1 = 30000.f; // (A)
  const float l3NominalCurrent2 = 12000.f; // (A)
  const float diNominalCurrent = 6000.f;   // (A)
 
  const float tolerance = 0.03; // relative current tolerance
  const float zero = 77.f;      // "zero" current (A)
 
  MagFieldParam::BMap_t map = MagFieldParam::k5kG;
  float sclL3, sclDip;
  float l3sav = l3, dipsav = dip;
 
  float l3Pol = l3 > 0 ? 1 : -1;
  float diPol = dip > 0 ? 1 : -1;
 
  l3 = TMath::Abs(l3);
  dip = TMath::Abs(dip);
  static bool overrideL3 = std::getenv("O2_OVERRIDE_L3_CURRENT") != nullptr;
  static bool overrideDIP = std::getenv("O2_OVERRIDE_DIPOLE_CURRENT") != nullptr;
  static bool warnL3Done = false, warnDipDone = false, warnPolarityDone = false;
 
  if (TMath::Abs((sclDip = dip / diNominalCurrent) - 1.) > tolerance && !uniform) {
    if (dip <= zero) {
      sclDip = 0.; // some small current.. -> Dipole OFF
    } else {
      if (!overrideDIP) {
        LOG(fatal) << "MagneticField::createFieldMap: Wrong dipole current (" << dipsav << " A)!";
      } else {
        if (!warnDipDone) {
          LOGP(error, "Dipole current was overridden to unsupported value {}", dipsav);
          warnDipDone = true;
        }
      }
    }
  }
  if (uniform) {
    // special treatment of special MC with uniform mag field (normalized to 0.5 T)
    // no check for scaling/polarities are done
    map = MagFieldParam::k5kGUniform;
    sclL3 = l3 / l3NominalCurrent1;
  } else {
    if (TMath::Abs((sclL3 = l3 / l3NominalCurrent1) - 1.) < tolerance) {
      map = MagFieldParam::k5kG;
    } else if (TMath::Abs((sclL3 = l3 / l3NominalCurrent2) - 1.) < tolerance) {
      map = MagFieldParam::k2kG;
    } else if (l3 <= zero && dip <= zero) {
      sclL3 = 0;
      sclDip = 0;
      map = MagFieldParam::k5kGUniform;
    } else {
      if (!overrideL3) {
        LOG(fatal) << "MagneticField::createFieldMap: Wrong L3 current (" << l3sav << "  A)!";
      } else {
        if (!warnL3Done) {
          LOGP(error, "L3 current was overridden to unsupported value {}", l3sav);
          warnL3Done = true;
        }
        map = MagFieldParam::k5kG;
        sclL3 = l3 / l3NominalCurrent1;
      }
    }
  }
  if (sclDip != 0 && map != MagFieldParam::k5kGUniform) {
    if ((l3 <= zero) ||
        ((convention == kConvLHC && l3Pol != diPol) || (convention == kConvDCS2008 && l3Pol == diPol))) {
      if (overrideL3 || overrideDIP) {
        if (!warnPolarityDone) {
          LOG(error) << "Overriden currents have wrong combination for L3/Dipole polarities ("
                     << (l3Pol > 0 ? '+' : '-') << "/" << (diPol > 0 ? '+' : '-') << ") for convention "
                     << getPolarityConvention();
          warnPolarityDone = true;
        }
      } else {
        LOG(fatal) << "MagneticField::createFieldMap: Wrong combination for L3/Dipole polarities ("
                   << (l3Pol > 0 ? '+' : '-') << "/" << (diPol > 0 ? '+' : '-') << ") for convention "
                   << getPolarityConvention();
      }
    }
  }
  l3 = (l3Pol < 0) ? -sclL3 : sclL3;
  dip = (diPol < 0) ? -sclDip : sclDip;
  return map;
}
 
MagneticField* MagneticField::createFieldMap(float l3Cur, float diCur, Int_t convention, Bool_t uniform,
                                             float beamenergy, const Char_t* beamtype, const std::string path)
{
  float sclL3 = l3Cur, sclDip = diCur;
  MagFieldParam::BMap_t map = getFieldMapScale(sclL3, sclDip, uniform);
  MagFieldParam::BeamType_t btype = MagFieldParam::kNoBeamField;
  TString btypestr = beamtype;
  btypestr.ToLower();
  TPRegexp protonBeam(R"((proton|p)\s*-?\s*\1)");
  TPRegexp ionBeam(R"((lead|pb|ion|a|A)\s*-?\s*\1)");
  TPRegexp protonionBeam(R"((proton|p)\s*-?\s*(lead|pb|ion|a|A))");
  TPRegexp ionprotonBeam(R"((lead|pb|ion|a|A)\s*-?\s*(proton|p))");
  if (btypestr.Contains(ionBeam)) {
    btype = MagFieldParam::kBeamTypeAA;
  } else if (btypestr.Contains(protonBeam)) {
    btype = MagFieldParam::kBeamTypepp;
  } else if (btypestr.Contains(protonionBeam)) {
    btype = MagFieldParam::kBeamTypepA;
  } else if (btypestr.Contains(ionprotonBeam)) {
    btype = MagFieldParam::kBeamTypeAp;
  } else {
    LOG(info) << "Assume no LHC magnet field for the beam type " << beamtype;
  }
  char ttl[80];
  snprintf(ttl, 79, "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 MagneticField("MagneticFieldMap", ttl, sclL3, sclDip, map, btype, beamenergy, 2, 10., path);
}
 
const char* MagneticField::getBeamTypeText() const
{
  const char* beamNA = "No Beam";
  const char* beamPP = "p-p";
  const char* beamPbPb = "A-A";
  const char* beamPPb = "p-A";
  const char* beamPbP = "A-p";
  switch (mBeamType) {
    case MagFieldParam::kBeamTypepp:
      return beamPP;
    case MagFieldParam::kBeamTypeAA:
      return beamPbPb;
    case MagFieldParam::kBeamTypepA:
      return beamPPb;
    case MagFieldParam::kBeamTypeAp:
      return beamPbP;
    case MagFieldParam::kNoBeamField:
    default:
      return beamNA;
  }
}
 
void MagneticField::Print(Option_t* opt) const
{
  TString opts = opt;
  opts.ToLower();
  LOG(info) << "MagneticField::Print: " << GetName() << ":" << GetTitle();
  LOG(info) << "MagneticField::Print: Solenoid (" << getFactorSolenoid() << "*)"
            << ((mMapType == MagFieldParam::k5kG || mMapType == MagFieldParam::k5kGUniform) ? 5. : 2) << " kG, Dipole "
            << (mDipoleOnOffFlag ? "OFF" : "ON") << " (" << getFactorDipole() << ") "
            << (mMapType == MagFieldParam::k5kGUniform ? " |Constant Field!" : "");
  if (opts.Contains("a")) {
    LOG(info) << "MagneticField::Print: Machine B fields for " << getBeamTypeText() << "  beam (" << mBeamEnergy
              << " GeV): QGrad: " << mQuadrupoleGradient << " Dipole: " << mDipoleField;
    LOG(info) << "MagneticField::Print: Uses " << getParameterName() << "  of " << getDataFileName();
  }
}
 
void MagneticField::FillParContainer()
{
  // fill field parameters
  FairRun* fRun = FairRun::Instance();
  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
  MagFieldParam* par = static_cast<MagFieldParam*>(rtdb->getContainer("MagFieldParam"));
  par->SetParam(this);
  par->setChanged();
}
 
//_____________________________________________________________________________
void MagneticField::AllowFastField(bool v)
{
  if (v) {
    if (!mFastField) {
      mFastField = std::make_unique<MagFieldFast>(getFactorSolenoid(), mMapType == MagFieldParam::k2kG ? 2 : 5);
    }
  } else {
    mFastField.reset(nullptr);
  }
}