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   return FieldCylSolBz(rphiz);
 300 }
 301
 302
 303 //__________________________________________________________________________________________
 304 void AliMagWrapCheb::Print(Option_t *) const
 305 {
 306   // print info
 307   printf("Alice magnetic field parameterized by Chebyshev polynomials\n");
 308   printf("Segmentation for Solenoid (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZSol,fMaxZSol,fMaxRSol);
 309   //
 310   if (fParamsSol) fParamsSol->Print();
 311   /*
 312   for (int iz=0;iz<fNSegZSol;iz++) {
 313     AliCheb3D* param = GetParamSol( fSegZIdSol[iz] );
 314     printf("*** Z Segment %2d (%+7.2f<Z<%+7.2f)\t***\n",iz,param->GetBoundMin(2),param->GetBoundMax(2));
 315     for (int ir=0;ir<fNSegRSol[iz];ir++) {
 316       param = GetParamSol( fSegZIdSol[iz]+ir );
 317       printf("    R Segment %2d (%+7.2f<R<%+7.2f, Precision: %.1e) (ID=%2d)\n",ir, param->GetBoundMin(0),
 318              param->GetBoundMax(0),param->GetPrecision(),fSegZIdSol[iz]+ir);
 319     }
 320   }
 321   */
 322   //
 323   printf("Segmentation for TPC field integral (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZTPCInt,fMaxZTPCInt,fMaxRTPCInt);
 324   //
 325   if (fParamsTPCInt) fParamsTPCInt->Print();
 326   /*
 327   for (int iz=0;iz<fNSegZTPCInt;iz++) {
 328     AliCheb3D* param = GetParamTPCInt( fSegZIdTPCInt[iz] );
 329     printf("*** Z Segment %2d (%+7.2f<Z<%+7.2f)\t***\n",iz,param->GetBoundMin(2),param->GetBoundMax(2));
 330     for (int ir=0;ir<fNSegRTPCInt[iz];ir++) {
 331       param = GetParamTPCInt( fSegZIdTPCInt[iz]+ir );
 332       printf("    R Segment %2d (%+7.2f<R<%+7.2f, Precision: %.1e) (ID=%2d)\n",ir, param->GetBoundMin(0),
 333              param->GetBoundMax(0),param->GetPrecision(),fSegZIdTPCInt[iz]+ir);
 334     }
 335   }
 336   */
 337   //
 338   printf("Segmentation for Dipole (%+.2f<Z<%+.2f cm)\n",fMinZDip,fMaxZDip);
 339   if (fParamsDip) fParamsDip->Print();
 340   //
 341
 342
 343 }
 344
 345
 346 //__________________________________________________________________________________________________
 347 Int_t    AliMagWrapCheb::FindDipSegment(const Double_t *xyz) const
 348 {
 349   // find the segment containing point xyz. If it is outside find the closest segment
 350   int xid,yid,zid = TMath::BinarySearch(fNZSegDip,fSegZDip,(Float_t)xyz[2]); // find zsegment
 351   int ysegBeg = fBegSegYDip[zid];
 352   //
 353   for (yid=0;yid<fNSegYDip[zid];yid++) if (xyz[1]<fSegYDip[ysegBeg+yid]) break;
 354   if ( --yid < 0 ) yid = 0;
 355   yid +=  ysegBeg;
 356   //
 357   int xsegBeg = fBegSegXDip[yid];
 358   for (xid=0;xid<fNSegXDip[yid];xid++) if (xyz[0]<fSegXDip[xsegBeg+xid]) break;
 359   if ( --xid < 0) xid = 0;
 360   xid +=  xsegBeg;
 361   //
 362   return fSegIDDip[xid];
 363 }
 364
 365 //__________________________________________________________________________________________
 366 void AliMagWrapCheb::GetTPCInt(const Double_t *xyz, Double_t *b) const
 367 {
 368   // compute TPC region field integral in cartesian coordinates.
 369   // If point is outside of the parameterized region get it at closeset valid point
 370   static Double_t rphiz[3];
 371   //
 372   // TPCInt region
 373   // convert coordinates to cyl system
 374   CartToCyl(xyz,rphiz);
 375 #ifndef _BRING_TO_BOUNDARY_
 376   if ( (rphiz[2]>GetMaxZTPCInt()||rphiz[2]<GetMinZTPCInt()) ||
 377        rphiz[0]>GetMaxRTPCInt()) {for (int i=3;i--;) b[i]=0; return;}
 378 #endif
 379   //
 380   GetTPCIntCyl(rphiz,b);
 381   //
 382   // convert field to cartesian system
 383   CylToCartCylB(rphiz, b,b);
 384   //
 385 }
 386
 387 //__________________________________________________________________________________________
 388 void AliMagWrapCheb::FieldCylSol(const Double_t *rphiz, Double_t *b) const
 389 {
 390   // compute Solenoid field in Cylindircal coordinates
 391   // note: if the point is outside the volume get the field in closest parameterized point
 392   int SolZId = 0;
 393   while (rphiz[2]>fSegZSol[SolZId] && SolZId<fNSegZSol-1) ++SolZId;    // find Z segment
 394   int SolRId = fSegZIdSol[SolZId];        // first R segment for this Z
 395   int SolRMax = SolRId + fNSegRSol[SolZId];
 396   while (rphiz[0]>fSegRSol[SolRId] && SolRId<SolRMax-1) ++SolRId;    // find R segment
 397   GetParamSol( SolRId )->Eval(rphiz,b);
 398   //
 399 }
 400
 401 //__________________________________________________________________________________________
 402 Double_t AliMagWrapCheb::FieldCylSolBz(const Double_t *rphiz) const
 403 {
 404   // compute Solenoid field in Cylindircal coordinates
 405   // note: if the point is outside the volume get the field in closest parameterized point
 406   int SolZId = 0;
 407   while (rphiz[2]>fSegZSol[SolZId] && SolZId<fNSegZSol-1) ++SolZId;    // find Z segment
 408   int SolRId = fSegZIdSol[SolZId];        // first R segment for this Z
 409   int SolRMax = SolRId + fNSegRSol[SolZId];
 410   while (rphiz[0]>fSegRSol[SolRId] && SolRId<SolRMax-1) ++SolRId;    // find R segment
 411   return GetParamSol( SolRId )->Eval(rphiz,2);
 412   //
 413 }
 414
 415 //__________________________________________________________________________________________
 416 void AliMagWrapCheb::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
 417 {
 418   // compute field integral in TPC region in Cylindircal coordinates
 419   // note: the check for the point being inside the parameterized region is done outside
 420   int tpcIntZId = 0;
 421   while (rphiz[2]>fSegZTPCInt[tpcIntZId] && tpcIntZId<fNSegZTPCInt) ++tpcIntZId;    // find Z segment
 422   int tpcIntRId = fSegZIdTPCInt[tpcIntZId];        // first R segment for this Z
 423   int tpcIntRIdMax = tpcIntRId + fNSegRTPCInt[tpcIntZId];
 424   while (rphiz[0]>fSegRTPCInt[tpcIntRId] && tpcIntRId<tpcIntRIdMax) ++tpcIntRId;    // find R segment
 425   GetParamTPCInt( tpcIntRId )->Eval(rphiz,b);
 426   //
 427 }
 428
 429
 430 #ifdef  _INC_CREATION_ALICHEB3D_
 431 //_______________________________________________
 432 void AliMagWrapCheb::LoadData(const char* inpfile)
 433 {
 434   // read coefficients data from the text file
 435   //
 436   TString strf = inpfile;
 437   gSystem->ExpandPathName(strf);
 438   FILE* stream = fopen(strf,"r");
 439   if (!stream) {
 440     printf("Did not find input file %s\n",strf.Data());
 441     return;
 442   }
 443   //
 444   TString buffs;
 445   AliCheb3DCalc::ReadLine(buffs,stream);
 446   if (!buffs.BeginsWith("START")) {
 447     Error("LoadData","Expected: \"START <name>\", found \"%s\"\nStop\n",buffs.Data());
 448     exit(1);
 449   }
 450   if (buffs.First(' ')>0) SetName(buffs.Data()+buffs.First(' ')+1);
 451   //
 452   // Solenoid part    -----------------------------------------------------------
 453   AliCheb3DCalc::ReadLine(buffs,stream);
 454   if (!buffs.BeginsWith("START SOLENOID")) {
 455     Error("LoadData","Expected: \"START SOLENOID\", found \"%s\"\nStop\n",buffs.Data());
 456     exit(1);
 457   }
 458   AliCheb3DCalc::ReadLine(buffs,stream); // nparam
 459   int nparSol = buffs.Atoi();
 460   //
 461   for (int ip=0;ip<nparSol;ip++) {
 462     AliCheb3D* cheb = new AliCheb3D();
 463     cheb->LoadData(stream);
 464     AddParamSol(cheb);
 465   }
 466   //
 467   AliCheb3DCalc::ReadLine(buffs,stream);
 468   if (!buffs.BeginsWith("END SOLENOID")) {
 469     Error("LoadData","Expected \"END SOLENOID\", found \"%s\"\nStop\n",buffs.Data());
 470     exit(1);
 471   }
 472   //
 473   // TPCInt part     -----------------------------------------------------------
 474   AliCheb3DCalc::ReadLine(buffs,stream);
 475   if (!buffs.BeginsWith("START TPCINT")) {
 476     Error("LoadData","Expected: \"START TPCINT\", found \"%s\"\nStop\n",buffs.Data());
 477     exit(1);
 478   }
 479   AliCheb3DCalc::ReadLine(buffs,stream); // nparam
 480   int nparTPCInt = buffs.Atoi();
 481   //
 482   for (int ip=0;ip<nparTPCInt;ip++) {
 483     AliCheb3D* cheb = new AliCheb3D();
 484     cheb->LoadData(stream);
 485     AddParamTPCInt(cheb);
 486   }
 487   //
 488   AliCheb3DCalc::ReadLine(buffs,stream);
 489   if (!buffs.BeginsWith("END TPCINT")) {
 490     Error("LoadData","Expected \"END TPCINT\", found \"%s\"\nStop\n",buffs.Data());
 491     exit(1);
 492   }
 493   //
 494   // Dipole part    -----------------------------------------------------------
 495   AliCheb3DCalc::ReadLine(buffs,stream);
 496   if (!buffs.BeginsWith("START DIPOLE")) {
 497     Error("LoadData","Expected: \"START DIPOLE\", found \"%s\"\nStop\n",buffs.Data());
 498     exit(1);
 499   }
 500   AliCheb3DCalc::ReadLine(buffs,stream); // nparam
 501   int nparDip = buffs.Atoi();
 502   //
 503   for (int ip=0;ip<nparDip;ip++) {
 504     AliCheb3D* cheb = new AliCheb3D();
 505     cheb->LoadData(stream);
 506     AddParamDip(cheb);
 507   }
 508   //
 509   AliCheb3DCalc::ReadLine(buffs,stream);
 510   if (!buffs.BeginsWith("END DIPOLE")) {
 511     Error("LoadData","Expected \"END DIPOLE\", found \"%s\"\nStop\n",GetName(),buffs.Data());
 512     exit(1);
 513   }
 514   //
 515   AliCheb3DCalc::ReadLine(buffs,stream);
 516   if (!buffs.BeginsWith("END") || !buffs.Contains(GetName())) {
 517     Error("LoadData","Expected: \"END %s\", found \"%s\"\nStop\n",GetName(),buffs.Data());
 518     exit(1);
 519   }
 520   //
 521   // ---------------------------------------------------------------------------
 522   fclose(stream);
 523   BuildTableSol();
 524   BuildTableTPCInt();
 525   BuildTableDip();
 526   printf("Loaded magnetic field \"%s\" from %s\n",GetName(),strf.Data());
 527   //
 528 }
 529 #endif
 530
 531 //_______________________________________________
 532 #ifdef  _INC_CREATION_ALICHEB3D_
 533
 534 //__________________________________________________________________________________________
 535 AliMagWrapCheb::AliMagWrapCheb(const char* inputFile) :
 536   fNParamsSol(0),
 537   fNSegZSol(0),
 538   fNParamsTPCInt(0),
 539   fNSegZTPCInt(0),
 540   fNParamsDip(0),
 541 //
 542   fNZSegDip(0),
 543   fNYSegDip(0),
 544   fNXSegDip(0),
 545 //
 546   fSegZSol(0),
 547   fSegRSol(0),
 548 //
 549   fSegZTPCInt(0),
 550   fSegRTPCInt(0),
 551 //
 552   fSegZDip(0),
 553   fSegYDip(0),
 554   fSegXDip(0),
 555 //
 556   fNSegRSol(0),
 557   fSegZIdSol(0),
 558 //
 559   fNSegRTPCInt(0),
 560   fSegZIdTPCInt(0),
 561 //
 562   fBegSegYDip(0),
 563   fNSegYDip(0),
 564   fBegSegXDip(0),
 565   fNSegXDip(0),
 566   fSegIDDip(0),
 567 //
 568   fMinZSol(0),
 569   fMaxZSol(0),
 570   fMaxRSol(0),
 571 //
 572   fMinZDip(0),
 573   fMaxZDip(0),
 574 //
 575   fMinZTPCInt(0),
 576   fMaxZTPCInt(0),
 577   fMaxRTPCInt(0),
 578 //
 579   fParamsSol(0),
 580   fParamsDip(0),
 581   fParamsTPCInt(0)
 582 //
 583 //
 584 {
 585   // construct from coeffs from the text file
 586   LoadData(inputFile);
 587 }
 588
 589 //__________________________________________________________________________________________
 590 void AliMagWrapCheb::AddParamSol(const AliCheb3D* param)
 591 {
 592   // adds new parameterization piece for Sol
 593   // NOTE: pieces must be added strictly in increasing R then increasing Z order
 594   //
 595   if (!fParamsSol) fParamsSol = new TObjArray();
 596   fParamsSol->Add( (AliCheb3D*)param );
 597   fNParamsSol++;
 598   //
 599 }
 600
 601 //__________________________________________________________________________________________
 602 void AliMagWrapCheb::AddParamTPCInt(const AliCheb3D* param)
 603 {
 604   // adds new parameterization piece for TPCInt
 605   // NOTE: pieces must be added strictly in increasing R then increasing Z order
 606   //
 607   if (!fParamsTPCInt) fParamsTPCInt = new TObjArray();
 608   fParamsTPCInt->Add( (AliCheb3D*)param);
 609   fNParamsTPCInt++;
 610   //
 611 }
 612
 613 //__________________________________________________________________________________________
 614 void AliMagWrapCheb::AddParamDip(const AliCheb3D* param)
 615 {
 616   // adds new parameterization piece for Dipole
 617   //
 618   if (!fParamsDip) fParamsDip = new TObjArray();
 619   fParamsDip->Add( (AliCheb3D*)param);
 620   fNParamsDip++;
 621   //
 622 }
 623
 624 //__________________________________________________________________________________________
 625 void AliMagWrapCheb::ResetTPCInt()
 626 {
 627   // clean TPC field integral (used for update)
 628   if (!fNParamsTPCInt) return;
 629   delete fParamsTPCInt;
 630   delete[] fSegZTPCInt;
 631   delete[] fSegRTPCInt;
 632   delete[] fNSegRTPCInt;
 633   delete[] fSegZIdTPCInt;
 634   //
 635   fNParamsTPCInt = 0;
 636   fNSegZTPCInt   = 0;
 637   fSegZTPCInt    = 0;
 638   fSegRTPCInt    = 0;
 639   fNSegRTPCInt   = 0;
 640   fSegZIdTPCInt  = 0;
 641   fMinZTPCInt    = 0;
 642   fMaxZTPCInt    = 0;
 643   fMaxRTPCInt    = 0;
 644   fParamsTPCInt  = 0;
 645   //
 646 }
 647
 648 //__________________________________________________
 649 void AliMagWrapCheb::BuildTableDip()
 650 {
 651   // build lookup table for dipole
 652   //
 653   TArrayF segY,segX;
 654   TArrayI begSegYDip,begSegXDip;
 655   TArrayI nsegYDip,nsegXDip;
 656   TArrayI segID;
 657   float *tmpSegZ,*tmpSegY,*tmpSegX;
 658   //
 659   // create segmentation in Z
 660   fNZSegDip = SegmentDipDimension(&tmpSegZ, fParamsDip, fNParamsDip, 2, 1,-1, 1,-1, 1,-1) - 1;
 661   fNYSegDip = 0;
 662   fNXSegDip = 0;
 663   //
 664   // for each Z slice create segmentation in Y
 665   begSegYDip.Set(fNZSegDip);
 666   nsegYDip.Set(fNZSegDip);
 667   float xyz[3];
 668   for (int iz=0;iz<fNZSegDip;iz++) {
 669     printf("\nZSegment#%d  %+e : %+e\n",iz,tmpSegZ[iz],tmpSegZ[iz+1]);
 670     int ny = SegmentDipDimension(&tmpSegY, fParamsDip, fNParamsDip, 1,
 671                                  1,-1, 1,-1, tmpSegZ[iz],tmpSegZ[iz+1]) - 1;
 672     segY.Set(ny + fNYSegDip);
 673     for (int iy=0;iy<ny;iy++) segY[fNYSegDip+iy] = tmpSegY[iy];
 674     begSegYDip[iz] = fNYSegDip;
 675     nsegYDip[iz] = ny;
 676     printf(" Found %d YSegments, to start from %d\n",ny, begSegYDip[iz]);
 677     //
 678     // for each slice in Z and Y create segmentation in X
 679     begSegXDip.Set(fNYSegDip+ny);
 680     nsegXDip.Set(fNYSegDip+ny);
 681     xyz[2] = (tmpSegZ[iz]+tmpSegZ[iz+1])/2.; // mean Z of this segment
 682     //
 683     for (int iy=0;iy<ny;iy++) {
 684       int isg = fNYSegDip+iy;
 685       printf("\n   YSegment#%d  %+e : %+e\n",iy, tmpSegY[iy],tmpSegY[iy+1]);
 686       int nx = SegmentDipDimension(&tmpSegX, fParamsDip, fNParamsDip, 0,
 687                                    1,-1, tmpSegY[iy],tmpSegY[iy+1], tmpSegZ[iz],tmpSegZ[iz+1]) - 1;
 688       //
 689       segX.Set(nx + fNXSegDip);
 690       for (int ix=0;ix<nx;ix++) segX[fNXSegDip+ix] = tmpSegX[ix];
 691       begSegXDip[isg] = fNXSegDip;
 692       nsegXDip[isg] = nx;
 693       printf("   Found %d XSegments, to start from %d\n",nx, begSegXDip[isg]);
 694       //
 695       segID.Set(fNXSegDip+nx);
 696       //
 697       // find corresponding params
 698       xyz[1] = (tmpSegY[iy]+tmpSegY[iy+1])/2.; // mean Y of this segment
 699       //
 700       for (int ix=0;ix<nx;ix++) {
 701         xyz[0] = (tmpSegX[ix]+tmpSegX[ix+1])/2.; // mean X of this segment
 702         for (int ipar=0;ipar<fNParamsDip;ipar++) {
 703           AliCheb3D* cheb = (AliCheb3D*) fParamsDip->At(ipar);
 704           if (!cheb->IsInside(xyz)) continue;
 705           segID[fNXSegDip+ix] = ipar;
 706           break;
 707         }
 708       }
 709       fNXSegDip += nx;
 710       //
 711       delete[] tmpSegX;
 712     }
 713     delete[] tmpSegY;
 714     fNYSegDip += ny;
 715   }
 716   //
 717   fMinZDip = tmpSegZ[0];
 718   fMaxZDip = tmpSegZ[fNZSegDip];
 719   fSegZDip    = new Float_t[fNZSegDip];
 720   for (int i=fNZSegDip;i--;) fSegZDip[i] = tmpSegZ[i];
 721   delete[] tmpSegZ;
 722   //
 723   fSegYDip    = new Float_t[fNYSegDip];
 724   fSegXDip    = new Float_t[fNXSegDip];
 725   fBegSegYDip = new Int_t[fNZSegDip];
 726   fNSegYDip   = new Int_t[fNZSegDip];
 727   fBegSegXDip = new Int_t[fNYSegDip];
 728   fNSegXDip   = new Int_t[fNYSegDip];
 729   fSegIDDip   = new Int_t[fNXSegDip];
 730   //
 731   for (int i=fNYSegDip;i--;) fSegYDip[i] = segY[i];
 732   for (int i=fNXSegDip;i--;) fSegXDip[i] = segX[i];
 733   for (int i=fNZSegDip;i--;) {fBegSegYDip[i] = begSegYDip[i]; fNSegYDip[i] = nsegYDip[i];}
 734   for (int i=fNYSegDip;i--;) {fBegSegXDip[i] = begSegXDip[i]; fNSegXDip[i] = nsegXDip[i];}
 735   for (int i=fNXSegDip;i--;) {fSegIDDip[i]   = segID[i];}
 736   //
 737 }
 738
 739 //__________________________________________________________________________________________
 740 void AliMagWrapCheb::BuildTableSol()
 741 {
 742   // build the indexes for each parameterization of Solenoid
 743   //
 744   const float kSafety=0.001;
 745   //
 746   if (fNParamsSol<1) return;
 747   fNSegZSol = 0;
 748   fMaxRSol = 0;
 749   fSegRSol   = new Float_t[fNParamsSol];
 750   float *tmpbufF  = new float[fNParamsSol+1];
 751   int   *tmpbufI  = new int[fNParamsSol+1];
 752   int   *tmpbufI1 = new int[fNParamsSol+1];
 753   //
 754   // count number of Z slices and number of R slices in each Z slice
 755   for (int ip=0;ip<fNParamsSol;ip++) {
 756     if (ip==0 || (GetParamSol(ip)->GetBoundMax(2)-GetParamSol(ip-1)->GetBoundMax(2))>kSafety) { // new Z slice
 757       tmpbufF[fNSegZSol] = GetParamSol(ip)->GetBoundMax(2);                                     //
 758       tmpbufI[fNSegZSol] = 0;
 759       tmpbufI1[fNSegZSol++] = ip;
 760     }
 761     fSegRSol[ip] = GetParamSol(ip)->GetBoundMax(0);  // upper R
 762     tmpbufI[fNSegZSol-1]++;
 763     if (fMaxRSol<fSegRSol[ip]) fMaxRSol = fSegRSol[ip];
 764   }
 765   //
 766   fSegZSol   = new Float_t[fNSegZSol];
 767   fSegZIdSol = new Int_t[fNSegZSol];
 768   fNSegRSol  = new Int_t[fNSegZSol];
 769   for (int iz=0;iz<fNSegZSol;iz++) {
 770     fSegZSol[iz]   = tmpbufF[iz];
 771     fNSegRSol[iz]  = tmpbufI[iz];
 772     fSegZIdSol[iz] = tmpbufI1[iz];
 773   }
 774   //
 775   fMinZSol = GetParamSol(0)->GetBoundMin(2);
 776   fMaxZSol = GetParamSol(fNParamsSol-1)->GetBoundMax(2);
 777   //
 778   delete[] tmpbufF;
 779   delete[] tmpbufI;
 780   delete[] tmpbufI1;
 781   //
 782   //
 783 }
 784
 785 //__________________________________________________________________________________________
 786 void AliMagWrapCheb::BuildTableTPCInt()
 787 {
 788   // build the indexes for each parameterization of TPC field integral
 789   //
 790   const float kSafety=0.001;
 791   //
 792   if (fNParamsTPCInt<1) return;
 793   fNSegZTPCInt = 0;
 794   fSegRTPCInt     = new Float_t[fNParamsTPCInt];
 795   float *tmpbufF  = new float[fNParamsTPCInt+1];
 796   int   *tmpbufI  = new int[fNParamsTPCInt+1];
 797   int   *tmpbufI1 = new int[fNParamsTPCInt+1];
 798   //
 799   // count number of Z slices and number of R slices in each Z slice
 800   for (int ip=0;ip<fNParamsTPCInt;ip++) {
 801     if (ip==0 || (GetParamTPCInt(ip)->GetBoundMax(2)-GetParamTPCInt(ip-1)->GetBoundMax(2))>kSafety) { // new Z slice
 802       tmpbufF[fNSegZTPCInt] = GetParamTPCInt(ip)->GetBoundMax(2);                                     //
 803       tmpbufI[fNSegZTPCInt] = 0;
 804       tmpbufI1[fNSegZTPCInt++] = ip;
 805     }
 806     fSegRTPCInt[ip] = GetParamTPCInt(ip)->GetBoundMax(0);  // upper R
 807     tmpbufI[fNSegZTPCInt-1]++;
 808   }
 809   //
 810   fSegZTPCInt   = new Float_t[fNSegZTPCInt];
 811   fSegZIdTPCInt = new Int_t[fNSegZTPCInt];
 812   fNSegRTPCInt  = new Int_t[fNSegZTPCInt];
 813   for (int iz=0;iz<fNSegZTPCInt;iz++) {
 814     fSegZTPCInt[iz]   = tmpbufF[iz];
 815     fNSegRTPCInt[iz]  = tmpbufI[iz];
 816     fSegZIdTPCInt[iz] = tmpbufI1[iz];
 817   }
 818   //
 819   fMinZTPCInt = GetParamTPCInt(0)->GetBoundMin(2);
 820   fMaxZTPCInt = GetParamTPCInt(fNParamsTPCInt-1)->GetBoundMax(2);
 821   fMaxRTPCInt = GetParamTPCInt(fNParamsTPCInt-1)->GetBoundMax(0);
 822   //
 823   delete[] tmpbufF;
 824   delete[] tmpbufI;
 825   delete[] tmpbufI1;
 826   //
 827   //
 828 }
 829
 830 void AliMagWrapCheb::SaveData(const char* outfile) const
 831 {
 832   // writes coefficients data to output text file
 833   TString strf = outfile;
 834   gSystem->ExpandPathName(strf);
 835   FILE* stream = fopen(strf,"w+");
 836   //
 837   // Sol part    ---------------------------------------------------------
 838   fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName());
 839   fprintf(stream,"START SOLENOID\n#Number of pieces\n%d\n",fNParamsSol);
 840   for (int ip=0;ip<fNParamsSol;ip++) GetParamSol(ip)->SaveData(stream);
 841   fprintf(stream,"#\nEND SOLENOID\n");
 842   //
 843   // TPCInt part ---------------------------------------------------------
 844   fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName());
 845   fprintf(stream,"START TPCINT\n#Number of pieces\n%d\n",fNParamsTPCInt);
 846   for (int ip=0;ip<fNParamsTPCInt;ip++) GetParamTPCInt(ip)->SaveData(stream);
 847   fprintf(stream,"#\nEND TPCINT\n");
 848   //
 849   // Dip part   ---------------------------------------------------------
 850   fprintf(stream,"START DIPOLE\n#Number of pieces\n%d\n",fNParamsDip);
 851   for (int ip=0;ip<fNParamsDip;ip++) GetParamDip(ip)->SaveData(stream);
 852   fprintf(stream,"#\nEND DIPOLE\n");
 853   //
 854   fprintf(stream,"#\nEND %s\n",GetName());
 855   //
 856   fclose(stream);
 857   //
 858 }
 859
 860 Int_t AliMagWrapCheb::SegmentDipDimension(float** seg,const TObjArray* par,int npar, int dim,
 861                                        float xmn,float xmx,float ymn,float ymx,float zmn,float zmx)
 862 {
 863   // find all boundaries in deimension dim for boxes in given region.
 864   // if mn>mx for given projection the check is not done for it.
 865   float *tmpC = new float[2*npar];
 866   int *tmpInd = new int[2*npar];
 867   int nseg0 = 0;
 868   for (int ip=0;ip<npar;ip++) {
 869     AliCheb3D* cheb = (AliCheb3D*) par->At(ip);
 870     if (xmn<xmx && (cheb->GetBoundMin(0)>(xmx+xmn)/2 || cheb->GetBoundMax(0)<(xmn+xmx)/2)) continue;
 871     if (ymn<ymx && (cheb->GetBoundMin(1)>(ymx+ymn)/2 || cheb->GetBoundMax(1)<(ymn+ymx)/2)) continue;
 872     if (zmn<zmx && (cheb->GetBoundMin(2)>(zmx+zmn)/2 || cheb->GetBoundMax(2)<(zmn+zmx)/2)) continue;
 873     //
 874     tmpC[nseg0++] = cheb->GetBoundMin(dim);
 875     tmpC[nseg0++] = cheb->GetBoundMax(dim);
 876   }
 877   // range Dim's boundaries in increasing order
 878   TMath::Sort(nseg0,tmpC,tmpInd,kFALSE);
 879   // count number of really different Z's
 880   int nseg = 0;
 881   float cprev = -1e6;
 882   for (int ip=0;ip<nseg0;ip++) {
 883     if (TMath::Abs(cprev-tmpC[ tmpInd[ip] ])>1e-4) {
 884       cprev = tmpC[ tmpInd[ip] ];
 885       nseg++;
 886     }
 887     else tmpInd[ip] = -1; // supress redundant Z
 888   }
 889   //
 890   *seg  = new float[nseg]; // create final Z segmenations
 891   nseg = 0;
 892   for (int ip=0;ip<nseg0;ip++) if (tmpInd[ip]>=0) (*seg)[nseg++] = tmpC[ tmpInd[ip] ];
 893   //
 894   delete[] tmpC;
 895   delete[] tmpInd;
 896   return nseg;
 897 }
 898
 899 #endif
 900