CbmFieldMap.h

Go to the documentation of this file.

/* Copyright (C) 2004-2020 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Mohammad Al-Turany, Denis Bertini [committer], Volker Friese */
 
// -------------------------------------------------------------------------
// -----                      CbmFieldMap header file                  -----
// -----          Created 12/01/04  by M. Al/Turany (CbmField.h)       -----
// -----                Redesign 13/02/06  by V. Friese                -----
// -------------------------------------------------------------------------
 

 
 
#ifndef CBMFIELDMAP_H
#define CBMFIELDMAP_H 1
 
 
#include <FairField.h>  // for FairField
 
#include <Rtypes.h>      // for THashConsistencyHolder, ClassDef
#include <RtypesCore.h>  // for Double_t, Int_t, Bool_t, Option_t
#include <TArrayF.h>     // for TArrayF
#include <TString.h>     // for TString
 
class CbmFieldMapCreator;
class CbmFieldMapData;
class CbmFieldPar;
class TArrayF;
 
class CbmFieldMap : public FairField {
 
 
 public:
  CbmFieldMap();
 

  CbmFieldMap(const char* mapName, const char* fileType = "R");
 

  CbmFieldMap(CbmFieldPar* fieldPar);
 

  CbmFieldMap(CbmFieldMapCreator* creator);
 

  virtual ~CbmFieldMap();
 

  virtual void Init() override;
 

  virtual void Init(Int_t nX, Double_t xMin, Double_t xMax, Int_t nY, Double_t yMin, Double_t yMax, Int_t nZ,
                    Double_t zMin, Double_t zMax, TArrayF* bx, TArrayF* by, TArrayF* bz);
 

  virtual Double_t GetBx(Double_t x, Double_t y, Double_t z) override;
  virtual Double_t GetBy(Double_t x, Double_t y, Double_t z) override;
  virtual Double_t GetBz(Double_t x, Double_t y, Double_t z) override;

  virtual void GetFieldValue(const Double_t point[3], Double_t* bField) override;

  void WriteAsciiFile(const char* fileName);
 

  void WriteRootFile(const char* fileName, const char* mapName);
 

  virtual void SetPosition(Double_t x, Double_t y, Double_t z);
 

  virtual void SetScale(Double_t factor) { fScale = factor; }
 

  Double_t GetXmin() const { return fXmin; }
  Double_t GetYmin() const { return fYmin; }
  Double_t GetZmin() const { return fZmin; }
  Double_t GetXmax() const { return fXmax; }
  Double_t GetYmax() const { return fYmax; }
  Double_t GetZmax() const { return fZmax; }
  Double_t GetXstep() const { return fXstep; }
  Double_t GetYstep() const { return fYstep; }
  Double_t GetZstep() const { return fZstep; }
  Int_t GetNx() const { return fNx; }
  Int_t GetNy() const { return fNy; }
  Int_t GetNz() const { return fNz; }
 

  Double_t GetPositionX() const { return fPosX; }
  Double_t GetPositionY() const { return fPosY; }
  Double_t GetPositionZ() const { return fPosZ; }
 

  Double_t GetScale() const { return fScale; }
 

  TArrayF* GetBx() const { return fBx; }
  TArrayF* GetBy() const { return fBy; }
  TArrayF* GetBz() const { return fBz; }
 

  const char* GetFileName() { return fFileName.Data(); }
 

  virtual void Print(Option_t* = "") const override;
 
 
 protected:
  struct LocalCoordinates {
    Double_t x{0.};      // Local x coordinate [cm]
    Double_t y{0.};      // Local y coordinate [cm]
    Double_t z{0.};      // Local z coordinate [cm]
    Double_t hemiX{1.};  // Hemisphere in x direction (1 or -1)
    Double_t hemiY{1.};  // Hemisphere in y direction (1 or -1)
    Double_t hemiZ{1.};  // Hemisphere in z direction (1 or -1)
  };

  [[gnu::always_inline]] LocalCoordinates Global2Local(Double_t x, Double_t y, Double_t z) const
  {
    return {x - fPosX, y - fPosY, z - fPosZ, (x < fPosX ? -1. : 1.), (y < fPosY ? -1. : 1.), (z < fPosZ ? -1. : 1.)};
  }

  struct GridCoordinates {
    Bool_t isValid{kFALSE};  // kTRUE if point is inside the grid
    Int_t ix{0};             // Grid cell index in x
    Int_t iy{0};             // Grid cell index in y
    Int_t iz{0};             // Grid cell index in z
    Double_t dx{0.};         // Relative distance from grid point in x [units of cell size]
    Double_t dy{0.};         // Relative distance from grid point in y [units of cell size]
    Double_t dz{0.};         // Relative distance from grid point in z [units of cell size]
  };

  GridCoordinates Local2Grid(const LocalCoordinates& local) const;
 

  [[gnu::always_inline]] Double_t Interpolate(const TArrayF* B, const GridCoordinates& point) const;
 

  void Reset();
 

