Alice B-field derives from FairField + stores its parameters in MagFieldParam

Common/Field/CMakeLists.txt

@@ -4,12 +4,14 @@ O2_SETUP(NAME ${MODULE_NAME})
 
 set(SRCS
     src/MagneticWrapperChebyshev.cxx
+    src/MagFieldParam.cxx
     src/MagneticField.cxx
     )
 
 set(HEADERS
     include/${MODULE_NAME}/MagneticWrapperChebyshev.h
     include/${MODULE_NAME}/MagneticField.h
+    include/${MODULE_NAME}/MagFieldParam.h
     )
 set(LINKDEF src/FieldLinkDef.h)
 set(LIBRARY_NAME ${MODULE_NAME})

Common/Field/include/Field/MagFieldParam.h

@@ -0,0 +1,62 @@
+/// \file MagFieldParam.h
+/// \brief Definition of the MagFieldParam: container for ALICE mag. field parameters
+/// \author ruben.shahoyan@cern.ch
+
+#ifndef ALICEO2_FIELD_MAGFIELDPARAM_H_
+#define ALICEO2_FIELD_MAGFIELDPARAM_H_
+
+#include "FairParGenericSet.h"
+#include <TString.h>
+
+class FairParamList;
+
+
+namespace AliceO2 {
+  namespace Field {
+
+class MagneticField;
+
+class MagFieldParam : public FairParGenericSet
+{
+  public:
+    enum BMap_t
+    {
+        k2kG, k5kG, k5kGUniform
+    };
+    enum BeamType_t
+    {
+        kNoBeamField, kBeamTypepp, kBeamTypeAA, kBeamTypepA, kBeamTypeAp
+    };
+
+    MagFieldParam(const char* name="", const char* title="", const char* context="");
+
+    void SetParam(const MagneticField* field);
+
+    BMap_t     GetMapType()                   const {return mMapType;}
+    BeamType_t GetBeamType()                  const {return mBeamType;}
+    Int_t                     GetDefInt()     const {return mDefaultIntegration;}
+    Double_t                  GetFactorSol()  const {return mFactorSol;}
+    Double_t                  GetFactorDip()  const {return mFactorDip;}
+    Double_t                  GetBeamEnergy() const {return mBeamEnergy;}
+    Double_t                  GetMaxField()   const {return mMaxField;}
+    const char*               GetMapPath()    const {return mMapPath.Data();}
+
+    virtual void   putParams(FairParamList* list);
+    virtual Bool_t getParams(FairParamList* list);
+
+  protected:
+    BMap_t     mMapType;  ///< map type ID
+    BeamType_t mBeamType; ///< beam type ID
+    Int_t    mDefaultIntegration;        ///< field integration type for MC
+    Double_t mFactorSol;                 ///< solenoid current factor
+    Double_t mFactorDip;                 ///< dipole current factor
+    Double_t mBeamEnergy;                ///< beam energy
+    Double_t mMaxField;                  ///< max field for geant
+    TString  mMapPath;                   ///< path to map file
+
+    ClassDef(MagFieldParam,1)
+};
+}
+}
+
+#endif

Common/Field/include/Field/MagneticField.h

@@ -5,33 +5,30 @@
 #ifndef ALICEO2_FIELD_MAGNETICFIELD_H_
 #define ALICEO2_FIELD_MAGNETICFIELD_H_
 
-#include <TVirtualMagField.h>  // for TVirtualMagField
+#include "FairField.h"         // for FairField
+#include "Field/MagFieldParam.h"
+#include "Field/MagneticWrapperChebyshev.h" // for MagneticWrapperChebyshev
 #include "TSystem.h"
 #include "Rtypes.h"            // for Double_t, Char_t, Int_t, Float_t, etc
 #include "TNamed.h"            // for TNamed
+#include <memory>              // for str::unique_ptr
 
 class FairLogger;  // lines 14-14
