MagneticWrapperChebyshev.h

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.
 
 
#ifndef ALICEO2_FIELD_MAGNETICWRAPPERCHEBYSHEV_H_
#define ALICEO2_FIELD_MAGNETICWRAPPERCHEBYSHEV_H_
 
#include <TMath.h>                     // for ATan2, Cos, Sin, Sqrt
#include <TNamed.h>                    // for TNamed
#include <TObjArray.h>                 // for TObjArray
#include "MathUtils/Chebyshev3D.h"     // for Chebyshev3D
#include "MathUtils/Chebyshev3DCalc.h" // for _INC_CREATION_Chebyshev3D_
#include "Rtypes.h"                    // for Double_t, Int_t, Float_t, etc
 
namespace o2
{
namespace field
{
 
class MagneticWrapperChebyshev : public TNamed
{
 
 public:
  MagneticWrapperChebyshev();
 
  MagneticWrapperChebyshev(const MagneticWrapperChebyshev& src);
 
  ~MagneticWrapperChebyshev() override
  {
    Clear();
  }
 
  void copyFrom(const MagneticWrapperChebyshev& src);
 
  MagneticWrapperChebyshev& operator=(const MagneticWrapperChebyshev& rhs);
 
  void Clear(const Option_t* = "") override;
 
  Int_t getNumberOfParametersSol() const
  {
    return mNumberOfParameterizationSolenoid;
  }
 
  Int_t getNumberOmCoordinatesSegmentsZSolenoid() const
  {
    return mNumberOfDistinctZSegmentsSolenoid;
  }
 
  Float_t* getSegZSol() const
  {
    return mCoordinatesSegmentsZSolenoid;
  }
 
  Int_t getNumberOfParametersTPCIntegral() const
  {
    return mNumberOfParameterizationTPC;
  }
 
  Int_t getNumberOmCoordinatesSegmentsZTPCInt() const
  {
    return mNumberOfDistinctZSegmentsTPC;
  }
 
  Int_t getNumberOfParametersTPCRatIntegral() const
  {
    return mNumberOfParameterizationTPCRat;
  }
 
  Int_t getNumberOmCoordinatesSegmentsZTPCRatIntegral() const
  {
    return mNumberOfDistinctZSegmentsTPCRat;
  }
 
  Int_t getNumberOfParametersDip() const
  {
    return mNumberOfParameterizationDipole;
  }
 
  Int_t getNumberOmCoordinatesSegmentsZDipole() const
  {
    return mNumberOfDistinctZSegmentsDipole;
  }
 
  Float_t getMaxZ() const
  {
    return getMaxZSol();
  }
 
  Float_t getMinZ() const
  {
    return mParameterizationDipole ? getMinZDip() : getMinZSol();
  }
 
  Float_t getMinZSol() const
  {
    return mMinZSolenoid;
  }
 
  Float_t getMaxZSol() const
  {
    return mMaxZSolenoid;
  }
 
  Float_t getMaxRSol() const
  {
    return mMaxRadiusSolenoid;
  }
 
  Float_t getMinZDip() const
  {
    return mMinDipoleZ;
  }
 
  Float_t getMaxZDip() const
  {
    return mMaxDipoleZ;
  }
 
  Float_t getMinZTPCIntegral() const
  {
    return mMinZTPC;
  }
 
  Float_t getMaxZTPCIntegral() const
  {
    return mMaxZTPC;
  }
 
  Float_t getMaxRTPCIntegral() const
  {
    return mMaxRadiusTPC;
  }
 
  Float_t getMinZTPCRatIntegral() const
  {
    return mMinZTPCRat;
  }
 
  Float_t getMaxZTPCRatIntegral() const
  {
    return mMaxZTPCRat;
  }
 
  Float_t getMaxRTPCRatIntegral() const
  {
    return mMaxRadiusTPCRat;
  }
 
