STEER/STEERBase/AliMagWrapCheb.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 #include "AliMagWrapCheb.h"
  17 #include "AliLog.h"
  18 #include <TSystem.h>
  19 #include <TArrayF.h>
  20 #include <TArrayI.h>
  21
  22 ClassImp(AliMagWrapCheb)
  23
  24 //__________________________________________________________________________________________
  25 AliMagWrapCheb::AliMagWrapCheb() :
  26 fNParamsSol(0),fNZSegSol(0),fNPSegSol(0),fNRSegSol(0),
  27   fSegZSol(0),fSegPSol(0),fSegRSol(0),
  28   fBegSegPSol(0),fNSegPSol(0),fBegSegRSol(0),fNSegRSol(0),fSegIDSol(0),fMinZSol(1.e6),fMaxZSol(-1.e6),fParamsSol(0),fMaxRSol(0),
  29 //
  30   fNParamsTPC(0),fNZSegTPC(0),fNPSegTPC(0),fNRSegTPC(0),
  31   fSegZTPC(0),fSegPTPC(0),fSegRTPC(0),
  32   fBegSegPTPC(0),fNSegPTPC(0),fBegSegRTPC(0),fNSegRTPC(0),fSegIDTPC(0),fMinZTPC(1.e6),fMaxZTPC(-1.e6),fParamsTPC(0),fMaxRTPC(0),
  33 //
  34   fNParamsTPCRat(0),fNZSegTPCRat(0),fNPSegTPCRat(0),fNRSegTPCRat(0),
  35   fSegZTPCRat(0),fSegPTPCRat(0),fSegRTPCRat(0),
  36   fBegSegPTPCRat(0),fNSegPTPCRat(0),fBegSegRTPCRat(0),fNSegRTPCRat(0),fSegIDTPCRat(0),fMinZTPCRat(1.e6),fMaxZTPCRat(-1.e6),fParamsTPCRat(0),fMaxRTPCRat(0),
  37 //
  38   fNParamsDip(0),fNZSegDip(0),fNYSegDip(0),fNXSegDip(0),
  39   fSegZDip(0),fSegYDip(0),fSegXDip(0),
  40   fBegSegYDip(0),fNSegYDip(0),fBegSegXDip(0),fNSegXDip(0),fSegIDDip(0),fMinZDip(1.e6),fMaxZDip(-1.e6),fParamsDip(0)
  41 //
  42 #ifdef _MAGCHEB_CACHE_
  43   ,fCacheSol(0),fCacheDip(0),fCacheTPCInt(0),fCacheTPCRat(0)
  44 #endif
  45 //
  46 {
  47   // default constructor
  48 }
  49
  50 //__________________________________________________________________________________________
  51 AliMagWrapCheb::AliMagWrapCheb(const AliMagWrapCheb& src) :
  52   TNamed(src),
  53   fNParamsSol(0),fNZSegSol(0),fNPSegSol(0),fNRSegSol(0),
  54   fSegZSol(0),fSegPSol(0),fSegRSol(0),
  55   fBegSegPSol(0),fNSegPSol(0),fBegSegRSol(0),fNSegRSol(0),fSegIDSol(0),fMinZSol(1.e6),fMaxZSol(-1.e6),fParamsSol(0),fMaxRSol(0),
  56 //
  57   fNParamsTPC(0),fNZSegTPC(0),fNPSegTPC(0),fNRSegTPC(0),
  58   fSegZTPC(0),fSegPTPC(0),fSegRTPC(0),
  59   fBegSegPTPC(0),fNSegPTPC(0),fBegSegRTPC(0),fNSegRTPC(0),fSegIDTPC(0),fMinZTPC(1.e6),fMaxZTPC(-1.e6),fParamsTPC(0),fMaxRTPC(0),
  60 //
  61   fNParamsTPCRat(0),fNZSegTPCRat(0),fNPSegTPCRat(0),fNRSegTPCRat(0),
  62   fSegZTPCRat(0),fSegPTPCRat(0),fSegRTPCRat(0),
  63   fBegSegPTPCRat(0),fNSegPTPCRat(0),fBegSegRTPCRat(0),fNSegRTPCRat(0),fSegIDTPCRat(0),fMinZTPCRat(1.e6),fMaxZTPCRat(-1.e6),fParamsTPCRat(0),fMaxRTPCRat(0),
  64 //
  65   fNParamsDip(0),fNZSegDip(0),fNYSegDip(0),fNXSegDip(0),
  66   fSegZDip(0),fSegYDip(0),fSegXDip(0),
  67   fBegSegYDip(0),fNSegYDip(0),fBegSegXDip(0),fNSegXDip(0),fSegIDDip(0),fMinZDip(1.e6),fMaxZDip(-1.e6),fParamsDip(0)
  68 //
  69 #ifdef _MAGCHEB_CACHE_
  70   ,fCacheSol(0),fCacheDip(0),fCacheTPCInt(0),fCacheTPCRat(0)
  71 #endif
  72 {
  73   // copy constructor
  74   CopyFrom(src);
  75   //
  76 }
  77
  78 //__________________________________________________________________________________________
  79 void AliMagWrapCheb::CopyFrom(const AliMagWrapCheb& src)
  80 {
  81   // copy method
  82   Clear();
  83   SetName(src.GetName());
  84   SetTitle(src.GetTitle());
  85   //
  86   fNParamsSol    = src.fNParamsSol;
  87   fNZSegSol      = src.fNZSegSol;
  88   fNPSegSol      = src.fNPSegSol;
  89   fNRSegSol      = src.fNRSegSol;
  90   fMinZSol       = src.fMinZSol;
  91   fMaxZSol       = src.fMaxZSol;
  92   fMaxRSol       = src.fMaxRSol;
  93   if (src.fNParamsSol) {
  94     memcpy(fSegZSol   = new Float_t[fNZSegSol], src.fSegZSol, sizeof(Float_t)*fNZSegSol);
  95     memcpy(fSegPSol   = new Float_t[fNPSegSol], src.fSegPSol, sizeof(Float_t)*fNPSegSol);
  96     memcpy(fSegRSol   = new Float_t[fNRSegSol], src.fSegRSol, sizeof(Float_t)*fNRSegSol);
  97     memcpy(fBegSegPSol= new Int_t[fNZSegSol], src.fBegSegPSol, sizeof(Int_t)*fNZSegSol);
  98     memcpy(fNSegPSol  = new Int_t[fNZSegSol], src.fNSegPSol, sizeof(Int_t)*fNZSegSol);
  99     memcpy(fBegSegRSol= new Int_t[fNPSegSol], src.fBegSegRSol, sizeof(Int_t)*fNPSegSol);
 100     memcpy(fNSegRSol  = new Int_t[fNPSegSol], src.fNSegRSol, sizeof(Int_t)*fNPSegSol);
 101     memcpy(fSegIDSol  = new Int_t[fNRSegSol], src.fSegIDSol, sizeof(Int_t)*fNRSegSol);
 102     fParamsSol        = new TObjArray(fNParamsSol);
 103     for (int i=0;i<fNParamsSol;i++) fParamsSol->AddAtAndExpand(new AliCheb3D(*src.GetParamSol(i)),i);
 104   }
 105   //
 106   fNParamsTPC    = src.fNParamsTPC;
 107   fNZSegTPC      = src.fNZSegTPC;
 108   fNPSegTPC      = src.fNPSegTPC;
 109   fNRSegTPC      = src.fNRSegTPC;
 110   fMinZTPC       = src.fMinZTPC;
 111   fMaxZTPC       = src.fMaxZTPC;
 112   fMaxRTPC       = src.fMaxRTPC;
 113   if (src.fNParamsTPC) {
 114     memcpy(fSegZTPC   = new Float_t[fNZSegTPC], src.fSegZTPC, sizeof(Float_t)*fNZSegTPC);
 115     memcpy(fSegPTPC   = new Float_t[fNPSegTPC], src.fSegPTPC, sizeof(Float_t)*fNPSegTPC);
 116     memcpy(fSegRTPC   = new Float_t[fNRSegTPC], src.fSegRTPC, sizeof(Float_t)*fNRSegTPC);
 117     memcpy(fBegSegPTPC= new Int_t[fNZSegTPC], src.fBegSegPTPC, sizeof(Int_t)*fNZSegTPC);
 118     memcpy(fNSegPTPC  = new Int_t[fNZSegTPC], src.fNSegPTPC, sizeof(Int_t)*fNZSegTPC);
 119     memcpy(fBegSegRTPC= new Int_t[fNPSegTPC], src.fBegSegRTPC, sizeof(Int_t)*fNPSegTPC);
 120     memcpy(fNSegRTPC  = new Int_t[fNPSegTPC], src.fNSegRTPC, sizeof(Int_t)*fNPSegTPC);
 121     memcpy(fSegIDTPC  = new Int_t[fNRSegTPC], src.fSegIDTPC, sizeof(Int_t)*fNRSegTPC);
 122     fParamsTPC        = new TObjArray(fNParamsTPC);
 123     for (int i=0;i<fNParamsTPC;i++) fParamsTPC->AddAtAndExpand(new AliCheb3D(*src.GetParamTPCInt(i)),i);
 124   }
 125   //
 126   fNParamsTPCRat    = src.fNParamsTPCRat;
 127   fNZSegTPCRat      = src.fNZSegTPCRat;
 128   fNPSegTPCRat      = src.fNPSegTPCRat;
 129   fNRSegTPCRat      = src.fNRSegTPCRat;
 130   fMinZTPCRat       = src.fMinZTPCRat;
 131   fMaxZTPCRat       = src.fMaxZTPCRat;
 132   fMaxRTPCRat       = src.fMaxRTPCRat;
 133   if (src.fNParamsTPCRat) {
 134     memcpy(fSegZTPCRat   = new Float_t[fNZSegTPCRat], src.fSegZTPCRat, sizeof(Float_t)*fNZSegTPCRat);
 135     memcpy(fSegPTPCRat   = new Float_t[fNPSegTPCRat], src.fSegPTPCRat, sizeof(Float_t)*fNPSegTPCRat);
 136     memcpy(fSegRTPCRat   = new Float_t[fNRSegTPCRat], src.fSegRTPCRat, sizeof(Float_t)*fNRSegTPCRat);
 137     memcpy(fBegSegPTPCRat= new Int_t[fNZSegTPCRat], src.fBegSegPTPCRat, sizeof(Int_t)*fNZSegTPCRat);
 138     memcpy(fNSegPTPCRat  = new Int_t[fNZSegTPCRat], src.fNSegPTPCRat, sizeof(Int_t)*fNZSegTPCRat);
 139     memcpy(fBegSegRTPCRat= new Int_t[fNPSegTPCRat], src.fBegSegRTPCRat, sizeof(Int_t)*fNPSegTPCRat);
 140     memcpy(fNSegRTPCRat  = new Int_t[fNPSegTPCRat], src.fNSegRTPCRat, sizeof(Int_t)*fNPSegTPCRat);
 141     memcpy(fSegIDTPCRat  = new Int_t[fNRSegTPCRat], src.fSegIDTPCRat, sizeof(Int_t)*fNRSegTPCRat);
 142     fParamsTPCRat        = new TObjArray(fNParamsTPCRat);
 143     for (int i=0;i<fNParamsTPCRat;i++) fParamsTPCRat->AddAtAndExpand(new AliCheb3D(*src.GetParamTPCRatInt(i)),i);
 144   }
 145   //
 146   fNParamsDip    = src.fNParamsDip;
 147   fNZSegDip      = src.fNZSegDip;
 148   fNYSegDip      = src.fNYSegDip;
 149   fNXSegDip      = src.fNXSegDip;
 150   fMinZDip       = src.fMinZDip;
 151   fMaxZDip       = src.fMaxZDip;
 152   if (src.fNParamsDip) {
 153     memcpy(fSegZDip   = new Float_t[fNZSegDip], src.fSegZDip, sizeof(Float_t)*fNZSegDip);
 154     memcpy(fSegYDip   = new Float_t[fNYSegDip], src.fSegYDip, sizeof(Float_t)*fNYSegDip);
 155     memcpy(fSegXDip   = new Float_t[fNXSegDip], src.fSegXDip, sizeof(Float_t)*fNXSegDip);
 156     memcpy(fBegSegYDip= new Int_t[fNZSegDip], src.fBegSegYDip, sizeof(Int_t)*fNZSegDip);
 157     memcpy(fNSegYDip  = new Int_t[fNZSegDip], src.fNSegYDip, sizeof(Int_t)*fNZSegDip);
 158     memcpy(fBegSegXDip= new Int_t[fNYSegDip], src.fBegSegXDip, sizeof(Int_t)*fNYSegDip);
 159     memcpy(fNSegXDip  = new Int_t[fNYSegDip], src.fNSegXDip, sizeof(Int_t)*fNYSegDip);
 160     memcpy(fSegIDDip  = new Int_t[fNXSegDip], src.fSegIDDip, sizeof(Int_t)*fNXSegDip);
 161     fParamsDip        = new TObjArray(fNParamsDip);
 162     for (int i=0;i<fNParamsDip;i++) fParamsDip->AddAtAndExpand(new AliCheb3D(*src.GetParamDip(i)),i);
 163   }
 164   //
 165 }
 166
 167 //__________________________________________________________________________________________
 168 AliMagWrapCheb& AliMagWrapCheb::operator=(const AliMagWrapCheb& rhs)
 169 {
 170   // assignment
 171   if (this != &rhs) {
 172     Clear();
 173     CopyFrom(rhs);
 174   }
 175   return *this;
 176   //
 177 }
 178
 179 //__________________________________________________________________________________________
 180 void AliMagWrapCheb::Clear(const Option_t *)
 181 {
 182   // clear all dynamic parts
 183   if (fNParamsSol) {
 184     fParamsSol->SetOwner(kTRUE);
 185     delete   fParamsSol;  fParamsSol = 0;
 186     delete[] fSegZSol;    fSegZSol   = 0;
 187     delete[] fSegPSol;    fSegPSol   = 0;
 188     delete[] fSegRSol;    fSegRSol   = 0;
 189     delete[] fBegSegPSol; fBegSegPSol = 0;
 190     delete[] fNSegPSol;   fNSegPSol   = 0;
 191     delete[] fBegSegRSol; fBegSegRSol = 0;
 192     delete[] fNSegRSol;   fNSegRSol   = 0;
 193     delete[] fSegIDSol;   fSegIDSol   = 0;
 194   }
 195   fNParamsSol = fNZSegSol = fNPSegSol = fNRSegSol = 0;
 196   fMinZSol = 1e6;
 197   fMaxZSol = -1e6;
 198   fMaxRSol = 0;
 199   //
 200   if (fNParamsTPC) {
 201     fParamsTPC->SetOwner(kTRUE);
 202     delete   fParamsTPC;  fParamsTPC = 0;
 203     delete[] fSegZTPC;    fSegZTPC   = 0;
 204     delete[] fSegPTPC;    fSegPTPC   = 0;
 205     delete[] fSegRTPC;    fSegRTPC   = 0;
 206     delete[] fBegSegPTPC; fBegSegPTPC = 0;
 207     delete[] fNSegPTPC;   fNSegPTPC   = 0;
 208     delete[] fBegSegRTPC; fBegSegRTPC = 0;
 209     delete[] fNSegRTPC;   fNSegRTPC   = 0;
 210     delete[] fSegIDTPC;   fSegIDTPC   = 0;
 211   }
 212   fNParamsTPC = fNZSegTPC = fNPSegTPC = fNRSegTPC = 0;
 213   fMinZTPC = 1e6;
 214   fMaxZTPC = -1e6;
 215   fMaxRTPC = 0;
 216   //
 217   if (fNParamsTPCRat) {
 218     fParamsTPCRat->SetOwner(kTRUE);
 219     delete   fParamsTPCRat;  fParamsTPCRat = 0;
 220     delete[] fSegZTPCRat;    fSegZTPCRat   = 0;
 221     delete[] fSegPTPCRat;    fSegPTPCRat   = 0;
 222     delete[] fSegRTPCRat;    fSegRTPCRat   = 0;
 223     delete[] fBegSegPTPCRat; fBegSegPTPCRat = 0;
 224     delete[] fNSegPTPCRat;   fNSegPTPCRat   = 0;
 225     delete[] fBegSegRTPCRat; fBegSegRTPCRat = 0;
 226     delete[] fNSegRTPCRat;   fNSegRTPCRat   = 0;
 227     delete[] fSegIDTPCRat;   fSegIDTPCRat   = 0;
 228   }
 229   fNParamsTPCRat = fNZSegTPCRat = fNPSegTPCRat = fNRSegTPCRat = 0;
 230   fMinZTPCRat = 1e6;
 231   fMaxZTPCRat = -1e6;
 232   fMaxRTPCRat = 0;
 233   //
 234   if (fNParamsDip) {
 235     fParamsDip->SetOwner(kTRUE);
 236     delete   fParamsDip;  fParamsDip = 0;
 237     delete[] fSegZDip;    fSegZDip   = 0;
 238     delete[] fSegYDip;    fSegYDip   = 0;
 239     delete[] fSegXDip;    fSegXDip   = 0;
 240     delete[] fBegSegYDip; fBegSegYDip = 0;
 241     delete[] fNSegYDip;   fNSegYDip   = 0;
 242     delete[] fBegSegXDip; fBegSegXDip = 0;
 243     delete[] fNSegXDip;   fNSegXDip   = 0;
 244     delete[] fSegIDDip;   fSegIDDip   = 0;
 245   }
 246   fNParamsDip = fNZSegDip = fNYSegDip = fNXSegDip = 0;
 247   fMinZDip = 1e6;
 248   fMaxZDip = -1e6;
 249   //
 250 #ifdef _MAGCHEB_CACHE_
 251   fCacheSol = 0;
 252   fCacheDip = 0;
 253   fCacheTPCInt = 0;
 254   fCacheTPCRat = 0;
 255 #endif
 256   //
 257 }
 258
 259 //__________________________________________________________________________________________
 260 void AliMagWrapCheb::Field(const Double_t *xyz, Double_t *b) const
 261 {
 262   // compute field in cartesian coordinates. If point is outside of the parameterized region
 263   // get it at closest valid point
 264   Double_t rphiz[3];
 265   //
 266 #ifndef _BRING_TO_BOUNDARY_  // exact matching to fitted volume is requested
 267   b[0] = b[1] = b[2] = 0;
 268 #endif
 269   //
 270   if (xyz[2]>fMinZSol) {
 271     CartToCyl(xyz,rphiz);
 272     //
 273 #ifdef _MAGCHEB_CACHE_
 274     if (fCacheSol && fCacheSol->IsInside(rphiz))
 275       fCacheSol->Eval(rphiz,b);
 276     else
 277 #endif //_MAGCHEB_CACHE_
 278       FieldCylSol(rphiz,b);
 279     // convert field to cartesian system
 280     CylToCartCylB(rphiz, b,b);
 281     return;
 282   }
 283   //
 284 #ifdef _MAGCHEB_CACHE_
 285   if (fCacheDip && fCacheDip->IsInside(xyz)) {
 286     fCacheDip->Eval(xyz,b); // check the cache first
 287     return;
 288   }
 289 #else //_MAGCHEB_CACHE_
 290   AliCheb3D* fCacheDip = 0;
 291 #endif //_MAGCHEB_CACHE_
 292   int iddip = FindDipSegment(xyz);
 293   if (iddip>=0) {
 294     fCacheDip = GetParamDip(iddip);
 295     //
 296 #ifndef _BRING_TO_BOUNDARY_
 297     if (!fCacheDip->IsInside(xyz)) return;
 298 #endif //_BRING_TO_BOUNDARY_
 299     //
 300     fCacheDip->Eval(xyz,b);
 301   }
 302   //
 303 }
 304
 305 //__________________________________________________________________________________________
 306 Double_t AliMagWrapCheb::GetBz(const Double_t *xyz) const
 307 {
 308   // compute Bz for the point in cartesian coordinates. If point is outside of the parameterized region
 309   // get it at closest valid point
 310   Double_t rphiz[3];
 311   //
 312   if (xyz[2]>fMinZSol) {
 313     //
 314     CartToCyl(xyz,rphiz);
 315     //
 316 #ifdef _MAGCHEB_CACHE_
 317     if (fCacheSol && fCacheSol->IsInside(rphiz)) return fCacheSol->Eval(rphiz,2);
 318 #endif //_MAGCHEB_CACHE_
 319     return FieldCylSolBz(rphiz);
 320   }
 321   //
 322 #ifdef _MAGCHEB_CACHE_
 323   if (fCacheDip && fCacheDip->IsInside(xyz)) return fCacheDip->Eval(xyz,2); // check the cache first
 324   //
 325 #else //_MAGCHEB_CACHE_
 326   AliCheb3D* fCacheDip = 0;
 327 #endif //_MAGCHEB_CACHE_
 328   //
 329   int iddip = FindDipSegment(xyz);
 330   if (iddip>=0) {
 331     fCacheDip = GetParamDip(iddip);
 332     //
 333 #ifndef _BRING_TO_BOUNDARY_
 334     if (!fCacheDip->IsInside(xyz)) return 0.;
 335 #endif // _BRING_TO_BOUNDARY_
 336     //
 337     return fCacheDip->Eval(xyz,2);
 338   //
 339   }
 340   //
 341   return 0;
 342   //
 343 }
 344
 345
 346 //__________________________________________________________________________________________
 347 void AliMagWrapCheb::Print(Option_t *) const
 348 {
 349   // print info
 350   printf("Alice magnetic field parameterized by Chebyshev polynomials\n");
 351   printf("Segmentation for Solenoid (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZSol,fMaxZSol,fMaxRSol);
 352   //
 353   if (fParamsSol) {
 354     for (int i=0;i<fNParamsSol;i++) {
 355       printf("SOL%4d ",i);
 356       GetParamSol(i)->Print();
 357     }
 358   }
 359   //
 360   printf("Segmentation for TPC field integral (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZTPC,fMaxZTPC,fMaxRTPC);
 361   //
 362   if (fParamsTPC) {
 363     for (int i=0;i<fNParamsTPC;i++) {
 364       printf("TPC%4d ",i);
 365       GetParamTPCInt(i)->Print();
 366     }
 367   }
 368   //
 369   printf("Segmentation for TPC field ratios integral (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZTPCRat,fMaxZTPCRat,fMaxRTPCRat);
 370   //
 371   if (fParamsTPCRat) {
 372     for (int i=0;i<fNParamsTPCRat;i++) {
 373       printf("TPCRat%4d ",i);
 374       GetParamTPCRatInt(i)->Print();
 375     }
 376   }
 377   //
 378   printf("Segmentation for Dipole (%+.2f<Z<%+.2f cm)\n",fMinZDip,fMaxZDip);
 379   if (fParamsDip) {
 380     for (int i=0;i<fNParamsDip;i++) {
 381       printf("DIP%4d ",i);
 382       GetParamDip(i)->Print();
 383     }
 384   }
 385   //
 386 }
 387
 388 //__________________________________________________________________________________________________
 389 Int_t AliMagWrapCheb::FindDipSegment(const Double_t *xyz) const
 390 {
 391   // find the segment containing point xyz. If it is outside find the closest segment
 392   if (!fNParamsDip) return -1;
 393   int xid,yid,zid = TMath::BinarySearch(fNZSegDip,fSegZDip,(Float_t)xyz[2]); // find zsegment
 394   //
 395   Bool_t reCheck = kFALSE;
 396   while(1) {
 397     int ysegBeg = fBegSegYDip[zid];
 398     //
 399     for (yid=0;yid<fNSegYDip[zid];yid++) if (xyz[1]<fSegYDip[ysegBeg+yid]) break;
 400     if ( --yid < 0 ) yid = 0;
 401     yid +=  ysegBeg;
 402     //
 403     int xsegBeg = fBegSegXDip[yid];
 404     for (xid=0;xid<fNSegXDip[yid];xid++) if (xyz[0]<fSegXDip[xsegBeg+xid]) break;
 405     //
 406     if ( --xid < 0) xid = 0;
 407     xid +=  xsegBeg;
 408     //
 409     // to make sure that due to the precision problems we did not pick the next Zbin
 410     if (!reCheck && (xyz[2] - fSegZDip[zid] < 3.e-5) && zid &&
 411         !GetParamDip(fSegIDDip[xid])->IsInside(xyz)) {  // check the previous Z bin
 412       zid--;
 413       reCheck = kTRUE;
 414       continue;
 415     }
 416     break;
 417   }
 418   //  AliInfo(Form("%+.2f %+.2f %+.2f %d %d %d %4d",xyz[0],xyz[1],xyz[2],xid,yid,zid,fSegIDDip[xid]));
 419   return fSegIDDip[xid];
 420 }
 421
 422 //__________________________________________________________________________________________________
 423 Int_t AliMagWrapCheb::FindSolSegment(const Double_t *rpz) const
 424 {
 425   // find the segment containing point xyz. If it is outside find the closest segment
 426   if (!fNParamsSol) return -1;
 427   int rid,pid,zid = TMath::BinarySearch(fNZSegSol,fSegZSol,(Float_t)rpz[2]); // find zsegment
 428   //
 429   Bool_t reCheck = kFALSE;
 430   while(1) {
 431     int psegBeg = fBegSegPSol[zid];
 432     for (pid=0;pid<fNSegPSol[zid];pid++) if (rpz[1]<fSegPSol[psegBeg+pid]) break;
 433     if ( --pid < 0 ) pid = 0;
 434     pid +=  psegBeg;
 435     //
 436     int rsegBeg = fBegSegRSol[pid];
 437     for (rid=0;rid<fNSegRSol[pid];rid++) if (rpz[0]<fSegRSol[rsegBeg+rid]) break;
 438     if ( --rid < 0) rid = 0;
 439     rid +=  rsegBeg;
 440     //
 441     // to make sure that due to the precision problems we did not pick the next Zbin
 442     if (!reCheck && (rpz[2] - fSegZSol[zid] < 3.e-5) && zid &&
 443         !GetParamSol(fSegIDSol[rid])->IsInside(rpz)) {  // check the previous Z bin
 444       zid--;
 445       reCheck = kTRUE;
 446       continue;
 447     }
 448     break;
 449   }
 450   //  AliInfo(Form("%+.2f %+.4f %+.2f %d %d %d %4d",rpz[0],rpz[1],rpz[2],rid,pid,zid,fSegIDSol[rid]));
 451   return fSegIDSol[rid];
 452 }
 453
 454 //__________________________________________________________________________________________________
 455 Int_t AliMagWrapCheb::FindTPCSegment(const Double_t *rpz) const
 456 {
 457   // find the segment containing point xyz. If it is outside find the closest segment
 458   if (!fNParamsTPC) return -1;
 459   int rid,pid,zid = TMath::BinarySearch(fNZSegTPC,fSegZTPC,(Float_t)rpz[2]); // find zsegment
 460   //
 461   Bool_t reCheck = kFALSE;
 462   while(1) {
 463     int psegBeg = fBegSegPTPC[zid];
 464     //
 465     for (pid=0;pid<fNSegPTPC[zid];pid++) if (rpz[1]<fSegPTPC[psegBeg+pid]) break;
 466     if ( --pid < 0 ) pid = 0;
 467     pid +=  psegBeg;
 468     //
 469     int rsegBeg = fBegSegRTPC[pid];
 470     for (rid=0;rid<fNSegRTPC[pid];rid++) if (rpz[0]<fSegRTPC[rsegBeg+rid]) break;
 471     if ( --rid < 0) rid = 0;
 472     rid +=  rsegBeg;
 473     //
 474     // to make sure that due to the precision problems we did not pick the next Zbin
 475     if (!reCheck && (rpz[2] - fSegZTPC[zid] < 3.e-5) && zid &&
 476         !GetParamTPCInt(fSegIDTPC[rid])->IsInside(rpz)) {  // check the previous Z bin
 477       zid--;
 478       reCheck = kTRUE;
 479       continue;
 480     }
 481     break;
 482   }
 483   //  AliInfo(Form("%+.2f %+.4f %+.2f %d %d %d %4d",rpz[0],rpz[1],rpz[2],rid,pid,zid,fSegIDTPC[rid]));
 484   return fSegIDTPC[rid];
 485 }
 486
 487 //__________________________________________________________________________________________________
 488 Int_t AliMagWrapCheb::FindTPCRatSegment(const Double_t *rpz) const
 489 {
 490   // find the segment containing point xyz. If it is outside find the closest segment
 491   if (!fNParamsTPCRat) return -1;
 492   int rid,pid,zid = TMath::BinarySearch(fNZSegTPCRat,fSegZTPCRat,(Float_t)rpz[2]); // find zsegment
 493   //
 494   Bool_t reCheck = kFALSE;
 495   while(1) {
 496     int psegBeg = fBegSegPTPCRat[zid];
 497     //
 498     for (pid=0;pid<fNSegPTPCRat[zid];pid++) if (rpz[1]<fSegPTPCRat[psegBeg+pid]) break;
 499     if ( --pid < 0 ) pid = 0;
 500     pid +=  psegBeg;
 501     //
 502     int rsegBeg = fBegSegRTPCRat[pid];
 503     for (rid=0;rid<fNSegRTPCRat[pid];rid++) if (rpz[0]<fSegRTPCRat[rsegBeg+rid]) break;
 504     if ( --rid < 0) rid = 0;
 505     rid +=  rsegBeg;
 506     //
 507     // to make sure that due to the precision problems we did not pick the next Zbin
 508     if (!reCheck && (rpz[2] - fSegZTPCRat[zid] < 3.e-5) && zid &&
 509         !GetParamTPCRatInt(fSegIDTPCRat[rid])->IsInside(rpz)) {  // check the previous Z bin
 510       zid--;
 511       reCheck = kTRUE;
 512       continue;
 513     }
 514     break;
 515   }
 516   //  AliInfo(Form("%+.2f %+.4f %+.2f %d %d %d %4d",rpz[0],rpz[1],rpz[2],rid,pid,zid,fSegIDTPCRat[rid]));
 517   return fSegIDTPCRat[rid];
 518 }
 519
 520
 521 //__________________________________________________________________________________________
 522 void AliMagWrapCheb::GetTPCInt(const Double_t *xyz, Double_t *b) const
 523 {
 524   // compute TPC region field integral in cartesian coordinates.
 525   // If point is outside of the parameterized region get it at closeset valid point
 526   static Double_t rphiz[3];
 527   //
 528   // TPCInt region
 529   // convert coordinates to cyl system
 530   CartToCyl(xyz,rphiz);
 531 #ifndef _BRING_TO_BOUNDARY_
 532   if ( (rphiz[2]>GetMaxZTPCInt()||rphiz[2]<GetMinZTPCInt()) ||
 533        rphiz[0]>GetMaxRTPCInt()) {for (int i=3;i--;) b[i]=0; return;}
 534 #endif
 535   //
 536   GetTPCIntCyl(rphiz,b);
 537   //
 538   // convert field to cartesian system
 539   CylToCartCylB(rphiz, b,b);
 540   //
 541 }
 542
 543 //__________________________________________________________________________________________
 544 void AliMagWrapCheb::GetTPCRatInt(const Double_t *xyz, Double_t *b) const
 545 {
 546   // compute TPCRat region field integral in cartesian coordinates.
 547   // If point is outside of the parameterized region get it at closeset valid point
 548   static Double_t rphiz[3];
 549   //
 550   // TPCRatInt region
 551   // convert coordinates to cyl system
 552   CartToCyl(xyz,rphiz);
 553 #ifndef _BRING_TO_BOUNDARY_
 554   if ( (rphiz[2]>GetMaxZTPCRatInt()||rphiz[2]<GetMinZTPCRatInt()) ||
 555        rphiz[0]>GetMaxRTPCRatInt()) {for (int i=3;i--;) b[i]=0; return;}
 556 #endif
 557   //
 558   GetTPCRatIntCyl(rphiz,b);
 559   //
 560   // convert field to cartesian system
 561   CylToCartCylB(rphiz, b,b);
 562   //
 563 }
 564
 565 //__________________________________________________________________________________________
 566 void AliMagWrapCheb::FieldCylSol(const Double_t *rphiz, Double_t *b) const
 567 {
 568   // compute Solenoid field in Cylindircal coordinates
 569   // note: if the point is outside the volume get the field in closest parameterized point
 570   int id = FindSolSegment(rphiz);
 571   if (id>=0) {
 572 #ifndef _MAGCHEB_CACHE_
 573     AliCheb3D* fCacheSol = 0;
 574 #endif
 575     fCacheSol = GetParamSol(id);
 576     //
 577 #ifndef _BRING_TO_BOUNDARY_  // exact matching to fitted volume is requested
 578     if (!fCacheSol->IsInside(rphiz)) return;
 579 #endif
 580     fCacheSol->Eval(rphiz,b);
 581   }
 582   //
 583 }
 584
 585 //__________________________________________________________________________________________
 586 Double_t AliMagWrapCheb::FieldCylSolBz(const Double_t *rphiz) const
 587 {
 588   // compute Solenoid field in Cylindircal coordinates
 589   // note: if the point is outside the volume get the field in closest parameterized point
 590   int id = FindSolSegment(rphiz);
 591   if (id<0) return 0.;
 592   //
 593 #ifndef _MAGCHEB_CACHE_
 594   AliCheb3D* fCacheSol = 0;
 595 #endif
 596   //
 597   fCacheSol = GetParamSol(id);
 598 #ifndef _BRING_TO_BOUNDARY_
 599   return fCacheSol->IsInside(rphiz) ? fCacheSol->Eval(rphiz,2) : 0;
 600 #else
 601   return fCacheSol->Eval(rphiz,2);
 602 #endif
 603   //
 604 }
 605
 606 //__________________________________________________________________________________________
 607 void AliMagWrapCheb::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
 608 {
 609   // compute field integral in TPC region in Cylindircal coordinates
 610   // note: the check for the point being inside the parameterized region is done outside
 611   //
 612 #ifdef _MAGCHEB_CACHE_
 613   //
 614   if (fCacheTPCInt && fCacheTPCInt->IsInside(rphiz)) {
 615     fCacheTPCInt->Eval(rphiz,b);
 616     return;
 617   }
 618 #else //_MAGCHEB_CACHE_
 619   AliCheb3D* fCacheTPCInt = 0;
 620 #endif //_MAGCHEB_CACHE_
 621   //
 622   int id = FindTPCSegment(rphiz);
 623   if (id>=0) {
 624     //    if (id>=fNParamsTPC) {
 625     //      b[0] = b[1] = b[2] = 0;
 626     //      AliError(Form("Wrong TPCParam segment %d",id));
 627     //      b[0] = b[1] = b[2] = 0;
 628     //      return;
 629     //    }
 630     fCacheTPCInt = GetParamTPCInt(id);
 631     if (fCacheTPCInt->IsInside(rphiz)) {
 632       fCacheTPCInt->Eval(rphiz,b);
 633       return;
 634     }
 635   }
 636   //
 637   b[0] = b[1] = b[2] = 0;
 638   //
 639 }
 640
 641 //__________________________________________________________________________________________
 642 void AliMagWrapCheb::GetTPCRatIntCyl(const Double_t *rphiz, Double_t *b) const
 643 {
 644   // compute field integral in TPCRat region in Cylindircal coordinates
 645   // note: the check for the point being inside the parameterized region is done outside
 646   //
 647 #ifdef _MAGCHEB_CACHE_
 648   if (fCacheTPCRat && fCacheTPCRat->IsInside(rphiz)) {
 649     fCacheTPCRat->Eval(rphiz,b);
 650     return;
 651   }
 652 #else
 653   AliCheb3D* fCacheTPCRat = 0;
 654 #endif //_MAGCHEB_CACHE_
 655   //
 656   int id = FindTPCRatSegment(rphiz);
 657   if (id>=0) {
 658     //    if (id>=fNParamsTPCRat) {
 659     //      AliError(Form("Wrong TPCRatParam segment %d",id));
 660     //      b[0] = b[1] = b[2] = 0;
 661     //      return;
 662     //    }
 663     fCacheTPCRat = GetParamTPCRatInt(id);
 664     if (fCacheTPCRat->IsInside(rphiz)) {
 665       fCacheTPCRat->Eval(rphiz,b);
 666       return;
 667     }
 668   }
 669   //
 670   b[0] = b[1] = b[2] = 0;
 671   //
 672 }
 673
 674
 675 #ifdef  _INC_CREATION_ALICHEB3D_
 676 //_______________________________________________
 677 void AliMagWrapCheb::LoadData(const char* inpfile)
 678 {
 679   // read coefficients data from the text file
 680   //
 681   TString strf = inpfile;
 682   gSystem->ExpandPathName(strf);
 683   FILE* stream = fopen(strf,"r");
 684   if (!stream) {
 685     printf("Did not find input file %s\n",strf.Data());
 686     return;
 687   }
 688   //
 689   TString buffs;
 690   AliCheb3DCalc::ReadLine(buffs,stream);
 691   if (!buffs.BeginsWith("START")) {
 692     Error("LoadData","Expected: \"START <name>\", found \"%s\"\nStop\n",buffs.Data());
 693     exit(1);
 694   }
 695   if (buffs.First(' ')>0) SetName(buffs.Data()+buffs.First(' ')+1);
 696   //
 697   // Solenoid part    -----------------------------------------------------------
 698   AliCheb3DCalc::ReadLine(buffs,stream);
 699   if (!buffs.BeginsWith("START SOLENOID")) {
 700     Error("LoadData","Expected: \"START SOLENOID\", found \"%s\"\nStop\n",buffs.Data());
 701     exit(1);
 702   }
 703   AliCheb3DCalc::ReadLine(buffs,stream); // nparam
 704   int nparSol = buffs.Atoi();
 705   //
 706   for (int ip=0;ip<nparSol;ip++) {
 707     AliCheb3D* cheb = new AliCheb3D();
 708     cheb->LoadData(stream);
 709     AddParamSol(cheb);
 710   }
 711   //
 712   AliCheb3DCalc::ReadLine(buffs,stream);
 713   if (!buffs.BeginsWith("END SOLENOID")) {
 714     Error("LoadData","Expected \"END SOLENOID\", found \"%s\"\nStop\n",buffs.Data());
 715     exit(1);
 716   }
 717   //
 718   // TPCInt part     -----------------------------------------------------------
 719   AliCheb3DCalc::ReadLine(buffs,stream);
 720   if (!buffs.BeginsWith("START TPCINT")) {
 721     Error("LoadData","Expected: \"START TPCINT\", found \"%s\"\nStop\n",buffs.Data());
 722     exit(1);
 723   }
 724   AliCheb3DCalc::ReadLine(buffs,stream); // nparam
 725   int nparTPCInt = buffs.Atoi();
 726   //
 727   for (int ip=0;ip<nparTPCInt;ip++) {
 728     AliCheb3D* cheb = new AliCheb3D();
 729     cheb->LoadData(stream);
 730     AddParamTPCInt(cheb);
 731   }
 732   //
 733   AliCheb3DCalc::ReadLine(buffs,stream);
 734   if (!buffs.BeginsWith("END TPCINT")) {
 735     Error("LoadData","Expected \"END TPCINT\", found \"%s\"\nStop\n",buffs.Data());
 736     exit(1);
 737   }
 738   //
 739   // TPCRatInt part     -----------------------------------------------------------
 740   AliCheb3DCalc::ReadLine(buffs,stream);
 741   if (!buffs.BeginsWith("START TPCRatINT")) {
 742     Error("LoadData","Expected: \"START TPCRatINT\", found \"%s\"\nStop\n",buffs.Data());
 743     exit(1);
 744   }
 745   AliCheb3DCalc::ReadLine(buffs,stream); // nparam
 746   int nparTPCRatInt = buffs.Atoi();
 747   //
 748   for (int ip=0;ip<nparTPCRatInt;ip++) {
 749     AliCheb3D* cheb = new AliCheb3D();
 750     cheb->LoadData(stream);
 751     AddParamTPCRatInt(cheb);
 752   }
 753   //
 754   AliCheb3DCalc::ReadLine(buffs,stream);
 755   if (!buffs.BeginsWith("END TPCRatINT")) {
 756     Error("LoadData","Expected \"END TPCRatINT\", found \"%s\"\nStop\n",buffs.Data());
 757     exit(1);
 758   }
 759   //
 760   // Dipole part    -----------------------------------------------------------
 761   AliCheb3DCalc::ReadLine(buffs,stream);
 762   if (!buffs.BeginsWith("START DIPOLE")) {
 763     Error("LoadData","Expected: \"START DIPOLE\", found \"%s\"\nStop\n",buffs.Data());
 764     exit(1);
 765   }
 766   AliCheb3DCalc::ReadLine(buffs,stream); // nparam
 767   int nparDip = buffs.Atoi();
 768   //
 769   for (int ip=0;ip<nparDip;ip++) {
 770     AliCheb3D* cheb = new AliCheb3D();
 771     cheb->LoadData(stream);
 772     AddParamDip(cheb);
 773   }
 774   //
 775   AliCheb3DCalc::ReadLine(buffs,stream);
 776   if (!buffs.BeginsWith("END DIPOLE")) {
 777     Error("LoadData","Expected \"END DIPOLE\", found \"%s\"\nStop\n",buffs.Data());
 778     exit(1);
 779   }
 780   //
 781   AliCheb3DCalc::ReadLine(buffs,stream);
 782   if (!buffs.BeginsWith("END DIPOLE") || !buffs.Contains(GetName())) {
 783     Error("LoadData","Expected: \"END DIPOLE\", found \"%s\"\nStop\n",buffs.Data());
 784     exit(1);
 785   }
 786   //
 787   // ---------------------------------------------------------------------------
 788   fclose(stream);
 789   BuildTableSol();
 790   BuildTableDip();
 791   BuildTableTPCInt();
 792   BuildTableTPCRatInt();
 793   //
 794   printf("Loaded magnetic field \"%s\" from %s\n",GetName(),strf.Data());
 795   //
 796 }
 797
 798 //__________________________________________________________________________________________
 799 void AliMagWrapCheb::BuildTableSol()
 800 {
 801   // build lookup table
 802   BuildTable(fNParamsSol,fParamsSol,
 803              fNZSegSol,fNPSegSol,fNRSegSol,
 804              fMinZSol,fMaxZSol,
 805              &fSegZSol,&fSegPSol,&fSegRSol,
 806              &fBegSegPSol,&fNSegPSol,
 807              &fBegSegRSol,&fNSegRSol,
 808              &fSegIDSol);
 809 }
 810
 811 //__________________________________________________________________________________________
 812 void AliMagWrapCheb::BuildTableDip()
 813 {
 814   // build lookup table
 815   BuildTable(fNParamsDip,fParamsDip,
 816              fNZSegDip,fNYSegDip,fNXSegDip,
 817              fMinZDip,fMaxZDip,
 818              &fSegZDip,&fSegYDip,&fSegXDip,
 819              &fBegSegYDip,&fNSegYDip,
 820              &fBegSegXDip,&fNSegXDip,
 821              &fSegIDDip);
 822 }
 823
 824 //__________________________________________________________________________________________
 825 void AliMagWrapCheb::BuildTableTPCInt()
 826 {
 827   // build lookup table
 828   BuildTable(fNParamsTPC,fParamsTPC,
 829              fNZSegTPC,fNPSegTPC,fNRSegTPC,
 830              fMinZTPC,fMaxZTPC,
 831              &fSegZTPC,&fSegPTPC,&fSegRTPC,
 832              &fBegSegPTPC,&fNSegPTPC,
 833              &fBegSegRTPC,&fNSegRTPC,
 834              &fSegIDTPC);
 835 }
 836
 837 //__________________________________________________________________________________________
 838 void AliMagWrapCheb::BuildTableTPCRatInt()
 839 {
 840   // build lookup table
 841   BuildTable(fNParamsTPCRat,fParamsTPCRat,
 842              fNZSegTPCRat,fNPSegTPCRat,fNRSegTPCRat,
 843              fMinZTPCRat,fMaxZTPCRat,
 844              &fSegZTPCRat,&fSegPTPCRat,&fSegRTPCRat,
 845              &fBegSegPTPCRat,&fNSegPTPCRat,
 846              &fBegSegRTPCRat,&fNSegRTPCRat,
 847              &fSegIDTPCRat);
 848 }
 849
 850 #endif
 851
 852 //_______________________________________________
 853 #ifdef  _INC_CREATION_ALICHEB3D_
 854
 855 //__________________________________________________________________________________________
 856 AliMagWrapCheb::AliMagWrapCheb(const char* inputFile) :
 857   fNParamsSol(0),fNZSegSol(0),fNPSegSol(0),fNRSegSol(0),
 858   fSegZSol(0),fSegPSol(0),fSegRSol(0),
 859   fBegSegPSol(0),fNSegPSol(0),fBegSegRSol(0),fNSegRSol(0),fSegIDSol(0),fMinZSol(1.e6),fMaxZSol(-1.e6),fParamsSol(0),fMaxRSol(0),
 860 //
 861   fNParamsTPC(0),fNZSegTPC(0),fNPSegTPC(0),fNRSegTPC(0),
 862   fSegZTPC(0),fSegPTPC(0),fSegRTPC(0),
 863   fBegSegPTPC(0),fNSegPTPC(0),fBegSegRTPC(0),fNSegRTPC(0),fSegIDTPC(0),fMinZTPC(1.e6),fMaxZTPC(-1.e6),fParamsTPC(0),fMaxRTPC(0),
 864 //
 865   fNParamsTPCRat(0),fNZSegTPCRat(0),fNPSegTPCRat(0),fNRSegTPCRat(0),
 866   fSegZTPCRat(0),fSegPTPCRat(0),fSegRTPCRat(0),
 867   fBegSegPTPCRat(0),fNSegPTPCRat(0),fBegSegRTPCRat(0),fNSegRTPCRat(0),fSegIDTPCRat(0),fMinZTPCRat(1.e6),fMaxZTPCRat(-1.e6),fParamsTPCRat(0),fMaxRTPCRat(0),
 868 //
 869   fNParamsDip(0),fNZSegDip(0),fNYSegDip(0),fNXSegDip(0),
 870   fSegZDip(0),fSegYDip(0),fSegXDip(0),
 871   fBegSegYDip(0),fNSegYDip(0),fBegSegXDip(0),fNSegXDip(0),fSegIDDip(0),fMinZDip(1.e6),fMaxZDip(-1.e6),fParamsDip(0)
 872 //
 873 {
 874   // construct from coeffs from the text file
 875   LoadData(inputFile);
 876 }
 877
 878 //__________________________________________________________________________________________
 879 void AliMagWrapCheb::AddParamSol(const AliCheb3D* param)
 880 {
 881   // adds new parameterization piece for Sol
 882   // NOTE: pieces must be added strictly in increasing R then increasing Z order
 883   //
 884   if (!fParamsSol) fParamsSol = new TObjArray();
 885   fParamsSol->Add( (AliCheb3D*)param );
 886   fNParamsSol++;
 887   if (fMaxRSol<param->GetBoundMax(0)) fMaxRSol = param->GetBoundMax(0);
 888   //
 889 }
 890
 891 //__________________________________________________________________________________________
 892 void AliMagWrapCheb::AddParamTPCInt(const AliCheb3D* param)
 893 {
 894   // adds new parameterization piece for TPCInt
 895   // NOTE: pieces must be added strictly in increasing R then increasing Z order
 896   //
 897   if (!fParamsTPC) fParamsTPC = new TObjArray();
 898   fParamsTPC->Add( (AliCheb3D*)param);
 899   fNParamsTPC++;
 900   if (fMaxRTPC<param->GetBoundMax(0)) fMaxRTPC = param->GetBoundMax(0);
 901   //
 902 }
 903
 904 //__________________________________________________________________________________________
 905 void AliMagWrapCheb::AddParamTPCRatInt(const AliCheb3D* param)
 906 {
 907   // adds new parameterization piece for TPCRatInt
 908   // NOTE: pieces must be added strictly in increasing R then increasing Z order
 909   //
 910   if (!fParamsTPCRat) fParamsTPCRat = new TObjArray();
 911   fParamsTPCRat->Add( (AliCheb3D*)param);
 912   fNParamsTPCRat++;
 913   if (fMaxRTPCRat<param->GetBoundMax(0)) fMaxRTPCRat = param->GetBoundMax(0);
 914   //
 915 }
 916
 917 //__________________________________________________________________________________________
 918 void AliMagWrapCheb::AddParamDip(const AliCheb3D* param)
 919 {
 920   // adds new parameterization piece for Dipole
 921   //
 922   if (!fParamsDip) fParamsDip = new TObjArray();
 923   fParamsDip->Add( (AliCheb3D*)param);
 924   fNParamsDip++;
 925   //
 926 }
 927
 928 //__________________________________________________________________________________________
 929 void AliMagWrapCheb::ResetDip()
 930 {
 931   // clean Dip field (used for update)
 932   if (fNParamsDip) {
 933     delete   fParamsDip;  fParamsDip = 0;
 934     delete[] fSegZDip;    fSegZDip   = 0;
 935     delete[] fSegXDip;    fSegXDip   = 0;
 936     delete[] fSegYDip;    fSegYDip   = 0;
 937     delete[] fBegSegYDip; fBegSegYDip = 0;
 938     delete[] fNSegYDip;   fNSegYDip   = 0;
 939     delete[] fBegSegXDip; fBegSegXDip = 0;
 940     delete[] fNSegXDip;   fNSegXDip   = 0;
 941     delete[] fSegIDDip;   fSegIDDip   = 0;
 942   }
 943   fNParamsDip = fNZSegDip = fNXSegDip = fNYSegDip = 0;
 944   fMinZDip = 1e6;
 945   fMaxZDip = -1e6;
 946   //
 947 }
 948
 949 //__________________________________________________________________________________________
 950 void AliMagWrapCheb::ResetSol()
 951 {
 952   // clean Sol field (used for update)
 953   if (fNParamsSol) {
 954     delete   fParamsSol;  fParamsSol = 0;
 955     delete[] fSegZSol;    fSegZSol   = 0;
 956     delete[] fSegPSol;    fSegPSol   = 0;
 957     delete[] fSegRSol;    fSegRSol   = 0;
 958     delete[] fBegSegPSol; fBegSegPSol = 0;
 959     delete[] fNSegPSol;   fNSegPSol   = 0;
 960     delete[] fBegSegRSol; fBegSegRSol = 0;
 961     delete[] fNSegRSol;   fNSegRSol   = 0;
 962     delete[] fSegIDSol;   fSegIDSol   = 0;
 963   }
 964   fNParamsSol = fNZSegSol = fNPSegSol = fNRSegSol = 0;
 965   fMinZSol = 1e6;
 966   fMaxZSol = -1e6;
 967   fMaxRSol = 0;
 968   //
 969 }
 970
 971 //__________________________________________________________________________________________
 972 void AliMagWrapCheb::ResetTPCInt()
 973 {
 974   // clean TPC field integral (used for update)
 975   if (fNParamsTPC) {
 976     delete   fParamsTPC;  fParamsTPC = 0;
 977     delete[] fSegZTPC;    fSegZTPC   = 0;
 978     delete[] fSegPTPC;    fSegPTPC   = 0;
 979     delete[] fSegRTPC;    fSegRTPC   = 0;
 980     delete[] fBegSegPTPC; fBegSegPTPC = 0;
 981     delete[] fNSegPTPC;   fNSegPTPC   = 0;
 982     delete[] fBegSegRTPC; fBegSegRTPC = 0;
 983     delete[] fNSegRTPC;   fNSegRTPC   = 0;
 984     delete[] fSegIDTPC;   fSegIDTPC   = 0;
 985   }
 986   fNParamsTPC = fNZSegTPC = fNPSegTPC = fNRSegTPC = 0;
 987   fMinZTPC = 1e6;
 988   fMaxZTPC = -1e6;
 989   fMaxRTPC = 0;
 990   //
 991 }
 992
 993 //__________________________________________________________________________________________
 994 void AliMagWrapCheb::ResetTPCRatInt()
 995 {
 996   // clean TPCRat field integral (used for update)
 997   if (fNParamsTPCRat) {
 998     delete   fParamsTPCRat;  fParamsTPCRat = 0;
 999     delete[] fSegZTPCRat;    fSegZTPCRat   = 0;
1000     delete[] fSegPTPCRat;    fSegPTPCRat   = 0;
1001     delete[] fSegRTPCRat;    fSegRTPCRat   = 0;
1002     delete[] fBegSegPTPCRat; fBegSegPTPCRat = 0;
1003     delete[] fNSegPTPCRat;   fNSegPTPCRat   = 0;
1004     delete[] fBegSegRTPCRat; fBegSegRTPCRat = 0;
1005     delete[] fNSegRTPCRat;   fNSegRTPCRat   = 0;
1006     delete[] fSegIDTPCRat;   fSegIDTPCRat   = 0;
1007   }
1008   fNParamsTPCRat = fNZSegTPCRat = fNPSegTPCRat = fNRSegTPCRat = 0;
1009   fMinZTPCRat = 1e6;
1010   fMaxZTPCRat = -1e6;
1011   fMaxRTPCRat = 0;
1012   //
1013 }
1014
1015
1016 //__________________________________________________
1017 void AliMagWrapCheb::BuildTable(Int_t npar,TObjArray *parArr, Int_t &nZSeg, Int_t &nYSeg, Int_t &nXSeg,
1018                                 Float_t &minZ,Float_t &maxZ,
1019                                 Float_t **segZ,Float_t **segY,Float_t **segX,
1020                                 Int_t **begSegY,Int_t **nSegY,
1021                                 Int_t **begSegX,Int_t **nSegX,
1022                                 Int_t **segID)
1023 {
1024   // build lookup table for dipole
1025   //
1026   if (npar<1) return;
1027   TArrayF segYArr,segXArr;
1028   TArrayI begSegYDipArr,begSegXDipArr;
1029   TArrayI nSegYDipArr,nSegXDipArr;
1030   TArrayI segIDArr;
1031   float *tmpSegZ,*tmpSegY,*tmpSegX;
1032   //
1033   // create segmentation in Z
1034   nZSeg = SegmentDimension(&tmpSegZ, parArr, npar, 2, 1,-1, 1,-1, 1,-1) - 1;
1035   nYSeg = 0;
1036   nXSeg = 0;
1037   //
1038   // for each Z slice create segmentation in Y
1039   begSegYDipArr.Set(nZSeg);
1040   nSegYDipArr.Set(nZSeg);
1041   float xyz[3];
1042   for (int iz=0;iz<nZSeg;iz++) {
1043     printf("\nZSegment#%d  %+e : %+e\n",iz,tmpSegZ[iz],tmpSegZ[iz+1]);
1044     int ny = SegmentDimension(&tmpSegY, parArr, npar, 1,
1045                               1,-1, 1,-1, tmpSegZ[iz],tmpSegZ[iz+1]) - 1;
1046     segYArr.Set(ny + nYSeg);
1047     for (int iy=0;iy<ny;iy++) segYArr[nYSeg+iy] = tmpSegY[iy];
1048     begSegYDipArr[iz] = nYSeg;
1049     nSegYDipArr[iz] = ny;
1050     printf(" Found %d YSegments, to start from %d\n",ny, begSegYDipArr[iz]);
1051     //
1052     // for each slice in Z and Y create segmentation in X
1053     begSegXDipArr.Set(nYSeg+ny);
1054     nSegXDipArr.Set(nYSeg+ny);
1055     xyz[2] = (tmpSegZ[iz]+tmpSegZ[iz+1])/2.; // mean Z of this segment
1056     //
1057     for (int iy=0;iy<ny;iy++) {
1058       int isg = nYSeg+iy;
1059       printf("\n   YSegment#%d  %+e : %+e\n",iy, tmpSegY[iy],tmpSegY[iy+1]);
1060       int nx = SegmentDimension(&tmpSegX, parArr, npar, 0,
1061                                 1,-1, tmpSegY[iy],tmpSegY[iy+1], tmpSegZ[iz],tmpSegZ[iz+1]) - 1;
1062       //
1063       segXArr.Set(nx + nXSeg);
1064       for (int ix=0;ix<nx;ix++) segXArr[nXSeg+ix] = tmpSegX[ix];
1065       begSegXDipArr[isg] = nXSeg;
1066       nSegXDipArr[isg] = nx;
1067       printf("   Found %d XSegments, to start from %d\n",nx, begSegXDipArr[isg]);
1068       //
1069       segIDArr.Set(nXSeg+nx);
1070       //
1071       // find corresponding params
1072       xyz[1] = (tmpSegY[iy]+tmpSegY[iy+1])/2.; // mean Y of this segment
1073       //
1074       for (int ix=0;ix<nx;ix++) {
1075         xyz[0] = (tmpSegX[ix]+tmpSegX[ix+1])/2.; // mean X of this segment
1076         for (int ipar=0;ipar<npar;ipar++) {
1077           AliCheb3D* cheb = (AliCheb3D*) parArr->At(ipar);
1078           if (!cheb->IsInside(xyz)) continue;
1079           segIDArr[nXSeg+ix] = ipar;
1080           break;
1081         }
1082       }
1083       nXSeg += nx;
1084       //
1085       delete[] tmpSegX;
1086     }
1087     delete[] tmpSegY;
1088     nYSeg += ny;
1089   }
1090   //
1091   minZ = tmpSegZ[0];
1092   maxZ = tmpSegZ[nZSeg];
1093   (*segZ)  = new Float_t[nZSeg];
1094   for (int i=nZSeg;i--;) (*segZ)[i] = tmpSegZ[i];
1095   delete[] tmpSegZ;
1096   //
1097   (*segY)    = new Float_t[nYSeg];
1098   (*segX)    = new Float_t[nXSeg];
1099   (*begSegY) = new Int_t[nZSeg];
1100   (*nSegY)   = new Int_t[nZSeg];
1101   (*begSegX) = new Int_t[nYSeg];
1102   (*nSegX)   = new Int_t[nYSeg];
1103   (*segID)   = new Int_t[nXSeg];
1104   //
1105   for (int i=nYSeg;i--;) (*segY)[i] = segYArr[i];
1106   for (int i=nXSeg;i--;) (*segX)[i] = segXArr[i];
1107   for (int i=nZSeg;i--;) {(*begSegY)[i] = begSegYDipArr[i]; (*nSegY)[i] = nSegYDipArr[i];}
1108   for (int i=nYSeg;i--;) {(*begSegX)[i] = begSegXDipArr[i]; (*nSegX)[i] = nSegXDipArr[i];}
1109   for (int i=nXSeg;i--;) {(*segID)[i]   = segIDArr[i];}
1110   //
1111 }
1112
1113 /*
1114 //__________________________________________________
1115 void AliMagWrapCheb::BuildTableDip()
1116 {
1117   // build lookup table for dipole
1118   //
1119   if (fNParamsDip<1) return;
1120   TArrayF segY,segX;
1121   TArrayI begSegYDip,begSegXDip;
1122   TArrayI nsegYDip,nsegXDip;
1123   TArrayI segID;
1124   float *tmpSegZ,*tmpSegY,*tmpSegX;
1125   //
1126   // create segmentation in Z
1127   fNZSegDip = SegmentDimension(&tmpSegZ, fParamsDip, fNParamsDip, 2, 1,-1, 1,-1, 1,-1) - 1;
1128   fNYSegDip = 0;
1129   fNXSegDip = 0;
1130   //
1131   // for each Z slice create segmentation in Y
1132   begSegYDip.Set(fNZSegDip);
1133   nsegYDip.Set(fNZSegDip);
1134   float xyz[3];
1135   for (int iz=0;iz<fNZSegDip;iz++) {
1136     printf("\nZSegment#%d  %+e : %+e\n",iz,tmpSegZ[iz],tmpSegZ[iz+1]);
1137     int ny = SegmentDimension(&tmpSegY, fParamsDip, fNParamsDip, 1,
1138                                  1,-1, 1,-1, tmpSegZ[iz],tmpSegZ[iz+1]) - 1;
1139     segY.Set(ny + fNYSegDip);
1140     for (int iy=0;iy<ny;iy++) segY[fNYSegDip+iy] = tmpSegY[iy];
1141     begSegYDip[iz] = fNYSegDip;
1142     nsegYDip[iz] = ny;
1143     printf(" Found %d YSegments, to start from %d\n",ny, begSegYDip[iz]);
1144     //
1145     // for each slice in Z and Y create segmentation in X
1146     begSegXDip.Set(fNYSegDip+ny);
1147     nsegXDip.Set(fNYSegDip+ny);
1148     xyz[2] = (tmpSegZ[iz]+tmpSegZ[iz+1])/2.; // mean Z of this segment
1149     //
1150     for (int iy=0;iy<ny;iy++) {
1151       int isg = fNYSegDip+iy;
1152       printf("\n   YSegment#%d  %+e : %+e\n",iy, tmpSegY[iy],tmpSegY[iy+1]);
1153       int nx = SegmentDimension(&tmpSegX, fParamsDip, fNParamsDip, 0,
1154                                 1,-1, tmpSegY[iy],tmpSegY[iy+1], tmpSegZ[iz],tmpSegZ[iz+1]) - 1;
1155       //
1156       segX.Set(nx + fNXSegDip);
1157       for (int ix=0;ix<nx;ix++) segX[fNXSegDip+ix] = tmpSegX[ix];
1158       begSegXDip[isg] = fNXSegDip;
1159       nsegXDip[isg] = nx;
1160       printf("   Found %d XSegments, to start from %d\n",nx, begSegXDip[isg]);
1161       //
1162       segID.Set(fNXSegDip+nx);
1163       //
1164       // find corresponding params
1165       xyz[1] = (tmpSegY[iy]+tmpSegY[iy+1])/2.; // mean Y of this segment
1166       //
1167       for (int ix=0;ix<nx;ix++) {
1168         xyz[0] = (tmpSegX[ix]+tmpSegX[ix+1])/2.; // mean X of this segment
1169         for (int ipar=0;ipar<fNParamsDip;ipar++) {
1170           AliCheb3D* cheb = (AliCheb3D*) fParamsDip->At(ipar);
1171           if (!cheb->IsInside(xyz)) continue;
1172           segID[fNXSegDip+ix] = ipar;
1173           break;
1174         }
1175       }
1176       fNXSegDip += nx;
1177       //
1178       delete[] tmpSegX;
1179     }
1180     delete[] tmpSegY;
1181     fNYSegDip += ny;
1182   }
1183   //
1184   fMinZDip = tmpSegZ[0];
1185   fMaxZDip = tmpSegZ[fNZSegDip];
1186   fSegZDip    = new Float_t[fNZSegDip];
1187   for (int i=fNZSegDip;i--;) fSegZDip[i] = tmpSegZ[i];
1188   delete[] tmpSegZ;
1189   //
1190   fSegYDip    = new Float_t[fNYSegDip];
1191   fSegXDip    = new Float_t[fNXSegDip];
1192   fBegSegYDip = new Int_t[fNZSegDip];
1193   fNSegYDip   = new Int_t[fNZSegDip];
1194   fBegSegXDip = new Int_t[fNYSegDip];
1195   fNSegXDip   = new Int_t[fNYSegDip];
1196   fSegIDDip   = new Int_t[fNXSegDip];
1197   //
1198   for (int i=fNYSegDip;i--;) fSegYDip[i] = segY[i];
1199   for (int i=fNXSegDip;i--;) fSegXDip[i] = segX[i];
1200   for (int i=fNZSegDip;i--;) {fBegSegYDip[i] = begSegYDip[i]; fNSegYDip[i] = nsegYDip[i];}
1201   for (int i=fNYSegDip;i--;) {fBegSegXDip[i] = begSegXDip[i]; fNSegXDip[i] = nsegXDip[i];}
1202   for (int i=fNXSegDip;i--;) {fSegIDDip[i]   = segID[i];}
1203   //
1204 }
1205 */
1206
1207 //________________________________________________________________
1208 void AliMagWrapCheb::SaveData(const char* outfile) const
1209 {
1210   // writes coefficients data to output text file
1211   TString strf = outfile;
1212   gSystem->ExpandPathName(strf);
1213   FILE* stream = fopen(strf,"w+");
1214   //
1215   // Sol part    ---------------------------------------------------------
1216   fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName());
1217   fprintf(stream,"START SOLENOID\n#Number of pieces\n%d\n",fNParamsSol);
1218   for (int ip=0;ip<fNParamsSol;ip++) GetParamSol(ip)->SaveData(stream);
1219   fprintf(stream,"#\nEND SOLENOID\n");
1220   //
1221   // TPCInt part ---------------------------------------------------------
1222   fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName());
1223   fprintf(stream,"START TPCINT\n#Number of pieces\n%d\n",fNParamsTPC);
1224   for (int ip=0;ip<fNParamsTPC;ip++) GetParamTPCInt(ip)->SaveData(stream);
1225   fprintf(stream,"#\nEND TPCINT\n");
1226   //
1227   // TPCRatInt part ---------------------------------------------------------
1228   fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName());
1229   fprintf(stream,"START TPCRatINT\n#Number of pieces\n%d\n",fNParamsTPCRat);
1230   for (int ip=0;ip<fNParamsTPCRat;ip++) GetParamTPCRatInt(ip)->SaveData(stream);
1231   fprintf(stream,"#\nEND TPCRatINT\n");
1232   //
1233   // Dip part   ---------------------------------------------------------
1234   fprintf(stream,"START DIPOLE\n#Number of pieces\n%d\n",fNParamsDip);
1235   for (int ip=0;ip<fNParamsDip;ip++) GetParamDip(ip)->SaveData(stream);
1236   fprintf(stream,"#\nEND DIPOLE\n");
1237   //
1238   fprintf(stream,"#\nEND %s\n",GetName());
1239   //
1240   fclose(stream);
1241   //
1242 }
1243
1244 Int_t AliMagWrapCheb::SegmentDimension(float** seg,const TObjArray* par,int npar, int dim,
1245                                        float xmn,float xmx,float ymn,float ymx,float zmn,float zmx)
1246 {
1247   // find all boundaries in deimension dim for boxes in given region.
1248   // if mn>mx for given projection the check is not done for it.
1249   float *tmpC = new float[2*npar];
1250   int *tmpInd = new int[2*npar];
1251   int nseg0 = 0;
1252   for (int ip=0;ip<npar;ip++) {
1253     AliCheb3D* cheb = (AliCheb3D*) par->At(ip);
1254     if (xmn<xmx && (cheb->GetBoundMin(0)>(xmx+xmn)/2 || cheb->GetBoundMax(0)<(xmn+xmx)/2)) continue;
1255     if (ymn<ymx && (cheb->GetBoundMin(1)>(ymx+ymn)/2 || cheb->GetBoundMax(1)<(ymn+ymx)/2)) continue;
1256     if (zmn<zmx && (cheb->GetBoundMin(2)>(zmx+zmn)/2 || cheb->GetBoundMax(2)<(zmn+zmx)/2)) continue;
1257     //
1258     tmpC[nseg0++] = cheb->GetBoundMin(dim);
1259     tmpC[nseg0++] = cheb->GetBoundMax(dim);
1260   }
1261   // range Dim's boundaries in increasing order
1262   TMath::Sort(nseg0,tmpC,tmpInd,kFALSE);
1263   // count number of really different Z's
1264   int nseg = 0;
1265   float cprev = -1e6;
1266   for (int ip=0;ip<nseg0;ip++) {
1267     if (TMath::Abs(cprev-tmpC[ tmpInd[ip] ])>1e-4) {
1268       cprev = tmpC[ tmpInd[ip] ];
1269       nseg++;
1270     }
1271     else tmpInd[ip] = -1; // supress redundant Z
1272   }
1273   //
1274   *seg  = new float[nseg]; // create final Z segmenations
1275   nseg = 0;
1276   for (int ip=0;ip<nseg0;ip++) if (tmpInd[ip]>=0) (*seg)[nseg++] = tmpC[ tmpInd[ip] ];
1277   //
1278   delete[] tmpC;
1279   delete[] tmpInd;
1280   return nseg;
1281 }
1282
1283 #endif
1284