From: morsch Date: Wed, 25 Apr 2007 16:31:57 +0000 (+0000) Subject: Parameterisation of the measured mag. field map. (R. Shahoyan) X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=0eea9d4dec177e0534b23e1e83a39c40a622ff17;p=u%2Fmrichter%2FAliRoot.git Parameterisation of the measured mag. field map. (R. Shahoyan) --- diff --git a/STEER/AliCheb3D.cxx b/STEER/AliCheb3D.cxx new file mode 100644 index 00000000000..0eaa7b335f7 --- /dev/null +++ b/STEER/AliCheb3D.cxx @@ -0,0 +1,895 @@ +/************************************************************************** + * 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$ */ + +// Author: ruben.shahoyan@cern.ch 09/09/2006 +// +//////////////////////////////////////////////////////////////////////////////// +// // +// AliCheb3D produces the interpolation of the user 3D->NDimOut arbitrary // +// function supplied in "void (*fcn)(float* inp,float* out)" format // +// either in a separate macro file or as a function pointer. // +// Only coefficients needed to guarantee the requested precision are kept. // +// // +// The user-callable methods are: // +// To create the interpolation use: // +// AliCheb3D(const char* funName, // name of the file with user function // +// or // +// AliCheb3D(void (*ptr)(float*,float*),// pointer on the user function // +// Int_t DimOut, // dimensionality of the function's output // +// Float_t *bmin, // lower 3D bounds of interpolation domain // +// Float_t *bmax, // upper 3D bounds of interpolation domain // +// Int_t *npoints, // number of points in each of 3 input // +// // dimension, defining the interpolation grid // +// Float_t prec=1E-6); // requested max.absolute difference between // +// // the interpolation and any point on grid // +// // +// To test obtained parameterization use the method // +// TH1* TestRMS(int idim,int npoints = 1000,TH1* histo=0); // +// it will compare the user output of the user function and interpolation // +// for idim-th output dimension and fill the difference in the supplied // +// histogram. If no histogram is supplied, it will be created. // +// // +// To save the interpolation data: // +// SaveData(const char* filename, Bool_t append ) // +// write text file with data. If append is kTRUE and the output file already // +// exists, data will be added in the end of the file. // +// Alternatively, SaveData(FILE* stream) will write the data to // +// already existing stream. // +// // +// To read back already stored interpolation use either the constructor // +// AliCheb3D(const char* inpFile); // +// or the default constructor AliCheb3D() followed by // +// AliCheb3D::LoadData(const char* inpFile); // +// // +// To compute the interpolation use Eval(float* par,float *res) method, with // +// par being 3D vector of arguments (inside the validity region) and res is // +// the array of DimOut elements for the output. // +// // +// If only one component (say, idim-th) of the output is needed, use faster // +// Float_t Eval(Float_t *par,int idim) method. // +// // +// void Print(option="") will print the name, the ranges of validity and // +// the absolute precision of the parameterization. Option "l" will also print // +// the information about the number of coefficients for each output // +// dimension. // +// // +// NOTE: during the evaluation no check is done for parameter vector being // +// outside the interpolation region. If there is such a risk, use // +// Bool_t IsInside(float *par) method. Chebyshev parameterization is not // +// good for extrapolation! // +// // +// For the properties of Chebyshev parameterization see: // +// H.Wind, CERN EP Internal Report, 81-12/Rev. // +// // +//////////////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include +#include "AliCheb3D.h" + + + +ClassImp(AliCheb3DCalc) + + +AliCheb3DCalc::AliCheb3DCalc(): + TNamed("", ""), + fNCoefs(0), + fNRows(0), + fNCols(0), + fNElemBound2D(0), + fNColsAtRow(0), + fColAtRowBg(0), + fCoefBound2D0(0), + fCoefBound2D1(0), + fCoefs(0), + fTmpCf1(0), + fTmpCf0(0) +{ + // Default constructor + Init0(); +} + +AliCheb3DCalc::AliCheb3DCalc(FILE* stream): + TNamed("", ""), + fNCoefs(0), + fNRows(0), + fNCols(0), + fNElemBound2D(0), + fNColsAtRow(0), + fColAtRowBg(0), + fCoefBound2D0(0), + fCoefBound2D1(0), + fCoefs(0), + fTmpCf1(0), + fTmpCf0(0) +{ + // Default constructor + Init0(); + LoadData(stream); +} + +//__________________________________________________________________________________________ +void AliCheb3DCalc::Clear(Option_t*) +{ + // delete all dynamycally allocated structures + if (fTmpCf1) { delete[] fTmpCf1; fTmpCf1 = 0;} + if (fTmpCf0) { delete[] fTmpCf0; fTmpCf0 = 0;} + if (fCoefs) { delete[] fCoefs; fCoefs = 0;} + if (fCoefBound2D0) { delete[] fCoefBound2D0; fCoefBound2D0 = 0; } + if (fCoefBound2D1) { delete[] fCoefBound2D1; fCoefBound2D1 = 0; } + if (fNColsAtRow) { delete[] fNColsAtRow; fNColsAtRow = 0; } + if (fColAtRowBg) { delete[] fColAtRowBg; fColAtRowBg = 0; } + // +} + +//__________________________________________________________________________________________ +void AliCheb3DCalc::Init0() +{ + // reset everything to 0 + fNCoefs = fNRows = fNCols = fNElemBound2D = 0; + fCoefs = 0; + fCoefBound2D0 = fCoefBound2D1 = 0; + fNColsAtRow = fColAtRowBg = 0; + fTmpCf0 = fTmpCf1 = 0; +} + +//__________________________________________________________________________________________ +void AliCheb3DCalc::Print(Option_t* ) const +{ + printf("Chebyshev parameterization data %s for 3D->1 function.\n",GetName()); + int nmax3d = 0; + for (int i=fNElemBound2D;i--;) if (fCoefBound2D0[i]>nmax3d) nmax3d = fCoefBound2D0[i]; + printf("%d coefficients in %dx%dx%d matrix\n",fNCoefs,fNRows,fNCols,nmax3d); + // +} + +//__________________________________________________________________________________________ +Float_t AliCheb3DCalc::Eval(Float_t *par) const +{ + // evaluate Chebyshev parameterization for 3D function. + // VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval + Float_t &z = par[2]; + Float_t &y = par[1]; + Float_t &x = par[0]; + // + int ncfRC; + for (int id0=fNRows;id0--;) { + int nCLoc = fNColsAtRow[id0]; // number of significant coefs on this row + int Col0 = fColAtRowBg[id0]; // beginning of local column in the 2D boundary matrix + for (int id1=nCLoc;id1--;) { + int id = id1+Col0; + fTmpCf1[id1] = (ncfRC=fCoefBound2D0[id]) ? ChebEval1D(z,fCoefs + fCoefBound2D1[id], ncfRC) : 0.0; + } + fTmpCf0[id0] = nCLoc>0 ? ChebEval1D(y,fTmpCf1,nCLoc):0.0; + } + return ChebEval1D(x,fTmpCf0,fNRows); + // +} + +//_______________________________________________ +#ifdef _INC_CREATION_ALICHEB3D_ +void AliCheb3DCalc::SaveData(const char* outfile,Bool_t append) const +{ + // writes coefficients data to output text file, optionallt appending on the end of existing file + TString strf = outfile; + gSystem->ExpandPathName(strf); + FILE* stream = fopen(strf,append ? "a":"w"); + SaveData(stream); + fclose(stream); + // +} +#endif + +//_______________________________________________ +#ifdef _INC_CREATION_ALICHEB3D_ +void AliCheb3DCalc::SaveData(FILE* stream) const +{ + // writes coefficients data to existing output stream + // Note: fNCols, fNElemBound2D and fColAtRowBg is not stored, will be computed on fly during the loading of this file + fprintf(stream,"#\nSTART %s\n",GetName()); + fprintf(stream,"# Number of rows\n%d\n",fNRows); + // + fprintf(stream,"# Number of columns per row\n"); + for (int i=0;i\", found \"%s\"\nStop\n",buffs.Data());exit(1);} + if (buffs.First(' ')>0) SetName(buffs.Data()+buffs.First(' ')+1); + // + ReadLine(buffs,stream); // NRows + fNRows = buffs.Atoi(); + if (fNRows<1) {Error("LoadData","Expected: '', found \"%s\"\nStop\n",buffs.Data());exit(1);} + // + fNCols = 0; + fNElemBound2D = 0; + InitRows(fNRows); + // + for (int id0=0;id0%dD function. Precision: %e\n",GetName(),fDimOut,fPrec); + printf("Region of validity: [%+.5e:%+.5e] [%+.5e:%+.5e] [%+.5e:%+.5e]\n",fBMin[0],fBMax[0],fBMin[1],fBMax[1],fBMin[2],fBMax[2]); + TString opts = opt; opts.ToLower(); + if (opts.Contains("l")) for (int i=0;iPrint();} + // +} + +//__________________________________________________________________________________________ +void AliCheb3D::Init0() +{ + for (int i=3;i--;) fBMin[i] = fBMax[i] = fBScale[i] = fBOffset[i] = 0; + fMaxCoefs = 0; + fGrid = 0; + fResTmp = 0; + fUsrFunName = ""; + fUsrMacro = 0; +#ifdef _INC_CREATION_ALICHEB3D_ + gUsrFunAliCheb3D = 0; +#endif +} + +//__________________________________________________________________________________________ +void AliCheb3D::PrepareBoundaries(Float_t *bmin,Float_t *bmax) +{ + // Set and check boundaries defined by user, prepare coefficients for their conversion to [-1:1] interval + // + for (int i=3;i--;) { + fBMin[i] = bmin[i]; + fBMax[i] = bmax[i]; + fBScale[i] = bmax[i]-bmin[i]; + if (fBScale[i]<=0) { + Error("PrepareBoundaries","Boundaries for %d-th dimension are not increasing: %+.4e %+.4e\nStop\n",i,fBMin[i],fBMax[i]); + exit(1); + } + fBOffset[i] = bmin[i] + fBScale[i]/2.0; + fBScale[i] = 2./fBScale[i]; + } + // +} + +//__________________________________________________________________________________________ +#ifdef _INC_CREATION_ALICHEB3D_ +void AliCheb3D::SetUsrFunction(const char* name) +{ + // load user macro with function definition and compile it + gUsrFunAliCheb3D = 0; + fUsrFunName = name; + gSystem->ExpandPathName(fUsrFunName); + if (fUsrMacro) delete fUsrMacro; + TString tmpst = fUsrFunName; + tmpst += "+"; // prepare filename to compile + if (gROOT->LoadMacro(tmpst.Data())) {Error("SetUsrFunction","Failed to load user function from %s\nStop\n",name); exit(1);} + fUsrMacro = new TMethodCall(); + tmpst = tmpst.Data() + tmpst.Last('/')+1; //Strip away any path preceding the macro file name + int dot = tmpst.Last('.'); + if (dot>0) tmpst.Resize(dot); + fUsrMacro->InitWithPrototype(tmpst.Data(),"Float_t *,Float_t *"); + long args[2]; + args[0] = (long)fArgsTmp; + args[1] = (long)fResTmp; + fUsrMacro->SetParamPtrs(args); + // +} +#endif + +//__________________________________________________________________________________________ +#ifdef _INC_CREATION_ALICHEB3D_ +void AliCheb3D::SetUsrFunction(void (*ptr)(float*,float*)) +{ + if (fUsrMacro) delete fUsrMacro; + fUsrMacro = 0; + fUsrFunName = ""; + gUsrFunAliCheb3D = ptr; +} +#endif + +//__________________________________________________________________________________________ +#ifdef _INC_CREATION_ALICHEB3D_ +void AliCheb3D::EvalUsrFunction(Float_t *x, Float_t *res) { + for (int i=3;i--;) fArgsTmp[i] = x[i]; + if (gUsrFunAliCheb3D) gUsrFunAliCheb3D(fArgsTmp,fResTmp); + else fUsrMacro->Execute(); + for (int i=fDimOut;i--;) res[i] = fResTmp[i]; +} +#endif + +//__________________________________________________________________________________________ +#ifdef _INC_CREATION_ALICHEB3D_ +Int_t AliCheb3D::CalcChebCoefs(Float_t *funval,int np, Float_t *outCoefs, Float_t prec) +{ + // Calculate Chebyshev coeffs using precomputed function values at np roots. + // If prec>0, estimate the highest coeff number providing the needed precision + // + double sm; // do summations in double to minimize the roundoff error + for (int ic=0;ic=prec) break; + } + if (++cfMax==0) cfMax=1; + return cfMax; + // +} +#endif + +//__________________________________________________________________________________________ +#ifdef _INC_CREATION_ALICHEB3D_ +void AliCheb3D::DefineGrid(Int_t* npoints) +{ + // prepare the grid of Chebyshev roots in each dimension + const int kMinPoints = 1; + int ntot = 0; + fMaxCoefs = 1; + for (int id=3;id--;) { + fNPoints[id] = npoints[id]; + if (fNPoints[id]0 && tmpCoefSurf[(id1-1)+id0*fNPoints[1]]==0) id1--; + tmpCols[id0] = id1; + } + // find max significant row + for (int id0=NRows;id0--;) {if (tmpCols[id0]>0) break; NRows--;} + // find max significant column and fill the permanent storage for the max sigificant column of each row + cheb->InitRows(NRows); // create needed arrays; + int *NColsAtRow = cheb->GetNColsAtRow(); + int *ColAtRowBg = cheb->GetColAtRowBg(); + int NCols = 0; + int NElemBound2D = 0; + for (int id0=0;id0InitCols(NCols); + delete[] tmpCols; + // + // create the 2D matrix defining the boundary of significance for 3D coeffs.matrix + // and count the number of siginifacnt coefficients + // + cheb->InitElemBound2D(NElemBound2D); + int *CoefBound2D0 = cheb->GetCoefBound2D0(); + int *CoefBound2D1 = cheb->GetCoefBound2D1(); + fMaxCoefs = 0; // redefine number of coeffs + for (int id0=0;id0InitCoefs(fMaxCoefs); + Float_t *Coefs = cheb->GetCoefs(); + int count = 0; + for (int id0=0;id0ExpandPathName(strf); + FILE* stream = fopen(strf,append ? "a":"w"); + SaveData(stream); + fclose(stream); + // +} +#endif + +//_______________________________________________ +#ifdef _INC_CREATION_ALICHEB3D_ +void AliCheb3D::SaveData(FILE* stream) const +{ + // writes coefficients data to existing output stream + // + fprintf(stream,"\n# These are automatically generated data for the Chebyshev interpolation of 3D->%dD function\n",fDimOut); + fprintf(stream,"#\nSTART %s\n",GetName()); + fprintf(stream,"# Dimensionality of the output\n%d\n",fDimOut); + fprintf(stream,"# Interpolation abs. precision\n%+.8e\n",fPrec); + // + fprintf(stream,"# Lower boundaries of interpolation region\n"); + for (int i=0;i<3;i++) fprintf(stream,"%+.8e\n",fBMin[i]); + fprintf(stream,"# Upper boundaries of interpolation region\n"); + for (int i=0;i<3;i++) fprintf(stream,"%+.8e\n",fBMax[i]); + fprintf(stream,"# Parameterization for each output dimension follows:\n",GetName()); + // + for (int i=0;iSaveData(stream); + fprintf(stream,"#\nEND %s\n#\n",GetName()); + // +} +#endif + +//_______________________________________________ +void AliCheb3D::LoadData(const char* inpFile) +{ + TString strf = inpFile; + gSystem->ExpandPathName(strf); + FILE* stream = fopen(strf.Data(),"r"); + LoadData(stream); + fclose(stream); + // +} + +//_______________________________________________ +void AliCheb3D::LoadData(FILE* stream) +{ + if (!stream) {Error("LoadData","No stream provided.\nStop"); exit(1);} + TString buffs; + Clear(); + AliCheb3DCalc::ReadLine(buffs,stream); + if (!buffs.BeginsWith("START")) {Error("LoadData","Expected: \"START \", found \"%s\"\nStop\n",buffs.Data());exit(1);} + SetName(buffs.Data()+buffs.First(' ')+1); + // + AliCheb3DCalc::ReadLine(buffs,stream); // N output dimensions + fDimOut = buffs.Atoi(); + if (fDimOut<1) {Error("LoadData","Expected: '', found \"%s\"\nStop\n",buffs.Data());exit(1);} + // + SetDimOut(fDimOut); + // + AliCheb3DCalc::ReadLine(buffs,stream); // Interpolation abs. precision + fPrec = buffs.Atof(); + if (fPrec<=0) {Error("LoadData","Expected: '', found \"%s\"\nStop\n",buffs.Data());exit(1);} + // + for (int i=0;i<3;i++) { // Lower boundaries of interpolation region + AliCheb3DCalc::ReadLine(buffs,stream); + fBMin[i] = buffs.Atof(); + } + for (int i=0;i<3;i++) { // Upper boundaries of interpolation region + AliCheb3DCalc::ReadLine(buffs,stream); + fBMax[i] = buffs.Atof(); + } + PrepareBoundaries(fBMin,fBMax); + // + // data for each output dimension + for (int i=0;iLoadData(stream); + // + // check end_of_data record + 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); + } + // +} + +//_______________________________________________ +void AliCheb3D::SetDimOut(int d) +{ + fDimOut = d; + if (fResTmp) delete fResTmp; + fResTmp = new Float_t[fDimOut]; // RRR + fChebCalc.Delete(); + for (int i=0;i2) {printf("Maximum 3 dimensions are supported\n"); return;} + fBMin[id] += dif; + fBMax[id] += dif; + fBOffset[id] += dif; +} + +//_______________________________________________ +#ifdef _INC_CREATION_ALICHEB3D_ +TH1* AliCheb3D::TestRMS(int idim,int npoints,TH1* histo) +{ + // fills the difference between the original function and parameterization (for idim-th component of the output) + // to supplied histogram. Calculations are done in npoints random points. + // If the hostgram was not supplied, it will be created. It is up to the user to delete it! + if (!fUsrMacro) { + printf("No user function is set\n"); + return 0; + } + if (!histo) histo = new TH1D(GetName(),"Control: Function - Parametrization",100,-2*fPrec,2*fPrec); + for (int ip=npoints;ip--;) { + gRandom->RndmArray(3,(Float_t *)fArgsTmp); + for (int i=3;i--;) fArgsTmp[i] = fBMin[i] + fArgsTmp[i]*(fBMax[i]-fBMin[i]); + EvalUsrFunction(); + Float_t valFun = fResTmp[idim]; + Eval(fArgsTmp,fResTmp); + Float_t valPar = fResTmp[idim]; + histo->Fill(valFun - valPar); + } + return histo; + // +} +#endif diff --git a/STEER/AliCheb3D.h b/STEER/AliCheb3D.h new file mode 100644 index 00000000000..9151951418f --- /dev/null +++ b/STEER/AliCheb3D.h @@ -0,0 +1,294 @@ +#ifndef ALICHEB3D_H +#define ALICHEB3D_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 09/09/2006 +// +//////////////////////////////////////////////////////////////////////////////// +// // +// AliCheb3D produces the interpolation of the user 3D->NDimOut arbitrary // +// function supplied in "void (*fcn)(float* inp,float* out)" format // +// either in a separate macro file or as a function pointer. // +// Only coefficients needed to guarantee the requested precision are kept. // +// // +// The user-callable methods are: // +// To create the interpolation use: // +// AliCheb3D(const char* funName, // name of the file with user function // +// or // +// AliCheb3D(void (*ptr)(float*,float*),// pointer on the user function // +// Int_t DimOut, // dimensionality of the function's output // +// Float_t *bmin, // lower 3D bounds of interpolation domain // +// Float_t *bmax, // upper 3D bounds of interpolation domain // +// Int_t *npoints, // number of points in each of 3 input // +// // dimension, defining the interpolation grid // +// Float_t prec=1E-6); // requested max.absolute difference between // +// // the interpolation and any point on grid // +// // +// To test obtained parameterization use the method // +// TH1* TestRMS(int idim,int npoints = 1000,TH1* histo=0); // +// it will compare the user output of the user function and interpolation // +// for idim-th output dimension and fill the difference in the supplied // +// histogram. If no histogram is supplied, it will be created. // +// // +// To save the interpolation data: // +// SaveData(const char* filename, Bool_t append ) // +// write text file with data. If append is kTRUE and the output file already // +// exists, data will be added in the end of the file. // +// Alternatively, SaveData(FILE* stream) will write the data to // +// already existing stream. // +// // +// To read back already stored interpolation use either the constructor // +// AliCheb3D(const char* inpFile); // +// or the default constructor AliCheb3D() followed by // +// AliCheb3D::LoadData(const char* inpFile); // +// // +// To compute the interpolation use Eval(float* par,float *res) method, with // +// par being 3D vector of arguments (inside the validity region) and res is // +// the array of DimOut elements for the output. // +// // +// If only one component (say, idim-th) of the output is needed, use faster // +// Float_t Eval(Float_t *par,int idim) method. // +// // +// void Print(option="") will print the name, the ranges of validity and // +// the absolute precision of the parameterization. Option "l" will also print // +// the information about the number of coefficients for each output // +// dimension. // +// // +// NOTE: during the evaluation no check is done for parameter vector being // +// outside the interpolation region. If there is such a risk, use // +// Bool_t IsInside(float *par) method. Chebyshev parameterization is not // +// good for extrapolation! // +// // +// For the properties of Chebyshev parameterization see: // +// H.Wind, CERN EP Internal Report, 81-12/Rev. // +// // +//////////////////////////////////////////////////////////////////////////////// + + +#include +#include +#include +#include +#include +#include + +class TString; +class TSystem; +class TRandom; + + +// to decrease the compilable code size comment this define. This will exclude the routines +// used for the calculation and saving of the coefficients. +// #define _INC_CREATION_ALICHEB3D_ + +class AliCheb3DCalc: public TNamed +{ + public: + AliCheb3DCalc(); + AliCheb3DCalc(FILE* stream); // read coefs from text file + ~AliCheb3DCalc() {Clear();} + // + void Print(Option_t* opt="") const; + void LoadData(FILE* stream); + Float_t Eval(Float_t *par) const; + // +#ifdef _INC_CREATION_ALICHEB3D_ + void SaveData(const char* outfile,Bool_t append=kFALSE) const; + void SaveData(FILE* stream=stdout) const; +#endif + // + static void ReadLine(TString& str,FILE* stream); + // + protected: + // + void Clear(Option_t* option = ""); + void Init0(); + Float_t ChebEval1D(Float_t x, const Float_t * array, int ncf) const; + void InitRows(int nr); + void InitCols(int nc); + void InitElemBound2D(int ne); + void InitCoefs(int nc); + Int_t* GetNColsAtRow() const {return fNColsAtRow;} + Int_t* GetColAtRowBg() const {return fColAtRowBg;} + Int_t* GetCoefBound2D0() const {return fCoefBound2D0;} + Int_t* GetCoefBound2D1() const {return fCoefBound2D1;} + Float_t * GetCoefs() const {return fCoefs;} + // + protected: + Int_t fNCoefs; // total number of coeeficients + Int_t fNRows; // number of significant rows in the 3D coeffs matrix + Int_t fNCols; // max number of significant cols in the 3D coeffs matrix + Int_t fNElemBound2D; // number of elements (fNRows*fNCols) to store for the 2D boundary of significant coeffs + Int_t* fNColsAtRow; //[fNRows] number of sighificant columns (2nd dim) at each row of 3D coefs matrix + Int_t* fColAtRowBg; //[fNRows] beginnig of significant columns (2nd dim) for row in the 2D boundary matrix + Int_t* fCoefBound2D0; //[fNElemBound2D] 2D matrix defining the boundary of significance for 3D coeffs.matrix (Ncoefs for col/row) + Int_t* fCoefBound2D1; //[fNElemBound2D] 2D matrix defining the start beginnig of significant coeffs for col/row + Float_t * fCoefs; //[fNCoefs] array of Chebyshev coefficients + // + Float_t * fTmpCf1; //[fNCols] temp. coeffs for 2d summation + Float_t * fTmpCf0; //[fNRows] temp. coeffs for 1d summation + // + ClassDef(AliCheb3DCalc,1) // Class for interpolation of 3D->1 function by Chebyshev parametrization +}; + + +class AliCheb3D: public TNamed +{ + public: + AliCheb3D(); + AliCheb3D(const char* inpFile); // read coefs from text file + AliCheb3D(FILE*); // read coefs from stream + // +#ifdef _INC_CREATION_ALICHEB3D_ + AliCheb3D(const char* funName, Int_t DimOut, Float_t *bmin,Float_t *bmax, Int_t *npoints, Float_t prec=1E-6); + AliCheb3D(void (*ptr)(float*,float*), Int_t DimOut, Float_t *bmin,Float_t *bmax, Int_t *npoints, Float_t prec=1E-6); +#endif + // + ~AliCheb3D() {Clear();} + // + void Eval(Float_t *par,Float_t *res); + Float_t Eval(Float_t *par,int idim); + void Print(Option_t* opt="") const; + Bool_t IsInside(Float_t *par) const; + AliCheb3DCalc* GetChebCalc(int i) const {return (AliCheb3DCalc*)fChebCalc.UncheckedAt(i);} + Float_t GetBoundMin(int i) const {return fBMin[i];} + Float_t GetBoundMax(int i) const {return fBMax[i];} + Float_t GetPrecision() const {return fPrec;} + void ShiftBound(int id,float dif); + // + void LoadData(const char* inpFile); + void LoadData(FILE* stream); + // +#ifdef _INC_CREATION_ALICHEB3D_ + void SaveData(const char* outfile,Bool_t append=kFALSE) const; + void SaveData(FILE* stream=stdout) const; + // + void SetUsrFunction(const char* name); + void SetUsrFunction(void (*ptr)(float*,float*)); + void EvalUsrFunction(Float_t *x, Float_t *res); + TH1* TestRMS(int idim,int npoints = 1000,TH1* histo=0); +#endif + // + protected: + void Init0(); + void Clear(Option_t* option = ""); + void SetDimOut(int d); + void PrepareBoundaries(Float_t *bmin,Float_t *bmax); + // +#ifdef _INC_CREATION_ALICHEB3D_ + void EvalUsrFunction(); + void DefineGrid(Int_t* npoints); + Int_t ChebFit(); // fit all output dimensions + Int_t ChebFit(int dmOut); + Int_t CalcChebCoefs(Float_t *funval,int np, Float_t *outCoefs, Float_t prec=-1); +#endif + // + void Cyl2CartCyl(float *rphiz, float *b) const; + void Cart2Cyl(float *xyz,float *rphiz) const; + // + Float_t MapToInternal(Float_t x,Int_t d) const {return (x-fBOffset[d])*fBScale[d];} // map x to [-1:1] + Float_t MapToExternal(Float_t x,Int_t d) const {return x/fBScale[d]+fBOffset[d];} // map from [-1:1] to x + // + protected: + Int_t fDimOut; // dimension of the ouput array + Float_t fPrec; // requested precision + Float_t fBMin[3]; // min boundaries in each dimension + Float_t fBMax[3]; // max boundaries in each dimension + Float_t fBScale[3]; // scale for boundary mapping to [-1:1] interval + Float_t fBOffset[3]; // offset for boundary mapping to [-1:1] interval + TObjArray fChebCalc; // Chebyshev parameterization for each output dimension + // + Int_t fMaxCoefs; //! max possible number of coefs per parameterization + Int_t fNPoints[3]; //! number of used points in each dimension + Float_t fArgsTmp[3]; //! temporary vector for coefs caluclation + Float_t fBuff[6]; //! buffer for coordinate transformations + Float_t * fResTmp; //! temporary vector for results of user function caluclation + Float_t * fGrid; //! temporary buffer for Chebyshef roots grid + Int_t fGridOffs[3]; //! start of grid for each dimension + TString fUsrFunName; //! name of user macro containing the function of "void (*fcn)(float*,float*)" format + TMethodCall* fUsrMacro; //! Pointer to MethodCall for function from user macro + // + ClassDef(AliCheb3D,1) // Chebyshev parametrization for 3D->N function +}; + +// Pointer on user function (faster altrnative to TMethodCall) +#ifdef _INC_CREATION_ALICHEB3D_ +void (*gUsrFunAliCheb3D) (float* ,float* ); +#endif + +//__________________________________________________________________________________________ +#ifdef _INC_CREATION_ALICHEB3D_ +inline void AliCheb3D::EvalUsrFunction() +{ + // call user supplied function + if (gUsrFunAliCheb3D) gUsrFunAliCheb3D(fArgsTmp,fResTmp); + else fUsrMacro->Execute(); +} +#endif + +//__________________________________________________________________________________________ +inline Bool_t AliCheb3D::IsInside(Float_t *par) const +{ + // check if the point is inside of the fitted box + for (int i=3;i--;) if(par[i]fBMax[i]) return kFALSE; + return kTRUE; +} + +//__________________________________________________________________________________________ +inline Float_t AliCheb3DCalc::ChebEval1D(Float_t x, const Float_t * array, int ncf ) const +{ + // evaluate 1D Chebyshev parameterization. x is the argument mapped to [-1:1] interval + Float_t b0, b1, b2; + Float_t x2 = x+x; + b0 = array[--ncf]; + b1 = b2 = 0; + for (int i=ncf;i--;) { + b2 = b1; + b1 = b0; + b0 = array[i] + x2*b1 -b2; + } + return b0 - x*b1; +} + +//__________________________________________________________________________________________ +inline void AliCheb3D::Eval(Float_t *par, Float_t *res) +{ + // evaluate Chebyshev parameterization for 3d->DimOut function + for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i); + for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->Eval(fArgsTmp); + // +} + +//__________________________________________________________________________________________ +inline Float_t AliCheb3D::Eval(Float_t *par, int idim) +{ + // evaluate Chebyshev parameterization for idim-th output dimension of 3d->DimOut function + for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i); + return GetChebCalc(idim)->Eval(fArgsTmp); + // +} + +//__________________________________________________________________________________________________ +inline void AliCheb3D::Cyl2CartCyl(float *rphiz, float *b) const +{ + // convert field in cylindrical coordinates to cartesian system, point is in cyl.system + float btr = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]); + float ang = TMath::ATan2(b[1],b[0]) + rphiz[1]; + b[0] = btr*TMath::Cos(ang); + b[1] = btr*TMath::Sin(ang); + // +} + +//__________________________________________________________________________________________________ +inline void AliCheb3D::Cart2Cyl(float *xyz,float *rphiz) const +{ + // convert cartesian coordinate to cylindrical one + rphiz[0] = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); + rphiz[1] = TMath::ATan2(xyz[1],xyz[0]); + rphiz[2] = xyz[2]; +} + +#endif diff --git a/STEER/AliMagFCheb.cxx b/STEER/AliMagFCheb.cxx new file mode 100644 index 00000000000..a954393458d --- /dev/null +++ b/STEER/AliMagFCheb.cxx @@ -0,0 +1,306 @@ +/************************************************************************** + * 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 // +// To obtain the field in cartesian coordinates/components use // +// Field(float* xyz, float* bxyz); // +// 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, -550Add(param); + fNParamsSol++; + // +} + +//__________________________________________________________________________________________ +void AliMagFCheb::AddParamDip(AliCheb3D* param) +{ + // adds new parameterization piece for Dipole + // + if (!fParamsDip) fParamsDip = new TObjArray(); + fParamsDip->Add(param); + fNParamsDip++; + // +} + +//__________________________________________________________________________________________ +void AliMagFCheb::BuildTableSol() +{ + // build the indexes for each parameterization of Solenoid + // + const float kSafety=0.001; + // + 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]++; + } + // + fSegZSol = new Float_t[fNSegZSol]; + fSegZIdSol = new Int_t[fNSegZSol]; + fNSegRSol = new Int_t[fNSegZSol]; + for (int iz=0;izGetBoundMin(2); + fMaxZSol = GetParamSol(fNParamsSol-1)->GetBoundMax(2); + fMaxRSol = GetParamSol(fNParamsSol-1)->GetBoundMax(0); + // + delete[] tmpbufF; + delete[] tmpbufI; + delete[] tmpbufI1; + // + // +} + +//__________________________________________________________________________________________ +void AliMagFCheb::Field(Float_t *xyz, Float_t *b) const +{ + // compute field in cartesian coordinates + float rphiz[3]; + if (xyz[2]>GetMaxZSol() || xyz[2]GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;} + // + FieldCylSol(rphiz,b); + // + // 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); + // +} + +//__________________________________________________________________________________________ +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); + // +} + +//__________________________________________________________________________________________ +void AliMagFCheb::Print(Option_t *) const +{ + printf("Alice magnetic field parameterized by Chebyshev polynomials\n"); + printf("Segmentation for Solenoid (%+.2fGetBoundMin(2),param->GetBoundMax(2)); + for (int ir=0;irGetBoundMin(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 + 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"); + // + // 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); + // +} +#endif + +//_______________________________________________ +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 \", 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);} + // + // 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(); + // BuildDipTable(); + printf("Loaded magnetic field \"%s\" from %s\n",GetName(),strf.Data()); + // +} diff --git a/STEER/AliMagFCheb.h b/STEER/AliMagFCheb.h new file mode 100644 index 00000000000..157176254a6 --- /dev/null +++ b/STEER/AliMagFCheb.h @@ -0,0 +1,107 @@ +#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 // +// To obtain the field in cartesian coordinates/components use // +// Field(float* xyz, float* bxyz); // +// 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 +#include +#include "AliCheb3D.h" + +class AliMagFCheb: public TNamed +{ + public: + AliMagFCheb() {Init0();} + AliMagFCheb(const char* inputFile); + ~AliMagFCheb(); + // + void AddParamSol(AliCheb3D* param); + void AddParamDip(AliCheb3D* param); + void BuildTableSol(); + // + Int_t GetNParamsSol() const {return fNParamsSol;} + Int_t GetNSegZSol() const {return fNSegZSol;} + Int_t GetNSegRSol(int iz) const {return izUncheckedAt(ipar);} + AliCheb3D* GetParamDip(Int_t ipar) const {return (AliCheb3D*)fParamsDip->UncheckedAt(ipar);} + // + void LoadData(const char* inpfile); + // + virtual void Print(Option_t * = "") const; + // + virtual void Field(Float_t *xyz, Float_t *b) const; + virtual void FieldCyl(Float_t *rphiz, Float_t *b) const; + // + // +#ifdef _INC_CREATION_ALICHEB3D_ // see AliCheb3D.h for explanation + void SaveData(const char* outfile) const; +#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 fNParamsDip; // Total number of parameterization pieces for dipole + // + Float_t* fSegZSol; //[fNSegZSol] upper boundaries of Z segments + Float_t* fSegRSol; //[fNParamsSol] upper boundaries of R segments + // + 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... + // + 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) + // + TObjArray* fParamsSol; // Parameterization pieces for Solenoid field + TObjArray* fParamsDip; // Parameterization pieces for Dipole field + // + ClassDef(AliMagFCheb,1) // Wrapper class for the set of Chebishev parameterizations of Alice mag.field + // + }; + + +//__________________________________________________________________________________________ +inline void AliMagFCheb::FieldCyl(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); +} + +#endif