  o2::math_utils::Chebyshev3D* getParameterSolenoid(Int_t ipar) const
  {
    return (o2::math_utils::Chebyshev3D*)mParameterizationSolenoid->UncheckedAt(ipar);
  }
 
  o2::math_utils::Chebyshev3D* getParameterTPCRatIntegral(Int_t ipar) const
  {
    return (o2::math_utils::Chebyshev3D*)mParameterizationTPCRat->UncheckedAt(ipar);
  }
 
  o2::math_utils::Chebyshev3D* getParameterTPCIntegral(Int_t ipar) const
  {
    return (o2::math_utils::Chebyshev3D*)mParameterizationTPC->UncheckedAt(ipar);
  }
 
  o2::math_utils::Chebyshev3D* getParameterDipole(Int_t ipar) const
  {
    return (o2::math_utils::Chebyshev3D*)mParameterizationDipole->UncheckedAt(ipar);
  }
 
  void Print(Option_t* = "") const override;
 
  virtual void Field(const Double_t* xyz, Double_t* b) const;
 
  Double_t getBz(const Double_t* xyz) const;
 
  void fieldCylindrical(const Double_t* rphiz, Double_t* b) const;
 
  void getTPCIntegral(const Double_t* xyz, Double_t* b) const;
 
  // Computes field integral in TPC region in Cylindircal coordinates
  // note: the check for the point being inside the parameterized region is done outside
  void getTPCIntegralCylindrical(const Double_t* rphiz, Double_t* b) const;
 
  void getTPCRatIntegral(const Double_t* xyz, Double_t* b) const;
 
  // Computes field integral in TPCRat region in Cylindircal coordinates
  // note: the check for the point being inside the parameterized region is done outside
  void getTPCRatIntegralCylindrical(const Double_t* rphiz, Double_t* b) const;
 
  Int_t findSolenoidSegment(const Double_t* xyz) const;
 
  Int_t findTPCSegment(const Double_t* xyz) const;
 
  Int_t findTPCRatSegment(const Double_t* xyz) const;
 
  Int_t findDipoleSegment(const Double_t* xyz) const;
 
  static void cylindricalToCartesianCylB(const Double_t* rphiz, const Double_t* brphiz, Double_t* bxyz);
 
  static void cylindricalToCartesianCartB(const Double_t* xyz, const Double_t* brphiz, Double_t* bxyz);
 
  static void cartesianToCylindricalCartB(const Double_t* xyz, const Double_t* bxyz, Double_t* brphiz);
 
  static void cartesianToCylindricalCylB(const Double_t* rphiz, const Double_t* bxyz, Double_t* brphiz);
 
  static void cartesianToCylindrical(const Double_t* xyz, Double_t* rphiz);
 
  static void cylindricalToCartesian(const Double_t* rphiz, Double_t* xyz);
 
#ifdef _INC_CREATION_Chebyshev3D_ // see Cheb3D.h for explanation
  void loadData(const char* inpfile);
 
  MagneticWrapperChebyshev(const char* inputFile);
 
  void saveData(const char* outfile) const;
 
  Int_t segmentDimension(Float_t** seg, const TObjArray* par, int npar, int dim, Float_t xmn, Float_t xmx, Float_t ymn,
                         Float_t ymx, Float_t zmn, Float_t zmx);
 
  void addParameterSolenoid(const o2::math_utils::Chebyshev3D* param);
 
  // Adds new parameterization piece for TPCIntegral
  // NOTE: pieces must be added strictly in increasing R then increasing Z order
 
  void addParameterTPCIntegral(const o2::math_utils::Chebyshev3D* param);
  // NOTE: pieces must be added strictly in increasing R then increasing Z order
  void addParameterTPCRatIntegral(const o2::math_utils::Chebyshev3D* param);
 
  void addParameterDipole(const o2::math_utils::Chebyshev3D* param);
 
