#ifndef ALIMAGWRAPCHEB_H
#define ALIMAGWRAPCHEB_H
#include <TMath.h>
#include <TNamed.h>
#include <TObjArray.h>
#include <TStopwatch.h>
#include "AliCheb3D.h"
#ifndef _MAGCHEB_CACHE_
#define _MAGCHEB_CACHE_ // use to spead up, but then Field calls are not thread safe
#endif
class TSystem;
class TArrayF;
class TArrayI;
class AliMagWrapCheb: public TNamed
{
public:
AliMagWrapCheb();
AliMagWrapCheb(const AliMagWrapCheb& src);
~AliMagWrapCheb() {Clear();}
void CopyFrom(const AliMagWrapCheb& src);
AliMagWrapCheb& operator=(const AliMagWrapCheb& rhs);
virtual void Clear(const Option_t * = "");
Int_t GetNParamsSol() const {return fNParamsSol;}
Int_t GetNSegZSol() const {return fNZSegSol;}
Float_t* GetSegZSol() const {return fSegZSol;}
Int_t GetNParamsTPCInt() const {return fNParamsTPC;}
Int_t GetNSegZTPCInt() const {return fNZSegTPC;}
Int_t GetNParamsTPCRatInt() const {return fNParamsTPCRat;}
Int_t GetNSegZTPCRatInt() const {return fNZSegTPCRat;}
Int_t GetNParamsDip() const {return fNParamsDip;}
Int_t GetNSegZDip() const {return fNZSegDip;}
Float_t GetMaxZ() const {return GetMaxZSol();}
Float_t GetMinZ() const {return fParamsDip ? GetMinZDip() : GetMinZSol();}
Float_t GetMinZSol() const {return fMinZSol;}
Float_t GetMaxZSol() const {return fMaxZSol;}
Float_t GetMaxRSol() const {return fMaxRSol;}
Float_t GetMinZDip() const {return fMinZDip;}
Float_t GetMaxZDip() const {return fMaxZDip;}
Float_t GetMinZTPCInt() const {return fMinZTPC;}
Float_t GetMaxZTPCInt() const {return fMaxZTPC;}
Float_t GetMaxRTPCInt() const {return fMaxRTPC;}
Float_t GetMinZTPCRatInt() const {return fMinZTPCRat;}
Float_t GetMaxZTPCRatInt() const {return fMaxZTPCRat;}
Float_t GetMaxRTPCRatInt() const {return fMaxRTPCRat;}
AliCheb3D* GetParamSol(Int_t ipar) const {return (AliCheb3D*)fParamsSol->UncheckedAt(ipar);}
AliCheb3D* GetParamTPCRatInt(Int_t ipar) const {return (AliCheb3D*)fParamsTPCRat->UncheckedAt(ipar);}
AliCheb3D* GetParamTPCInt(Int_t ipar) const {return (AliCheb3D*)fParamsTPC->UncheckedAt(ipar);}
AliCheb3D* GetParamDip(Int_t ipar) const {return (AliCheb3D*)fParamsDip->UncheckedAt(ipar);}
virtual void Print(Option_t * = "") const;
virtual void Field(const Double_t *xyz, Double_t *b) const;
Double_t GetBz(const Double_t *xyz) const;
void FieldCyl(const Double_t *rphiz, Double_t *b) const;
void GetTPCInt(const Double_t *xyz, Double_t *b) const;
void GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const;
void GetTPCRatInt(const Double_t *xyz, Double_t *b) const;
void GetTPCRatIntCyl(const Double_t *rphiz, Double_t *b) const;
Int_t FindSolSegment(const Double_t *xyz) const;
Int_t FindTPCSegment(const Double_t *xyz) const;
Int_t FindTPCRatSegment(const Double_t *xyz) const;
Int_t FindDipSegment(const Double_t *xyz) const;
static void CylToCartCylB(const Double_t *rphiz, const Double_t *brphiz,Double_t *bxyz);
static void CylToCartCartB(const Double_t *xyz, const Double_t *brphiz,Double_t *bxyz);
static void CartToCylCartB(const Double_t *xyz, const Double_t *bxyz, Double_t *brphiz);
static void CartToCylCylB(const Double_t *rphiz, const Double_t *bxyz, Double_t *brphiz);
static void CartToCyl(const Double_t *xyz, Double_t *rphiz);
static void CylToCart(const Double_t *rphiz,Double_t *xyz);
#ifdef _INC_CREATION_ALICHEB3D_ // see AliCheb3D.h for explanation
void LoadData(const char* inpfile);
AliMagWrapCheb(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 AddParamSol(const AliCheb3D* param);
void AddParamTPCInt(const AliCheb3D* param);
void AddParamTPCRatInt(const AliCheb3D* param);
void AddParamDip(const AliCheb3D* 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 BuildTableSol();
void BuildTableDip();
void BuildTableTPCInt();
void BuildTableTPCRatInt();
void ResetTPCInt();
void ResetTPCRatInt();
void ResetSol();
void ResetDip();
#endif
protected:
void FieldCylSol(const Double_t *rphiz, Double_t *b) const;
Double_t FieldCylSolBz(const Double_t *rphiz) const;
protected:
Int_t fNParamsSol;
Int_t fNZSegSol;
Int_t fNPSegSol;
Int_t fNRSegSol;
Float_t* fSegZSol;
Float_t* fSegPSol;
Float_t* fSegRSol;
Int_t* fBegSegPSol;
Int_t* fNSegPSol;
Int_t* fBegSegRSol;
Int_t* fNSegRSol;
Int_t* fSegIDSol;
Float_t fMinZSol;
Float_t fMaxZSol;
TObjArray* fParamsSol;
Float_t fMaxRSol;
Int_t fNParamsTPC;
Int_t fNZSegTPC;
Int_t fNPSegTPC;
Int_t fNRSegTPC;
Float_t* fSegZTPC;
Float_t* fSegPTPC;
Float_t* fSegRTPC;
Int_t* fBegSegPTPC;
Int_t* fNSegPTPC;
Int_t* fBegSegRTPC;
Int_t* fNSegRTPC;
Int_t* fSegIDTPC;
Float_t fMinZTPC;
Float_t fMaxZTPC;
TObjArray* fParamsTPC;
Float_t fMaxRTPC;
Int_t fNParamsTPCRat;
Int_t fNZSegTPCRat;
Int_t fNPSegTPCRat;
Int_t fNRSegTPCRat;
Float_t* fSegZTPCRat;
Float_t* fSegPTPCRat;
Float_t* fSegRTPCRat;
Int_t* fBegSegPTPCRat;
Int_t* fNSegPTPCRat;
Int_t* fBegSegRTPCRat;
Int_t* fNSegRTPCRat;
Int_t* fSegIDTPCRat;
Float_t fMinZTPCRat;
Float_t fMaxZTPCRat;
TObjArray* fParamsTPCRat;
Float_t fMaxRTPCRat;
Int_t fNParamsDip;
Int_t fNZSegDip;
Int_t fNYSegDip;
Int_t fNXSegDip;
Float_t* fSegZDip;
Float_t* fSegYDip;
Float_t* fSegXDip;
Int_t* fBegSegYDip;
Int_t* fNSegYDip;
Int_t* fBegSegXDip;
Int_t* fNSegXDip;
Int_t* fSegIDDip;
Float_t fMinZDip;
Float_t fMaxZDip;
TObjArray* fParamsDip;
#ifdef _MAGCHEB_CACHE_
mutable AliCheb3D* fCacheSol;
mutable AliCheb3D* fCacheDip;
mutable AliCheb3D* fCacheTPCInt;
mutable AliCheb3D* fCacheTPCRat;
#endif
ClassDef(AliMagWrapCheb,8)
};
inline void AliMagWrapCheb::FieldCyl(const Double_t *rphiz, Double_t *b) const
{
b[0] = b[1] = b[2] = 0;
FieldCylSol(rphiz,b);
}
inline void AliMagWrapCheb::CylToCartCylB(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 AliMagWrapCheb::CylToCartCartB(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 AliMagWrapCheb::CartToCylCartB(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 AliMagWrapCheb::CartToCylCylB(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 AliMagWrapCheb::CartToCyl(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 AliMagWrapCheb::CylToCart(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];
}
#endif