+class FairParamList;
+
 namespace AliceO2 { namespace Field { class MagneticWrapperChebyshev; }}  // lines 19-19
 namespace AliceO2 {
 namespace Field {
 
 class MagneticWrapperChebyshev;
-
+ 
 /// Interface between the TVirtualMagField and MagneticWrapperChebyshev: wrapper to the set of magnetic field data +
 /// Tosca
 /// parametrization by Chebyshev polynomials
-class MagneticField : public TVirtualMagField
+class MagneticField : public FairField
 {
 
   public:
-    enum BMap_t
-    {
-        k2kG, k5kG, k5kGUniform
-    };
-    enum BeamType_t
-    {
-        kNoBeamField, kBeamTypepp, kBeamTypeAA, kBeamTypepA, kBeamTypeAp
-    };
     enum PolarityConvention_t
     {
         kConvLHC, kConvDCS2008, kConvMap2005
@@ -48,21 +45,55 @@ class MagneticField : public TVirtualMagField
     /// Impose scaling of parameterized L3 field by factorSol and of dipole by factorDip.
     /// The "be" is the energy of the beam in GeV/nucleon
     MagneticField(const char *name, const char *title, Double_t factorSol = 1., Double_t factorDip = 1.,
-                  BMap_t maptype = k5kG, BeamType_t btype = kBeamTypepp, Double_t benergy = -1, Int_t integ = 2,
+                  MagFieldParam::BMap_t maptype = MagFieldParam::k5kG,
+		  MagFieldParam::BeamType_t btype = MagFieldParam::kBeamTypepp,
+		  Double_t benergy = -1, Int_t integ = 2,
                   Double_t fmax = 15, const std::string path = std::string(gSystem->Getenv("VMCWORKDIR")) +
                                                                std::string("/Common/maps/mfchebKGI_sym.root")
     );
 
-    MagneticField(const MagneticField &src);
+
+    MagneticField(const MagFieldParam& param);
 
     MagneticField &operator=(const MagneticField &src);
 
     /// Default destructor
-    virtual ~MagneticField();
+    virtual ~MagneticField() {}
+
+    /// real field creation is here
+    void CreateField();
+
+    /// Virtual methods from FairField
+
+    /// X component, avoid using since slow
+    virtual Double_t GetBx(Double_t x, Double_t y, Double_t z) {
+      double xyz[3]={x,y,z},b[3];
+      GetFieldValue(xyz,b);
+      return b[0];
+    } 
+
+    /// Y component, avoid using since slow
+    virtual Double_t GetBy(Double_t x, Double_t y, Double_t z) {
+      double xyz[3]={x,y,z},b[3];
+      GetFieldValue(xyz,b);
+      return b[1];
+    }
+
+    /// Z component
+    virtual Double_t GetBz(Double_t x, Double_t y, Double_t z) {
+      double xyz[3]={x,y,z};
+      return getBz(xyz); 
+    } 
 
     /// Method to calculate the field at point xyz
-    virtual void Field(const Double_t *x, Double_t *b);
+    virtual void GetFieldValue(const Double_t point[3], Double_t* bField);
+
+    /// 3d field query alias for Alias Method to calculate the field at point xyz
+    virtual void GetBxyz(const Double_t p[3], Double_t* b) {Field(p,b);}
 
+    /// Fill Paramater
+    virtual void FillParContainer();
+
     /// Method to calculate the integral_0^z of br,bt,bz
     void getTPCIntegral(const Double_t *xyz, Double_t *b) const;
 
@@ -79,10 +110,7 @@ class MagneticField : public TVirtualMagField
     /// Method to calculate the field at point xyz
     Double_t getBz(const Double_t *xyz) const;
 
-    MagneticWrapperChebyshev *getMeasuredMap() const
-    {
-      return mMeasuredMap;
-    }
+    MagneticWrapperChebyshev *getMeasuredMap() const { return mMeasuredMap.get();}
 
     // Former MagF methods or their aliases
 
@@ -105,7 +133,7 @@ class MagneticField : public TVirtualMagField
 
     Double_t getCurrentSolenoid() const
     {
-      return getFactorSolenoid() * (mMapType == k2kG ? 12000 : 30000);
+      return getFactorSolenoid() * (mMapType == MagFieldParam::k2kG ? 12000 : 30000);
     }
 
     Double_t getCurrentDipole() const
@@ -115,17 +143,17 @@ class MagneticField : public TVirtualMagField
 
     Bool_t IsUniform() const
     {
-      return mMapType == k5kGUniform;
+      return mMapType == MagFieldParam::k5kGUniform;
     }
 
     void MachineField(const Double_t *x, Double_t *b) const;
 
-    BMap_t getMapType() const
+    MagFieldParam::BMap_t getMapType() const
     {
       return mMapType;
     }
 
-    BeamType_t getBeamType() const
+    MagFieldParam::BeamType_t getBeamType() const
     {
       return mBeamType;
     }
@@ -202,9 +230,9 @@ class MagneticField : public TVirtualMagField
 
   protected:
     // not supposed to be changed during the run, set only at the initialization via constructor
-    void initializeMachineField(BeamType_t btype, Double_t benergy);
+    void initializeMachineField(MagFieldParam::BeamType_t btype, Double_t benergy);
 
-    void setBeamType(BeamType_t type)
+    void setBeamType(MagFieldParam::BeamType_t type)
     {
       mBeamType = type;
     }
@@ -215,10 +243,10 @@ class MagneticField : public TVirtualMagField
     }
 
   protected:
-    MagneticWrapperChebyshev *mMeasuredMap; //! Measured part of the field map
-    BMap_t mMapType;                        ///< field map type
+    std::unique_ptr<MagneticWrapperChebyshev> mMeasuredMap; //! Measured part of the field map
+    MagFieldParam::BMap_t mMapType;         ///< field map type
     Double_t mSolenoid;                     ///< Solenoid field setting
-    BeamType_t mBeamType;                   ///< Beam type: A-A (mBeamType=0) or p-p (mBeamType=1)
+    MagFieldParam::BeamType_t mBeamType;    ///< Beam type: A-A (mBeamType=0) or p-p (mBeamType=1)
     Double_t mBeamEnergy;                   ///< Beam energy in GeV
 
     Int_t mDefaultIntegration;             ///< Default integration method as indicated in Geant
@@ -234,15 +262,20 @@ class MagneticField : public TVirtualMagField
     Double_t mCompensatorField1A; ///< Side A 1st compensator field
     Double_t mCompensatorField2A; ///< Side A 2nd compensator field
 
-    TNamed mParameterNames; ///< file and parameterization loadad
+    TNamed mParameterNames; ///< file and parameterization loaded
 
     static const Double_t sSolenoidToDipoleZ;  ///< conventional Z of transition from L3 to Dipole field
     static const UShort_t sPolarityConvention; ///< convention for the mapping of the curr.sign on main component sign
 
     FairLogger *mLogger;
 
+   private:
+    MagneticField(const MagneticField &src);
+
+
+
     ClassDef(AliceO2::Field::MagneticField,
-    2) // Class for all Alice MagField wrapper for measured data + Tosca parameterization
+    3) // Class for all Alice MagField wrapper for measured data + Tosca parameterization
 };
 }
 }

