}
//__________________________________________________________________________________________
-void AliMagWrapCheb::Field(const float *xyz, float *b) const
+void AliMagWrapCheb::Field(const Double_t *xyz, Double_t *b) const
{
// compute field in cartesian coordinates. If point is outside of the parameterized region
// get it at closest valid point
- static float rphiz[3];
+ static Double_t rphiz[3];
//
#ifndef _BRING_TO_BOUNDARY_ // exact matching to fitted volume is requested
if ( !(xyz[2]>=GetMinZSol()&&xyz[2]<=GetMaxZSol()) &&
}
//__________________________________________________________________________________________
-void AliMagWrapCheb::Field(const Double_t *xyz, Double_t *b) const
+Double_t AliMagWrapCheb::GetBz(const Double_t *xyz) const
{
- // compute field in cartesian coordinates. If point is outside of the parameterized region
+ // compute Bz for the point in cartesian coordinates. If point is outside of the parameterized region
// get it at closest valid point
static Double_t 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;}
+ !(xyz[2]>=GetMinZDip()&&xyz[2]<=GetMaxZDip()) ) return 0;
#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;
+ if (par->IsInside(xyz)) return par->Eval(xyz,2);
+ else return 0;
#else
- GetParamDip(FindDipSegment(xyz))->Eval(xyz,b); return;
+ return GetParamDip(FindDipSegment(xyz))->Eval(xyz,2);
#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
+ return FieldCylSolBz(rphiz);
+}
+
+
+//__________________________________________________________________________________________
+void AliMagWrapCheb::Print(Option_t *) const
+{
+ // print info
+ printf("Alice magnetic field parameterized by Chebyshev polynomials\n");
+ printf("Segmentation for Solenoid (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZSol,fMaxZSol,fMaxRSol);
//
- FieldCylSol(rphiz,b);
+ 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);
+ }
+ }
+ */
//
- // convert field to cartesian system
- CylToCartCylB(rphiz, b,b);
+ 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();
+ //
+
+
+}
+
+
+//__________________________________________________________________________________________________
+Int_t AliMagWrapCheb::FindDipSegment(const Double_t *xyz) const
+{
+ // find the segment containing point xyz. If it is outside find the closest segment
+ int xid,yid,zid = TMath::BinarySearch(fNZSegDip,fSegZDip,(Float_t)xyz[2]); // find zsegment
+ int ysegBeg = fBegSegYDip[zid];
+ //
+ for (yid=0;yid<fNSegYDip[zid];yid++) if (xyz[1]<fSegYDip[ysegBeg+yid]) break;
+ if ( --yid < 0 ) yid = 0;
+ yid += ysegBeg;
+ //
+ 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];
}
//__________________________________________________________________________________________
-void AliMagWrapCheb::GetTPCInt(const Float_t *xyz, Float_t *b) const
+void AliMagWrapCheb::GetTPCInt(const Double_t *xyz, Double_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];
+ static Double_t rphiz[3];
//
// TPCInt region
// convert coordinates to cyl system
}
//__________________________________________________________________________________________
-void AliMagWrapCheb::FieldCylSol(const float *rphiz, float *b) const
+void AliMagWrapCheb::FieldCylSol(const Double_t *rphiz, Double_t *b) const
{
// compute Solenoid field in Cylindircal coordinates
// note: if the point is outside the volume get the field in closest parameterized point
}
//__________________________________________________________________________________________
-void AliMagWrapCheb::FieldCylSol(const Double_t *rphiz, Double_t *b) const
+Double_t AliMagWrapCheb::FieldCylSolBz(const Double_t *rphiz) const
{
// compute Solenoid field in Cylindircal coordinates
// note: if the point is outside the volume get the field in closest parameterized point
int SolRId = fSegZIdSol[SolZId]; // first R segment for this Z
int SolRMax = SolRId + fNSegRSol[SolZId];
while (rphiz[0]>fSegRSol[SolRId] && SolRId<SolRMax-1) ++SolRId; // find R segment
- GetParamSol( SolRId )->Eval(rphiz,b);
+ return GetParamSol( SolRId )->Eval(rphiz,2);
//
}
//__________________________________________________________________________________________
-void AliMagWrapCheb::GetTPCIntCyl(const Float_t *rphiz, Float_t *b) const
+void AliMagWrapCheb::GetTPCIntCyl(const Double_t *rphiz, Double_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
}
-//__________________________________________________________________________________________
-void AliMagWrapCheb::Print(Option_t *) const
-{
- // print info
- 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();
- //
-
-
-}
#ifdef _INC_CREATION_ALICHEB3D_
//_______________________________________________
void AliMagWrapCheb::LoadData(const char* inpfile)
//
if (fNParamsTPCInt<1) return;
fNSegZTPCInt = 0;
- fSegRTPCInt = new Float_t[fNParamsTPCInt];
+ fSegRTPCInt = new Float_t[fNParamsTPCInt];
float *tmpbufF = new float[fNParamsTPCInt+1];
int *tmpbufI = new int[fNParamsTPCInt+1];
int *tmpbufI1 = new int[fNParamsTPCInt+1];
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff