STEER/AliCheb3DCalc.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 #include <cstdlib>
  17 #include "AliCheb3DCalc.h"
  18 #include <TSystem.h>
  19
  20 ClassImp(AliCheb3DCalc)
  21
  22 //__________________________________________________________________________________________
  23 AliCheb3DCalc::AliCheb3DCalc() :
  24   fNCoefs(0),
  25   fNRows(0),
  26   fNCols(0),
  27   fNElemBound2D(0),
  28   fNColsAtRow(0),
  29   fColAtRowBg(0),
  30   fCoefBound2D0(0),
  31   fCoefBound2D1(0),
  32   fCoefs(0),
  33   fTmpCf1(0),
  34   fTmpCf0(0)
  35 {
  36   // default constructor
  37 }
  38
  39 //__________________________________________________________________________________________
  40 AliCheb3DCalc::AliCheb3DCalc(const AliCheb3DCalc& src) :
  41   TNamed(src),
  42   fNCoefs(src.fNCoefs),
  43   fNRows(src.fNRows),
  44   fNCols(src.fNCols),
  45   fNElemBound2D(src.fNElemBound2D),
  46   fNColsAtRow(0),
  47   fColAtRowBg(0),
  48   fCoefBound2D0(0),
  49   fCoefBound2D1(0),
  50   fCoefs(0),
  51   fTmpCf1(0),
  52   fTmpCf0(0)
  53 {
  54   // copy constructor
  55   //
  56   if (src.fNColsAtRow) {
  57     fNColsAtRow = new Int_t[fNRows];
  58     for (int i=fNRows;i--;) fNColsAtRow[i] = src.fNColsAtRow[i];
  59   }
  60   if (src.fColAtRowBg) {
  61     fColAtRowBg = new Int_t[fNRows];
  62     for (int i=fNRows;i--;) fColAtRowBg[i] = src.fColAtRowBg[i];
  63   }
  64   if (src.fCoefBound2D0) {
  65     fCoefBound2D0 = new Int_t[fNElemBound2D];
  66     for (int i=fNElemBound2D;i--;) fCoefBound2D0[i] = src.fCoefBound2D0[i];
  67   }
  68   if (src.fCoefBound2D1) {
  69     fCoefBound2D1 = new Int_t[fNElemBound2D];
  70     for (int i=fNElemBound2D;i--;) fCoefBound2D1[i] = src.fCoefBound2D1[i];
  71   }
  72   if (src.fCoefs) {
  73     fCoefs = new Float_t[fNCoefs];
  74     for (int i=fNCoefs;i--;) fCoefs[i] = src.fCoefs[i];
  75   }
  76   if (src.fTmpCf1) fTmpCf1 = new Float_t[fNCols];
  77   if (src.fTmpCf0) fTmpCf0 = new Float_t[fNRows];
  78 }
  79
  80 //__________________________________________________________________________________________
  81 AliCheb3DCalc::AliCheb3DCalc(FILE* stream) :
  82   fNCoefs(0),
  83   fNRows(0),
  84   fNCols(0),
  85   fNElemBound2D(0),
  86   fNColsAtRow(0),
  87   fColAtRowBg(0),
  88   fCoefBound2D0(0),
  89   fCoefBound2D1(0),
  90   fCoefs(0),
  91   fTmpCf1(0),
  92   fTmpCf0(0)
  93 {
  94   // constructor from coeffs. streem
  95   LoadData(stream);
  96 }
  97
  98 //__________________________________________________________________________________________
  99 AliCheb3DCalc& AliCheb3DCalc::operator=(const AliCheb3DCalc& rhs)
 100 {
 101   // assignment operator
 102   if (this != &rhs) {
 103     Clear();
 104     SetName(rhs.GetName());
 105     SetTitle(rhs.GetTitle());
 106     fNCoefs = rhs.fNCoefs;
 107     fNRows  = rhs.fNRows;
 108     fNCols  = rhs.fNCols;
 109     if (rhs.fNColsAtRow) {
 110       fNColsAtRow = new Int_t[fNRows];
 111       for (int i=fNRows;i--;) fNColsAtRow[i] = rhs.fNColsAtRow[i];
 112     }
 113     if (rhs.fColAtRowBg) {
 114       fColAtRowBg = new Int_t[fNRows];
 115       for (int i=fNRows;i--;) fColAtRowBg[i] = rhs.fColAtRowBg[i];
 116     }
 117     if (rhs.fCoefBound2D0) {
 118       fCoefBound2D0 = new Int_t[fNElemBound2D];
 119       for (int i=fNElemBound2D;i--;) fCoefBound2D0[i] = rhs.fCoefBound2D0[i];
 120     }
 121     if (rhs.fCoefBound2D1) {
 122       fCoefBound2D1 = new Int_t[fNElemBound2D];
 123       for (int i=fNElemBound2D;i--;) fCoefBound2D1[i] = rhs.fCoefBound2D1[i];
 124     }
 125     if (rhs.fCoefs) {
 126       fCoefs = new Float_t[fNCoefs];
 127       for (int i=fNCoefs;i--;) fCoefs[i] = rhs.fCoefs[i];
 128     }
 129     if (rhs.fTmpCf1) fTmpCf1 = new Float_t[fNCols];
 130     if (rhs.fTmpCf0) fTmpCf0 = new Float_t[fNRows];
 131   }
 132   return *this;
 133 }
 134
 135 //__________________________________________________________________________________________
 136 void AliCheb3DCalc::Clear(const Option_t*)
 137 {
 138   // delete all dynamycally allocated structures
 139   if (fTmpCf1)       { delete[] fTmpCf1;  fTmpCf1 = 0;}
 140   if (fTmpCf0)       { delete[] fTmpCf0;  fTmpCf0 = 0;}
 141   if (fCoefs)        { delete[] fCoefs;   fCoefs  = 0;}
 142   if (fCoefBound2D0) { delete[] fCoefBound2D0; fCoefBound2D0 = 0; }
 143   if (fCoefBound2D1) { delete[] fCoefBound2D1; fCoefBound2D1 = 0; }
 144   if (fNColsAtRow)   { delete[] fNColsAtRow;   fNColsAtRow = 0; }
 145   if (fColAtRowBg)   { delete[] fColAtRowBg;   fColAtRowBg = 0; }
 146   //
 147 }
 148
 149 //__________________________________________________________________________________________
 150 void AliCheb3DCalc::Print(const Option_t* ) const
 151 {
 152   // print info
 153   printf("Chebyshev parameterization data %s for 3D->1 function.\n",GetName());
 154   int nmax3d = 0;
 155   for (int i=fNElemBound2D;i--;) if (fCoefBound2D0[i]>nmax3d) nmax3d = fCoefBound2D0[i];
 156   printf("%d coefficients in %dx%dx%d matrix\n",fNCoefs,fNRows,fNCols,nmax3d);
 157   //
 158 }
 159
 160 //__________________________________________________________________________________________
 161 Float_t  AliCheb3DCalc::Eval(const Float_t  *par) const
 162 {
 163   // evaluate Chebyshev parameterization for 3D function.
 164   // VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval
 165   const Float_t  &z = par[2];
 166   const Float_t  &y = par[1];
 167   const Float_t  &x = par[0];
 168   //
 169   int ncfRC;
 170   for (int id0=fNRows;id0--;) {
 171     int nCLoc = fNColsAtRow[id0];                   // number of significant coefs on this row
 172     int col0  = fColAtRowBg[id0];                   // beginning of local column in the 2D boundary matrix
 173     for (int id1=nCLoc;id1--;) {
 174       int id = id1+col0;
 175       fTmpCf1[id1] = (ncfRC=fCoefBound2D0[id]) ? ChebEval1D(z,fCoefs + fCoefBound2D1[id], ncfRC) : 0.0;
 176     }
 177     fTmpCf0[id0] = nCLoc>0 ? ChebEval1D(y,fTmpCf1,nCLoc):0.0;
 178   }
 179   return ChebEval1D(x,fTmpCf0,fNRows);
 180   //
 181 }
 182
 183 //__________________________________________________________________________________________
 184 Float_t  AliCheb3DCalc::EvalDeriv(int dim, const Float_t  *par) const
 185 {
 186   // evaluate Chebyshev parameterization derivative in given dimension  for 3D function.
 187   // VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval
 188   const Float_t  &z = par[2];
 189   const Float_t  &y = par[1];
 190   const Float_t  &x = par[0];
 191   //
 192   int ncfRC;
 193   for (int id0=fNRows;id0--;) {
 194     int nCLoc = fNColsAtRow[id0];                   // number of significant coefs on this row
 195     if (!nCLoc) {fTmpCf0[id0]=0; continue;}
 196     //
 197     int col0  = fColAtRowBg[id0];                   // beginning of local column in the 2D boundary matrix
 198     for (int id1=nCLoc;id1--;) {
 199       int id = id1+col0;
 200       if (!(ncfRC=fCoefBound2D0[id])) { fTmpCf1[id1]=0; continue;}
 201       if (dim==2) fTmpCf1[id1] = ChebEval1Deriv(z,fCoefs + fCoefBound2D1[id], ncfRC);
 202       else        fTmpCf1[id1] = ChebEval1D(z,fCoefs + fCoefBound2D1[id], ncfRC);
 203     }
 204     if (dim==1)   fTmpCf0[id0] = ChebEval1Deriv(y,fTmpCf1,nCLoc);
 205     else          fTmpCf0[id0] = ChebEval1D(y,fTmpCf1,nCLoc);
 206   }
 207   return (dim==0) ? ChebEval1Deriv(x,fTmpCf0,fNRows) : ChebEval1D(x,fTmpCf0,fNRows);
 208   //
 209 }
 210
 211 //__________________________________________________________________________________________
 212 Float_t  AliCheb3DCalc::EvalDeriv2(int dim1,int dim2, const Float_t  *par) const
 213 {
 214   // evaluate Chebyshev parameterization 2n derivative in given dimensions  for 3D function.
 215   // VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval
 216   const Float_t  &z = par[2];
 217   const Float_t  &y = par[1];
 218   const Float_t  &x = par[0];
 219   //
 220   Bool_t same = dim1==dim2;
 221   int ncfRC;
 222   for (int id0=fNRows;id0--;) {
 223     int nCLoc = fNColsAtRow[id0];                   // number of significant coefs on this row
 224     if (!nCLoc) {fTmpCf0[id0]=0; continue;}
 225     //
 226     int col0  = fColAtRowBg[id0];                   // beginning of local column in the 2D boundary matrix
 227     for (int id1=nCLoc;id1--;) {
 228       int id = id1+col0;
 229       if (!(ncfRC=fCoefBound2D0[id])) { fTmpCf1[id1]=0; continue;}
 230       if (dim1==2||dim2==2) fTmpCf1[id1] = same ? ChebEval1Deriv2(z,fCoefs + fCoefBound2D1[id], ncfRC)
 231                               :                   ChebEval1Deriv(z,fCoefs + fCoefBound2D1[id], ncfRC);
 232       else        fTmpCf1[id1] = ChebEval1D(z,fCoefs + fCoefBound2D1[id], ncfRC);
 233     }
 234     if (dim1==1||dim2==1) fTmpCf0[id0] = same ? ChebEval1Deriv2(y,fTmpCf1,nCLoc):ChebEval1Deriv(y,fTmpCf1,nCLoc);
 235     else                  fTmpCf0[id0] = ChebEval1D(y,fTmpCf1,nCLoc);
 236   }
 237   return (dim1==0||dim2==0) ? (same ? ChebEval1Deriv2(x,fTmpCf0,fNRows):ChebEval1Deriv(x,fTmpCf0,fNRows)) : ChebEval1D(x,fTmpCf0,fNRows);
 238   //
 239 }
 240
 241 //_______________________________________________
 242 #ifdef _INC_CREATION_ALICHEB3D_
 243 void AliCheb3DCalc::SaveData(const char* outfile,Bool_t append) const
 244 {
 245   // writes coefficients data to output text file, optionallt appending on the end of existing file
 246   TString strf = outfile;
 247   gSystem->ExpandPathName(strf);
 248   FILE* stream = fopen(strf,append ? "a":"w");
 249   SaveData(stream);
 250   fclose(stream);
 251   //
 252 }
 253 #endif
 254
 255 //_______________________________________________
 256 #ifdef _INC_CREATION_ALICHEB3D_
 257 void AliCheb3DCalc::SaveData(FILE* stream) const
 258 {
 259   // writes coefficients data to existing output stream
 260   // Note: fNCols, fNElemBound2D and fColAtRowBg is not stored, will be computed on fly during the loading of this file
 261   fprintf(stream,"#\nSTART %s\n",GetName());
 262   fprintf(stream,"# Number of rows\n%d\n",fNRows);
 263   //
 264   fprintf(stream,"# Number of columns per row\n");
 265   for (int i=0;i<fNRows;i++) fprintf(stream,"%d\n",fNColsAtRow[i]);
 266   //
 267   fprintf(stream,"# Number of Coefs in each significant block of third dimension\n");
 268   for (int i=0;i<fNElemBound2D;i++) fprintf(stream,"%d\n",fCoefBound2D0[i]);
 269   //
 270   fprintf(stream,"# Coefficients\n");
 271   for (int i=0;i<fNCoefs;i++) fprintf(stream,"%+.8e\n",fCoefs[i]);
 272   fprintf(stream,"END %s\n",GetName());
 273   //
 274 }
 275 #endif
 276
 277 //_______________________________________________
 278 void AliCheb3DCalc::LoadData(FILE* stream)
 279 {
 280   // Load coefs. from the stream
 281   if (!stream) {Error("LoadData","No stream provided.\nStop"); exit(1);}
 282   TString buffs;
 283   Clear();
 284   ReadLine(buffs,stream);
 285   if (!buffs.BeginsWith("START")) {Error("LoadData","Expected: \"START <fit_name>\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
 286   if (buffs.First(' ')>0) SetName(buffs.Data()+buffs.First(' ')+1);
 287   //
 288   ReadLine(buffs,stream); // NRows
 289   fNRows = buffs.Atoi();
 290   if (fNRows<1) {Error("LoadData","Expected: '<number_of_rows>', found \"%s\"\nStop\n",buffs.Data());exit(1);}
 291   //
 292   fNCols = 0;
 293   fNElemBound2D = 0;
 294   InitRows(fNRows);
 295   //
 296   for (int id0=0;id0<fNRows;id0++) {
 297     ReadLine(buffs,stream);               // n.cols at this row
 298     fNColsAtRow[id0] = buffs.Atoi();
 299     fColAtRowBg[id0] = fNElemBound2D;     // begining of this row in 2D boundary surface
 300     fNElemBound2D   += fNColsAtRow[id0];
 301     if (fNCols<fNColsAtRow[id0]) fNCols = fNColsAtRow[id0];
 302   }
 303   InitCols(fNCols);
 304   //
 305   fNCoefs = 0;
 306   InitElemBound2D(fNElemBound2D);
 307   //
 308   for (int i=0;i<fNElemBound2D;i++) {
 309     ReadLine(buffs,stream);               // n.coeffs at 3-d dimension for the given column/row
 310     fCoefBound2D0[i] = buffs.Atoi();
 311     fCoefBound2D1[i] = fNCoefs;
 312     fNCoefs += fCoefBound2D0[i];
 313   }
 314   //
 315   if (fNCoefs<=0) {Error("LoadData","Negtive (%d) number of Chebychef coeffs. is obtained.\nStop\n",fNCoefs);exit(1);}
 316   //
 317   InitCoefs(fNCoefs);
 318   for (int i=0;i<fNCoefs;i++) {
 319     ReadLine(buffs,stream);
 320     fCoefs[i] = buffs.Atof();
 321   }
 322   // check end_of_data record
 323   ReadLine(buffs,stream);
 324   if (!buffs.BeginsWith("END") || !buffs.Contains(GetName())) {
 325     Error("LoadData","Expected \"END %s\", found \"%s\".\nStop\n",GetName(),buffs.Data());
 326     exit(1);
 327   }
 328   //
 329 }
 330
 331 //_______________________________________________
 332 void AliCheb3DCalc::ReadLine(TString& str,FILE* stream)
 333 {
 334   // read single line from the stream, skipping empty and commented lines. EOF is not expected
 335   while (str.Gets(stream)) {
 336     str = str.Strip(TString::kBoth,' ');
 337     if (str.IsNull()||str.BeginsWith("#")) continue;
 338     return;
 339   }
 340   fprintf(stderr,"AliCheb3D::ReadLine: Failed to read from stream.\nStop");exit(1); // normally, should not reach here
 341 }
 342
 343 //_______________________________________________
 344 void AliCheb3DCalc::InitCols(int nc)
 345 {
 346   // Set max.number of significant columns in the coefs matrix
 347   fNCols = nc;
 348   if (fTmpCf1) delete[] fTmpCf1;
 349   fTmpCf1 = new Float_t [fNCols];
 350 }
 351
 352 //_______________________________________________
 353 void AliCheb3DCalc::InitRows(int nr)
 354 {
 355   // Set max.number of significant rows in the coefs matrix
 356   if (fNColsAtRow) delete[] fNColsAtRow;
 357   if (fColAtRowBg) delete[] fColAtRowBg;
 358   if (fTmpCf0)     delete[] fTmpCf0;
 359   fNRows = nr;
 360   fNColsAtRow = new Int_t[fNRows];
 361   fTmpCf0     = new Float_t [fNRows];
 362   fColAtRowBg = new Int_t[fNRows];
 363   for (int i=fNRows;i--;) fNColsAtRow[i] = fColAtRowBg[i] = 0;
 364 }
 365
 366 //_______________________________________________
 367 void AliCheb3DCalc::InitElemBound2D(int ne)
 368 {
 369   // Set max number of significant coefs for given row/column of coefs 3D matrix
 370   if (fCoefBound2D0) delete[] fCoefBound2D0;
 371   if (fCoefBound2D1) delete[] fCoefBound2D1;
 372   fNElemBound2D = ne;
 373   fCoefBound2D0 = new Int_t[fNElemBound2D];
 374   fCoefBound2D1 = new Int_t[fNElemBound2D];
 375   for (int i=fNElemBound2D;i--;) fCoefBound2D0[i] = fCoefBound2D1[i] = 0;
 376 }
 377
 378 //_______________________________________________
 379 void AliCheb3DCalc::InitCoefs(int nc)
 380 {
 381   // Set total number of significant coefs
 382   if (fCoefs) delete[] fCoefs;
 383   fNCoefs = nc;
 384   fCoefs = new Float_t [fNCoefs];
 385   for (int i=fNCoefs;i--;) fCoefs[i] = 0.0;
 386 }
 387
 388 //__________________________________________________________________________________________
 389 Float_t AliCheb3DCalc::ChebEval1Deriv(Float_t  x, const Float_t * array, int ncf )
 390 {
 391   // evaluate 1D Chebyshev parameterization's derivative. x is the argument mapped to [-1:1] interval
 392   if (--ncf<1) return 0;
 393   Float_t b0, b1, b2;
 394   Float_t x2 = x+x;
 395   b1 = b2 = 0;
 396   float dcf0=0,dcf1,dcf2=0;
 397   b0 = dcf1 = 2*ncf*array[ncf];
 398   if (!(--ncf)) return b0/2;
 399   //
 400   for (int i=ncf;i--;) {
 401     b2 = b1;
 402     b1 = b0;
 403     dcf0 = dcf2 + 2*(i+1)*array[i+1];
 404     b0 = dcf0 + x2*b1 -b2;
 405     dcf2 = dcf1;
 406     dcf1 = dcf0;
 407   }
 408   //
 409   return b0 - x*b1 - dcf0/2;
 410 }
 411
 412 //__________________________________________________________________________________________
 413 Float_t AliCheb3DCalc::ChebEval1Deriv2(Float_t  x, const Float_t * array, int ncf )
 414 {
 415   // evaluate 1D Chebyshev parameterization's 2nd derivative. x is the argument mapped to [-1:1] interval
 416   if (--ncf<2) return 0;
 417   Float_t b0, b1, b2;
 418   Float_t x2 = x+x;
 419   b1 = b2 = 0;
 420   float dcf0=0,dcf1=0,dcf2=0;
 421   float ddcf0=0,ddcf1,ddcf2=0;
 422   //
 423   dcf2 = 2*ncf*array[ncf];
 424   --ncf;
 425
 426   dcf1 = 2*ncf*array[ncf];
 427   b0 = ddcf1 = 2*ncf*dcf2;
 428   //
 429   if (!(--ncf)) return b0/2;
 430   //
 431   for (int i=ncf;i--;) {
 432     b2 = b1;
 433     b1 = b0;
 434     dcf0  = dcf2 + 2*(i+1)*array[i+1];
 435     ddcf0 = ddcf2 + 2*(i+1)*dcf1;
 436     b0 = ddcf0 + x2*b1 -b2;
 437     //
 438     ddcf2 = ddcf1;
 439     ddcf1 = ddcf0;
 440     //
 441     dcf2 = dcf1;
 442     dcf1 = dcf0;
 443     //
 444   }
 445   //
 446   return b0 - x*b1 - ddcf0/2;
 447 }
 448