STEER/AliMagWrapCheb.h

   1
   2 // Author: ruben.shahoyan@cern.ch   20/03/2007
   3
   4 ///////////////////////////////////////////////////////////////////////////////////
   5 //                                                                               //
   6 //  Wrapper for the set of mag.field parameterizations by Chebyshev polinomials  //
   7 //  To obtain the field in cartesian coordinates/components use                  //
   8 //    Field(float* xyz, float* bxyz);                                            //
   9 //  For cylindrical coordinates/components:                                      //
  10 //    FieldCyl(float* rphiz, float* brphiz)                                      //
  11 //                                                                               //
  12 //  The solenoid part is parameterized in the volume  R<500, -550<Z<550 cm       //
  13 //                                                                               //
  14 //  The region R<423 cm,  -343.3<Z<481.3 for 30kA and -343.3<Z<481.3 for 12kA    //
  15 //  is parameterized using measured data while outside the Tosca calculation     //
  16 //  is used (matched to data on the boundary of the measurements)                //
  17 //                                                                               //
  18 //  Two options are possible:                                                    //
  19 //  1) _BRING_TO_BOUNDARY_ is defined in the AliCheb3D:                          //
  20 //     If the querried point is outside of the validity region then the field    //
  21 //     at the closest point on the fitted surface is returned.                   //
  22 //  2) _BRING_TO_BOUNDARY_ is not defined in the AliCheb3D:                      //
  23 //     If the querried point is outside of the validity region the return        //
  24 //     value for the field components are set to 0.                              //
  25 //                                                                               //
  26 //  To obtain the field integral in the TPC region from given point to nearest   //
  27 //  cathod plane (+- 250 cm) use:                                                //
  28 //  GetTPCInt(float* xyz, float* bxyz);  for Cartesian frame                     //
  29 //  or                                                                           //
  30 //  GetTPCIntCyl(Float_t *rphiz, Float_t *b); for Cylindrical frame              //
  31 //                                                                               //
  32 //                                                                               //
  33 //  The units are kiloGauss and cm.                                              //
  34 //                                                                               //
  35 ///////////////////////////////////////////////////////////////////////////////////
  36
  37 #ifndef ALIMAGWRAPCHEB_H
  38 #define ALIMAGWRAPCHEB_H
  39
  40 #include <TMath.h>
  41 #include <TNamed.h>
  42 #include "AliCheb3D.h"
  43
  44 class TSystem;
  45 class TArrayF;
  46 class TArrayI;
  47
  48 class AliMagWrapCheb: public TNamed
  49 {
  50  public:
  51   AliMagWrapCheb();
  52   AliMagWrapCheb(const AliMagWrapCheb& src);
  53   ~AliMagWrapCheb() {Clear();}
  54   //
  55   void       CopyFrom(const AliMagWrapCheb& src);
  56   AliMagWrapCheb& operator=(const AliMagWrapCheb& rhs);
  57   virtual void Clear(const Option_t * = "");
  58   //
  59   Int_t      GetNParamsSol()                              const {return fNParamsSol;}
  60   Int_t      GetNSegZSol()                                const {return fNSegZSol;}
  61   Float_t*     GetSegZSol() const {return fSegZSol;}
  62   //
  63   Int_t      GetNParamsTPCInt()                           const {return fNParamsTPCInt;}
  64   Int_t      GetNSegZTPCInt()                             const {return fNSegZTPCInt;}
  65   //
  66   Int_t      GetNParamsDip()                              const {return fNParamsDip;}
  67   Int_t      GetNSegZDip()                                const {return fNZSegDip;}
  68   //
  69   //
  70   Float_t    GetMinZSol()                                 const {return fMinZSol;}
  71   Float_t    GetMaxZSol()                                 const {return fMaxZSol;}
  72   Float_t    GetMaxRSol()                                 const {return fMaxRSol;}
  73   //
  74   Float_t    GetMinZDip()                                 const {return fMinZDip;}
  75   Float_t    GetMaxZDip()                                 const {return fMaxZDip;}
  76   //
  77   Float_t    GetMinZTPCInt()                              const {return fMinZTPCInt;}
  78   Float_t    GetMaxZTPCInt()                              const {return fMaxZTPCInt;}
  79   Float_t    GetMaxRTPCInt()                              const {return fMaxRTPCInt;}
  80   //
  81   AliCheb3D* GetParamSol(Int_t ipar)                      const {return (AliCheb3D*)fParamsSol->UncheckedAt(ipar);}
  82   AliCheb3D* GetParamTPCInt(Int_t ipar)                   const {return (AliCheb3D*)fParamsTPCInt->UncheckedAt(ipar);}
  83   AliCheb3D* GetParamDip(Int_t ipar)                      const {return (AliCheb3D*)fParamsDip->UncheckedAt(ipar);}
  84   //
  85   virtual void Print(Option_t * = "")                     const;
  86   //
  87   virtual void Field(const Float_t *xyz, Float_t *b)                const;
  88   virtual void Field(const Double_t *xyz, Double_t *b)              const;
  89   //
  90   virtual void FieldCyl(const Float_t *rphiz, Float_t *b)     const;
  91   virtual void FieldCyl(const Double_t *rphiz, Double_t *b)   const;
  92   //
  93   virtual void GetTPCInt(const Float_t *xyz, Float_t *b)        const;
  94   virtual void GetTPCIntCyl(const Float_t *rphiz, Float_t *b)   const;
  95   //
  96   template <class T>
  97     Int_t      FindDipSegment(const T *xyz)               const;
  98   //
  99   template <class T>
 100     static void CylToCartCylB(const T *rphiz, const T *brphiz,T *bxyz);
 101   template <class T>
 102     static void CylToCartCartB(const T *xyz,  const T *brphiz,T *bxyz);
 103   template <class T>
 104     static void CartToCylCartB(const T *xyz,  const T *bxyz,  T *brphiz);
 105   template <class T>
 106     static void CartToCylCylB(const T *rphiz, const T *bxyz,  T *brphiz);
 107   template <class T>
 108     static void CartToCyl(const T *xyz,  T *rphiz);
 109   template <class T>
 110     static void CylToCart(const T *rphiz,T *xyz);
 111   //
 112 #ifdef  _INC_CREATION_ALICHEB3D_                          // see AliCheb3D.h for explanation
 113   void         LoadData(const char* inpfile);
 114   //
 115   AliMagWrapCheb(const char* inputFile);
 116   void       SaveData(const char* outfile)              const;
 117   Int_t      SegmentDipDimension(Float_t** seg,const TObjArray* par,int npar, int dim,
 118                                  Float_t xmn,Float_t xmx,Float_t ymn,Float_t ymx,Float_t zmn,Float_t zmx);
 119   //
 120   void       AddParamSol(const AliCheb3D* param);
 121   void       AddParamTPCInt(const AliCheb3D* param);
 122   void       AddParamDip(const AliCheb3D* param);
 123   void       BuildTableDip();
 124   void       BuildTableSol();
 125   void       BuildTableTPCInt();
 126   void       ResetTPCInt();
 127   //
 128   //
 129 #endif
 130   //
 131  protected:
 132     virtual void FieldCylSol(const Float_t *rphiz, Float_t *b)      const;
 133     virtual void FieldCylSol(const Double_t *rphiz, Double_t *b)    const;
 134   //
 135  protected:
 136   //
 137   Int_t      fNParamsSol;            // Total number of parameterization pieces for Sol
 138   Int_t      fNSegZSol;              // Number of segments in Z for Solenoid field
 139   //
 140   Int_t      fNParamsTPCInt;         // Total number of parameterization pieces for TPC field integral
 141   Int_t      fNSegZTPCInt;           // Number of segments in Z for TPC field integral
 142   //
 143   Int_t      fNParamsDip;            // Total number of parameterization pieces for dipole
 144   Int_t      fNZSegDip;              // number of distinct Z segments in Dipole
 145   Int_t      fNYSegDip;              // number of distinct Y segments in Dipole
 146   Int_t      fNXSegDip;              // number of distinct X segments in Dipole
 147   //
 148   Float_t*   fSegZSol;               //[fNSegZSol]      upper boundaries of Z segments
 149   Float_t*   fSegRSol;               //[fNParamsSol]    upper boundaries of R segments
 150   //
 151   Float_t*   fSegZTPCInt;            //[fNSegZTPCInt]    upper boundaries of Z segments
 152   Float_t*   fSegRTPCInt;            //[fNParamsTPCInt]  upper boundaries of R segments
 153   //
 154   Float_t*   fSegZDip;               //[fNZSegDip] coordinates of distinct Z segments in Dipole
 155   Float_t*   fSegYDip;               //[fNYSegDip] coordinated of Y segments for each Zsegment in Dipole
 156   Float_t*   fSegXDip;               //[fNXSegDip] coordinated of X segments for each Ysegment in Dipole
 157   //
 158   Int_t*     fNSegRSol;              //[fNSegZSol]      number of R segments for each Z segment
 159   Int_t*     fSegZIdSol;             //[fNSegZSol]      Id of the first R segment of each Z segment in the fSegRSol...
 160   //
 161   Int_t*     fNSegRTPCInt;           //[fNSegZTPCInt]   number of R segments for each Z segment
 162   Int_t*     fSegZIdTPCInt;          //[fNSegZTPCInt]   Id of the first R segment of each Z segment in the fSegRTPCInt...
 163   //
 164   Int_t*     fBegSegYDip;            //[fNZSegDip] beginning of Y segments array for each Z segment
 165   Int_t*     fNSegYDip;              //[fNZSegDip] number of Y segments for each Z segment
 166   Int_t*     fBegSegXDip;            //[fNYSegDip] beginning of X segments array for each Y segment
 167   Int_t*     fNSegXDip;              //[fNYSegDip] number of X segments for each Y segment
 168   Int_t*     fSegIDDip;              //[fNXSegDip] ID of the dipole parameterization for given XYZ segment
 169   //
 170   Float_t    fMinZSol;               // Min Z of Sol parameterization (in CYL. coordinates)
 171   Float_t    fMaxZSol;               // Max Z of Sol parameterization (in CYL. coordinates)
 172   Float_t    fMaxRSol;               // Max R of Sol parameterization (in CYL. coordinates)
 173   //
 174   Float_t    fMinZDip;               // Min Z of Dipole parameterization
 175   Float_t    fMaxZDip;               // Max Z of Dipole parameterization
 176   //
 177   Float_t    fMinZTPCInt;            // Min Z of TPCInt parameterization (in CYL. coordinates)
 178   Float_t    fMaxZTPCInt;            // Max Z of TPCInt parameterization (in CYL. coordinates)
 179   Float_t    fMaxRTPCInt;            // Max R of TPCInt parameterization (in CYL. coordinates)
 180   //
 181   TObjArray* fParamsSol;             // Parameterization pieces for Solenoid field
 182   TObjArray* fParamsDip;             // Parameterization pieces for Dipole field
 183   TObjArray* fParamsTPCInt;          // Parameterization pieces for Solenoid field integrals in TPC region
 184   //
 185   ClassDef(AliMagWrapCheb,3)            // Wrapper class for the set of Chebishev parameterizations of Alice mag.field
 186   //
 187  };
 188
 189
 190 //__________________________________________________________________________________________
 191 inline void AliMagWrapCheb::FieldCyl(const Float_t *rphiz, Float_t *b) const
 192 {
 193   // compute field in Cylindircal coordinates
 194   //  if (rphiz[2]<GetMinZSol() || rphiz[2]>GetMaxZSol() || rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;}
 195   FieldCylSol(rphiz,b);
 196 }
 197
 198
 199 //__________________________________________________________________________________________
 200 inline void AliMagWrapCheb::FieldCyl(const Double_t *rphiz, Double_t *b) const
 201 {
 202   // compute field in Cylindircal coordinates
 203   //  if (rphiz[2]<GetMinZSol() || rphiz[2]>GetMaxZSol() || rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;}
 204   FieldCylSol(rphiz,b);
 205 }
 206
 207 //__________________________________________________________________________________________________
 208 template <class T>
 209 inline void AliMagWrapCheb::CylToCartCylB(const T *rphiz, const T *brphiz,T *bxyz)
 210 {
 211   // convert field in cylindrical coordinates to cartesian system, point is in cyl.system
 212   T btr = TMath::Sqrt(brphiz[0]*brphiz[0]+brphiz[1]*brphiz[1]);
 213   T psiPLUSphi = TMath::ATan2(brphiz[1],brphiz[0]) + rphiz[1];
 214   bxyz[0] = btr*TMath::Cos(psiPLUSphi);
 215   bxyz[1] = btr*TMath::Sin(psiPLUSphi);
 216   bxyz[2] = brphiz[2];
 217   //
 218 }
 219
 220 //__________________________________________________________________________________________________
 221 template <class T>
 222 inline void AliMagWrapCheb::CylToCartCartB(const T *xyz, const T *brphiz, T *bxyz)
 223 {
 224   // convert field in cylindrical coordinates to cartesian system, point is in cart.system
 225   T btr = TMath::Sqrt(brphiz[0]*brphiz[0]+brphiz[1]*brphiz[1]);
 226   T phiPLUSpsi = TMath::ATan2(xyz[1],xyz[0]) +  TMath::ATan2(brphiz[1],brphiz[0]);
 227   bxyz[0] = btr*TMath::Cos(phiPLUSpsi);
 228   bxyz[1] = btr*TMath::Sin(phiPLUSpsi);
 229   bxyz[2] = brphiz[2];
 230   //
 231 }
 232
 233 //__________________________________________________________________________________________________
 234 template <class T>
 235 inline void AliMagWrapCheb::CartToCylCartB(const T *xyz, const T *bxyz, T *brphiz)
 236 {
 237   // convert field in cylindrical coordinates to cartesian system, poin is in cart.system
 238   T btr = TMath::Sqrt(bxyz[0]*bxyz[0]+bxyz[1]*bxyz[1]);
 239   T psiMINphi = TMath::ATan2(bxyz[1],bxyz[0]) - TMath::ATan2(xyz[1],xyz[0]);
 240   //
 241   brphiz[0] = btr*TMath::Cos(psiMINphi);
 242   brphiz[1] = btr*TMath::Sin(psiMINphi);
 243   brphiz[2] = bxyz[2];
 244   //
 245 }
 246
 247 //__________________________________________________________________________________________________
 248 template <class T>
 249 inline void AliMagWrapCheb::CartToCylCylB(const T *rphiz, const T *bxyz, T *brphiz)
 250 {
 251   // convert field in cylindrical coordinates to cartesian system, point is in cyl.system
 252   T btr = TMath::Sqrt(bxyz[0]*bxyz[0]+bxyz[1]*bxyz[1]);
 253   T psiMINphi =  TMath::ATan2(bxyz[1],bxyz[0]) - rphiz[1];
 254   brphiz[0] = btr*TMath::Cos(psiMINphi);
 255   brphiz[1] = btr*TMath::Sin(psiMINphi);
 256   brphiz[2] = bxyz[2];
 257   //
 258 }
 259
 260 //__________________________________________________________________________________________________
 261 template <class T>
 262 inline void AliMagWrapCheb::CartToCyl(const T *xyz,T *rphiz)
 263 {
 264   rphiz[0] = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
 265   rphiz[1] = TMath::ATan2(xyz[1],xyz[0]);
 266   rphiz[2] = xyz[2];
 267 }
 268
 269 //__________________________________________________________________________________________________
 270 template <class T>
 271 inline void AliMagWrapCheb::CylToCart(const T *rphiz, T *xyz)
 272 {
 273   xyz[0] = rphiz[0]*TMath::Cos(rphiz[1]);
 274   xyz[1] = rphiz[0]*TMath::Sin(rphiz[1]);
 275   xyz[2] = rphiz[2];
 276 }
 277
 278 //__________________________________________________________________________________________________
 279 template <class T>
 280 Int_t    AliMagWrapCheb::FindDipSegment(const T *xyz) const
 281 {
 282   // find the segment containing point xyz. If it is outside find the closest segment
 283   int xid,yid,zid = TMath::BinarySearch(fNZSegDip,fSegZDip,(Float_t)xyz[2]); // find zsegment
 284   int ysegBeg = fBegSegYDip[zid];
 285   //
 286   for (yid=0;yid<fNSegYDip[zid];yid++) if (xyz[1]<fSegYDip[ysegBeg+yid]) break;
 287   if ( --yid < 0 ) yid = 0;
 288   yid +=  ysegBeg;
 289   //
 290   int xsegBeg = fBegSegXDip[yid];
 291   for (xid=0;xid<fNSegXDip[yid];xid++) if (xyz[0]<fSegXDip[xsegBeg+xid]) break;
 292   if ( --xid < 0) xid = 0;
 293   xid +=  xsegBeg;
 294   //
 295   return fSegIDDip[xid];
 296 }
 297
 298 #endif