From: shahoian Date: Wed, 17 Dec 2008 15:37:09 +0000 (+0000) Subject: Swapped the names AliMagFCheb and AliMagWrapCheb. The former should be used X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=da7cd221ccf4b0444be94b77343b843510063e88 Swapped the names AliMagFCheb and AliMagWrapCheb. The former should be used to construct the field. --- diff --git a/STEER/AliCheb3D.cxx b/STEER/AliCheb3D.cxx index 65eb5186f03..6d93b1f3053 100644 --- a/STEER/AliCheb3D.cxx +++ b/STEER/AliCheb3D.cxx @@ -357,7 +357,7 @@ void AliCheb3D::SetUsrFunction(void (*ptr)(float*,float*)) //__________________________________________________________________________________________ #ifdef _INC_CREATION_ALICHEB3D_ -void AliCheb3D::EvalUsrFunction(const Float_t *x, const Float_t *res) +void AliCheb3D::EvalUsrFunction(const Float_t *x, Float_t *res) { // evaluate user function value // diff --git a/STEER/AliCheb3D.h b/STEER/AliCheb3D.h index f147129a5a4..513995dfce3 100644 --- a/STEER/AliCheb3D.h +++ b/STEER/AliCheb3D.h @@ -131,7 +131,7 @@ class AliCheb3D: public TNamed // void SetUsrFunction(const char* name); void SetUsrFunction(void (*ptr)(float*,float*)); - void EvalUsrFunction(const Float_t *x, const Float_t *res); + void EvalUsrFunction(const Float_t *x, Float_t *res); TH1* TestRMS(int idim,int npoints = 1000,TH1* histo=0); static Int_t CalcChebCoefs(const Float_t *funval,int np, Float_t *outCoefs, Float_t prec=-1); #endif diff --git a/STEER/AliMagFCheb.cxx b/STEER/AliMagFCheb.cxx index 7a4e7263290..34a31f4b4cf 100644 --- a/STEER/AliMagFCheb.cxx +++ b/STEER/AliMagFCheb.cxx @@ -13,873 +13,148 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -#include "AliMagFCheb.h" -#include - -ClassImp(AliMagFCheb) -//__________________________________________________________________________________________ -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) -// -{ - // default constructor -} - -//__________________________________________________________________________________________ -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) -{ - // copy constructor - CopyFrom(src); -} - -//__________________________________________________________________________________________ -void AliMagFCheb::CopyFrom(const AliMagFCheb& src) -{ - // copy method - 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;iAddAtAndExpand(new AliCheb3D(*src.GetParamSol(i)),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;iAddAtAndExpand(new AliCheb3D(*src.GetParamDip(i)),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;iAddAtAndExpand(new AliCheb3D(*src.GetParamTPCInt(i)),i); - } - // -} - -//__________________________________________________________________________________________ -AliMagFCheb& AliMagFCheb::operator=(const AliMagFCheb& rhs) -{ - // assignment - if (this != &rhs) { - Clear(); - CopyFrom(rhs); - } - return *this; - // -} - -//__________________________________________________________________________________________ -void AliMagFCheb::Clear(const Option_t *) -{ - // clear all dynamic parts - if (fNParamsSol) { - delete fParamsSol; - delete[] fSegZSol; - delete[] fSegRSol; - delete[] fNSegRSol; - delete[] fSegZIdSol; - } - // - if (fNParamsTPCInt) { - delete fParamsTPCInt; - delete[] fSegZTPCInt; - delete[] fSegRTPCInt; - delete[] fNSegRTPCInt; - delete[] fSegZIdTPCInt; - } - // - if (fNParamsDip) { - 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::Field(const float *xyz, float *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]; - // -#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]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::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 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;} -#endif - // - if (xyz[2]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(const 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]GetMaxRTPCInt()) {for (int i=3;i--;) b[i]=0; return;} -#endif - // - GetTPCIntCyl(rphiz,b); - // - // convert field to cartesian system - CylToCartCylB(rphiz, b,b); - // -} +#include +#include +#include -//__________________________________________________________________________________________ -void AliMagFCheb::FieldCylSol(const float *rphiz, float *b) const -{ - // compute Solenoid field in Cylindircal coordinates - // note: if the point is outside the volume get the field in closest parameterized point - int SolZId = 0; - while (rphiz[2]>fSegZSol[SolZId] && SolZIdfSegRSol[SolRId] && SolRIdEval(rphiz,b); - // -} +#include "AliMagFCheb.h" +#include "AliMagWrapCheb.h" +#include "AliLog.h" -//__________________________________________________________________________________________ -void AliMagFCheb::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 - int SolZId = 0; - while (rphiz[2]>fSegZSol[SolZId] && SolZIdfSegRSol[SolRId] && SolRIdEval(rphiz,b); - // -} +ClassImp(AliMagFCheb) + -//__________________________________________________________________________________________ -void AliMagFCheb::GetTPCIntCyl(const Float_t *rphiz, Float_t *b) const +//_______________________________________________________________________ +AliMagFCheb::AliMagFCheb(): + AliMagFC(), + fMeasuredMap(0), + fSolenoid(5.) { - // compute field integral in TPC region in Cylindircal coordinates - // note: the check for the point being inside the parameterized region is done outside - int tpcIntZId = 0; - while (rphiz[2]>fSegZTPCInt[tpcIntZId] && tpcIntZIdfSegRTPCInt[tpcIntRId] && tpcIntRIdEval(rphiz,b); + // Default constructor // } - -//__________________________________________________________________________________________ -void AliMagFCheb::Print(Option_t *) const +//_______________________________________________________________________ +AliMagFCheb::AliMagFCheb(const char *name, const char *title, Int_t integ, + Float_t factor, Float_t fmax, Int_t map, + Bool_t dipoleON,const char* path): + AliMagFC(name, title, integ, factor, fmax), + fMeasuredMap(0), + fSolenoid(5.) { - // print info - printf("Alice magnetic field parameterized by Chebyshev polynomials\n"); - printf("Segmentation for Solenoid (%+.2fPrint(); - /* - for (int iz=0;izGetBoundMin(2),param->GetBoundMax(2)); - for (int ir=0;irGetBoundMin(0), - param->GetBoundMax(0),param->GetPrecision(),fSegZIdSol[iz]+ir); - } - } - */ - // - printf("Segmentation for TPC field integral (%+.2fPrint(); - /* - for (int iz=0;izGetBoundMin(2),param->GetBoundMax(2)); - for (int ir=0;irGetBoundMin(0), - param->GetBoundMax(0),param->GetPrecision(),fSegZIdTPCInt[iz]+ir); - } - } - */ // - printf("Segmentation for Dipole (%+.2fPrint(); - // - - -} -#ifdef _INC_CREATION_ALICHEB3D_ -//_______________________________________________ -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()); + fMap = map; + char* fname = gSystem->ExpandPathName(path); + TFile* file = TFile::Open(fname); + if (!file) { + AliError(Form("Failed to open magnetic field data file %s\n",fname)); return; } - // - TString buffs; - AliCheb3DCalc::ReadLine(buffs,stream); - if (!buffs.BeginsWith("START")) { - Error("LoadData","Expected: \"START \", 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;ipLoadData(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;ipLoadData(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;ipLoadData(stream); - AddParamDip(cheb); - } - // - AliCheb3DCalc::ReadLine(buffs,stream); - if (!buffs.BeginsWith("END DIPOLE")) { - Error("LoadData","Expected \"END DIPOLE\", found \"%s\"\nStop\n",GetName(),buffs.Data()); - exit(1); + const char* parname = 0; + if (fMap == k2kG) { + fSolenoid = 2.; + parname = dipoleON ? "Sol12_Dip6_Hole":"Sol12_Dip0_Hole"; + } else if (fMap == k5kG) { + fSolenoid = 5.; + parname = dipoleON ? "Sol30_Dip6_Hole":"Sol30_Dip0_Hole"; + } else { + AliError(Form("Unknown field identifier %d is requested\n",fMap)); + return; } // - 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); + fMeasuredMap = dynamic_cast(file->Get(parname)); + if (!fMeasuredMap) { + AliError(Form("Did not find field %s in %s\n",parname,fname)); + return; } - // - // --------------------------------------------------------------------------- - fclose(stream); - BuildTableSol(); - BuildTableTPCInt(); - BuildTableDip(); - printf("Loaded magnetic field \"%s\" from %s\n",GetName(),strf.Data()); - // + file->Close(); + delete file; } -#endif -//_______________________________________________ -#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) -// -// -{ - // construct from coeffs from the text file - LoadData(inputFile); -} -//__________________________________________________________________________________________ -void AliMagFCheb::AddParamSol(const AliCheb3D* param) +//_______________________________________________________________________ +AliMagFCheb::AliMagFCheb(const AliMagFCheb &src): + AliMagFC(src), + fMeasuredMap(0), + fSolenoid(src.fSolenoid) { - // adds new parameterization piece for Sol - // NOTE: pieces must be added strictly in increasing R then increasing Z order - // - if (!fParamsSol) fParamsSol = new TObjArray(); - fParamsSol->Add(param); - fNParamsSol++; - // + if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap); } -//__________________________________________________________________________________________ -void AliMagFCheb::AddParamTPCInt(const AliCheb3D* param) +//_______________________________________________________________________ +AliMagFCheb::~AliMagFCheb() { - // 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++; - // + delete fMeasuredMap; } -//__________________________________________________________________________________________ -void AliMagFCheb::AddParamDip(const AliCheb3D* param) +//_______________________________________________________________________ +void AliMagFCheb::GetTPCInt(const Float_t *xyz, Float_t *b) const { - // adds new parameterization piece for Dipole - // - if (!fParamsDip) fParamsDip = new TObjArray(); - fParamsDip->Add(param); - fNParamsDip++; + // Method to calculate the integral of magnetic integral from xyz to nearest cathode plane // + b[0]=b[1]=b[2]=0.0; + if (fMeasuredMap) fMeasuredMap->GetTPCInt(xyz,b); + for (int i=3;i--;) b[i] *= fFactor; } -//__________________________________________________________________________________________ -void AliMagFCheb::ResetTPCInt() +//_______________________________________________________________________ +void AliMagFCheb::GetTPCIntCyl(const Float_t *rphiz, Float_t *b) const { - // 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; - // + // Method to calculate the integral of magnetic integral from point to nearest cathode plane + // in cylindrical coordiates ( -piGetTPCIntCyl(rphiz,b); + for (int i=3;i--;) b[i] *= fFactor; } -//__________________________________________________ -void AliMagFCheb::BuildTableDip() +//_______________________________________________________________________ +void AliMagFCheb::Field(const Float_t *xyz, Float_t *b) const { - // build lookup table for dipole + // Method to calculate the field at point xyz // - 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;izAt(ipar); - if (!cheb->IsInside(xyz)) continue; - segID[fNXSegDip+ix] = ipar; - break; - } - } - fNXSegDip += nx; - // - delete[] tmpSegX; + b[0]=b[1]=b[2]=0.0; + if (xyz[2] > 919. || xyz[2] < -1972.) { + ZDCField(xyz, b); + } else { + if (fMeasuredMap && fFactor !=0.) { + fMeasuredMap->Field(xyz,b); + for (int i=3;i--;) b[i] *= fFactor; } - 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() +//_______________________________________________________________________ +void AliMagFCheb::Field(const Double_t *xyz, Double_t *b) const { - // build the indexes for each parameterization of Solenoid - // - const float kSafety=0.001; + // Method to calculate the field at point xyz // - if (fNParamsSol<1) return; - fNSegZSol = 0; - fMaxRSol = 0; - fSegRSol = new Float_t[fNParamsSol]; - float *tmpbufF = new float[fNParamsSol+1]; - int *tmpbufI = new int[fNParamsSol+1]; - int *tmpbufI1 = new int[fNParamsSol+1]; - // - // count number of Z slices and number of R slices in each Z slice - for (int ip=0;ipGetBoundMax(2)-GetParamSol(ip-1)->GetBoundMax(2))>kSafety) { // new Z slice - tmpbufF[fNSegZSol] = GetParamSol(ip)->GetBoundMax(2); // - tmpbufI[fNSegZSol] = 0; - tmpbufI1[fNSegZSol++] = ip; + b[0]=b[1]=b[2]=0.0; + if (xyz[2] > 919. || xyz[2] < -1972.) { + ZDCField(xyz, b); + } + else { + if (fMeasuredMap && fFactor !=0.) { + fMeasuredMap->Field(xyz,b); + for (int i=3;i--;) b[i] *= fFactor; } - fSegRSol[ip] = GetParamSol(ip)->GetBoundMax(0); // upper R - tmpbufI[fNSegZSol-1]++; - if (fMaxRSolGetBoundMin(2); - fMaxZSol = GetParamSol(fNParamsSol-1)->GetBoundMax(2); - // - delete[] tmpbufF; - delete[] tmpbufI; - delete[] tmpbufI1; - // - // } -//__________________________________________________________________________________________ -void AliMagFCheb::BuildTableTPCInt() +//_______________________________________________________________________ +AliMagFCheb& AliMagFCheb::operator=(const AliMagFCheb& maps) { - // build the indexes for each parameterization of TPC field integral - // - const float kSafety=0.001; - // - if (fNParamsTPCInt<1) return; - fNSegZTPCInt = 0; - fSegRTPCInt = new Float_t[fNParamsTPCInt]; - float *tmpbufF = new float[fNParamsTPCInt+1]; - int *tmpbufI = new int[fNParamsTPCInt+1]; - int *tmpbufI1 = new int[fNParamsTPCInt+1]; - // - // count number of Z slices and number of R slices in each Z slice - for (int ip=0;ipGetBoundMax(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]++; - } - // - fSegZTPCInt = new Float_t[fNSegZTPCInt]; - fSegZIdTPCInt = new Int_t[fNSegZTPCInt]; - fNSegRTPCInt = new Int_t[fNSegZTPCInt]; - for (int iz=0;izGetBoundMin(2); - fMaxZTPCInt = GetParamTPCInt(fNParamsTPCInt-1)->GetBoundMax(2); - fMaxRTPCInt = GetParamTPCInt(fNParamsTPCInt-1)->GetBoundMax(0); - // - delete[] tmpbufF; - delete[] tmpbufI; - delete[] tmpbufI1; - // - // -} - -void AliMagFCheb::SaveData(const char* outfile) const -{ - // writes coefficients data to output text file - TString strf = outfile; - gSystem->ExpandPathName(strf); - FILE* stream = fopen(strf,"w+"); - // - // 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;ipSaveData(stream); - fprintf(stream,"#\nEND SOLENOID\n"); - // - // 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;ipSaveData(stream); - fprintf(stream,"#\nEND TPCINT\n"); - // - // Dip part --------------------------------------------------------- - fprintf(stream,"START DIPOLE\n#Number of pieces\n%d\n",fNParamsDip); - for (int ip=0;ipSaveData(stream); - fprintf(stream,"#\nEND DIPOLE\n"); - // - fprintf(stream,"#\nEND %s\n",GetName()); - // - fclose(stream); - // + return *this; } -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) +//_______________________________________________________________________ +void AliMagFCheb::SetMeasuredMap(AliMagWrapCheb* parm) { - // 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;ipAt(ip); - if (xmnGetBoundMin(0)>(xmx+xmn)/2 || cheb->GetBoundMax(0)<(xmn+xmx)/2)) continue; - if (ymnGetBoundMin(1)>(ymx+ymn)/2 || cheb->GetBoundMax(1)<(ymn+ymx)/2)) continue; - if (zmnGetBoundMin(2)>(zmx+zmn)/2 || cheb->GetBoundMax(2)<(zmn+zmx)/2)) continue; - // - tmpC[nseg0++] = cheb->GetBoundMin(dim); - tmpC[nseg0++] = cheb->GetBoundMax(dim); - } - // 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;ip1e-4) { - cprev = tmpC[ tmpInd[ip] ]; - nseg++; - } - else tmpInd[ip] = -1; // supress redundant Z - } - // - *seg = new float[nseg]; // create final Z segmenations - nseg = 0; - for (int ip=0;ip=0) (*seg)[nseg++] = tmpC[ tmpInd[ip] ]; - // - delete[] tmpC; - delete[] tmpInd; - return nseg; + if (fMeasuredMap) delete fMeasuredMap; + fMeasuredMap = parm; } - -#endif - diff --git a/STEER/AliMagFCheb.h b/STEER/AliMagFCheb.h index 724f0e5585a..74e706fed89 100644 --- a/STEER/AliMagFCheb.h +++ b/STEER/AliMagFCheb.h @@ -1,296 +1,47 @@ - -// Author: ruben.shahoyan@cern.ch 20/03/2007 - -/////////////////////////////////////////////////////////////////////////////////// -// // -// Wrapper for the set of mag.field parameterizations by Chebyshev polinomials // -// To obtain the field in cartesian coordinates/components use // -// Field(float* xyz, float* bxyz); // -// For cylindrical coordinates/components: // -// FieldCyl(float* rphiz, float* brphiz) // -// // -// The solenoid part is parameterized in the volume R<500, -550 -#include -#include "AliCheb3D.h" +// +// Interface between the AliMagWrapCheb and AliMagF: set of magnetic field data + Tosca +// parameterization by Chebyshev polynomials +// +// Author: ruben.shahoyan@cern.ch +// -class TSystem; +#include "AliMagFC.h" +class AliMagWrapCheb; -class AliMagFCheb: public TNamed + +class AliMagFCheb : public AliMagFC { - public: +public: + enum constants {k2kG, k4kG, k5kG}; AliMagFCheb(); - AliMagFCheb(const AliMagFCheb& src); - ~AliMagFCheb() {Clear();} - // - void CopyFrom(const AliMagFCheb& src); - AliMagFCheb& operator=(const AliMagFCheb& rhs); - virtual void Clear(const Option_t * = ""); - // - Int_t GetNParamsSol() const {return fNParamsSol;} - Int_t GetNSegZSol() const {return fNSegZSol;} - Float_t* GetSegZSol() const {return fSegZSol;} - // - 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;} - // - 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);} - // - virtual void Print(Option_t * = "") const; - // - virtual void Field(const Float_t *xyz, Float_t *b) const; - virtual void Field(const Double_t *xyz, Double_t *b) const; - // - virtual void FieldCyl(const Float_t *rphiz, Float_t *b) const; - virtual void FieldCyl(const Double_t *rphiz, Double_t *b) const; - // + AliMagFCheb(const char *name, const char *title, Int_t integ, + Float_t factor=1, Float_t fmax=15, Int_t map = k2kG, + Bool_t dipoleON = kTRUE, + const char* path="$(ALICE_ROOT)/data/maps/mfchebKGI_sym.root"); + AliMagFCheb(const AliMagFCheb& maps); + AliMagFCheb& operator=(const AliMagFCheb& maps); + virtual ~AliMagFCheb(); + // + virtual void Field(const Float_t *x, Float_t *b) const; + virtual void Field(const Double_t *x, Double_t *b) const; virtual void GetTPCInt(const Float_t *xyz, Float_t *b) const; virtual void GetTPCIntCyl(const Float_t *rphiz, Float_t *b) const; // - template - Int_t FindDipSegment(const T *xyz) const; - // - template - static void CylToCartCylB(const T *rphiz, const T *brphiz,T *bxyz); - template - static void CylToCartCartB(const T *xyz, const T *brphiz,T *bxyz); - template - static void CartToCylCartB(const T *xyz, const T *bxyz, T *brphiz); - template - static void CartToCylCylB(const T *rphiz, const T *bxyz, T *brphiz); - template - static void CartToCyl(const T *xyz, T *rphiz); - template - static void CylToCart(const T *rphiz,T *xyz); - // -#ifdef _INC_CREATION_ALICHEB3D_ // see AliCheb3D.h for explanation - void LoadData(const char* inpfile); - // - AliMagFCheb(const char* inputFile); - void SaveData(const char* outfile) const; - Int_t SegmentDipDimension(Float_t** seg,const TObjArray* par,int npar, int dim, - Float_t xmn,Float_t xmx,Float_t ymn,Float_t ymx,Float_t zmn,Float_t zmx); - // - void AddParamSol(const AliCheb3D* param); - void AddParamTPCInt(const AliCheb3D* param); - void AddParamDip(const AliCheb3D* param); - void BuildTableDip(); - void BuildTableSol(); - void BuildTableTPCInt(); - void ResetTPCInt(); - // - // -#endif - // - protected: - virtual void FieldCylSol(const Float_t *rphiz, Float_t *b) const; - virtual void FieldCylSol(const Double_t *rphiz, Double_t *b) const; + AliMagWrapCheb* GetMeasuredMap() const {return fMeasuredMap;} + void SetMeasuredMap(AliMagWrapCheb* parm); + virtual Float_t SolenoidField() const {return -Factor()*fSolenoid;} // protected: - // - Int_t fNParamsSol; // Total number of parameterization pieces for Sol - 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* 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* 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,3) // Wrapper class for the set of Chebishev parameterizations of Alice mag.field - // - }; - - -//__________________________________________________________________________________________ -inline void AliMagFCheb::FieldCyl(const Float_t *rphiz, Float_t *b) const -{ - // compute field in Cylindircal coordinates - // if (rphiz[2]GetMaxZSol() || rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;} - FieldCylSol(rphiz,b); -} + AliMagWrapCheb* fMeasuredMap; // Measured part of the field map + Float_t fSolenoid; // Solenoid field setting + // + ClassDef(AliMagFCheb, 2) // Class for all Alice MagField wrapper for measured data + Tosca parameterization +}; -//__________________________________________________________________________________________ -inline void AliMagFCheb::FieldCyl(const Double_t *rphiz, Double_t *b) const -{ - // compute field in Cylindircal coordinates - // if (rphiz[2]GetMaxZSol() || rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;} - FieldCylSol(rphiz,b); -} - -//__________________________________________________________________________________________________ -template -inline void AliMagFCheb::CylToCartCylB(const T *rphiz, const T *brphiz,T *bxyz) -{ - // convert field in cylindrical coordinates to cartesian system, point is in cyl.system - T btr = TMath::Sqrt(brphiz[0]*brphiz[0]+brphiz[1]*brphiz[1]); - T 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]; - // -} - -//__________________________________________________________________________________________________ -template -inline void AliMagFCheb::CylToCartCartB(const T *xyz, const T *brphiz, T *bxyz) -{ - // convert field in cylindrical coordinates to cartesian system, point is in cart.system - T btr = TMath::Sqrt(brphiz[0]*brphiz[0]+brphiz[1]*brphiz[1]); - T 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]; - // -} - -//__________________________________________________________________________________________________ -template -inline void AliMagFCheb::CartToCylCartB(const T *xyz, const T *bxyz, T *brphiz) -{ - // convert field in cylindrical coordinates to cartesian system, poin is in cart.system - T btr = TMath::Sqrt(bxyz[0]*bxyz[0]+bxyz[1]*bxyz[1]); - T 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]; - // -} - -//__________________________________________________________________________________________________ -template -inline void AliMagFCheb::CartToCylCylB(const T *rphiz, const T *bxyz, T *brphiz) -{ - // convert field in cylindrical coordinates to cartesian system, point is in cyl.system - T btr = TMath::Sqrt(bxyz[0]*bxyz[0]+bxyz[1]*bxyz[1]); - T 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]; - // -} - -//__________________________________________________________________________________________________ -template -inline void AliMagFCheb::CartToCyl(const T *xyz,T *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]; -} - -//__________________________________________________________________________________________________ -template -inline void AliMagFCheb::CylToCart(const T *rphiz, T *xyz) -{ - xyz[0] = rphiz[0]*TMath::Cos(rphiz[1]); - xyz[1] = rphiz[0]*TMath::Sin(rphiz[1]); - xyz[2] = rphiz[2]; -} - -//__________________________________________________________________________________________________ -template -Int_t AliMagFCheb::FindDipSegment(const 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 -#include -#include - #include "AliMagWrapCheb.h" -#include "AliLog.h" +#include +#include +#include ClassImp(AliMagWrapCheb) - -//_______________________________________________________________________ -AliMagWrapCheb::AliMagWrapCheb(): - AliMagFC(), - fMeasuredMap(0), - fSolenoid(5.) +//__________________________________________________________________________________________ +AliMagWrapCheb::AliMagWrapCheb() : + 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) +// { - // Default constructor - // + // default constructor } -//_______________________________________________________________________ -AliMagWrapCheb::AliMagWrapCheb(const char *name, const char *title, Int_t integ, - Float_t factor, Float_t fmax, Int_t map, - Bool_t dipoleON,const char* path): - AliMagFC(name, title, integ, factor, fmax), - fMeasuredMap(0), - fSolenoid(5.) +//__________________________________________________________________________________________ +AliMagWrapCheb::AliMagWrapCheb(const AliMagWrapCheb& 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) { + // copy constructor + CopyFrom(src); +} + +//__________________________________________________________________________________________ +void AliMagWrapCheb::CopyFrom(const AliMagWrapCheb& src) +{ + // copy method + Clear(); + SetName(src.GetName()); + SetTitle(src.GetTitle()); + fNParamsSol = src.fNParamsSol; + fNSegZSol = src.fNSegZSol; + fNParamsTPCInt = src.fNParamsTPCInt; + fNSegZTPCInt = src.fNSegZTPCInt; + fNParamsDip = src.fNParamsDip; // - fMap = map; - char* fname = gSystem->ExpandPathName(path); - TFile* file = TFile::Open(fname); - if (!file) { - AliError(Form("Failed to open magnetic field data file %s\n",fname)); - return; + 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;iAddAtAndExpand(new AliCheb3D(*src.GetParamSol(i)),i); } - const char* parname = 0; - if (fMap == k2kG) { - fSolenoid = 2.; - parname = dipoleON ? "Sol12_Dip6_Hole":"Sol12_Dip0_Hole"; - } else if (fMap == k5kG) { - fSolenoid = 5.; - parname = dipoleON ? "Sol30_Dip6_Hole":"Sol30_Dip0_Hole"; - } else { - AliError(Form("Unknown field identifier %d is requested\n",fMap)); - return; + // + 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;iAddAtAndExpand(new AliCheb3D(*src.GetParamDip(i)),i); } // - fMeasuredMap = dynamic_cast(file->Get(parname)); - if (!fMeasuredMap) { - AliError(Form("Did not find field %s in %s\n",parname,fname)); - return; + 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;iAddAtAndExpand(new AliCheb3D(*src.GetParamTPCInt(i)),i); + } + // +} + +//__________________________________________________________________________________________ +AliMagWrapCheb& AliMagWrapCheb::operator=(const AliMagWrapCheb& rhs) +{ + // assignment + if (this != &rhs) { + Clear(); + CopyFrom(rhs); } - file->Close(); - delete file; + return *this; + // } +//__________________________________________________________________________________________ +void AliMagWrapCheb::Clear(const Option_t *) +{ + // clear all dynamic parts + if (fNParamsSol) { + delete fParamsSol; + delete[] fSegZSol; + delete[] fSegRSol; + delete[] fNSegRSol; + delete[] fSegZIdSol; + } + // + if (fNParamsTPCInt) { + delete fParamsTPCInt; + delete[] fSegZTPCInt; + delete[] fSegRTPCInt; + delete[] fNSegRTPCInt; + delete[] fSegZIdTPCInt; + } + // + if (fNParamsDip) { + 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; + // +} -//_______________________________________________________________________ -AliMagWrapCheb::AliMagWrapCheb(const AliMagWrapCheb &src): - AliMagFC(src), - fMeasuredMap(0), - fSolenoid(src.fSolenoid) +//__________________________________________________________________________________________ +void AliMagWrapCheb::Field(const float *xyz, float *b) const { - if (src.fMeasuredMap) fMeasuredMap = new AliMagFCheb(*src.fMeasuredMap); + // 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]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); + // } -//_______________________________________________________________________ -AliMagWrapCheb::~AliMagWrapCheb() +//__________________________________________________________________________________________ +void AliMagWrapCheb::Field(const Double_t *xyz, Double_t *b) const { - delete fMeasuredMap; + // compute field 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;} +#endif + // + if (xyz[2]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 AliMagWrapCheb::GetTPCInt(const Float_t *xyz, Float_t *b) const { - // Method to calculate the integral of magnetic integral from xyz to nearest cathode plane + // 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]GetMaxRTPCInt()) {for (int i=3;i--;) b[i]=0; return;} +#endif + // + GetTPCIntCyl(rphiz,b); + // + // convert field to cartesian system + CylToCartCylB(rphiz, b,b); + // +} + +//__________________________________________________________________________________________ +void AliMagWrapCheb::FieldCylSol(const float *rphiz, float *b) const +{ + // compute Solenoid field in Cylindircal coordinates + // note: if the point is outside the volume get the field in closest parameterized point + int SolZId = 0; + while (rphiz[2]>fSegZSol[SolZId] && SolZIdfSegRSol[SolRId] && SolRIdEval(rphiz,b); + // +} + +//__________________________________________________________________________________________ +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 + int SolZId = 0; + while (rphiz[2]>fSegZSol[SolZId] && SolZIdfSegRSol[SolRId] && SolRIdEval(rphiz,b); // - b[0]=b[1]=b[2]=0.0; - if (fMeasuredMap) fMeasuredMap->GetTPCInt(xyz,b); - for (int i=3;i--;) b[i] *= fFactor; } -//_______________________________________________________________________ +//__________________________________________________________________________________________ void AliMagWrapCheb::GetTPCIntCyl(const Float_t *rphiz, Float_t *b) const { - // Method to calculate the integral of magnetic integral from point to nearest cathode plane - // in cylindrical coordiates ( -piGetTPCIntCyl(rphiz,b); - for (int i=3;i--;) b[i] *= fFactor; + // compute field integral in TPC region in Cylindircal coordinates + // note: the check for the point being inside the parameterized region is done outside + int tpcIntZId = 0; + while (rphiz[2]>fSegZTPCInt[tpcIntZId] && tpcIntZIdfSegRTPCInt[tpcIntRId] && tpcIntRIdEval(rphiz,b); + // } -//_______________________________________________________________________ -void AliMagWrapCheb::Field(const Float_t *xyz, Float_t *b) const + +//__________________________________________________________________________________________ +void AliMagWrapCheb::Print(Option_t *) const { - // Method to calculate the field at point xyz + // print info + printf("Alice magnetic field parameterized by Chebyshev polynomials\n"); + printf("Segmentation for Solenoid (%+.2f 919. || xyz[2] < -1972.) { - ZDCField(xyz, b); - } else { - if (fMeasuredMap && fFactor !=0.) { - fMeasuredMap->Field(xyz,b); - for (int i=3;i--;) b[i] *= fFactor; - } + if (fParamsSol) fParamsSol->Print(); + /* + for (int iz=0;izGetBoundMin(2),param->GetBoundMax(2)); + for (int ir=0;irGetBoundMin(0), + param->GetBoundMax(0),param->GetPrecision(),fSegZIdSol[iz]+ir); + } + } + */ + // + printf("Segmentation for TPC field integral (%+.2fPrint(); + /* + for (int iz=0;izGetBoundMin(2),param->GetBoundMax(2)); + for (int ir=0;irGetBoundMin(0), + param->GetBoundMax(0),param->GetPrecision(),fSegZIdTPCInt[iz]+ir); } + } + */ + // + printf("Segmentation for Dipole (%+.2fPrint(); + // + + } +#ifdef _INC_CREATION_ALICHEB3D_ +//_______________________________________________ +void AliMagWrapCheb::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 \", 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;ipLoadData(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;ipLoadData(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;ipLoadData(stream); + AddParamDip(cheb); + } + // + 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()); + // +} +#endif -//_______________________________________________________________________ -void AliMagWrapCheb::Field(const Double_t *xyz, Double_t *b) const +//_______________________________________________ +#ifdef _INC_CREATION_ALICHEB3D_ + +//__________________________________________________________________________________________ +AliMagWrapCheb::AliMagWrapCheb(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) +// +// +{ + // construct from coeffs from the text file + LoadData(inputFile); +} + +//__________________________________________________________________________________________ +void AliMagWrapCheb::AddParamSol(const AliCheb3D* param) +{ + // adds new parameterization piece for Sol + // NOTE: pieces must be added strictly in increasing R then increasing Z order + // + if (!fParamsSol) fParamsSol = new TObjArray(); + fParamsSol->Add( (AliCheb3D*)param ); + fNParamsSol++; + // +} + +//__________________________________________________________________________________________ +void AliMagWrapCheb::AddParamTPCInt(const 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( (AliCheb3D*)param); + fNParamsTPCInt++; + // +} + +//__________________________________________________________________________________________ +void AliMagWrapCheb::AddParamDip(const AliCheb3D* param) { - // Method to calculate the field at point xyz + // adds new parameterization piece for Dipole // - b[0]=b[1]=b[2]=0.0; - if (xyz[2] > 919. || xyz[2] < -1972.) { - ZDCField(xyz, b); - } else { - if (fMeasuredMap && fFactor !=0.) { - fMeasuredMap->Field(xyz,b); - for (int i=3;i--;) b[i] *= fFactor; + if (!fParamsDip) fParamsDip = new TObjArray(); + fParamsDip->Add( (AliCheb3D*)param); + fNParamsDip++; + // +} + +//__________________________________________________________________________________________ +void AliMagWrapCheb::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 AliMagWrapCheb::BuildTableDip() +{ + // build lookup table for dipole + // + 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;izAt(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];} + // } -//_______________________________________________________________________ -AliMagWrapCheb& AliMagWrapCheb::operator=(const AliMagWrapCheb& maps) +//__________________________________________________________________________________________ +void AliMagWrapCheb::BuildTableSol() { - fSolenoid=maps.fSolenoid; - if (this != &maps && maps.fMeasuredMap) { - if (fMeasuredMap) delete fMeasuredMap; - fMeasuredMap = new AliMagFCheb(*maps.fMeasuredMap); + // build the indexes for each parameterization of Solenoid + // + const float kSafety=0.001; + // + if (fNParamsSol<1) return; + fNSegZSol = 0; + fMaxRSol = 0; + fSegRSol = new Float_t[fNParamsSol]; + float *tmpbufF = new float[fNParamsSol+1]; + int *tmpbufI = new int[fNParamsSol+1]; + int *tmpbufI1 = new int[fNParamsSol+1]; + // + // count number of Z slices and number of R slices in each Z slice + for (int ip=0;ipGetBoundMax(2)-GetParamSol(ip-1)->GetBoundMax(2))>kSafety) { // new Z slice + tmpbufF[fNSegZSol] = GetParamSol(ip)->GetBoundMax(2); // + tmpbufI[fNSegZSol] = 0; + tmpbufI1[fNSegZSol++] = ip; + } + fSegRSol[ip] = GetParamSol(ip)->GetBoundMax(0); // upper R + tmpbufI[fNSegZSol-1]++; + if (fMaxRSolGetBoundMin(2); + fMaxZSol = GetParamSol(fNParamsSol-1)->GetBoundMax(2); + // + delete[] tmpbufF; + delete[] tmpbufI; + delete[] tmpbufI1; + // + // } + +//__________________________________________________________________________________________ +void AliMagWrapCheb::BuildTableTPCInt() +{ + // build the indexes for each parameterization of TPC field integral + // + const float kSafety=0.001; + // + if (fNParamsTPCInt<1) return; + fNSegZTPCInt = 0; + fSegRTPCInt = new Float_t[fNParamsTPCInt]; + float *tmpbufF = new float[fNParamsTPCInt+1]; + int *tmpbufI = new int[fNParamsTPCInt+1]; + int *tmpbufI1 = new int[fNParamsTPCInt+1]; + // + // count number of Z slices and number of R slices in each Z slice + for (int ip=0;ipGetBoundMax(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]++; + } + // + fSegZTPCInt = new Float_t[fNSegZTPCInt]; + fSegZIdTPCInt = new Int_t[fNSegZTPCInt]; + fNSegRTPCInt = new Int_t[fNSegZTPCInt]; + for (int iz=0;izGetBoundMin(2); + fMaxZTPCInt = GetParamTPCInt(fNParamsTPCInt-1)->GetBoundMax(2); + fMaxRTPCInt = GetParamTPCInt(fNParamsTPCInt-1)->GetBoundMax(0); + // + delete[] tmpbufF; + delete[] tmpbufI; + delete[] tmpbufI1; + // + // +} + +void AliMagWrapCheb::SaveData(const char* outfile) const +{ + // writes coefficients data to output text file + TString strf = outfile; + gSystem->ExpandPathName(strf); + FILE* stream = fopen(strf,"w+"); + // + // 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;ipSaveData(stream); + fprintf(stream,"#\nEND SOLENOID\n"); + // + // 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;ipSaveData(stream); + fprintf(stream,"#\nEND TPCINT\n"); + // + // Dip part --------------------------------------------------------- + fprintf(stream,"START DIPOLE\n#Number of pieces\n%d\n",fNParamsDip); + for (int ip=0;ipSaveData(stream); + fprintf(stream,"#\nEND DIPOLE\n"); + // + fprintf(stream,"#\nEND %s\n",GetName()); + // + fclose(stream); + // +} + +Int_t AliMagWrapCheb::SegmentDipDimension(float** seg,const TObjArray* par,int npar, int dim, + float xmn,float xmx,float ymn,float ymx,float zmn,float zmx) +{ + // 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;ipAt(ip); + if (xmnGetBoundMin(0)>(xmx+xmn)/2 || cheb->GetBoundMax(0)<(xmn+xmx)/2)) continue; + if (ymnGetBoundMin(1)>(ymx+ymn)/2 || cheb->GetBoundMax(1)<(ymn+ymx)/2)) continue; + if (zmnGetBoundMin(2)>(zmx+zmn)/2 || cheb->GetBoundMax(2)<(zmn+zmx)/2)) continue; + // + tmpC[nseg0++] = cheb->GetBoundMin(dim); + tmpC[nseg0++] = cheb->GetBoundMax(dim); + } + // 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;ip1e-4) { + cprev = tmpC[ tmpInd[ip] ]; + nseg++; + } + else tmpInd[ip] = -1; // supress redundant Z + } + // + *seg = new float[nseg]; // create final Z segmenations + nseg = 0; + for (int ip=0;ip=0) (*seg)[nseg++] = tmpC[ tmpInd[ip] ]; + // + delete[] tmpC; + delete[] tmpInd; + return nseg; +} + +#endif + diff --git a/STEER/AliMagWrapCheb.h b/STEER/AliMagWrapCheb.h index 189be0e735b..195dc50d71b 100644 --- a/STEER/AliMagWrapCheb.h +++ b/STEER/AliMagWrapCheb.h @@ -1,46 +1,298 @@ + +// Author: ruben.shahoyan@cern.ch 20/03/2007 + +/////////////////////////////////////////////////////////////////////////////////// +// // +// Wrapper for the set of mag.field parameterizations by Chebyshev polinomials // +// To obtain the field in cartesian coordinates/components use // +// Field(float* xyz, float* bxyz); // +// For cylindrical coordinates/components: // +// FieldCyl(float* rphiz, float* brphiz) // +// // +// The solenoid part is parameterized in the volume R<500, -550 +#include +#include "AliCheb3D.h" -#include "AliMagFC.h" -#include "AliMagFCheb.h" +class TSystem; +class TArrayF; +class TArrayI; - -class AliMagWrapCheb : public AliMagFC +class AliMagWrapCheb: public TNamed { -public: - enum constants {k2kG, k4kG, k5kG}; + public: AliMagWrapCheb(); - AliMagWrapCheb(const char *name, const char *title, Int_t integ, - Float_t factor=1, Float_t fmax=15, Int_t map = k2kG, - Bool_t dipoleON = kTRUE, - const char* path="$(ALICE_ROOT)/data/maps/mfchebKGI_sym.root"); - AliMagWrapCheb(const AliMagWrapCheb& maps); - AliMagWrapCheb& operator=(const AliMagWrapCheb& maps); - virtual ~AliMagWrapCheb(); - // - virtual void Field(const Float_t *x, Float_t *b) const; - virtual void Field(const Double_t *x, Double_t *b) const; - virtual void GetTPCInt(const Float_t *xyz, Float_t *b) const; - virtual void GetTPCIntCyl(const Float_t *rphiz, Float_t *b) const; - // - AliMagFCheb* GetMeasuredMap() const {return fMeasuredMap;} - void SetMeasuredMap(AliMagFCheb* parm) {if (fMeasuredMap) delete fMeasuredMap; fMeasuredMap = parm;} - virtual Float_t SolenoidField() const {return -Factor()*fSolenoid;} + AliMagWrapCheb(const AliMagWrapCheb& src); + ~AliMagWrapCheb() {Clear();} + // + void CopyFrom(const AliMagWrapCheb& src); + AliMagWrapCheb& operator=(const AliMagWrapCheb& rhs); + virtual void Clear(const Option_t * = ""); + // + Int_t GetNParamsSol() const {return fNParamsSol;} + Int_t GetNSegZSol() const {return fNSegZSol;} + Float_t* GetSegZSol() const {return fSegZSol;} + // + 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;} + // + 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);} + // + virtual void Print(Option_t * = "") const; + // + virtual void Field(const Float_t *xyz, Float_t *b) const; + virtual void Field(const Double_t *xyz, Double_t *b) const; + // + virtual void FieldCyl(const Float_t *rphiz, Float_t *b) const; + virtual void FieldCyl(const Double_t *rphiz, Double_t *b) const; + // + virtual void GetTPCInt(const Float_t *xyz, Float_t *b) const; + virtual void GetTPCIntCyl(const Float_t *rphiz, Float_t *b) const; + // + template + Int_t FindDipSegment(const T *xyz) const; + // + template + static void CylToCartCylB(const T *rphiz, const T *brphiz,T *bxyz); + template + static void CylToCartCartB(const T *xyz, const T *brphiz,T *bxyz); + template + static void CartToCylCartB(const T *xyz, const T *bxyz, T *brphiz); + template + static void CartToCylCylB(const T *rphiz, const T *bxyz, T *brphiz); + template + static void CartToCyl(const T *xyz, T *rphiz); + template + static void CylToCart(const T *rphiz,T *xyz); + // +#ifdef _INC_CREATION_ALICHEB3D_ // see AliCheb3D.h for explanation + void LoadData(const char* inpfile); + // + AliMagWrapCheb(const char* inputFile); + void SaveData(const char* outfile) const; + Int_t SegmentDipDimension(Float_t** seg,const TObjArray* par,int npar, int dim, + Float_t xmn,Float_t xmx,Float_t ymn,Float_t ymx,Float_t zmn,Float_t zmx); + // + void AddParamSol(const AliCheb3D* param); + void AddParamTPCInt(const AliCheb3D* param); + void AddParamDip(const AliCheb3D* param); + void BuildTableDip(); + void BuildTableSol(); + void BuildTableTPCInt(); + void ResetTPCInt(); + // + // +#endif // protected: - AliMagFCheb* fMeasuredMap; // Measured part of the field map - Float_t fSolenoid; // Solenoid field setting - // - ClassDef(AliMagWrapCheb, 2) // Class for all Alice MagField wrapper for measured data + Tosca parameterization -}; + virtual void FieldCylSol(const Float_t *rphiz, Float_t *b) const; + virtual void FieldCylSol(const Double_t *rphiz, Double_t *b) const; + // + protected: + // + Int_t fNParamsSol; // Total number of parameterization pieces for Sol + 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* 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* 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(AliMagWrapCheb,3) // Wrapper class for the set of Chebishev parameterizations of Alice mag.field + // + }; + + +//__________________________________________________________________________________________ +inline void AliMagWrapCheb::FieldCyl(const Float_t *rphiz, Float_t *b) const +{ + // compute field in Cylindircal coordinates + // if (rphiz[2]GetMaxZSol() || rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;} + FieldCylSol(rphiz,b); +} + + +//__________________________________________________________________________________________ +inline void AliMagWrapCheb::FieldCyl(const Double_t *rphiz, Double_t *b) const +{ + // compute field in Cylindircal coordinates + // if (rphiz[2]GetMaxZSol() || rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;} + FieldCylSol(rphiz,b); +} + +//__________________________________________________________________________________________________ +template +inline void AliMagWrapCheb::CylToCartCylB(const T *rphiz, const T *brphiz,T *bxyz) +{ + // convert field in cylindrical coordinates to cartesian system, point is in cyl.system + T btr = TMath::Sqrt(brphiz[0]*brphiz[0]+brphiz[1]*brphiz[1]); + T 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]; + // +} + +//__________________________________________________________________________________________________ +template +inline void AliMagWrapCheb::CylToCartCartB(const T *xyz, const T *brphiz, T *bxyz) +{ + // convert field in cylindrical coordinates to cartesian system, point is in cart.system + T btr = TMath::Sqrt(brphiz[0]*brphiz[0]+brphiz[1]*brphiz[1]); + T 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]; + // +} + +//__________________________________________________________________________________________________ +template +inline void AliMagWrapCheb::CartToCylCartB(const T *xyz, const T *bxyz, T *brphiz) +{ + // convert field in cylindrical coordinates to cartesian system, poin is in cart.system + T btr = TMath::Sqrt(bxyz[0]*bxyz[0]+bxyz[1]*bxyz[1]); + T 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]; + // +} + +//__________________________________________________________________________________________________ +template +inline void AliMagWrapCheb::CartToCylCylB(const T *rphiz, const T *bxyz, T *brphiz) +{ + // convert field in cylindrical coordinates to cartesian system, point is in cyl.system + T btr = TMath::Sqrt(bxyz[0]*bxyz[0]+bxyz[1]*bxyz[1]); + T 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]; + // +} +//__________________________________________________________________________________________________ +template +inline void AliMagWrapCheb::CartToCyl(const T *xyz,T *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]; +} + +//__________________________________________________________________________________________________ +template +inline void AliMagWrapCheb::CylToCart(const T *rphiz, T *xyz) +{ + xyz[0] = rphiz[0]*TMath::Cos(rphiz[1]); + xyz[1] = rphiz[0]*TMath::Sin(rphiz[1]); + xyz[2] = rphiz[2]; +} + +//__________________________________________________________________________________________________ +template +Int_t AliMagWrapCheb::FindDipSegment(const 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