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