Common/Field/src/FieldLinkDef.h

@@ -6,5 +6,6 @@
 
 #pragma link C++ class AliceO2::Field::MagneticField+;
 #pragma link C++ class AliceO2::Field::MagneticWrapperChebyshev+;
+#pragma link C++ class AliceO2::Field::MagFieldParam+;
 
 #endif

Common/Field/src/MagFieldParam.cxx

@@ -0,0 +1,86 @@
+/// \file MagFieldParam.cxx
+/// \brief Implementation of the MagFieldParam class
+/// \author ruben.shahoyan@cern.ch
+
+#include "Field/MagFieldParam.h"
+#include "Field/MagneticField.h"
+#include "FairParamList.h"
+#include "FairRun.h"
+#include "FairRuntimeDb.h"
+
+using namespace AliceO2::Field;
+
+
+//========================================
+ClassImp(MagFieldParam);
+
+MagFieldParam::MagFieldParam(const char* name, const char* title, const char* context)
+  :FairParGenericSet(name, title, context)
+  ,mMapType(k5kG)
+  ,mBeamType(kNoBeamField)
+  ,mDefaultIntegration(0)
+  ,mFactorSol(0.)
+  ,mFactorDip(0.)
+  ,mBeamEnergy(0.)
+  ,mMaxField(0.)
+  ,mMapPath()
+{
+  /// create param for alice mag. field
+}
+
+void MagFieldParam::SetParam(const MagneticField* field)
+{
+  /// fill parameters from the initialized field
+  //  SetName(field->GetName()); ? is this needed
+  //  SetTitle(field->GetTitle());
+  mMapType = field->getMapType();
+  mBeamType = field->getBeamType();
+  mDefaultIntegration = field->Integral();
+  mFactorSol = field->getFactorSolenoid();
+  mFactorDip = field->getFactorDipole();
+  mBeamEnergy = field->getBeamEnergy();
+  mMaxField = field->Max();
+  mMapPath = field->getDataFileName();
+  //
+}
+
+void MagFieldParam::putParams(FairParamList* list)
+{
+  /// store parameters in the list
+  if (!list) return;
+  list->add("Map  Type  ID", int(mMapType));
+  list->add("Beam Type  ID", int(mBeamType));
+  list->add("Integral Type", mDefaultIntegration);
+  list->add("Fact.Solenoid", mFactorSol);
+  list->add("Fact.Dipole  ", mFactorDip);
+  list->add("Beam Energy  ", mBeamEnergy);
+  list->add("Max. Field   ", mMaxField);
+  list->add("Path to map  ", mMapPath.Data());
+  //
+}
+
+Bool_t MagFieldParam::getParams(FairParamList* list)
+{
+  /// retried parameters
+  int int2enum=0;
+  if (!list->fill("Map  Type  ID", &int2enum)) return kFALSE;
+  mMapType = static_cast<BMap_t>(int2enum);
+  if (!list->fill("Beam Type  ID", &int2enum)) return kFALSE;
+  mBeamType = static_cast<BeamType_t>(int2enum);
+  //
+  if (!list->fill("Integral Type", &mDefaultIntegration)) return kFALSE;
+  if (!list->fill("Fact.Solenoid", &mFactorSol)) return kFALSE;
+  if (!list->fill("Fact.Dipole  ", &mFactorDip)) return kFALSE;
+  if (!list->fill("Beam Energy  ", &mBeamEnergy)) return kFALSE;
+  if (!list->fill("Max. Field   ", &mMaxField)) return kFALSE;
+  FairParamObj* parpath = list->find("Path to map  ");
+  if (!parpath) return kFALSE;
+  int lgt = parpath->getLength();
+  // RS: is there a bug in FairParamList::fill(const Text_t* name,Text_t* value,const Int_t length)?
+  // I think the "if (l<length-1)" should be "if (l<length)"
+  char cbuff[lgt+2];
+  memset(cbuff,0,sizeof(char)*(lgt+2));
+  if (!list->fill("Path to map  ", cbuff, lgt+2)) return kFALSE;
+  mMapPath = cbuff;
+  return kTRUE;
+}