StRoot: StarRoot/TMDFParameters.cxx Source File

00001 #include "Riostream.h"
00002 #include "TMDFParameters.h"
00003 #include "TMath.h"
00004 #include <assert.h>
00005 TMDFParameters* TMDFParameters::fgTMDFParameters = 0;
00006 ClassImp(TMDFParameters);
00007 //________________________________________________________________________________
00008 TArrayD *TMDFParameters::GetTerms(Double_t *x) {
00009   TArrayD *TV = new TArrayD[fNvar];
00010   for (Int_t i = 0; i < fNvar; i++) {
00011     Int_t N = fMaxPower[i];
00012     Double_t xx = x[i];
00013     if (xx < fVmin[i]) xx = fVmin[i];
00014     if (xx > fVmax[i]) xx = fVmax[i];
00015     Double_t X = 1 + 2./(fVmax[i] - fVmin[i])*(xx - fVmax[i]);
00016     Int_t nn = N+1;
00017     if (nn < 3) nn = 3;
00018     TV[i] = TArrayD(nn);
00019     Tcheb(X,N,TV[i].GetArray());
00020   }
00021   return TV;
00022 }
00023 //________________________________________________________________________________
00024 Double_t TMDFParameters::Eval(Double_t *x) {
00025   TArrayD *TV = GetTerms(x);
00026   Double_t value = fMean;
00027   for (Int_t i = 0; i < fNcoef; i++) {
00028     Double_t term = fCoef[i];
00029     Int_t p = fCode[i];
00030     for (Int_t j = fNvar - 1; j >= 0; j--) {
00031       Int_t k = p%10; assert(k >= 0);
00032       term *= TV[j][k-1]; 
00033       p /= 10;
00034     }
00035     value += term;
00036   }
00037   delete [] TV;
00038   return value;
00039 }
00040 //________________________________________________________________________________
00041 Double_t TMDFParameters::dEval(Double_t *x) {
00042   TArrayD *TV = GetTerms(x);
00043   Double_t value = 0;
00044   for (Int_t i = 0; i < fNcoef; i++) {
00045     Double_t term = fdCoef[i];
00046     Int_t p = fCode[i];
00047     for (Int_t j = fNvar - 1; j >= 0; j--) {
00048       Int_t k = p%10; assert(k >= 0);
00049       term *= TV[j][k-1]; 
00050       p /= 10;
00051     }
00052     value += term*term;
00053   }
00054   delete [] TV;
00055   return TMath::Sqrt(value);
00056 }
00057 //________________________________________________________________________________
00058 Double_t *TMDFParameters::Tcheb(Double_t x, Int_t N, Double_t *T) {
00059   T[0] = 1; T[1] = T[2] = 0;
00060   for (Int_t j = 1; j <= N; j++) {
00061     if (j == 1) T[j] = x; 
00062     else        T[j] = 2 * x * T[j-1] - T[j-2];
00063   }
00064   return &T[0];
00065 }
00066 //________________________________________________________________________________
00067 Double_t TMDFParameters::Func(Double_t *x, Double_t *p) {
00068   if (! p) return Instance()->Eval(x);
00069   TArrayD X(Instance()->Nvar(),x);
00070   Int_t j = 0;
00071   for (Int_t i = 0; i < Instance()->Nvar(); i++) {
00072     if (p[i] > -999.0) X[i] = p[i];
00073     else               X[i] = x[j++];
00074   }
00075   return Instance()->Eval(X.GetArray());
00076 }
00077 //________________________________________________________________________________
00078 TF1 *TMDFParameters::ProjectionX(Int_t code) {
00079   if (code < 0 || code >= Instance()->Nvar()) return 0;
00080   TF1 *f = new TF1(Form("Func%i",code),Func,Vmin(code),Vmax(code),4);
00081   Double_t params[4]; 
00082   for (Int_t i = 0; i < Instance()->Nvar(); i++) {
00083     params[i] = -999.;
00084     if (i != code) params[i] = Vmax(i) - (Vmean(i) - 1)*(Vmin(i) - Vmax(i))/2.;
00085   }
00086   f->SetParameters(params);
00087   return f;
00088 }
00089 //________________________________________________________________________________
00090 TF2 *TMDFParameters::ProjectionXY(Int_t code1, Int_t code2) {
00091   if (code1 == code2) return 0;
00092   if (code1 < 0 || code1 >= Instance()->Nvar()) return 0;
00093   if (code2 < 0 || code2 >= Instance()->Nvar()) return 0;
00094   TF2 *f = new TF2(Form("Func%i_%i",code1,code2),Func,Vmin(code1),Vmax(code1),Vmin(code2),Vmax(code2),4);
00095   Double_t params[4]; 
00096   for (Int_t i = 0; i < Instance()->Nvar(); i++) {
00097     params[i] = -999.;
00098     if (i != code1 && i != code2) params[i] = Vmax(i) - (Vmean(i) - 1)*(Vmin(i) - Vmax(i))/2.;
00099   }
00100   f->SetParameters(params);
00101   return f;
00102 }
00103 //________________________________________________________________________________
00104 void TMDFParameters::Print(Option_t */* option */) const {
00105   cout << "Sample statistics:" << endl
00106        << "------------------" << endl
00107        << "          ";          
00108   for (Int_t i = 0; i < fNvar; i++)
00109     cout << " " << setw(10) << i+1 << flush;
00110   cout << endl << " Max:   " << flush;
00111   for (Int_t i = 0; i < fNvar; i++)
00112     cout << " " << setw(10) << setprecision(4)
00113          << fVmax[i] << flush;
00114   cout << endl << " Min:   " << flush;
00115   for (Int_t i = 0; i < fNvar; i++)
00116     cout << " " << setw(10) << setprecision(4)
00117          << fVmin[i] << flush;
00118   cout << endl << " Mean:  " << flush;
00119   for (Int_t i = 0; i < fNvar; i++)
00120     cout << " " << setw(10) << setprecision(4)
00121          << fVmean[i] << flush;
00122   cout << endl;
00123   cout << "Coefficients:" << endl
00124        << "-------------" << endl
00125        << "   #         Value        Error   Powers" << endl
00126        << " ---------------------------------------" << endl;
00127   for (Int_t i = 0; i < fNcoef; i++) {
00128     cout << " " << setw(3) << i << "  "
00129          << setw(12) << fCoef[i] << "  "
00130          << setw(12) << fdCoef[i] << "  " << flush;
00131     Int_t p = fCode[i];
00132     TArrayI Power(fNvar);
00133     for (Int_t j = fNvar - 1; j >= 0; j--) {Power[j] = p%10; p /= 10;}
00134     for (Int_t j = 0; j < fNvar; j++) cout << " " << setw(3) << Power[j] - 1 << flush;
00135     cout << endl;
00136   }
00137 }