  void ReadAsciiFile(const char* fileName);

  void ReadAsciiFile2018(const char* fileName);

  void ReadRootFile(const char* fileName, const char* mapName);
 

  void SetField(const CbmFieldMapData* data);

  TString fFileName;
 

  Double_t fScale;
 

  Double_t fPosX, fPosY, fPosZ;
 

  Double_t fXmin, fXmax, fXstep;
  Double_t fYmin, fYmax, fYstep;
  Double_t fZmin, fZmax, fZstep;
 

  Int_t fNx, fNy, fNz;
 

  TArrayF* fBx;
  TArrayF* fBy;
  TArrayF* fBz;
 

 
  Double_t fBxOrigin;  
  Double_t fByOrigin;  
  Double_t fBzOrigin;  
 
 private:
  CbmFieldMap(const CbmFieldMap&);
  CbmFieldMap& operator=(const CbmFieldMap&);

  [[gnu::always_inline]] GridCoordinates Global2Grid(Double_t x, Double_t y, Double_t z) const
  {
    return Local2Grid(Global2Local(x, y, z));
  }
 
  ClassDefOverride(CbmFieldMap, 2)
};
 
 
// -----------   Check whether a point is inside the map   ----------------
[[gnu::always_inline]] inline CbmFieldMap::GridCoordinates
CbmFieldMap::Local2Grid(const CbmFieldMap::LocalCoordinates& loc) const
{
  // --- Determine grid cell
  Double_t x = (loc.x - fXmin) / fXstep;
  Double_t y = (loc.y - fYmin) / fYstep;
  Double_t z = (loc.z - fZmin) / fZstep;
 
  Int_t ix = (Int_t) std::floor(x);
  Int_t iy = (Int_t) std::floor(y);
  Int_t iz = (Int_t) std::floor(z);
 
  bool isValid = (ix >= 0 && ix < fNx - 1 && iy >= 0 && iy < fNy - 1 && iz >= 0 && iz < fNz - 1);
 
  // Relative distance from grid point (in units of cell size)
  Double_t dx = x - (Double_t) ix;
  Double_t dy = y - (Double_t) iy;
  Double_t dz = z - (Double_t) iz;
 
  return {isValid, ix, iy, iz, dx, dy, dz};
}
// ------------------------------------------------------------------------
 
 
// ------------   Interpolation in a grid cell (private)  -----------------
[[gnu::always_inline]] inline Double_t
CbmFieldMap::Interpolate(const TArrayF* B,  // Array with field values - fBx, fBy, or fBz
                         const CbmFieldMap::GridCoordinates& p) const
{
  if (!p.isValid) {
    return 0.;
  }
 
  int Nyz    = fNy * fNz;  // Number of grid points in y-z plane
  int ind000 = p.ix * Nyz + p.iy * fNz + p.iz;
  int ind111 = ind000 + Nyz + fNz + 1;
 
  if (ind000 < 0 || ind111 >= B->GetSize()) {
    LOG(error) << "CbmFieldMap::Interpolate: Index out of bounds!";
    return 0.;
  }
 
  const auto* arr = B->GetArray();
 
  // Get field values at grid cell corners
 
  Double_t H000 = arr[ind000];              // Bx, By, or Bz at (ix, iy, iz)
  Double_t H001 = arr[ind000 + 1];          // Bx, By, or Bz at (ix, iy, iz+1)
  Double_t H010 = arr[ind000 + fNz];        // Bx, By, or Bz at (ix, iy+1, iz)
  Double_t H011 = arr[ind000 + fNz + 1];    // Bx, By, or Bz at (ix, iy+1, iz+1)
  Double_t H100 = arr[ind000 + Nyz];        // Bx, By, or Bz at (ix+1, iy, iz)
  Double_t H101 = arr[ind000 + Nyz + 1];    // Bx, By, or Bz at (ix+1, iy, iz+1)
  Double_t H110 = arr[ind000 + Nyz + fNz];  // Bx, By, or Bz at (ix+1, iy+1, iz)
  Double_t H111 = arr[ind111];              // Bx, By, or Bz at (ix+1, iy+1, iz+1)
 
  // Interpolate in x coordinate. Interpolated field in yz (2-dim)
 
  Double_t Hyz00 = H000 + (H100 - H000) * p.dx;
  Double_t Hyz10 = H010 + (H110 - H010) * p.dx;
  Double_t Hyz01 = H001 + (H101 - H001) * p.dx;
  Double_t Hyz11 = H011 + (H111 - H011) * p.dx;
 
  // Interpolate in y coordinate. Interpolated field in z (1-dim)
  Double_t Hz0 = Hyz00 + (Hyz10 - Hyz00) * p.dy;
  Double_t Hz1 = Hyz01 + (Hyz11 - Hyz01) * p.dy;
 
  // Interpolate in z coordinate
  return Hz0 + (Hz1 - Hz0) * p.dz;
}
// ------------------------------------------------------------------------
 
 
#endif