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