Updates for full field map (R. Shahoyan)
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 7 Apr 2008 15:53:15 +0000 (15:53 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 7 Apr 2008 15:53:15 +0000 (15:53 +0000)
STEER/AliMagFCheb.cxx
STEER/AliMagFCheb.h

index fcc42b3..b486249 100644 (file)
@@ -1,20 +1,3 @@
-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- *                                                                        *
- * Author: The ALICE Off-line Project.                                    *
- * Contributors are mentioned in the code where appropriate.              *
- *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
- * documentation strictly for non-commercial purposes is hereby granted   *
- * without fee, provided that the above copyright notice appears in all   *
- * copies and that both the copyright notice and this permission notice   *
- * appear in the supporting documentation. The authors make no claims     *
- * about the suitability of this software for any purpose. It is          *
- * provided "as is" without express or implied warranty.                  *
- **************************************************************************/
-
-/* $Id$ */
-
 ///////////////////////////////////////////////////////////////////////////////////
 //                                                                               //
 //  Wrapper for the set of mag.field parameterizations by Chebyshev polinomials  //
 //  For cylindrical coordinates/components:                                      //
 //    FieldCyl(float* rphiz, float* brphiz)                                      //
 //                                                                               //
-//  For the moment only the solenoid part is parameterized in the volume defined //
-//  by R<500, -550<Z<550 cm                                                      //
+//  The solenoid part is parameterized in the volume  R<500, -550<Z<550 cm       //
 //                                                                               //
 //  The region R<423 cm,  -343.3<Z<481.3 for 30kA and -343.3<Z<481.3 for 12kA    //
 //  is parameterized using measured data while outside the Tosca calculation     //
 //  is used (matched to data on the boundary of the measurements)                //
 //                                                                               //
-//  If the querried point is outside the validity region no the return values    //
-//  for the field components are set to 0.                                       //
+//  Two options are possible:                                                    //
+//  1) _BRING_TO_BOUNDARY_ is defined in the AliCheb3D:                          //
+//     If the querried point is outside of the validity region then the field    //
+//     at the closest point on the fitted surface is returned.                   //
+//  2) _BRING_TO_BOUNDARY_ is not defined in the AliCheb3D:                      //
+//     If the querried point is outside of the validity region the return        //
+//     value for the field components are set to 0.                              //
+//                                                                               //
+//  To obtain the field integral in the TPC region from given point to nearest   //
+//  cathod plane (+- 250 cm) use:                                                //
+//  GetTPCInt(float* xyz, float* bxyz);  for Cartesian frame                     //
+//  or                                                                           //
+//  GetTPCIntCyl(Float_t *rphiz, Float_t *b); for Cylindrical frame              //
+//                                                                               //
+//                                                                               //
+//  The units are kiloGauss and cm.                                              //
 //                                                                               //
 ///////////////////////////////////////////////////////////////////////////////////
 
-#include "AliLog.h"
 #include "AliMagFCheb.h"
 
 ClassImp(AliMagFCheb)
 
-
-
 //__________________________________________________________________________________________
-AliMagFCheb::AliMagFCheb() :
-    TNamed(),
-    fNParamsSol(0),
-    fNSegZSol(0),
-    fNParamsDip(0),
-    fSegZSol(0),
-    fSegRSol(0),
-    fNSegRSol(0),
-    fSegZIdSol(0),
-    fMinZSol(0.),
-    fMaxZSol(0.),
-    fMaxRSol(0.),
-    fParamsSol(0),
-    fParamsDip(0)
-{
-  Init0();
-}
+AliMagFCheb::AliMagFCheb() : 
+  fNParamsSol(0),
+  fNSegZSol(0),
+  fNParamsTPCInt(0),
+  fNSegZTPCInt(0),
+  fNParamsDip(0),
+//
+  fNZSegDip(0),
+  fNYSegDip(0),
+  fNXSegDip(0),
+//
+  fSegZSol(0),
+  fSegRSol(0),
+//
+  fSegZTPCInt(0),
+  fSegRTPCInt(0),
+//
+  fSegZDip(0),
+  fSegYDip(0),
+  fSegXDip(0),
+//
+  fNSegRSol(0),
+  fSegZIdSol(0),
+//
+  fNSegRTPCInt(0),
+  fSegZIdTPCInt(0),
+//
+  fBegSegYDip(0),
+  fNSegYDip(0),
+  fBegSegXDip(0),
+  fNSegXDip(0),
+  fSegIDDip(0),
+//
+  fMinZSol(1e6),
+  fMaxZSol(-1e6),
+  fMaxRSol(-1e6), 
+//
+  fMinZDip(1e6),
+  fMaxZDip(-1e6),
+//
+  fMinZTPCInt(1e6),
+  fMaxZTPCInt(-1e6),
+  fMaxRTPCInt(-1e6), 
+//
+  fParamsSol(0),
+  fParamsDip(0),
+  fParamsTPCInt(0)
+//
+{}
 
-AliMagFCheb::AliMagFCheb(const char* inputFile) :
-    TNamed("Field Map", inputFile),
-    fNParamsSol(0),
-    fNSegZSol(0),
-    fNParamsDip(0),
-    fSegZSol(0),
-    fSegRSol(0),
-    fNSegRSol(0),
-    fSegZIdSol(0),
-    fMinZSol(0.),
-    fMaxZSol(0.),
-    fMaxRSol(0.),
-    fParamsSol(0),
-    fParamsDip(0)
+//__________________________________________________________________________________________
+AliMagFCheb::AliMagFCheb(const AliMagFCheb& src) : 
+  TNamed(src),
+  fNParamsSol(0),
+  fNSegZSol(0),
+  fNParamsTPCInt(0),
+  fNSegZTPCInt(0),
+  fNParamsDip(0),
+//
+  fNZSegDip(0),
+  fNYSegDip(0),
+  fNXSegDip(0),
+//
+  fSegZSol(0),
+  fSegRSol(0),
+//
+  fSegZTPCInt(0),
+  fSegRTPCInt(0),
+//
+  fSegZDip(0),
+  fSegYDip(0),
+  fSegXDip(0),
+//
+  fNSegRSol(0),
+  fSegZIdSol(0),
+//
+  fNSegRTPCInt(0),
+  fSegZIdTPCInt(0),
+//
+  fBegSegYDip(0),
+  fNSegYDip(0),
+  fBegSegXDip(0),
+  fNSegXDip(0),
+  fSegIDDip(0),
+//
+  fMinZSol(1e6),
+  fMaxZSol(-1e6),
+  fMaxRSol(-1e6), 
+//
+  fMinZDip(1e6),
+  fMaxZDip(-1e6),
+//
+  fMinZTPCInt(1e6),
+  fMaxZTPCInt(-1e6),
+  fMaxRTPCInt(-1e6), 
+//
+  fParamsSol(0),
+  fParamsDip(0),
+  fParamsTPCInt(0)
 {
-  Init0();
-  LoadData(inputFile);
+  CopyFrom(src);
+  //
 }
 
 //__________________________________________________________________________________________
-AliMagFCheb::AliMagFCheb(const AliMagFCheb& src)
-    : TNamed(src), 
-      fNParamsSol(src.fNParamsSol),
-      fNSegZSol(src.fNSegZSol), 
-      fNParamsDip(src.fNParamsDip), 
-      fSegZSol(0), 
-      fSegRSol(0), 
-      fNSegRSol(0), 
-      fSegZIdSol(0),
-      fMinZSol(src.fMinZSol), 
-      fMaxZSol(src.fMaxZSol), 
-      fMaxRSol(src.fMaxRSol), 
-      fParamsSol(0), 
-      fParamsDip(0)
-{
-// Copy constructor
-  if (src.fSegZSol) {
-    fSegZSol = new Float_t[fNSegZSol];
-    for (int i=fNSegZSol;i--;) fSegZSol[i] = src.fSegZSol[i];
-  }
-  if (src.fSegRSol) {
-    fSegRSol = new Float_t[fNParamsSol];
-    for (int i=fNParamsSol;i--;) fSegRSol[i] = src.fSegRSol[i];
-  }
-  if (src.fNSegRSol) {
-    fNSegRSol = new Int_t[fNSegZSol];
-    for (int i=fNSegZSol;i--;) fNSegRSol[i] = src.fNSegRSol[i];
+void AliMagFCheb::CopyFrom(const AliMagFCheb& src) 
+{ 
+  Clear();
+  SetName(src.GetName());
+  SetTitle(src.GetTitle());
+  fNParamsSol    = src.fNParamsSol;
+  fNSegZSol      = src.fNSegZSol;
+  fNParamsTPCInt = src.fNParamsTPCInt;
+  fNSegZTPCInt   = src.fNSegZTPCInt; 
+  fNParamsDip    = src.fNParamsDip;
+  //
+  fNZSegDip      = src.fNZSegDip;
+  fNYSegDip      = src.fNYSegDip;
+  fNXSegDip      = src.fNXSegDip;  
+  //
+  fMinZSol       = src.fMinZSol; 
+  fMaxZSol       = src.fMaxZSol;
+  fMaxRSol       = src.fMaxRSol; 
+  //
+  fMinZDip       = src.fMinZDip;
+  fMaxZDip       = src.fMaxZDip;
+  //
+  fMinZTPCInt    = src.fMinZTPCInt;
+  fMaxZTPCInt    = src.fMaxZTPCInt;
+  fMaxRTPCInt    = src.fMaxRTPCInt; 
+  // 
+  if (src.fNParamsSol) {
+    memcpy(fSegZSol  = new Float_t[fNSegZSol], src.fSegZSol, sizeof(Float_t)*fNSegZSol);
+    memcpy(fSegRSol  = new Float_t[fNParamsSol], src.fSegRSol, sizeof(Float_t)*fNParamsSol);
+    memcpy(fNSegRSol = new Int_t[fNSegZSol], src.fNSegRSol, sizeof(Int_t)*fNSegZSol);
+    memcpy(fSegZIdSol= new Int_t[fNSegZSol], src.fSegZIdSol, sizeof(Int_t)*fNSegZSol);
+    fParamsSol       = new TObjArray(fNParamsSol);
+    for (int i=0;i<fNParamsSol;i++) fParamsSol->AddAtAndExpand(new AliCheb3D(*src.GetParamSol(i)),i);
   }
-  if (src.fSegZIdSol) {
-    fSegZIdSol = new Int_t[fNSegZSol];
-    for (int i=fNSegZSol;i--;) fSegZIdSol[i] = src.fSegZIdSol[i];
-  }
-  if (src.fParamsSol) {
-    fParamsSol = new TObjArray(fNParamsSol);
-    for (int i=0;i<fNParamsSol;i++) {
-      AliCheb3D* pr = src.GetParamSol(i);
-      if (pr) fParamsSol->AddAtAndExpand(new AliCheb3D(*pr),i);
-    }
+  //
+  if (src.fNParamsDip) {
+    memcpy(fSegZDip   = new Float_t[fNZSegDip], src.fSegZDip, sizeof(Float_t)*fNZSegDip);
+    memcpy(fSegYDip   = new Float_t[fNYSegDip], src.fSegYDip, sizeof(Float_t)*fNYSegDip);
+    memcpy(fSegXDip   = new Float_t[fNXSegDip], src.fSegZDip, sizeof(Float_t)*fNXSegDip);
+    memcpy(fBegSegYDip= new Int_t[fNZSegDip], src.fBegSegYDip, sizeof(Int_t)*fNZSegDip);
+    memcpy(fNSegYDip  = new Int_t[fNZSegDip], src.fNSegYDip, sizeof(Int_t)*fNZSegDip);
+    memcpy(fBegSegXDip= new Int_t[fNYSegDip], src.fBegSegXDip, sizeof(Int_t)*fNYSegDip);
+    memcpy(fNSegXDip  = new Int_t[fNYSegDip], src.fNSegXDip, sizeof(Int_t)*fNYSegDip);
+    memcpy(fSegIDDip  = new Int_t[fNXSegDip], src.fSegIDDip, sizeof(Int_t)*fNXSegDip);
+    fParamsDip        = new TObjArray(fNParamsDip);
+    for (int i=0;i<fNParamsDip;i++) fParamsDip->AddAtAndExpand(new AliCheb3D(*src.GetParamDip(i)),i);
   }
-  if (src.fParamsDip) {
-    fParamsDip = new TObjArray(fNParamsDip);
-    for (int i=0;i<fNParamsDip;i++) {
-      AliCheb3D* pr = src.GetParamDip(i);
-      if (pr) fParamsDip->AddAtAndExpand(new AliCheb3D(*pr),i);
-    }
+  //
+  if (src.fNParamsTPCInt) {
+    memcpy(fSegZTPCInt  = new Float_t[fNSegZTPCInt], src.fSegZTPCInt, sizeof(Float_t)*fNSegZTPCInt);
+    memcpy(fSegRTPCInt  = new Float_t[fNParamsTPCInt], src.fSegRTPCInt, sizeof(Float_t)*fNParamsTPCInt);
+    memcpy(fNSegRTPCInt = new Int_t[fNSegZTPCInt], src.fNSegRTPCInt, sizeof(Int_t)*fNSegZTPCInt);
+    memcpy(fSegZIdTPCInt= new Int_t[fNSegZTPCInt], src.fSegZIdTPCInt, sizeof(Int_t)*fNSegZTPCInt);
+    fParamsTPCInt       = new TObjArray(fNParamsTPCInt);
+    for (int i=0;i<fNParamsTPCInt;i++) fParamsTPCInt->AddAtAndExpand(new AliCheb3D(*src.GetParamTPCInt(i)),i);
   }
   //
 }
 
-
+//__________________________________________________________________________________________
 AliMagFCheb& AliMagFCheb::operator=(const AliMagFCheb& rhs)
 {
-// Assignment operator
-    if (this != &rhs) {  
-       Clear();
-       SetName(rhs.GetName());
-       SetTitle(rhs.GetTitle());
-       fNParamsSol = rhs.fNParamsSol;
-       fNSegZSol   = rhs.fNSegZSol;
-       fMinZSol    = rhs.fMinZSol;
-       fMaxZSol    = rhs.fMaxZSol;
-       fMaxRSol    = rhs.fMaxRSol;
-       fNParamsDip = rhs.fNParamsDip;
-       fSegZSol    = fSegRSol    = 0; 
-       fNSegRSol   = fSegZIdSol  = 0;
-       fParamsSol  = fParamsDip  = 0;  
-       //
-       if (rhs.fSegZSol) {
-           fSegZSol = new Float_t[fNSegZSol];
-           for (int i=fNSegZSol;i--;) fSegZSol[i] = rhs.fSegZSol[i];
-       }
-       if (rhs.fSegRSol) {
-           fSegRSol = new Float_t[fNParamsSol];
-           for (int i=fNParamsSol;i--;) fSegRSol[i] = rhs.fSegRSol[i];
-       }
-       if (rhs.fNSegRSol) {
-           fNSegRSol = new Int_t[fNSegZSol];
-           for (int i=fNSegZSol;i--;) fNSegRSol[i] = rhs.fNSegRSol[i];
-       }
-       if (rhs.fSegZIdSol) {
-           fSegZIdSol = new Int_t[fNSegZSol];
-           for (int i=fNSegZSol;i--;) fSegZIdSol[i] = rhs.fSegZIdSol[i];
-       }
-       if (rhs.fParamsSol) {
-           fParamsSol = new TObjArray(fNParamsSol);
-           for (int i=0;i<fNParamsSol;i++) {
-               AliCheb3D* pr = rhs.GetParamSol(i);
-               if (pr) fParamsSol->AddAtAndExpand(new AliCheb3D(*pr),i);
-           }
-       }
-       if (rhs.fParamsDip) {
-           fParamsDip = new TObjArray(fNParamsDip);
-           for (int i=0;i<fNParamsDip;i++) {
-               AliCheb3D* pr = rhs.GetParamDip(i);
-               if (pr) fParamsDip->AddAtAndExpand(new AliCheb3D(*pr),i);
-           }
-       }
-    }
-    return *this;  
-    //
+  if (this != &rhs) {  
+    Clear();
+    CopyFrom(rhs);
+  }
+  return *this;  
+  //
 }
 
-
 //__________________________________________________________________________________________
-AliMagFCheb::~AliMagFCheb()
+void AliMagFCheb::Clear(Option_t *)
 {
   if (fNParamsSol) {
-    delete fParamsSol;
+    delete   fParamsSol;
     delete[] fSegZSol;
     delete[] fSegRSol;
     delete[] fNSegRSol;
     delete[] fSegZIdSol;
   }
   //
-  // Dipole part ...
+  if (fNParamsTPCInt) {
+    delete   fParamsTPCInt;
+    delete[] fSegZTPCInt;
+    delete[] fSegRTPCInt;
+    delete[] fNSegRTPCInt;
+    delete[] fSegZIdTPCInt;
+  }
+  // 
   if (fNParamsDip) {
-    delete fParamsDip;
+    delete   fParamsDip;
+    delete[] fSegZDip;
+    delete[] fSegYDip;
+    delete[] fSegXDip;
+    delete[] fBegSegYDip;
+    delete[] fNSegYDip;
+    delete[] fBegSegXDip;
+    delete[] fNSegXDip;
+    delete[] fSegIDDip;
   }
+  fNParamsSol = fNParamsTPCInt = fNParamsDip = fNZSegDip = fNYSegDip = fNXSegDip = 0;
+  fNSegZSol = fNSegZTPCInt = 0;
+  fMinZSol = fMinZDip = fMinZTPCInt = 1e6;
+  fMaxZSol = fMaxZDip = fMaxZTPCInt = fMaxRSol = fMaxRTPCInt = -1e6;
+  //
 }
 
 //__________________________________________________________________________________________
-void AliMagFCheb::Init0()
+void AliMagFCheb::Field(Float_t *xyz, Float_t *b) const
 {
-  // Solenoid part
-  fNParamsSol = 0;
-  fNSegZSol   = 0;
+  // compute field in cartesian coordinates. If point is outside of the parameterized region
+  // get it at closest valid point
+  static float rphiz[3];
+  //
+#ifndef _BRING_TO_BOUNDARY_  // exact matching to fitted volume is requested
+  if ( !(xyz[2]>=GetMinZSol()&&xyz[2]<=GetMaxZSol()) && 
+       !(xyz[2]>=GetMinZDip()&&xyz[2]<=GetMaxZDip())  ) {for (int i=3;i--;) b[i]=0; return;}
+#endif
+  //
+  if (xyz[2]<fMaxZDip) {    // dipole part?
+#ifndef _BRING_TO_BOUNDARY_
+    AliCheb3D* par = GetParamDip(FindDipSegment(xyz));
+    if (par->IsInside(xyz)) {par->Eval(xyz,b); return;}
+    for (int i=3;i--;) b[i]=0; return;
+#else
+    GetParamDip(FindDipSegment(xyz))->Eval(xyz,b); return;  
+#endif
+  }
+  //
+  // Sol region: convert coordinates to cyl system
+  CartToCyl(xyz,rphiz);
+#ifndef _BRING_TO_BOUNDARY_
+  if (rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;}
+#endif
+  //
+  FieldCylSol(rphiz,b);
+  //
+  // convert field to cartesian system
+  CylToCartCylB(rphiz, b,b);
+  //
+}
+
+//__________________________________________________________________________________________
+void AliMagFCheb::GetTPCInt(Float_t *xyz, Float_t *b) const
+{
+  // compute TPC region field integral in cartesian coordinates.
+  // If point is outside of the parameterized region get it at closeset valid point
+  static float rphiz[3];
+  //
+  // TPCInt region
+  // convert coordinates to cyl system
+  CartToCyl(xyz,rphiz);
+#ifndef _BRING_TO_BOUNDARY_
+  if ( (rphiz[2]>GetMaxZTPCInt()||rphiz[2]<GetMinZTPCInt()) ||
+       rphiz[0]>GetMaxRTPCInt()) {for (int i=3;i--;) b[i]=0; return;}
+#endif
+  //
+  GetTPCIntCyl(rphiz,b);
+  //
+  // convert field to cartesian system
+  CylToCartCylB(rphiz, b,b);
+  //
+}
+
+//__________________________________________________________________________________________
+void AliMagFCheb::FieldCylSol(Float_t *rphiz, Float_t *b) const
+{
+  // compute Solenoid field in Cylindircal coordinates
+  // note: if the point is outside the volume get the field in closest parameterized point
+  float &r = rphiz[0];
+  float &z = rphiz[2];
+  int SolZId = 0;
+  while (z>fSegZSol[SolZId] && SolZId<fNSegZSol) ++SolZId;    // find Z segment
+  int SolRId = fSegZIdSol[SolZId];        // first R segment for this Z
+  int SolRMax = SolRId + fNSegRSol[SolZId];
+  while (r>fSegRSol[SolRId] && SolRId<SolRMax) ++SolRId;    // find R segment
+  GetParamSol( SolRId )->Eval(rphiz,b);
+  //
+}
+
+//__________________________________________________________________________________________
+void AliMagFCheb::GetTPCIntCyl(Float_t *rphiz, Float_t *b) const
+{
+  // compute field integral in TPC region in Cylindircal coordinates
+  // note: the check for the point being inside the parameterized region is done outside
+  float &r = rphiz[0];
+  float &z = rphiz[2];
+  int TPCIntZId = 0;
+  while (z>fSegZTPCInt[TPCIntZId] && TPCIntZId<fNSegZTPCInt) ++TPCIntZId;    // find Z segment
+  int TPCIntRId = fSegZIdTPCInt[TPCIntZId];        // first R segment for this Z
+  int TPCIntRIdMax = TPCIntRId + fNSegRTPCInt[TPCIntZId];
+  while (r>fSegRTPCInt[TPCIntRId] && TPCIntRId<TPCIntRIdMax) ++TPCIntRId;    // find R segment
+  GetParamTPCInt( TPCIntRId )->Eval(rphiz,b);
+  //
+}
+
+
+//__________________________________________________________________________________________
+void AliMagFCheb::Print(Option_t *) const
+{
+  printf("Alice magnetic field parameterized by Chebyshev polynomials\n");
+  printf("Segmentation for Solenoid (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZSol,fMaxZSol,fMaxRSol);
+  //
+  if (fParamsSol) fParamsSol->Print();
+  /*
+  for (int iz=0;iz<fNSegZSol;iz++) {
+    AliCheb3D* param = GetParamSol( fSegZIdSol[iz] );
+    printf("*** Z Segment %2d (%+7.2f<Z<%+7.2f)\t***\n",iz,param->GetBoundMin(2),param->GetBoundMax(2));
+    for (int ir=0;ir<fNSegRSol[iz];ir++) {
+      param = GetParamSol( fSegZIdSol[iz]+ir );
+      printf("    R Segment %2d (%+7.2f<R<%+7.2f, Precision: %.1e) (ID=%2d)\n",ir, param->GetBoundMin(0),
+            param->GetBoundMax(0),param->GetPrecision(),fSegZIdSol[iz]+ir);
+    }
+  }
+  */
+  //
+  printf("Segmentation for TPC field integral (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZTPCInt,fMaxZTPCInt,fMaxRTPCInt);
+  //
+  if (fParamsTPCInt) fParamsTPCInt->Print();
+  /*
+  for (int iz=0;iz<fNSegZTPCInt;iz++) {
+    AliCheb3D* param = GetParamTPCInt( fSegZIdTPCInt[iz] );
+    printf("*** Z Segment %2d (%+7.2f<Z<%+7.2f)\t***\n",iz,param->GetBoundMin(2),param->GetBoundMax(2));
+    for (int ir=0;ir<fNSegRTPCInt[iz];ir++) {
+      param = GetParamTPCInt( fSegZIdTPCInt[iz]+ir );
+      printf("    R Segment %2d (%+7.2f<R<%+7.2f, Precision: %.1e) (ID=%2d)\n",ir, param->GetBoundMin(0),
+            param->GetBoundMax(0),param->GetPrecision(),fSegZIdTPCInt[iz]+ir);
+    }
+  }
+  */
+  //
+  printf("Segmentation for Dipole (%+.2f<Z<%+.2f cm)\n",fMinZDip,fMaxZDip);
+  if (fParamsDip) fParamsDip->Print();
+  //
+  
+
+}
+
+//_______________________________________________
+void AliMagFCheb::LoadData(const char* inpfile)
+{
+  // read coefficients data from the text file
+  //
+  TString strf = inpfile;
+  gSystem->ExpandPathName(strf);
+  FILE* stream = fopen(strf,"r");
+  if (!stream) {
+    printf("Did not find input file %s\n",strf.Data());
+    return;
+  }
+  //
+  TString buffs;
+  AliCheb3DCalc::ReadLine(buffs,stream);
+  if (!buffs.BeginsWith("START")) {
+    Error("LoadData","Expected: \"START <name>\", found \"%s\"\nStop\n",buffs.Data());
+    exit(1);
+  }
+  if (buffs.First(' ')>0) SetName(buffs.Data()+buffs.First(' ')+1);
+  //
+  // Solenoid part    -----------------------------------------------------------
+  AliCheb3DCalc::ReadLine(buffs,stream);
+  if (!buffs.BeginsWith("START SOLENOID")) {
+    Error("LoadData","Expected: \"START SOLENOID\", found \"%s\"\nStop\n",buffs.Data());
+    exit(1);
+  }
+  AliCheb3DCalc::ReadLine(buffs,stream); // nparam
+  int nparSol = buffs.Atoi(); 
+  //
+  for (int ip=0;ip<nparSol;ip++) {
+    AliCheb3D* cheb = new AliCheb3D();
+    cheb->LoadData(stream);
+    AddParamSol(cheb);
+  }
+  //
+  AliCheb3DCalc::ReadLine(buffs,stream);
+  if (!buffs.BeginsWith("END SOLENOID")) {
+    Error("LoadData","Expected \"END SOLENOID\", found \"%s\"\nStop\n",buffs.Data());
+    exit(1);
+  }
+  //
+  // TPCInt part     -----------------------------------------------------------
+  AliCheb3DCalc::ReadLine(buffs,stream);
+  if (!buffs.BeginsWith("START TPCINT")) {
+    Error("LoadData","Expected: \"START TPCINT\", found \"%s\"\nStop\n",buffs.Data());
+    exit(1);
+  }
+  AliCheb3DCalc::ReadLine(buffs,stream); // nparam
+  int nparTPCInt = buffs.Atoi(); 
+  //
+  for (int ip=0;ip<nparTPCInt;ip++) {
+    AliCheb3D* cheb = new AliCheb3D();
+    cheb->LoadData(stream);
+    AddParamTPCInt(cheb);
+  }
+  //
+  AliCheb3DCalc::ReadLine(buffs,stream);
+  if (!buffs.BeginsWith("END TPCINT")) {
+    Error("LoadData","Expected \"END TPCINT\", found \"%s\"\nStop\n",buffs.Data());
+    exit(1);
+  }
+  //
+  // Dipole part    -----------------------------------------------------------
+  AliCheb3DCalc::ReadLine(buffs,stream);
+  if (!buffs.BeginsWith("START DIPOLE")) {
+    Error("LoadData","Expected: \"START DIPOLE\", found \"%s\"\nStop\n",buffs.Data());
+    exit(1);
+  }
+  AliCheb3DCalc::ReadLine(buffs,stream); // nparam
+  int nparDip = buffs.Atoi();  
+  //
+  for (int ip=0;ip<nparDip;ip++) {
+    AliCheb3D* cheb = new AliCheb3D();
+    cheb->LoadData(stream);
+    AddParamDip(cheb);
+  }
   //
-  fSegZSol    = 0;
-  fSegRSol    = 0;
+  AliCheb3DCalc::ReadLine(buffs,stream);
+  if (!buffs.BeginsWith("END DIPOLE")) {
+    Error("LoadData","Expected \"END DIPOLE\", found \"%s\"\nStop\n",GetName(),buffs.Data());
+    exit(1);
+  }
+  //
+  AliCheb3DCalc::ReadLine(buffs,stream);
+  if (!buffs.BeginsWith("END") || !buffs.Contains(GetName())) {
+    Error("LoadData","Expected: \"END %s\", found \"%s\"\nStop\n",GetName(),buffs.Data());
+    exit(1);
+  }
+  //
+  // ---------------------------------------------------------------------------
+  fclose(stream);
+  BuildTableSol();
+  BuildTableTPCInt();
+  BuildTableDip();
+  printf("Loaded magnetic field \"%s\" from %s\n",GetName(),strf.Data());
   //
-  fNSegRSol   = 0;
-  fSegZIdSol  = 0;
+}
+
+//_________________________________________________________________________
+Int_t AliMagFCheb::FindDipSegment(float *xyz) const
+{
+  // find the segment containing point xyz. If it is outside find the closest segment 
+  int xid,yid,zid = TMath::BinarySearch(fNZSegDip,fSegZDip,xyz[2]); // find zsegment
+  int ysegBeg = fBegSegYDip[zid];
   //
-  fMinZSol = fMaxZSol = fMaxRSol = fMaxRSol = 0;
-  fParamsSol = 0;
+  for (yid=0;yid<fNSegYDip[zid];yid++) if (xyz[1]<fSegYDip[ysegBeg+yid]) break;
+  if ( --yid < 0 ) yid = 0;
+  yid +=  ysegBeg;
   //
-  // Dipole part ...
-  fNParamsDip = 0;
-  fParamsDip  = 0;
+  int xsegBeg = fBegSegXDip[yid];
+  for (xid=0;xid<fNSegXDip[yid];xid++) if (xyz[0]<fSegXDip[xsegBeg+xid]) break;
+  if ( --xid < 0) xid = 0;
+  xid +=  xsegBeg;
   //
+  return fSegIDDip[xid];
+}
+
+//_______________________________________________
+#ifdef  _INC_CREATION_ALICHEB3D_
+
+
+//__________________________________________________________________________________________
+AliMagFCheb::AliMagFCheb(const char* inputFile) : 
+  fNParamsSol(0),
+  fNSegZSol(0),
+  fNParamsTPCInt(0),
+  fNSegZTPCInt(0),
+  fNParamsDip(0),
+//
+  fNZSegDip(0),
+  fNYSegDip(0),
+  fNXSegDip(0),
+//
+  fSegZSol(0),
+  fSegRSol(0),
+//
+  fSegZTPCInt(0),
+  fSegRTPCInt(0),
+//
+  fSegZDip(0),
+  fSegYDip(0),
+  fSegXDip(0),
+//
+  fNSegRSol(0),
+  fSegZIdSol(0),
+//
+  fNSegRTPCInt(0),
+  fSegZIdTPCInt(0),
+//
+  fBegSegYDip(0),
+  fNSegYDip(0),
+  fBegSegXDip(0),
+  fNSegXDip(0),
+  fSegIDDip(0),
+//
+  fMinZSol(0),
+  fMaxZSol(0),
+  fMaxRSol(0), 
+//
+  fMinZDip(0),
+  fMaxZDip(0),
+//
+  fMinZTPCInt(0),
+  fMaxZTPCInt(0),
+  fMaxRTPCInt(0), 
+//
+  fParamsSol(0),
+  fParamsDip(0),
+  fParamsTPCInt(0)
+//
+//
+{
+  LoadData(inputFile);
 }
 
 //__________________________________________________________________________________________
@@ -236,6 +568,18 @@ void AliMagFCheb::AddParamSol(AliCheb3D* param)
 }
 
 //__________________________________________________________________________________________
+void AliMagFCheb::AddParamTPCInt(AliCheb3D* param)
+{
+  // adds new parameterization piece for TPCInt
+  // NOTE: pieces must be added strictly in increasing R then increasing Z order
+  //
+  if (!fParamsTPCInt) fParamsTPCInt = new TObjArray();
+  fParamsTPCInt->Add(param);
+  fNParamsTPCInt++;
+  //
+}
+
+//__________________________________________________________________________________________
 void AliMagFCheb::AddParamDip(AliCheb3D* param)
 {
   // adds new parameterization piece for Dipole
@@ -247,12 +591,127 @@ void AliMagFCheb::AddParamDip(AliCheb3D* param)
 }
 
 //__________________________________________________________________________________________
+void AliMagFCheb::ResetTPCInt()
+{
+  // clean TPC field integral (used for update)
+  if (!fNParamsTPCInt) return;
+  delete fParamsTPCInt;  
+  delete[] fSegZTPCInt;
+  delete[] fSegRTPCInt;
+  delete[] fNSegRTPCInt;
+  delete[] fSegZIdTPCInt;
+  //
+  fNParamsTPCInt = 0; 
+  fNSegZTPCInt   = 0;
+  fSegZTPCInt    = 0; 
+  fSegRTPCInt    = 0; 
+  fNSegRTPCInt   = 0; 
+  fSegZIdTPCInt  = 0;
+  fMinZTPCInt    = 0; 
+  fMaxZTPCInt    = 0; 
+  fMaxRTPCInt    = 0; 
+  fParamsTPCInt  = 0;
+  //
+}
+
+//__________________________________________________
+void AliMagFCheb::BuildTableDip()
+{
+  //
+  TArrayF segY,segX;
+  TArrayI begSegYDip,begSegXDip;
+  TArrayI nsegYDip,nsegXDip;
+  TArrayI segID;
+  float *tmpSegZ,*tmpSegY,*tmpSegX;
+  //
+  // create segmentation in Z
+  fNZSegDip = SegmentDipDimension(&tmpSegZ, fParamsDip, fNParamsDip, 2, 1,-1, 1,-1, 1,-1) - 1;
+  fNYSegDip = 0;
+  fNXSegDip = 0;
+  //
+  // for each Z slice create segmentation in Y
+  begSegYDip.Set(fNZSegDip);
+  nsegYDip.Set(fNZSegDip);
+  float xyz[3];
+  for (int iz=0;iz<fNZSegDip;iz++) {
+    printf("\nZSegment#%d  %+e : %+e\n",iz,tmpSegZ[iz],tmpSegZ[iz+1]);
+    int ny = SegmentDipDimension(&tmpSegY, fParamsDip, fNParamsDip, 1, 
+                                1,-1, 1,-1, tmpSegZ[iz],tmpSegZ[iz+1]) - 1;
+    segY.Set(ny + fNYSegDip);
+    for (int iy=0;iy<ny;iy++) segY[fNYSegDip+iy] = tmpSegY[iy];
+    begSegYDip[iz] = fNYSegDip;
+    nsegYDip[iz] = ny;
+    printf(" Found %d YSegments, to start from %d\n",ny, begSegYDip[iz]);
+    //
+    // for each slice in Z and Y create segmentation in X
+    begSegXDip.Set(fNYSegDip+ny);
+    nsegXDip.Set(fNYSegDip+ny);
+    xyz[2] = (tmpSegZ[iz]+tmpSegZ[iz+1])/2.; // mean Z of this segment
+    //
+    for (int iy=0;iy<ny;iy++) {
+      int isg = fNYSegDip+iy;
+      printf("\n   YSegment#%d  %+e : %+e\n",iy, tmpSegY[iy],tmpSegY[iy+1]);
+      int nx = SegmentDipDimension(&tmpSegX, fParamsDip, fNParamsDip, 0, 
+                                  1,-1, tmpSegY[iy],tmpSegY[iy+1], tmpSegZ[iz],tmpSegZ[iz+1]) - 1;
+      //
+      segX.Set(nx + fNXSegDip);
+      for (int ix=0;ix<nx;ix++) segX[fNXSegDip+ix] = tmpSegX[ix];
+      begSegXDip[isg] = fNXSegDip;
+      nsegXDip[isg] = nx;
+      printf("   Found %d XSegments, to start from %d\n",nx, begSegXDip[isg]);
+      //
+      segID.Set(fNXSegDip+nx);
+      //
+      // find corresponding params
+      xyz[1] = (tmpSegY[iy]+tmpSegY[iy+1])/2.; // mean Y of this segment
+      //
+      for (int ix=0;ix<nx;ix++) {
+       xyz[0] = (tmpSegX[ix]+tmpSegX[ix+1])/2.; // mean X of this segment
+       for (int ipar=0;ipar<fNParamsDip;ipar++) {
+         AliCheb3D* cheb = (AliCheb3D*) fParamsDip->At(ipar);
+         if (!cheb->IsInside(xyz)) continue;
+         segID[fNXSegDip+ix] = ipar;
+         break;
+       }
+      }
+      fNXSegDip += nx;
+      //
+      delete[] tmpSegX;
+    }
+    delete[] tmpSegY;
+    fNYSegDip += ny;
+  }
+  //
+  fMinZDip = tmpSegZ[0];
+  fMaxZDip = tmpSegZ[fNZSegDip];
+  fSegZDip    = new Float_t[fNZSegDip];
+  for (int i=fNZSegDip;i--;) fSegZDip[i] = tmpSegZ[i];
+  delete[] tmpSegZ;
+  //
+  fSegYDip    = new Float_t[fNYSegDip];
+  fSegXDip    = new Float_t[fNXSegDip];
+  fBegSegYDip = new Int_t[fNZSegDip];
+  fNSegYDip   = new Int_t[fNZSegDip];
+  fBegSegXDip = new Int_t[fNYSegDip];
+  fNSegXDip   = new Int_t[fNYSegDip];
+  fSegIDDip   = new Int_t[fNXSegDip];
+  //
+  for (int i=fNYSegDip;i--;) fSegYDip[i] = segY[i];
+  for (int i=fNXSegDip;i--;) fSegXDip[i] = segX[i];
+  for (int i=fNZSegDip;i--;) {fBegSegYDip[i] = begSegYDip[i]; fNSegYDip[i] = nsegYDip[i];}
+  for (int i=fNYSegDip;i--;) {fBegSegXDip[i] = begSegXDip[i]; fNSegXDip[i] = nsegXDip[i];}
+  for (int i=fNXSegDip;i--;) {fSegIDDip[i]   = segID[i];}
+  //
+}
+
+//__________________________________________________________________________________________
 void AliMagFCheb::BuildTableSol()
 {
   // build the indexes for each parameterization of Solenoid
   //
   const float kSafety=0.001;
   //
+  if (fNParamsSol<1) return;
   fSegRSol   = new Float_t[fNParamsSol];
   float *tmpbufF  = new float[fNParamsSol+1];
   int   *tmpbufI  = new int[fNParamsSol+1];
@@ -290,63 +749,49 @@ void AliMagFCheb::BuildTableSol()
 }
 
 //__________________________________________________________________________________________
-void AliMagFCheb::Field(Float_t *xyz, Float_t *b) const
+void AliMagFCheb::BuildTableTPCInt()
 {
-  // compute field in cartesian coordinates
-  float rphiz[3];
-  if (xyz[2]>GetMaxZSol() || xyz[2]<GetMinZSol()) {for (int i=3;i--;) b[i]=0; return;}
+  // build the indexes for each parameterization of TPC field integral
   //
-  // Sol region
-  // convert coordinates to cyl system
-  rphiz[0] = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
-  rphiz[1] = TMath::ATan2(xyz[1],xyz[0]);
-  rphiz[2] = xyz[2];
-  if (rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;}
+  const float kSafety=0.001;
   //
-  FieldCylSol(rphiz,b);
+  if (fNParamsTPCInt<1) return;
+  fSegRTPCInt   = new Float_t[fNParamsTPCInt];
+  float *tmpbufF  = new float[fNParamsTPCInt+1];
+  int   *tmpbufI  = new int[fNParamsTPCInt+1];
+  int   *tmpbufI1 = new int[fNParamsTPCInt+1];
   //
-  // convert field to cartesian system
-  float btr = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
-  float psiPLUSphi = rphiz[1] + TMath::ATan2(b[1],b[0]);
-  b[0] = btr*TMath::Cos(psiPLUSphi);
-  b[1] = btr*TMath::Sin(psiPLUSphi);
+  // count number of Z slices and number of R slices in each Z slice
+  for (int ip=0;ip<fNParamsTPCInt;ip++) {
+    if (ip==0 || (GetParamTPCInt(ip)->GetBoundMax(2)-GetParamTPCInt(ip-1)->GetBoundMax(2))>kSafety) { // new Z slice
+      tmpbufF[fNSegZTPCInt] = GetParamTPCInt(ip)->GetBoundMax(2);                                     // 
+      tmpbufI[fNSegZTPCInt] = 0;
+      tmpbufI1[fNSegZTPCInt++] = ip;
+    }
+    fSegRTPCInt[ip] = GetParamTPCInt(ip)->GetBoundMax(0);  // upper R
+    tmpbufI[fNSegZTPCInt-1]++;
+  }
   //
-}
-
-//__________________________________________________________________________________________
-void AliMagFCheb::FieldCylSol(Float_t *rphiz, Float_t *b) const
-{
-  // compute Solenoid field in Cylindircal coordinates
-  // note: the check for the point being inside the parameterized region is done outside
-  float &r = rphiz[0];
-  float &z = rphiz[2];
-  int SolZId = 0;
-  while (z>fSegZSol[SolZId]) ++SolZId;    // find Z segment
-  int SolRId = fSegZIdSol[SolZId];        // first R segment for this Z
-  while (r>fSegRSol[SolRId]) ++SolRId;    // find R segment
-  GetParamSol( SolRId )->Eval(rphiz,b);
+  fSegZTPCInt   = new Float_t[fNSegZTPCInt];
+  fSegZIdTPCInt = new Int_t[fNSegZTPCInt];
+  fNSegRTPCInt  = new Int_t[fNSegZTPCInt];
+  for (int iz=0;iz<fNSegZTPCInt;iz++) {
+    fSegZTPCInt[iz]   = tmpbufF[iz];
+    fNSegRTPCInt[iz]  = tmpbufI[iz];
+    fSegZIdTPCInt[iz] = tmpbufI1[iz];
+  }
+  //
+  fMinZTPCInt = GetParamTPCInt(0)->GetBoundMin(2);
+  fMaxZTPCInt = GetParamTPCInt(fNParamsTPCInt-1)->GetBoundMax(2);
+  fMaxRTPCInt = GetParamTPCInt(fNParamsTPCInt-1)->GetBoundMax(0);
+  //
+  delete[] tmpbufF;
+  delete[] tmpbufI;
+  delete[] tmpbufI1;
   //
-}
-
-//__________________________________________________________________________________________
-void AliMagFCheb::Print(Option_t *) const
-{
-  printf("Alice magnetic field parameterized by Chebyshev polynomials\n");
-  printf("Segmentation for Solenoid (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZSol,fMaxZSol,fMaxRSol);
   //
-  for (int iz=0;iz<fNSegZSol;iz++) {
-    AliCheb3D* param = GetParamSol( fSegZIdSol[iz] );
-    printf("*** Z Segment %2d (%+7.2f<Z<%+7.2f)\t***\n",iz,param->GetBoundMin(2),param->GetBoundMax(2));
-    for (int ir=0;ir<fNSegRSol[iz];ir++) {
-      param = GetParamSol( fSegZIdSol[iz]+ir );
-      printf("    R Segment %2d (%+7.2f<R<%+7.2f, Precision: %.1e) (ID=%2d)\n",ir, param->GetBoundMin(0),param->GetBoundMax(0),
-            param->GetPrecision(),fSegZIdSol[iz]+ir);
-    }
-  }
 }
 
-//_______________________________________________
-#ifdef  _INC_CREATION_ALICHEB3D_
 void AliMagFCheb::SaveData(const char* outfile) const
 {
   // writes coefficients data to output text file
@@ -354,13 +799,19 @@ void AliMagFCheb::SaveData(const char* outfile) const
   gSystem->ExpandPathName(strf);
   FILE* stream = fopen(strf,"w+");
   //
-  // Sol part
+  // Sol part    ---------------------------------------------------------
   fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName());
   fprintf(stream,"START SOLENOID\n#Number of pieces\n%d\n",fNParamsSol);
   for (int ip=0;ip<fNParamsSol;ip++) GetParamSol(ip)->SaveData(stream);
   fprintf(stream,"#\nEND SOLENOID\n");
   //
-  // Dip part
+  // TPCInt part ---------------------------------------------------------
+  fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName());
+  fprintf(stream,"START TPCINT\n#Number of pieces\n%d\n",fNParamsTPCInt);
+  for (int ip=0;ip<fNParamsTPCInt;ip++) GetParamTPCInt(ip)->SaveData(stream);
+  fprintf(stream,"#\nEND TPCINT\n");
+  //
+  // Dip part   ---------------------------------------------------------
   fprintf(stream,"START DIPOLE\n#Number of pieces\n%d\n",fNParamsDip);
   for (int ip=0;ip<fNParamsDip;ip++) GetParamDip(ip)->SaveData(stream);
   fprintf(stream,"#\nEND DIPOLE\n");
@@ -370,62 +821,45 @@ void AliMagFCheb::SaveData(const char* outfile) const
   fclose(stream);
   //
 }
-#endif
 
-//_______________________________________________
-void AliMagFCheb::LoadData(const char* inpfile)
+Int_t AliMagFCheb::SegmentDipDimension(float** seg,const TObjArray* par,int npar, int dim, 
+                                      float xmn,float xmx,float ymn,float ymx,float zmn,float zmx)
 {
-  // read coefficients data from the text file
-  //
-  TString strf = inpfile;
-  gSystem->ExpandPathName(strf);
-  FILE* stream = fopen(strf,"r");
-  if (!stream) {
-    printf("Did not find input file %s\n",strf.Data());
-    return;
-  }
-  //
-  TString buffs;
-  AliCheb3DCalc::ReadLine(buffs,stream);
-  if (!buffs.BeginsWith("START")) {Error("LoadData","Expected: \"START <name>\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
-  if (buffs.First(' ')>0) SetName(buffs.Data()+buffs.First(' ')+1);
-  //
-  // Solenoid part
-  AliCheb3DCalc::ReadLine(buffs,stream);
-  if (!buffs.BeginsWith("START SOLENOID")) {Error("LoadData","Expected: \"START SOLENOID\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
-  AliCheb3DCalc::ReadLine(buffs,stream); // nparam
-  int nparSol = buffs.Atoi(); 
-  //
-  for (int ip=0;ip<nparSol;ip++) {
-    AliCheb3D* cheb = new AliCheb3D();
-    cheb->LoadData(stream);
-    AddParamSol(cheb);
+  // find all boundaries in deimension dim for boxes in given region.
+  // if mn>mx for given projection the check is not done for it.
+  float *tmpC = new float[2*npar];
+  int *tmpInd = new int[2*npar];
+  int nseg0 = 0;
+  for (int ip=0;ip<npar;ip++) {
+    AliCheb3D* cheb = (AliCheb3D*) par->At(ip);
+    if (xmn<xmx && (cheb->GetBoundMin(0)>(xmx+xmn)/2 || cheb->GetBoundMax(0)<(xmn+xmx)/2)) continue;
+    if (ymn<ymx && (cheb->GetBoundMin(1)>(ymx+ymn)/2 || cheb->GetBoundMax(1)<(ymn+ymx)/2)) continue;
+    if (zmn<zmx && (cheb->GetBoundMin(2)>(zmx+zmn)/2 || cheb->GetBoundMax(2)<(zmn+zmx)/2)) continue;
+    //
+    tmpC[nseg0++] = cheb->GetBoundMin(dim);
+    tmpC[nseg0++] = cheb->GetBoundMax(dim);    
   }
-  //
-  AliCheb3DCalc::ReadLine(buffs,stream);
-  if (!buffs.BeginsWith("END SOLENOID")) {Error("LoadData","Expected \"END SOLENOID\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
-  //
-  // Dipole part
-  AliCheb3DCalc::ReadLine(buffs,stream);
-  if (!buffs.BeginsWith("START DIPOLE")) {Error("LoadData","Expected: \"START DIPOLE\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
-  AliCheb3DCalc::ReadLine(buffs,stream); // nparam
-  int nparDip = buffs.Atoi();  
-  //
-  for (int ip=0;ip<nparDip;ip++) {
-    AliCheb3D* cheb = new AliCheb3D();
-    cheb->LoadData(stream);
-    AddParamDip(cheb);
+  // range Dim's boundaries in increasing order
+  TMath::Sort(nseg0,tmpC,tmpInd,kFALSE);
+  // count number of really different Z's
+  int nseg = 0;
+  float cprev = -1e6;
+  for (int ip=0;ip<nseg0;ip++) {
+    if (TMath::Abs(cprev-tmpC[ tmpInd[ip] ])>1e-4) {
+      cprev = tmpC[ tmpInd[ip] ];
+      nseg++;
+    }
+    else tmpInd[ip] = -1; // supress redundant Z
   }
-  //
-  AliCheb3DCalc::ReadLine(buffs,stream);
-  if (!buffs.BeginsWith("END DIPOLE")) {Error("LoadData","Expected \"END DIPOLE\", found \"%s\"\nStop\n",GetName(),buffs.Data());exit(1);}
-  //
-  AliCheb3DCalc::ReadLine(buffs,stream);
-  if (!buffs.BeginsWith("END") || !buffs.Contains(GetName())) {Error("LoadData","Expected: \"END %s\", found \"%s\"\nStop\n",GetName(),buffs.Data());exit(1);}
-  //
-  fclose(stream);
-  BuildTableSol();
-  //  BuildDipTable();
-  printf("Loaded magnetic field \"%s\" from %s\n",GetName(),strf.Data());
-  //
+  // 
+  *seg  = new float[nseg]; // create final Z segmenations
+  nseg = 0;
+  for (int ip=0;ip<nseg0;ip++) if (tmpInd[ip]>=0) (*seg)[nseg++] = tmpC[ tmpInd[ip] ];
+  //
+  delete[] tmpC;
+  delete[] tmpInd;
+  return nseg;
 }
+
+#endif
+
index 46e9410..3d6b918 100644 (file)
@@ -1,13 +1,6 @@
-#ifndef ALIMAGFCHEB_H
-#define ALIMAGFCHEB_H
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * See cxx source for full Copyright notice                               */
-
-/* $Id$ */
-
 
 // Author: ruben.shahoyan@cern.ch   20/03/2007
-//
+
 ///////////////////////////////////////////////////////////////////////////////////
 //                                                                               //
 //  Wrapper for the set of mag.field parameterizations by Chebyshev polinomials  //
@@ -16,8 +9,7 @@
 //  For cylindrical coordinates/components:                                      //
 //    FieldCyl(float* rphiz, float* brphiz)                                      //
 //                                                                               //
-//  For the moment only the solenoid part is parameterized in the volume defined //
-//  by R<500, -550<Z<550 cm                                                      //
+//  The solenoid part is parameterized in the volume  R<500, -550<Z<550 cm       //
 //                                                                               //
 //  The region R<423 cm,  -343.3<Z<481.3 for 30kA and -343.3<Z<481.3 for 12kA    //
 //  is parameterized using measured data while outside the Tosca calculation     //
 //  If the querried point is outside the validity region no the return values    //
 //  for the field components are set to 0.                                       //
 //                                                                               //
+//  To obtain the field integral in the TPC region from given point to nearest   //
+//  cathod plane (+- 250 cm) use:                                                //
+//  GetTPCInt(float* xyz, float* bxyz);  for Cartesian frame                     //
+//  or                                                                           //
+//  GetTPCIntCyl(Float_t *rphiz, Float_t *b); for Cylindrical frame              //
+//                                                                               //
+//                                                                               //
+//  The units are kiloGauss and cm.                                              //
+//                                                                               //
 ///////////////////////////////////////////////////////////////////////////////////
+#ifndef _ALIMAGFCHEB_
+#define _ALIMAGFCHEB_
 
-
-#include <TSystem.h>
 #include <TNamed.h>
+#include <TSystem.h>
 #include "AliCheb3D.h"
-#include "AliCheb3DCalc.h"
 
 class AliMagFCheb: public TNamed
 {
  public:
-    AliMagFCheb();
-    AliMagFCheb(const char* inputFile);
-    AliMagFCheb(const AliMagFCheb &src);
-    AliMagFCheb& operator= (const AliMagFCheb &rhs);
-    
-   ~AliMagFCheb();
+  AliMagFCheb();
+  AliMagFCheb(const AliMagFCheb& src);
+  ~AliMagFCheb() {Clear();}
   //
-  void       AddParamSol(AliCheb3D* param);
-  void       AddParamDip(AliCheb3D* param);
-  void       BuildTableSol();
+  void       CopyFrom(const AliMagFCheb& src);
+  AliMagFCheb& operator=(const AliMagFCheb& rhs);
+  virtual void Clear(Option_t * = "");
   //
   Int_t      GetNParamsSol()                              const {return fNParamsSol;}
   Int_t      GetNSegZSol()                                const {return fNSegZSol;}
-  Int_t      GetNSegRSol(int iz)                          const {return iz<fNParamsSol ? fNSegRSol[iz]:0;}
-  Int_t      GetSegIDSol(int iz,int ir)                   const {return iz<fNParamsSol&&ir<fNSegRSol[iz] ? fSegZIdSol[iz]+ir:-1;}
+  //
+  Int_t      GetNParamsTPCInt()                           const {return fNParamsTPCInt;}
+  Int_t      GetNSegZTPCInt()                             const {return fNSegZTPCInt;}
+  //
+  Int_t      GetNParamsDip()                              const {return fNParamsDip;}
+  Int_t      GetNSegZDip()                                const {return fNZSegDip;}
+  //
   //
   Float_t    GetMinZSol()                                 const {return fMinZSol;}
   Float_t    GetMaxZSol()                                 const {return fMaxZSol;}
   Float_t    GetMaxRSol()                                 const {return fMaxRSol;}
+  //
+  Float_t    GetMinZDip()                                 const {return fMinZDip;}
+  Float_t    GetMaxZDip()                                 const {return fMaxZDip;}
+  //
+  Float_t    GetMinZTPCInt()                              const {return fMinZTPCInt;}
+  Float_t    GetMaxZTPCInt()                              const {return fMaxZTPCInt;}
+  Float_t    GetMaxRTPCInt()                              const {return fMaxRTPCInt;}
+  //
+  Int_t      FindDipSegment(float *xyz)                   const;
   AliCheb3D* GetParamSol(Int_t ipar)                      const {return (AliCheb3D*)fParamsSol->UncheckedAt(ipar);}
+  AliCheb3D* GetParamTPCInt(Int_t ipar)                   const {return (AliCheb3D*)fParamsTPCInt->UncheckedAt(ipar);}
   AliCheb3D* GetParamDip(Int_t ipar)                      const {return (AliCheb3D*)fParamsDip->UncheckedAt(ipar);}
   //
   void         LoadData(const char* inpfile);
@@ -65,40 +78,87 @@ class AliMagFCheb: public TNamed
   //
   virtual void Field(Float_t *xyz, Float_t *b)            const;
   virtual void FieldCyl(Float_t *rphiz, Float_t *b)       const;
-  virtual void GetTPCInt(Float_t */*xyz*/,     Float_t */*b*/)    const {;}
-  virtual void GetTPCIntCyl(Float_t */*rphiz*/, Float_t */*b*/)   const {;}
   //
+  virtual void GetTPCInt(Float_t *xyz, Float_t *b)        const;
+  virtual void GetTPCIntCyl(Float_t *rphiz, Float_t *b)   const;
   //
+  static void CylToCartCylB(float *rphiz,  float *brphiz,float *bxyz);
+  static void CylToCartCartB(float *xyz,   float *brphiz,float *bxyz);
+  static void CartToCylCartB(float *xyz,   float *bxyz, float *brphiz);
+  static void CartToCylCylB(float *rphiz,  float *bxyz, float *brphiz);
+  static void CartToCyl(float *xyz,  float *rphiz);
+  static void CylToCart(float *rphiz,float *xyz);
   //
 #ifdef  _INC_CREATION_ALICHEB3D_                          // see AliCheb3D.h for explanation
-  void         SaveData(const char* outfile)              const;
+  AliMagFCheb(const char* inputFile);
+  void       SaveData(const char* outfile)              const;
+  Int_t      SegmentDipDimension(float** seg,const TObjArray* par,int npar, int dim, 
+                                float xmn,float xmx,float ymn,float ymx,float zmn,float zmx);
+  //
+  void       AddParamSol(AliCheb3D* param);
+  void       AddParamTPCInt(AliCheb3D* param);
+  void       AddParamDip(AliCheb3D* param);
+  void       BuildTableDip();
+  void       BuildTableSol();
+  void       BuildTableTPCInt();
+  void       ResetTPCInt();
+
 #endif
   //
  protected:
-  void         Init0();
   virtual void FieldCylSol(Float_t *rphiz, Float_t *b)    const;
   //
  protected:
   //
   Int_t      fNParamsSol;            // Total number of parameterization pieces for Sol 
-  Int_t      fNSegZSol;              // Number of segments is Z
+  Int_t      fNSegZSol;              // Number of segments in Z for Solenoid field
+  //
+  Int_t      fNParamsTPCInt;         // Total number of parameterization pieces for TPC field integral 
+  Int_t      fNSegZTPCInt;           // Number of segments in Z for TPC field integral
   //
   Int_t      fNParamsDip;            // Total number of parameterization pieces for dipole 
+  Int_t      fNZSegDip;              // number of distinct Z segments in Dipole
+  Int_t      fNYSegDip;              // number of distinct Y segments in Dipole
+  Int_t      fNXSegDip;              // number of distinct X segments in Dipole
+  //
+  Float_t*   fSegZSol;               //[fNSegZSol]      upper boundaries of Z segments
+  Float_t*   fSegRSol;               //[fNParamsSol]    upper boundaries of R segments
+  //
+  Float_t*   fSegZTPCInt;            //[fNSegZTPCInt]    upper boundaries of Z segments
+  Float_t*   fSegRTPCInt;            //[fNParamsTPCInt]  upper boundaries of R segments
   //
-  Float_t*   fSegZSol;               //[fNSegZSol]       upper boundaries of Z segments
-  Float_t*   fSegRSol;               //[fNParamsSol]     upper boundaries of R segments
+  Float_t*   fSegZDip;               //[fNZSegDip] coordinates of distinct Z segments in Dipole
+  Float_t*   fSegYDip;               //[fNYSegDip] coordinated of Y segments for each Zsegment in Dipole
+  Float_t*   fSegXDip;               //[fNXSegDip] coordinated of X segments for each Ysegment in Dipole
   //
-  Int_t*     fNSegRSol;              //[fNSegZSol]       number of R segments for each Z segment
-  Int_t*     fSegZIdSol;             //[fNSegZSol]       Id of the first R segment of each Z segment in the fSegRSol...
+  Int_t*     fNSegRSol;              //[fNSegZSol]      number of R segments for each Z segment
+  Int_t*     fSegZIdSol;             //[fNSegZSol]      Id of the first R segment of each Z segment in the fSegRSol...
+  //
+  Int_t*     fNSegRTPCInt;           //[fNSegZTPCInt]   number of R segments for each Z segment
+  Int_t*     fSegZIdTPCInt;          //[fNSegZTPCInt]   Id of the first R segment of each Z segment in the fSegRTPCInt...
+  //
+  Int_t*     fBegSegYDip;            //[fNZSegDip] beginning of Y segments array for each Z segment
+  Int_t*     fNSegYDip;              //[fNZSegDip] number of Y segments for each Z segment
+  Int_t*     fBegSegXDip;            //[fNYSegDip] beginning of X segments array for each Y segment
+  Int_t*     fNSegXDip;              //[fNYSegDip] number of X segments for each Y segment
+  Int_t*     fSegIDDip;              //[fNXSegDip] ID of the dipole parameterization for given XYZ segment
   //
   Float_t    fMinZSol;               // Min Z of Sol parameterization (in CYL. coordinates)
   Float_t    fMaxZSol;               // Max Z of Sol parameterization (in CYL. coordinates)
   Float_t    fMaxRSol;               // Max R of Sol parameterization (in CYL. coordinates)
   //
+  Float_t    fMinZDip;               // Min Z of Dipole parameterization
+  Float_t    fMaxZDip;               // Max Z of Dipole parameterization
+  //
+  Float_t    fMinZTPCInt;            // Min Z of TPCInt parameterization (in CYL. coordinates)
+  Float_t    fMaxZTPCInt;            // Max Z of TPCInt parameterization (in CYL. coordinates)
+  Float_t    fMaxRTPCInt;            // Max R of TPCInt parameterization (in CYL. coordinates)
+  // 
   TObjArray* fParamsSol;             // Parameterization pieces for Solenoid field
   TObjArray* fParamsDip;             // Parameterization pieces for Dipole field
+  TObjArray* fParamsTPCInt;          // Parameterization pieces for Solenoid field integrals in TPC region
   //
-  ClassDef(AliMagFCheb,1)            // Wrapper class for the set of Chebishev parameterizations of Alice mag.field
+  ClassDef(AliMagFCheb,3)            // Wrapper class for the set of Chebishev parameterizations of Alice mag.field
   //
  };
 
@@ -107,8 +167,73 @@ class AliMagFCheb: public TNamed
 inline void AliMagFCheb::FieldCyl(Float_t *rphiz, Float_t *b) const
 {
   // compute field in Cylindircal coordinates
-  if (rphiz[2]<GetMinZSol() || rphiz[2]>GetMaxZSol() || rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;}
+  //  if (rphiz[2]<GetMinZSol() || rphiz[2]>GetMaxZSol() || rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;}
   FieldCylSol(rphiz,b);
 }
 
+//__________________________________________________________________________________________________
+inline void AliMagFCheb::CylToCartCylB(float *rphiz, float *brphiz,float *bxyz)
+{
+  // convert field in cylindrical coordinates to cartesian system, point is in cyl.system
+  float btr = TMath::Sqrt(brphiz[0]*brphiz[0]+brphiz[1]*brphiz[1]);
+  float psiPLUSphi = TMath::ATan2(brphiz[1],brphiz[0]) + rphiz[1];
+  bxyz[0] = btr*TMath::Cos(psiPLUSphi);
+  bxyz[1] = btr*TMath::Sin(psiPLUSphi);
+  bxyz[2] = brphiz[2];
+  //
+}
+
+//__________________________________________________________________________________________________
+inline void AliMagFCheb::CylToCartCartB(float *xyz, float *brphiz,float *bxyz)
+{
+  // convert field in cylindrical coordinates to cartesian system, point is in cart.system
+  float btr = TMath::Sqrt(brphiz[0]*brphiz[0]+brphiz[1]*brphiz[1]);
+  float phiPLUSpsi = TMath::ATan2(xyz[1],xyz[0]) +  TMath::ATan2(brphiz[1],brphiz[0]);
+  bxyz[0] = btr*TMath::Cos(phiPLUSpsi);
+  bxyz[1] = btr*TMath::Sin(phiPLUSpsi);
+  bxyz[2] = brphiz[2];
+  //
+}
+
+//__________________________________________________________________________________________________
+inline void AliMagFCheb::CartToCylCartB(float *xyz ,float *bxyz, float *brphiz)
+{
+  // convert field in cylindrical coordinates to cartesian system, poin is in cart.system
+  float btr = TMath::Sqrt(bxyz[0]*bxyz[0]+bxyz[1]*bxyz[1]);
+  float psiMINphi = TMath::ATan2(bxyz[1],bxyz[0]) - TMath::ATan2(xyz[1],xyz[0]);
+  //
+  brphiz[0] = btr*TMath::Cos(psiMINphi);
+  brphiz[1] = btr*TMath::Sin(psiMINphi);
+  brphiz[2] = bxyz[2];
+  //
+}
+
+//__________________________________________________________________________________________________
+inline void AliMagFCheb::CartToCylCylB(float *rphiz,float *bxyz, float *brphiz)
+{
+  // convert field in cylindrical coordinates to cartesian system, point is in cyl.system
+  float btr = TMath::Sqrt(bxyz[0]*bxyz[0]+bxyz[1]*bxyz[1]);
+  float psiMINphi =  TMath::ATan2(bxyz[1],bxyz[0]) - rphiz[1];
+  brphiz[0] = btr*TMath::Cos(psiMINphi);
+  brphiz[1] = btr*TMath::Sin(psiMINphi);
+  brphiz[2] = bxyz[2];
+  //
+}
+
+//__________________________________________________________________________________________________
+inline void AliMagFCheb::CartToCyl(float *xyz,float *rphiz)
+{
+  rphiz[0] = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
+  rphiz[1] = TMath::ATan2(xyz[1],xyz[0]);
+  rphiz[2] = xyz[2];
+}
+
+//__________________________________________________________________________________________________
+inline void AliMagFCheb::CylToCart(float *rphiz,float *xyz)
+{
+  xyz[0] = rphiz[0]*TMath::Cos(rphiz[1]);
+  xyz[1] = rphiz[0]*TMath::Sin(rphiz[1]);
+  xyz[2] = rphiz[2];
+}
+
 #endif