  void buildTable(Int_t npar, TObjArray* parArr, Int_t& nZSeg, Int_t& nYSeg, Int_t& nXSeg, Float_t& minZ, Float_t& maxZ,
                  Float_t** segZ, Float_t** segY, Float_t** segX, Int_t** begSegY, Int_t** nSegY, Int_t** begSegX,
                  Int_t** nSegX, Int_t** segID);
 
  void buildTableSolenoid();
 
  void buildTableDipole();
 
  void buildTableTPCIntegral();
 
  void buildTableTPCRatIntegral();
 
  void resetTPCIntegral();
 
  void resetTPCRatIntegral();
 
  void resetSolenoid();
 
  void resetDipole();
#endif
 
 protected:
  void fieldCylindricalSolenoid(const Double_t* rphiz, Double_t* b) const;
 
  Double_t fieldCylindricalSolenoidBz(const Double_t* rphiz) const;
 
 private:
  Int_t mNumberOfParameterizationSolenoid;  
  Int_t mNumberOfDistinctZSegmentsSolenoid; 
  Int_t mNumberOfDistinctPSegmentsSolenoid; 
  Int_t mNumberOfDistinctRSegmentsSolenoid; 
  Float_t*
    mCoordinatesSegmentsZSolenoid;        //[mNumberOfDistinctZSegmentsSolenoid] coordinates of distinct Z segments in Solenoid
  Float_t* mCoordinatesSegmentsPSolenoid; //[mNumberOfDistinctPSegmentsSolenoid] coordinates of P segments for each
  // Zsegment in Solenoid
  Float_t* mCoordinatesSegmentsRSolenoid; //[mNumberOfDistinctRSegmentsSolenoid] coordinates of R segments for each
  // Psegment in Solenoid
  Int_t* mBeginningOfSegmentsPSolenoid; //[mNumberOfDistinctPSegmentsSolenoid] beginning of P segments array for each Z
  // segment
  Int_t* mNumberOfSegmentsPSolenoid;    //[mNumberOfDistinctZSegmentsSolenoid] number of P segments for each Z segment
  Int_t* mBeginningOfSegmentsRSolenoid; //[mNumberOfDistinctPSegmentsSolenoid] beginning of R segments array for each P
  // segment
  Int_t* mNumberOfRSegmentsSolenoid; //[mNumberOfDistinctPSegmentsSolenoid] number of R segments for each P segment
  Int_t*
    mSegmentIdSolenoid;                 //[mNumberOfDistinctRSegmentsSolenoid] ID of the solenoid parameterization for given RPZ segment
  Float_t mMinZSolenoid;                
  Float_t mMaxZSolenoid;                
  TObjArray* mParameterizationSolenoid; 
  Float_t mMaxRadiusSolenoid;           
 
  Int_t mNumberOfParameterizationTPC;  
  Int_t mNumberOfDistinctZSegmentsTPC; 
  Int_t mNumberOfDistinctPSegmentsTPC; 
  Int_t mNumberOfDistinctRSegmentsTPC; 
  Float_t* mCoordinatesSegmentsZTPC;   //[mNumberOfDistinctZSegmentsTPC] coordinates of distinct Z segments in TPCint
  Float_t*
    mCoordinatesSegmentsPTPC; //[mNumberOfDistinctPSegmentsTPC] coordinates of P segments for each Zsegment in TPCint
  Float_t*
    mCoordinatesSegmentsRTPC;      //[mNumberOfDistinctRSegmentsTPC] coordinates of R segments for each Psegment in TPCint
  Int_t* mBeginningOfSegmentsPTPC; //[mNumberOfDistinctPSegmentsTPC] beginning of P segments array for each Z segment
  Int_t* mNumberOfSegmentsPTPC;    //[mNumberOfDistinctZSegmentsTPC] number of P segments for each Z segment
  Int_t* mBeginningOfSegmentsRTPC; //[mNumberOfDistinctPSegmentsTPC] beginning of R segments array for each P segment
  Int_t* mNumberOfRSegmentsTPC;    //[mNumberOfDistinctPSegmentsTPC] number of R segments for each P segment
  Int_t* mSegmentIdTPC;            //[mNumberOfDistinctRSegmentsTPC] ID of the TPCint parameterization for given RPZ segment
  Float_t mMinZTPC;                
  Float_t mMaxZTPC;                
  TObjArray* mParameterizationTPC; 
  Float_t mMaxRadiusTPC;           
 
