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