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