  Int_t
    mNumberOfParameterizationTPCRat; 
  Int_t mNumberOfDistinctZSegmentsTPCRat; 
  Int_t mNumberOfDistinctPSegmentsTPCRat; 
  Int_t mNumberOfDistinctRSegmentsTPCRat; 
  Float_t*
    mCoordinatesSegmentsZTPCRat;        //[mNumberOfDistinctZSegmentsTPCRat] coordinates of distinct Z segments in TpcRatInt
  Float_t* mCoordinatesSegmentsPTPCRat; //[mNumberOfDistinctPSegmentsTPCRat] coordinates of P segments for each Zsegment
  // in TpcRatInt
  Float_t* mCoordinatesSegmentsRTPCRat; //[mNumberOfDistinctRSegmentsTPCRat] coordinates of R segments for each Psegment
  // in TpcRatInt
  Int_t*
    mBeginningOfSegmentsPTPCRat;   //[mNumberOfDistinctPSegmentsTPCRat] beginning of P segments array for each Z segment
  Int_t* mNumberOfSegmentsPTPCRat; //[mNumberOfDistinctZSegmentsTPCRat] number of P segments for each Z segment
  Int_t*
    mBeginningOfSegmentsRTPCRat;   //[mNumberOfDistinctPSegmentsTPCRat] beginning of R segments array for each P segment
  Int_t* mNumberOfRSegmentsTPCRat; //[mNumberOfDistinctPSegmentsTPCRat] number of R segments for each P segment
  Int_t*
    mSegmentIdTPCRat;                 //[mNumberOfDistinctRSegmentsTPCRat] ID of the TpcRatInt parameterization for given RPZ segment
  Float_t mMinZTPCRat;                
  Float_t mMaxZTPCRat;                
  TObjArray* mParameterizationTPCRat; 
  Float_t mMaxRadiusTPCRat;           
 
  Int_t mNumberOfParameterizationDipole;  
  Int_t mNumberOfDistinctZSegmentsDipole; 
  Int_t mNumberOfDistinctYSegmentsDipole; 
  Int_t mNumberOfDistinctXSegmentsDipole; 
  Float_t*
    mCoordinatesSegmentsZDipole;        //[mNumberOfDistinctZSegmentsDipole] coordinates of distinct Z segments in Dipole
  Float_t* mCoordinatesSegmentsYDipole; //[mNumberOfDistinctYSegmentsDipole] coordinates of Y segments for each Zsegment
  // in Dipole
  Float_t* mCoordinatesSegmentsXDipole; //[mNumberOfDistinctXSegmentsDipole] coordinates of X segments for each Ysegment
  // in Dipole
  Int_t*
    mBeginningOfSegmentsYDipole;   //[mNumberOfDistinctZSegmentsDipole] beginning of Y segments array for each Z segment
  Int_t* mNumberOfSegmentsYDipole; //[mNumberOfDistinctZSegmentsDipole] number of Y segments for each Z segment
  Int_t*
    mBeginningOfSegmentsXDipole;      //[mNumberOfDistinctYSegmentsDipole] beginning of X segments array for each Y segment
  Int_t* mNumberOfSegmentsXDipole;    //[mNumberOfDistinctYSegmentsDipole] number of X segments for each Y segment
  Int_t* mSegmentIdDipole;            //[mNumberOfDistinctXSegmentsDipole] ID of the dipole parameterization for given XYZ segment
  Float_t mMinDipoleZ;                
  Float_t mMaxDipoleZ;                
  TObjArray* mParameterizationDipole; 
 
