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),
  26   fNSegZSol(0),
  27   fNParamsTPCInt(0),
  28   fNSegZTPCInt(0),
  29   fNParamsDip(0),
  30 //
  31   fNZSegDip(0),
  32   fNYSegDip(0),
  33   fNXSegDip(0),
  34 //
  35   fSegZSol(0),
  36   fSegRSol(0),
  37 //
  38   fSegZTPCInt(0),
  39   fSegRTPCInt(0),
  40 //
  41   fSegZDip(0),
  42   fSegYDip(0),
  43   fSegXDip(0),
  44 //
  45   fNSegRSol(0),
  46   fSegZIdSol(0),
  47 //
  48   fNSegRTPCInt(0),
  49   fSegZIdTPCInt(0),
  50 //
  51   fBegSegYDip(0),
  52   fNSegYDip(0),
  53   fBegSegXDip(0),
  54   fNSegXDip(0),
  55   fSegIDDip(0),
  56 //
  57   fMinZSol(1e6),
  58   fMaxZSol(-1e6),
  59   fMaxRSol(-1e6),
  60 //
  61   fMinZDip(1e6),
  62   fMaxZDip(-1e6),
  63 //
  64   fMinZTPCInt(1e6),
  65   fMaxZTPCInt(-1e6),
  66   fMaxRTPCInt(-1e6),
  67 //
  68   fParamsSol(0),
  69   fParamsDip(0),
  70   fParamsTPCInt(0)
  71 //
  72 {
  73   // default constructor
  74 }
  75
  76 //__________________________________________________________________________________________
  77 AliMagWrapCheb::AliMagWrapCheb(const AliMagWrapCheb& src) :
  78   TNamed(src),
  79   fNParamsSol(0),
  80   fNSegZSol(0),
  81   fNParamsTPCInt(0),
  82   fNSegZTPCInt(0),
  83   fNParamsDip(0),
  84 //
  85   fNZSegDip(0),
  86   fNYSegDip(0),
  87   fNXSegDip(0),
  88 //
  89   fSegZSol(0),
  90   fSegRSol(0),
  91 //
  92   fSegZTPCInt(0),
  93   fSegRTPCInt(0),
  94 //
  95   fSegZDip(0),
  96   fSegYDip(0),
  97   fSegXDip(0),
  98 //
  99   fNSegRSol(0),
 100   fSegZIdSol(0),
 101 //
 102   fNSegRTPCInt(0),
 103   fSegZIdTPCInt(0),
 104 //
 105   fBegSegYDip(0),
 106   fNSegYDip(0),
 107   fBegSegXDip(0),
 108   fNSegXDip(0),
 109   fSegIDDip(0),
 110 //
 111   fMinZSol(1e6),
 112   fMaxZSol(-1e6),
 113   fMaxRSol(-1e6),
 114 //
 115   fMinZDip(1e6),
 116   fMaxZDip(-1e6),
 117 //
 118   fMinZTPCInt(1e6),
 119   fMaxZTPCInt(-1e6),
 120   fMaxRTPCInt(-1e6),
 121 //
 122   fParamsSol(0),
 123   fParamsDip(0),
 124   fParamsTPCInt(0)
 125 {
 126   // copy constructor
 127   CopyFrom(src);
 128 }
 129
 130 //__________________________________________________________________________________________
 131 void AliMagWrapCheb::CopyFrom(const AliMagWrapCheb& src)
 132 {
 133   // copy method
 134   Clear();
 135   SetName(src.GetName());
 136   SetTitle(src.GetTitle());
 137   fNParamsSol    = src.fNParamsSol;
 138   fNSegZSol      = src.fNSegZSol;
 139   fNParamsTPCInt = src.fNParamsTPCInt;
 140   fNSegZTPCInt   = src.fNSegZTPCInt;
 141   fNParamsDip    = src.fNParamsDip;
 142   //
 143   fNZSegDip      = src.fNZSegDip;
 144   fNYSegDip      = src.fNYSegDip;
 145   fNXSegDip      = src.fNXSegDip;
 146   //
 147   fMinZSol       = src.fMinZSol;
 148   fMaxZSol       = src.fMaxZSol;
 149   fMaxRSol       = src.fMaxRSol;
 150   //
 151   fMinZDip       = src.fMinZDip;
 152   fMaxZDip       = src.fMaxZDip;
 153   //
 154   fMinZTPCInt    = src.fMinZTPCInt;
 155   fMaxZTPCInt    = src.fMaxZTPCInt;
 156   fMaxRTPCInt    = src.fMaxRTPCInt;
 157   //
 158   if (src.fNParamsSol) {
 159     memcpy(fSegZSol  = new Float_t[fNSegZSol], src.fSegZSol, sizeof(Float_t)*fNSegZSol);
 160     memcpy(fSegRSol  = new Float_t[fNParamsSol], src.fSegRSol, sizeof(Float_t)*fNParamsSol);
 161     memcpy(fNSegRSol = new Int_t[fNSegZSol], src.fNSegRSol, sizeof(Int_t)*fNSegZSol);
 162     memcpy(fSegZIdSol= new Int_t[fNSegZSol], src.fSegZIdSol, sizeof(Int_t)*fNSegZSol);
 163     fParamsSol       = new TObjArray(fNParamsSol);
 164     for (int i=0;i<fNParamsSol;i++) fParamsSol->AddAtAndExpand(new AliCheb3D(*src.GetParamSol(i)),i);
 165   }
 166   //
 167   if (src.fNParamsDip) {
 168     memcpy(fSegZDip   = new Float_t[fNZSegDip], src.fSegZDip, sizeof(Float_t)*fNZSegDip);
 169     memcpy(fSegYDip   = new Float_t[fNYSegDip], src.fSegYDip, sizeof(Float_t)*fNYSegDip);
 170     memcpy(fSegXDip   = new Float_t[fNXSegDip], src.fSegZDip, sizeof(Float_t)*fNXSegDip);
 171     memcpy(fBegSegYDip= new Int_t[fNZSegDip], src.fBegSegYDip, sizeof(Int_t)*fNZSegDip);
 172     memcpy(fNSegYDip  = new Int_t[fNZSegDip], src.fNSegYDip, sizeof(Int_t)*fNZSegDip);
 173     memcpy(fBegSegXDip= new Int_t[fNYSegDip], src.fBegSegXDip, sizeof(Int_t)*fNYSegDip);
 174     memcpy(fNSegXDip  = new Int_t[fNYSegDip], src.fNSegXDip, sizeof(Int_t)*fNYSegDip);
 175     memcpy(fSegIDDip  = new Int_t[fNXSegDip], src.fSegIDDip, sizeof(Int_t)*fNXSegDip);
 176     fParamsDip        = new TObjArray(fNParamsDip);
 177     for (int i=0;i<fNParamsDip;i++) fParamsDip->AddAtAndExpand(new AliCheb3D(*src.GetParamDip(i)),i);
 178   }
 179   //
 180   if (src.fNParamsTPCInt) {
 181     memcpy(fSegZTPCInt  = new Float_t[fNSegZTPCInt], src.fSegZTPCInt, sizeof(Float_t)*fNSegZTPCInt);
 182     memcpy(fSegRTPCInt  = new Float_t[fNParamsTPCInt], src.fSegRTPCInt, sizeof(Float_t)*fNParamsTPCInt);
 183     memcpy(fNSegRTPCInt = new Int_t[fNSegZTPCInt], src.fNSegRTPCInt, sizeof(Int_t)*fNSegZTPCInt);
 184     memcpy(fSegZIdTPCInt= new Int_t[fNSegZTPCInt], src.fSegZIdTPCInt, sizeof(Int_t)*fNSegZTPCInt);
 185     fParamsTPCInt       = new TObjArray(fNParamsTPCInt);
 186     for (int i=0;i<fNParamsTPCInt;i++) fParamsTPCInt->AddAtAndExpand(new AliCheb3D(*src.GetParamTPCInt(i)),i);
 187   }
 188   //
 189 }
 190
 191 //__________________________________________________________________________________________
 192 AliMagWrapCheb& AliMagWrapCheb::operator=(const AliMagWrapCheb& rhs)
 193 {
 194   // assignment
 195   if (this != &rhs) {
 196     Clear();
 197     CopyFrom(rhs);
 198   }
 199   return *this;
 200   //
 201 }
 202
 203 //__________________________________________________________________________________________
 204 void AliMagWrapCheb::Clear(const Option_t *)
 205 {
 206   // clear all dynamic parts
 207   if (fNParamsSol) {
 208     delete   fParamsSol;
 209     delete[] fSegZSol;
 210     delete[] fSegRSol;
 211     delete[] fNSegRSol;
 212     delete[] fSegZIdSol;
 213   }
 214   //
 215   if (fNParamsTPCInt) {
 216     delete   fParamsTPCInt;
 217     delete[] fSegZTPCInt;
 218     delete[] fSegRTPCInt;
 219     delete[] fNSegRTPCInt;
 220     delete[] fSegZIdTPCInt;
 221   }
 222   //
 223   if (fNParamsDip) {
 224     delete   fParamsDip;
 225     delete[] fSegZDip;
 226     delete[] fSegYDip;
 227     delete[] fSegXDip;
 228     delete[] fBegSegYDip;
 229     delete[] fNSegYDip;
 230     delete[] fBegSegXDip;
 231     delete[] fNSegXDip;
 232     delete[] fSegIDDip;
 233   }
 234   fNParamsSol = fNParamsTPCInt = fNParamsDip = fNZSegDip = fNYSegDip = fNXSegDip = 0;
 235   fNSegZSol = fNSegZTPCInt = 0;
 236   fMinZSol = fMinZDip = fMinZTPCInt = 1e6;
 237   fMaxZSol = fMaxZDip = fMaxZTPCInt = fMaxRSol = fMaxRTPCInt = -1e6;
 238   //
 239 }
 240
 241 //__________________________________________________________________________________________
 242 void AliMagWrapCheb::Field(const Double_t *xyz, Double_t *b) const
 243 {
 244   // compute field in cartesian coordinates. If point is outside of the parameterized region
 245   // get it at closest valid point
 246   static Double_t rphiz[3];
 247   //
 248 #ifndef _BRING_TO_BOUNDARY_  // exact matching to fitted volume is requested
 249   if ( !(xyz[2]>=GetMinZSol()&&xyz[2]<=GetMaxZSol()) &&
 250        !(xyz[2]>=GetMinZDip()&&xyz[2]<=GetMaxZDip())  ) {for (int i=3;i--;) b[i]=0; return;}
 251 #endif
 252   //
 253   if (xyz[2]<fMaxZDip) {    // dipole part?
 254 #ifndef _BRING_TO_BOUNDARY_
 255     AliCheb3D* par = GetParamDip(FindDipSegment(xyz));
 256     if (par->IsInside(xyz)) {par->Eval(xyz,b); return;}
 257     for (int i=3;i--;) b[i]=0; return;
 258 #else
 259     GetParamDip(FindDipSegment(xyz))->Eval(xyz,b); return;
 260 #endif
 261   }
 262   //
 263   // Sol region: convert coordinates to cyl system
 264   CartToCyl(xyz,rphiz);
 265 #ifndef _BRING_TO_BOUNDARY_
 266   if (rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;}
 267 #endif
 268   //
 269   FieldCylSol(rphiz,b);
 270   //
 271   // convert field to cartesian system
 272   CylToCartCylB(rphiz, b,b);
 273   //
 274 }
 275
 276 //__________________________________________________________________________________________
 277 Double_t AliMagWrapCheb::GetBz(const Double_t *xyz) const
 278 {
 279   // compute Bz for the point in cartesian coordinates. If point is outside of the parameterized region
 280   // get it at closest valid point
 281   static Double_t rphiz[3];
 282   //
 283 #ifndef _BRING_TO_BOUNDARY_  // exact matching to fitted volume is requested
 284   if ( !(xyz[2]>=GetMinZSol()&&xyz[2]<=GetMaxZSol()) &&
 285        !(xyz[2]>=GetMinZDip()&&xyz[2]<=GetMaxZDip())  ) return 0.;
 286 #endif
 287   //
 288   if (xyz[2]<fMaxZDip) {    // dipole part?
 289 #ifndef _BRING_TO_BOUNDARY_
 290     AliCheb3D* par = GetParamDip(FindDipSegment(xyz));
 291     if (par->IsInside(xyz)) return par->Eval(xyz,2);
 292     else return 0.;
 293 #else
 294     return GetParamDip(FindDipSegment(xyz))->Eval(xyz,2);
 295 #endif
 296   }
 297   // Sol region: convert coordinates to cyl system
 298   CartToCyl(xyz,rphiz);
 299 #ifndef _BRING_TO_BOUNDARY_
 300   if (rphiz[0]>GetMaxRSol()) return 0.;
 301 #endif
 302   //
 303   return FieldCylSolBz(rphiz);
 304 }
 305
 306
 307 //__________________________________________________________________________________________
 308 void AliMagWrapCheb::Print(Option_t *) const
 309 {
 310   // print info
 311   printf("Alice magnetic field parameterized by Chebyshev polynomials\n");
 312   printf("Segmentation for Solenoid (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZSol,fMaxZSol,fMaxRSol);
 313   //
 314   if (fParamsSol) fParamsSol->Print();
 315   /*
 316   for (int iz=0;iz<fNSegZSol;iz++) {
 317     AliCheb3D* param = GetParamSol( fSegZIdSol[iz] );
 318     printf("*** Z Segment %2d (%+7.2f<Z<%+7.2f)\t***\n",iz,param->GetBoundMin(2),param->GetBoundMax(2));
 319     for (int ir=0;ir<fNSegRSol[iz];ir++) {
 320       param = GetParamSol( fSegZIdSol[iz]+ir );
 321       printf("    R Segment %2d (%+7.2f<R<%+7.2f, Precision: %.1e) (ID=%2d)\n",ir, param->GetBoundMin(0),
 322              param->GetBoundMax(0),param->GetPrecision(),fSegZIdSol[iz]+ir);
 323     }
 324   }
 325   */
 326   //
 327   printf("Segmentation for TPC field integral (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZTPCInt,fMaxZTPCInt,fMaxRTPCInt);
 328   //
 329   if (fParamsTPCInt) fParamsTPCInt->Print();
 330   /*
 331   for (int iz=0;iz<fNSegZTPCInt;iz++) {
 332     AliCheb3D* param = GetParamTPCInt( fSegZIdTPCInt[iz] );
 333     printf("*** Z Segment %2d (%+7.2f<Z<%+7.2f)\t***\n",iz,param->GetBoundMin(2),param->GetBoundMax(2));
 334     for (int ir=0;ir<fNSegRTPCInt[iz];ir++) {
 335       param = GetParamTPCInt( fSegZIdTPCInt[iz]+ir );
 336       printf("    R Segment %2d (%+7.2f<R<%+7.2f, Precision: %.1e) (ID=%2d)\n",ir, param->GetBoundMin(0),
 337              param->GetBoundMax(0),param->GetPrecision(),fSegZIdTPCInt[iz]+ir);
 338     }
 339   }
 340   */
 341   //
 342   printf("Segmentation for Dipole (%+.2f<Z<%+.2f cm)\n",fMinZDip,fMaxZDip);
 343   if (fParamsDip) fParamsDip->Print();
 344   //
 345
 346
 347 }
 348
 349
 350 //__________________________________________________________________________________________________
 351 Int_t    AliMagWrapCheb::FindDipSegment(const Double_t *xyz) const
 352 {
 353   // find the segment containing point xyz. If it is outside find the closest segment
 354   int xid,yid,zid = TMath::BinarySearch(fNZSegDip,fSegZDip,(Float_t)xyz[2]); // find zsegment
 355   int ysegBeg = fBegSegYDip[zid];
 356   //
 357   for (yid=0;yid<fNSegYDip[zid];yid++) if (xyz[1]<fSegYDip[ysegBeg+yid]) break;
 358   if ( --yid < 0 ) yid = 0;
 359   yid +=  ysegBeg;
 360   //
 361   int xsegBeg = fBegSegXDip[yid];
 362   for (xid=0;xid<fNSegXDip[yid];xid++) if (xyz[0]<fSegXDip[xsegBeg+xid]) break;
 363   if ( --xid < 0) xid = 0;
 364   xid +=  xsegBeg;
 365   //
 366   return fSegIDDip[xid];
 367 }
 368
 369 //__________________________________________________________________________________________
 370 void AliMagWrapCheb::GetTPCInt(const Double_t *xyz, Double_t *b) const
 371 {
 372   // compute TPC region field integral in cartesian coordinates.
 373   // If point is outside of the parameterized region get it at closeset valid point
 374   static Double_t rphiz[3];
 375   //
 376   // TPCInt region
 377   // convert coordinates to cyl system
 378   CartToCyl(xyz,rphiz);
 379 #ifndef _BRING_TO_BOUNDARY_
 380   if ( (rphiz[2]>GetMaxZTPCInt()||rphiz[2]<GetMinZTPCInt()) ||
 381        rphiz[0]>GetMaxRTPCInt()) {for (int i=3;i--;) b[i]=0; return;}
 382 #endif
 383   //
 384   GetTPCIntCyl(rphiz,b);
 385   //
 386   // convert field to cartesian system
 387   CylToCartCylB(rphiz, b,b);
 388   //
 389 }
 390
 391 //__________________________________________________________________________________________
 392 void AliMagWrapCheb::FieldCylSol(const Double_t *rphiz, Double_t *b) const
 393 {
 394   // compute Solenoid field in Cylindircal coordinates
 395   // note: if the point is outside the volume get the field in closest parameterized point
 396   int SolZId = 0;
 397   while (rphiz[2]>fSegZSol[SolZId] && SolZId<fNSegZSol-1) ++SolZId;    // find Z segment
 398   int SolRId = fSegZIdSol[SolZId];        // first R segment for this Z
 399   int SolRMax = SolRId + fNSegRSol[SolZId];
 400   while (rphiz[0]>fSegRSol[SolRId] && SolRId<SolRMax-1) ++SolRId;    // find R segment
 401   GetParamSol( SolRId )->Eval(rphiz,b);
 402   //
 403 }
 404
 405 //__________________________________________________________________________________________
 406 Double_t AliMagWrapCheb::FieldCylSolBz(const Double_t *rphiz) const
 407 {
 408   // compute Solenoid field in Cylindircal coordinates
 409   // note: if the point is outside the volume get the field in closest parameterized point
 410   int SolZId = 0;
 411   while (rphiz[2]>fSegZSol[SolZId] && SolZId<fNSegZSol-1) ++SolZId;    // find Z segment
 412   int SolRId = fSegZIdSol[SolZId];        // first R segment for this Z
 413   int SolRMax = SolRId + fNSegRSol[SolZId];
 414   while (rphiz[0]>fSegRSol[SolRId] && SolRId<SolRMax-1) ++SolRId;    // find R segment
 415   return GetParamSol( SolRId )->Eval(rphiz,2);
 416   //
 417 }
 418
 419 //__________________________________________________________________________________________
 420 void AliMagWrapCheb::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
 421 {
 422   // compute field integral in TPC region in Cylindircal coordinates
 423   // note: the check for the point being inside the parameterized region is done outside
 424   int tpcIntZId = 0;
 425   while (rphiz[2]>fSegZTPCInt[tpcIntZId] && tpcIntZId<fNSegZTPCInt) ++tpcIntZId;    // find Z segment
 426   int tpcIntRId = fSegZIdTPCInt[tpcIntZId];        // first R segment for this Z
 427   int tpcIntRIdMax = tpcIntRId + fNSegRTPCInt[tpcIntZId];
 428   while (rphiz[0]>fSegRTPCInt[tpcIntRId] && tpcIntRId<tpcIntRIdMax) ++tpcIntRId;    // find R segment
 429   GetParamTPCInt( tpcIntRId )->Eval(rphiz,b);
 430   //
 431 }
 432
 433
 434 #ifdef  _INC_CREATION_ALICHEB3D_
 435 //_______________________________________________
 436 void AliMagWrapCheb::LoadData(const char* inpfile)
 437 {
 438   // read coefficients data from the text file
 439   //
 440   TString strf = inpfile;
 441   gSystem->ExpandPathName(strf);
 442   FILE* stream = fopen(strf,"r");
 443   if (!stream) {
 444     printf("Did not find input file %s\n",strf.Data());
 445     return;
 446   }
 447   //
 448   TString buffs;
 449   AliCheb3DCalc::ReadLine(buffs,stream);
 450   if (!buffs.BeginsWith("START")) {
 451     Error("LoadData","Expected: \"START <name>\", found \"%s\"\nStop\n",buffs.Data());
 452     exit(1);
 453   }
 454   if (buffs.First(' ')>0) SetName(buffs.Data()+buffs.First(' ')+1);
 455   //
 456   // Solenoid part    -----------------------------------------------------------
 457   AliCheb3DCalc::ReadLine(buffs,stream);
 458   if (!buffs.BeginsWith("START SOLENOID")) {
 459     Error("LoadData","Expected: \"START SOLENOID\", found \"%s\"\nStop\n",buffs.Data());
 460     exit(1);
 461   }
 462   AliCheb3DCalc::ReadLine(buffs,stream); // nparam
 463   int nparSol = buffs.Atoi();
 464   //
 465   for (int ip=0;ip<nparSol;ip++) {
 466     AliCheb3D* cheb = new AliCheb3D();
 467     cheb->LoadData(stream);
 468     AddParamSol(cheb);
 469   }
 470   //
 471   AliCheb3DCalc::ReadLine(buffs,stream);
 472   if (!buffs.BeginsWith("END SOLENOID")) {
 473     Error("LoadData","Expected \"END SOLENOID\", found \"%s\"\nStop\n",buffs.Data());
 474     exit(1);
 475   }
 476   //
 477   // TPCInt part     -----------------------------------------------------------
 478   AliCheb3DCalc::ReadLine(buffs,stream);
 479   if (!buffs.BeginsWith("START TPCINT")) {
 480     Error("LoadData","Expected: \"START TPCINT\", found \"%s\"\nStop\n",buffs.Data());
 481     exit(1);
 482   }
 483   AliCheb3DCalc::ReadLine(buffs,stream); // nparam
 484   int nparTPCInt = buffs.Atoi();
 485   //
 486   for (int ip=0;ip<nparTPCInt;ip++) {
 487     AliCheb3D* cheb = new AliCheb3D();
 488     cheb->LoadData(stream);
 489     AddParamTPCInt(cheb);
 490   }
 491   //
 492   AliCheb3DCalc::ReadLine(buffs,stream);
 493   if (!buffs.BeginsWith("END TPCINT")) {
 494     Error("LoadData","Expected \"END TPCINT\", found \"%s\"\nStop\n",buffs.Data());
 495     exit(1);
 496   }
 497   //
 498   // Dipole part    -----------------------------------------------------------
 499   AliCheb3DCalc::ReadLine(buffs,stream);
 500   if (!buffs.BeginsWith("START DIPOLE")) {
 501     Error("LoadData","Expected: \"START DIPOLE\", found \"%s\"\nStop\n",buffs.Data());
 502     exit(1);
 503   }
 504   AliCheb3DCalc::ReadLine(buffs,stream); // nparam
 505   int nparDip = buffs.Atoi();
 506   //
 507   for (int ip=0;ip<nparDip;ip++) {
 508     AliCheb3D* cheb = new AliCheb3D();
 509     cheb->LoadData(stream);
 510     AddParamDip(cheb);
 511   }
 512   //
 513   AliCheb3DCalc::ReadLine(buffs,stream);
 514   if (!buffs.BeginsWith("END DIPOLE")) {
 515     Error("LoadData","Expected \"END DIPOLE\", found \"%s\"\nStop\n",GetName(),buffs.Data());
 516     exit(1);
 517   }
 518   //
 519   AliCheb3DCalc::ReadLine(buffs,stream);
 520   if (!buffs.BeginsWith("END") || !buffs.Contains(GetName())) {
 521     Error("LoadData","Expected: \"END %s\", found \"%s\"\nStop\n",GetName(),buffs.Data());
 522     exit(1);
 523   }
 524   //
 525   // ---------------------------------------------------------------------------
 526   fclose(stream);
 527   BuildTableSol();
 528   BuildTableTPCInt();
 529   BuildTableDip();
 530   printf("Loaded magnetic field \"%s\" from %s\n",GetName(),strf.Data());
 531   //
 532 }
 533 #endif
 534
 535 //_______________________________________________
 536 #ifdef  _INC_CREATION_ALICHEB3D_
 537
 538 //__________________________________________________________________________________________
 539 AliMagWrapCheb::AliMagWrapCheb(const char* inputFile) :
 540   fNParamsSol(0),
 541   fNSegZSol(0),
 542   fNParamsTPCInt(0),
 543   fNSegZTPCInt(0),
 544   fNParamsDip(0),
 545 //
 546   fNZSegDip(0),
 547   fNYSegDip(0),
 548   fNXSegDip(0),
 549 //
 550   fSegZSol(0),
 551   fSegRSol(0),
 552 //
 553   fSegZTPCInt(0),
 554   fSegRTPCInt(0),
 555 //
 556   fSegZDip(0),
 557   fSegYDip(0),
 558   fSegXDip(0),
 559 //
 560   fNSegRSol(0),
 561   fSegZIdSol(0),
 562 //
 563   fNSegRTPCInt(0),
 564   fSegZIdTPCInt(0),
 565 //
 566   fBegSegYDip(0),
 567   fNSegYDip(0),
 568   fBegSegXDip(0),
 569   fNSegXDip(0),
 570   fSegIDDip(0),
 571 //
 572   fMinZSol(0),
 573   fMaxZSol(0),
 574   fMaxRSol(0),
 575 //
 576   fMinZDip(0),
 577   fMaxZDip(0),
 578 //
 579   fMinZTPCInt(0),
 580   fMaxZTPCInt(0),
 581   fMaxRTPCInt(0),
 582 //
 583   fParamsSol(0),
 584   fParamsDip(0),
 585   fParamsTPCInt(0)
 586 //
 587 //
 588 {
 589   // construct from coeffs from the text file
 590   LoadData(inputFile);
 591 }
 592
 593 //__________________________________________________________________________________________
 594 void AliMagWrapCheb::AddParamSol(const AliCheb3D* param)
 595 {
 596   // adds new parameterization piece for Sol
 597   // NOTE: pieces must be added strictly in increasing R then increasing Z order
 598   //
 599   if (!fParamsSol) fParamsSol = new TObjArray();
 600   fParamsSol->Add( (AliCheb3D*)param );
 601   fNParamsSol++;
 602   //
 603 }
 604
 605 //__________________________________________________________________________________________
 606 void AliMagWrapCheb::AddParamTPCInt(const AliCheb3D* param)
 607 {
 608   // adds new parameterization piece for TPCInt
 609   // NOTE: pieces must be added strictly in increasing R then increasing Z order
 610   //
 611   if (!fParamsTPCInt) fParamsTPCInt = new TObjArray();
 612   fParamsTPCInt->Add( (AliCheb3D*)param);
 613   fNParamsTPCInt++;
 614   //
 615 }
 616
 617 //__________________________________________________________________________________________
 618 void AliMagWrapCheb::AddParamDip(const AliCheb3D* param)
 619 {
 620   // adds new parameterization piece for Dipole
 621   //
 622   if (!fParamsDip) fParamsDip = new TObjArray();
 623   fParamsDip->Add( (AliCheb3D*)param);
 624   fNParamsDip++;
 625   //
 626 }
 627
 628 //__________________________________________________________________________________________
 629 void AliMagWrapCheb::ResetTPCInt()
 630 {
 631   // clean TPC field integral (used for update)
 632   if (!fNParamsTPCInt) return;
 633   delete fParamsTPCInt;
 634   delete[] fSegZTPCInt;
 635   delete[] fSegRTPCInt;
 636   delete[] fNSegRTPCInt;
 637   delete[] fSegZIdTPCInt;
 638   //
 639   fNParamsTPCInt = 0;
 640   fNSegZTPCInt   = 0;
 641   fSegZTPCInt    = 0;
 642   fSegRTPCInt    = 0;
 643   fNSegRTPCInt   = 0;
 644   fSegZIdTPCInt  = 0;
 645   fMinZTPCInt    = 0;
 646   fMaxZTPCInt    = 0;
 647   fMaxRTPCInt    = 0;
 648   fParamsTPCInt  = 0;
 649   //
 650 }
 651
 652 //__________________________________________________
 653 void AliMagWrapCheb::BuildTableDip()
 654 {
 655   // build lookup table for dipole
 656   //
 657   TArrayF segY,segX;
 658   TArrayI begSegYDip,begSegXDip;
 659   TArrayI nsegYDip,nsegXDip;
 660   TArrayI segID;
 661   float *tmpSegZ,*tmpSegY,*tmpSegX;
 662   //
 663   // create segmentation in Z
 664   fNZSegDip = SegmentDipDimension(&tmpSegZ, fParamsDip, fNParamsDip, 2, 1,-1, 1,-1, 1,-1) - 1;
 665   fNYSegDip = 0;
 666   fNXSegDip = 0;
 667   //
 668   // for each Z slice create segmentation in Y
 669   begSegYDip.Set(fNZSegDip);
 670   nsegYDip.Set(fNZSegDip);
 671   float xyz[3];
 672   for (int iz=0;iz<fNZSegDip;iz++) {
 673     printf("\nZSegment#%d  %+e : %+e\n",iz,tmpSegZ[iz],tmpSegZ[iz+1]);
 674     int ny = SegmentDipDimension(&tmpSegY, fParamsDip, fNParamsDip, 1,
 675                                  1,-1, 1,-1, tmpSegZ[iz],tmpSegZ[iz+1]) - 1;
 676     segY.Set(ny + fNYSegDip);
 677     for (int iy=0;iy<ny;iy++) segY[fNYSegDip+iy] = tmpSegY[iy];
 678     begSegYDip[iz] = fNYSegDip;
 679     nsegYDip[iz] = ny;
 680     printf(" Found %d YSegments, to start from %d\n",ny, begSegYDip[iz]);
 681     //
 682     // for each slice in Z and Y create segmentation in X
 683     begSegXDip.Set(fNYSegDip+ny);
 684     nsegXDip.Set(fNYSegDip+ny);
 685     xyz[2] = (tmpSegZ[iz]+tmpSegZ[iz+1])/2.; // mean Z of this segment
 686     //
 687     for (int iy=0;iy<ny;iy++) {
 688       int isg = fNYSegDip+iy;
 689       printf("\n   YSegment#%d  %+e : %+e\n",iy, tmpSegY[iy],tmpSegY[iy+1]);
 690       int nx = SegmentDipDimension(&tmpSegX, fParamsDip, fNParamsDip, 0,
 691                                    1,-1, tmpSegY[iy],tmpSegY[iy+1], tmpSegZ[iz],tmpSegZ[iz+1]) - 1;
 692       //
 693       segX.Set(nx + fNXSegDip);
 694       for (int ix=0;ix<nx;ix++) segX[fNXSegDip+ix] = tmpSegX[ix];
 695       begSegXDip[isg] = fNXSegDip;
 696       nsegXDip[isg] = nx;
 697       printf("   Found %d XSegments, to start from %d\n",nx, begSegXDip[isg]);
 698       //
 699       segID.Set(fNXSegDip+nx);
 700       //
 701       // find corresponding params
 702       xyz[1] = (tmpSegY[iy]+tmpSegY[iy+1])/2.; // mean Y of this segment
 703       //
 704       for (int ix=0;ix<nx;ix++) {
 705         xyz[0] = (tmpSegX[ix]+tmpSegX[ix+1])/2.; // mean X of this segment
 706         for (int ipar=0;ipar<fNParamsDip;ipar++) {
 707           AliCheb3D* cheb = (AliCheb3D*) fParamsDip->At(ipar);
 708           if (!cheb->IsInside(xyz)) continue;
 709           segID[fNXSegDip+ix] = ipar;
 710           break;
 711         }
 712       }
 713       fNXSegDip += nx;
 714       //
 715       delete[] tmpSegX;
 716     }
 717     delete[] tmpSegY;
 718     fNYSegDip += ny;
 719   }
 720   //
 721   fMinZDip = tmpSegZ[0];
 722   fMaxZDip = tmpSegZ[fNZSegDip];
 723   fSegZDip    = new Float_t[fNZSegDip];
 724   for (int i=fNZSegDip;i--;) fSegZDip[i] = tmpSegZ[i];
 725   delete[] tmpSegZ;
 726   //
 727   fSegYDip    = new Float_t[fNYSegDip];
 728   fSegXDip    = new Float_t[fNXSegDip];
 729   fBegSegYDip = new Int_t[fNZSegDip];
 730   fNSegYDip   = new Int_t[fNZSegDip];
 731   fBegSegXDip = new Int_t[fNYSegDip];
 732   fNSegXDip   = new Int_t[fNYSegDip];
 733   fSegIDDip   = new Int_t[fNXSegDip];
 734   //
 735   for (int i=fNYSegDip;i--;) fSegYDip[i] = segY[i];
 736   for (int i=fNXSegDip;i--;) fSegXDip[i] = segX[i];
 737   for (int i=fNZSegDip;i--;) {fBegSegYDip[i] = begSegYDip[i]; fNSegYDip[i] = nsegYDip[i];}
 738   for (int i=fNYSegDip;i--;) {fBegSegXDip[i] = begSegXDip[i]; fNSegXDip[i] = nsegXDip[i];}
 739   for (int i=fNXSegDip;i--;) {fSegIDDip[i]   = segID[i];}
 740   //
 741 }
 742
 743 //__________________________________________________________________________________________
 744 void AliMagWrapCheb::BuildTableSol()
 745 {
 746   // build the indexes for each parameterization of Solenoid
 747   //
 748   const float kSafety=0.001;
 749   //
 750   if (fNParamsSol<1) return;
 751   fNSegZSol = 0;
 752   fMaxRSol = 0;
 753   fSegRSol   = new Float_t[fNParamsSol];
 754   float *tmpbufF  = new float[fNParamsSol+1];
 755   int   *tmpbufI  = new int[fNParamsSol+1];
 756   int   *tmpbufI1 = new int[fNParamsSol+1];
 757   //
 758   // count number of Z slices and number of R slices in each Z slice
 759   for (int ip=0;ip<fNParamsSol;ip++) {
 760     if (ip==0 || (GetParamSol(ip)->GetBoundMax(2)-GetParamSol(ip-1)->GetBoundMax(2))>kSafety) { // new Z slice
 761       tmpbufF[fNSegZSol] = GetParamSol(ip)->GetBoundMax(2);                                     //
 762       tmpbufI[fNSegZSol] = 0;
 763       tmpbufI1[fNSegZSol++] = ip;
 764     }
 765     fSegRSol[ip] = GetParamSol(ip)->GetBoundMax(0);  // upper R
 766     tmpbufI[fNSegZSol-1]++;
 767     if (fMaxRSol<fSegRSol[ip]) fMaxRSol = fSegRSol[ip];
 768   }
 769   //
 770   fSegZSol   = new Float_t[fNSegZSol];
 771   fSegZIdSol = new Int_t[fNSegZSol];
 772   fNSegRSol  = new Int_t[fNSegZSol];
 773   for (int iz=0;iz<fNSegZSol;iz++) {
 774     fSegZSol[iz]   = tmpbufF[iz];
 775     fNSegRSol[iz]  = tmpbufI[iz];
 776     fSegZIdSol[iz] = tmpbufI1[iz];
 777   }
 778   //
 779   fMinZSol = GetParamSol(0)->GetBoundMin(2);
 780   fMaxZSol = GetParamSol(fNParamsSol-1)->GetBoundMax(2);
 781   //
 782   delete[] tmpbufF;
 783   delete[] tmpbufI;
 784   delete[] tmpbufI1;
 785   //
 786   //
 787 }
 788
 789 //__________________________________________________________________________________________
 790 void AliMagWrapCheb::BuildTableTPCInt()
 791 {
 792   // build the indexes for each parameterization of TPC field integral
 793   //
 794   const float kSafety=0.001;
 795   //
 796   if (fNParamsTPCInt<1) return;
 797   fNSegZTPCInt = 0;
 798   fSegRTPCInt     = new Float_t[fNParamsTPCInt];
 799   float *tmpbufF  = new float[fNParamsTPCInt+1];
 800   int   *tmpbufI  = new int[fNParamsTPCInt+1];
 801   int   *tmpbufI1 = new int[fNParamsTPCInt+1];
 802   //
 803   // count number of Z slices and number of R slices in each Z slice
 804   for (int ip=0;ip<fNParamsTPCInt;ip++) {
 805     if (ip==0 || (GetParamTPCInt(ip)->GetBoundMax(2)-GetParamTPCInt(ip-1)->GetBoundMax(2))>kSafety) { // new Z slice
 806       tmpbufF[fNSegZTPCInt] = GetParamTPCInt(ip)->GetBoundMax(2);                                     //
 807       tmpbufI[fNSegZTPCInt] = 0;
 808       tmpbufI1[fNSegZTPCInt++] = ip;
 809     }
 810     fSegRTPCInt[ip] = GetParamTPCInt(ip)->GetBoundMax(0);  // upper R
 811     tmpbufI[fNSegZTPCInt-1]++;
 812   }
 813   //
 814   fSegZTPCInt   = new Float_t[fNSegZTPCInt];
 815   fSegZIdTPCInt = new Int_t[fNSegZTPCInt];
 816   fNSegRTPCInt  = new Int_t[fNSegZTPCInt];
 817   for (int iz=0;iz<fNSegZTPCInt;iz++) {
 818     fSegZTPCInt[iz]   = tmpbufF[iz];
 819     fNSegRTPCInt[iz]  = tmpbufI[iz];
 820     fSegZIdTPCInt[iz] = tmpbufI1[iz];
 821   }
 822   //
 823   fMinZTPCInt = GetParamTPCInt(0)->GetBoundMin(2);
 824   fMaxZTPCInt = GetParamTPCInt(fNParamsTPCInt-1)->GetBoundMax(2);
 825   fMaxRTPCInt = GetParamTPCInt(fNParamsTPCInt-1)->GetBoundMax(0);
 826   //
 827   delete[] tmpbufF;
 828   delete[] tmpbufI;
 829   delete[] tmpbufI1;
 830   //
 831   //
 832 }
 833
 834 void AliMagWrapCheb::SaveData(const char* outfile) const
 835 {
 836   // writes coefficients data to output text file
 837   TString strf = outfile;
 838   gSystem->ExpandPathName(strf);
 839   FILE* stream = fopen(strf,"w+");
 840   //
 841   // Sol part    ---------------------------------------------------------
 842   fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName());
 843   fprintf(stream,"START SOLENOID\n#Number of pieces\n%d\n",fNParamsSol);
 844   for (int ip=0;ip<fNParamsSol;ip++) GetParamSol(ip)->SaveData(stream);
 845   fprintf(stream,"#\nEND SOLENOID\n");
 846   //
 847   // TPCInt part ---------------------------------------------------------
 848   fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName());
 849   fprintf(stream,"START TPCINT\n#Number of pieces\n%d\n",fNParamsTPCInt);
 850   for (int ip=0;ip<fNParamsTPCInt;ip++) GetParamTPCInt(ip)->SaveData(stream);
 851   fprintf(stream,"#\nEND TPCINT\n");
 852   //
 853   // Dip part   ---------------------------------------------------------
 854   fprintf(stream,"START DIPOLE\n#Number of pieces\n%d\n",fNParamsDip);
 855   for (int ip=0;ip<fNParamsDip;ip++) GetParamDip(ip)->SaveData(stream);
 856   fprintf(stream,"#\nEND DIPOLE\n");
 857   //
 858   fprintf(stream,"#\nEND %s\n",GetName());
 859   //
 860   fclose(stream);
 861   //
 862 }
 863
 864 Int_t AliMagWrapCheb::SegmentDipDimension(float** seg,const TObjArray* par,int npar, int dim,
 865                                        float xmn,float xmx,float ymn,float ymx,float zmn,float zmx)
 866 {
 867   // find all boundaries in deimension dim for boxes in given region.
 868   // if mn>mx for given projection the check is not done for it.
 869   float *tmpC = new float[2*npar];
 870   int *tmpInd = new int[2*npar];
 871   int nseg0 = 0;
 872   for (int ip=0;ip<npar;ip++) {
 873     AliCheb3D* cheb = (AliCheb3D*) par->At(ip);
 874     if (xmn<xmx && (cheb->GetBoundMin(0)>(xmx+xmn)/2 || cheb->GetBoundMax(0)<(xmn+xmx)/2)) continue;
 875     if (ymn<ymx && (cheb->GetBoundMin(1)>(ymx+ymn)/2 || cheb->GetBoundMax(1)<(ymn+ymx)/2)) continue;
 876     if (zmn<zmx && (cheb->GetBoundMin(2)>(zmx+zmn)/2 || cheb->GetBoundMax(2)<(zmn+zmx)/2)) continue;
 877     //
 878     tmpC[nseg0++] = cheb->GetBoundMin(dim);
 879     tmpC[nseg0++] = cheb->GetBoundMax(dim);
 880   }
 881   // range Dim's boundaries in increasing order
 882   TMath::Sort(nseg0,tmpC,tmpInd,kFALSE);
 883   // count number of really different Z's
 884   int nseg = 0;
 885   float cprev = -1e6;
 886   for (int ip=0;ip<nseg0;ip++) {
 887     if (TMath::Abs(cprev-tmpC[ tmpInd[ip] ])>1e-4) {
 888       cprev = tmpC[ tmpInd[ip] ];
 889       nseg++;
 890     }
 891     else tmpInd[ip] = -1; // supress redundant Z
 892   }
 893   //
 894   *seg  = new float[nseg]; // create final Z segmenations
 895   nseg = 0;
 896   for (int ip=0;ip<nseg0;ip++) if (tmpInd[ip]>=0) (*seg)[nseg++] = tmpC[ tmpInd[ip] ];
 897   //
 898   delete[] tmpC;
 899   delete[] tmpInd;
 900   return nseg;
 901 }
 902
 903 #endif
 904