  ClassDefOverride(o2::field::MagneticWrapperChebyshev,
                   2) // Wrapper class for the set of Chebishev parameterizations of Alice mag.field
};
 
inline void MagneticWrapperChebyshev::fieldCylindrical(const Double_t* rphiz, Double_t* b) const
{
  //  if (rphiz[2]<GetMinZSol() || rphiz[2]>GetMaxZSol() || rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;}
  b[0] = b[1] = b[2] = 0;
  fieldCylindricalSolenoid(rphiz, b);
}
 
inline void MagneticWrapperChebyshev::cylindricalToCartesianCylB(const Double_t* rphiz, const Double_t* brphiz,
                                                                 Double_t* bxyz)
{
  Double_t btr = TMath::Sqrt(brphiz[0] * brphiz[0] + brphiz[1] * brphiz[1]);
  Double_t psiPLUSphi = TMath::ATan2(brphiz[1], brphiz[0]) + rphiz[1];
  bxyz[0] = btr * TMath::Cos(psiPLUSphi);
  bxyz[1] = btr * TMath::Sin(psiPLUSphi);
  bxyz[2] = brphiz[2];
}
 
inline void MagneticWrapperChebyshev::cylindricalToCartesianCartB(const Double_t* xyz, const Double_t* brphiz,
                                                                  Double_t* bxyz)
{
  Double_t btr = TMath::Sqrt(brphiz[0] * brphiz[0] + brphiz[1] * brphiz[1]);
  Double_t phiPLUSpsi = TMath::ATan2(xyz[1], xyz[0]) + TMath::ATan2(brphiz[1], brphiz[0]);
  bxyz[0] = btr * TMath::Cos(phiPLUSpsi);
  bxyz[1] = btr * TMath::Sin(phiPLUSpsi);
  bxyz[2] = brphiz[2];
}
 
inline void MagneticWrapperChebyshev::cartesianToCylindricalCartB(const Double_t* xyz, const Double_t* bxyz,
                                                                  Double_t* brphiz)
{
  Double_t btr = TMath::Sqrt(bxyz[0] * bxyz[0] + bxyz[1] * bxyz[1]);
  Double_t psiMINphi = TMath::ATan2(bxyz[1], bxyz[0]) - TMath::ATan2(xyz[1], xyz[0]);
 
  brphiz[0] = btr * TMath::Cos(psiMINphi);
  brphiz[1] = btr * TMath::Sin(psiMINphi);
  brphiz[2] = bxyz[2];
}
 
inline void MagneticWrapperChebyshev::cartesianToCylindricalCylB(const Double_t* rphiz, const Double_t* bxyz,
                                                                 Double_t* brphiz)
{
  Double_t btr = TMath::Sqrt(bxyz[0] * bxyz[0] + bxyz[1] * bxyz[1]);
  Double_t psiMINphi = TMath::ATan2(bxyz[1], bxyz[0]) - rphiz[1];
  brphiz[0] = btr * TMath::Cos(psiMINphi);
  brphiz[1] = btr * TMath::Sin(psiMINphi);
  brphiz[2] = bxyz[2];
}
 
inline void MagneticWrapperChebyshev::cartesianToCylindrical(const Double_t* xyz, Double_t* rphiz)
{
  rphiz[0] = TMath::Sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]);
  rphiz[1] = TMath::ATan2(xyz[1], xyz[0]);
  rphiz[2] = xyz[2];
}
 
inline void MagneticWrapperChebyshev::cylindricalToCartesian(const Double_t* rphiz, Double_t* xyz)
{
  xyz[0] = rphiz[0] * TMath::Cos(rphiz[1]);
  xyz[1] = rphiz[0] * TMath::Sin(rphiz[1]);
  xyz[2] = rphiz[2];
}
} // namespace field
} // namespace o2
 
#endif