%%cpp -d
#include "Riostream.h"
#include "TROOT.h"
#include "TApplication.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TSystem.h"
#include "TBrowser.h"
#include "TFile.h"
#include "TRandom.h"
#include "TMultiDimFit.h"
#include "TVectorD.h"
#include "TMath.h"

%%cpp -d
void makeData(double* x, double& d, double& e)
{
  // Make data points
  double upp[5] = { 10, 10, 10, 10,  1 };
  double low[5] = {  0,  0,  0,  0, .1 };
  for (int i = 0; i < 4; i++)
    x[i] = (upp[i] - low[i]) * gRandom->Rndm() + low[i];

  d = x[0] * TMath::Sqrt(x[1] * x[1] + x[2] * x[2] + x[3] * x[3]);

  e = gRandom->Gaus(upp[4],low[4]);
}

%%cpp -d
int CompareResults(TMultiDimFit *fit, bool doFit)
{
      //Compare results with reference run


      // the right coefficients (before fit)
   double GoodCoeffsNoFit[] = {
   -4.37056,
   43.1468,
   13.432,
   13.4632,
   13.3964,
   13.328,
   13.3016,
   13.3519,
   4.49724,
   4.63876,
   4.89036,
   -3.69982,
   -3.98618,
   -3.86195,
   4.36054,
   -4.02597,
   4.57037,
   4.69845,
   2.83819,
   -3.48855,
   -3.97612
   };

      // the right coefficients (after fit)
   double GoodCoeffs[] = {
      -4.399,
      43.15,
      13.41,
      13.49,
      13.4,
      13.23,
      13.34,
      13.29,
      4.523,
      4.659,
      4.948,
      -4.026,
      -4.045,
      -3.939,
      4.421,
      -4.006,
      4.626,
      4.378,
      3.516,
      -4.111,
      -3.823,
   };

   // Good Powers
   int GoodPower[] = {
   1,  1,  1,  1,
   2,  1,  1,  1,
   1,  1,  1,  2,
   1,  1,  2,  1,
   1,  2,  1,  1,
   2,  2,  1,  1,
   2,  1,  1,  2,
   2,  1,  2,  1,
   1,  1,  1,  3,
   1,  3,  1,  1,
   1,  1,  5,  1,
   1,  1,  2,  2,
   1,  2,  1,  2,
   1,  2,  2,  1,
   2,  1,  1,  3,
   2,  2,  1,  2,
   2,  1,  3,  1,
   2,  3,  1,  1,
   1,  2,  2,  2,
   2,  1,  2,  2,
   2,  2,  2,  1
   };

   int nc = fit->GetNCoefficients();
   int nv = fit->GetNVariables();
   const int *powers = fit->GetPowers();
   const int *pindex = fit->GetPowerIndex();
   if (nc != 21) return 1;
   const TVectorD *coeffs = fit->GetCoefficients();
   int k = 0;
   for (int i=0;i<nc;i++) {
      if (doFit) {
         if (!TMath::AreEqualRel((*coeffs)[i],GoodCoeffs[i],1e-3)) return 2;
      }
      else {
         if (TMath::Abs((*coeffs)[i] - GoodCoeffsNoFit[i]) > 5e-5) return 2;
      }
      for (int j=0;j<nv;j++) {
         if (powers[pindex[i]*nv+j] != GoodPower[k]) return 3;
         k++;
      }
   }

   // now test the result of the generated function
   gROOT->ProcessLine(".L MDF.C");

   double refMDF = (doFit) ? 43.95 : 43.98;
   // this does not work in CLing since the function is not defined
   //double x[]    = {5,5,5,5};
   //double rMDF   = MDF(x);
   //LM:  need to return the address of the result since it is casted to a long (this should not be in a tutorial !)
   std::intptr_t iret = gROOT->ProcessLine(" double xvalues[] = {5,5,5,5}; double result=MDF(xvalues); &result;");
   double rMDF = * ( (double*)iret);
   //printf("%f\n",rMDF);
   if (TMath::Abs(rMDF -refMDF) > 1e-2) return 4;
   return 0;
}

bool doFit = true;

std::cout << "*************************************************" << std::endl;
std::cout << "*             Multidimensional Fit              *" << std::endl;
std::cout << "*                                               *" << std::endl;
std::cout << "* By Christian Holm <cholm@nbi.dk> 14/10/00     *" << std::endl;
std::cout << "*************************************************" << std::endl;
std::cout << std::endl;

gRandom = new TRandom();

TFile* output = new TFile("mdf.root", "RECREATE");

int nVars       = 4;
int nData       = 500;
double x[4];

TMultiDimFit* fit = new TMultiDimFit(nVars, TMultiDimFit::kMonomials,"v");

int mPowers[]   = { 6 , 6, 6, 6 };
fit->SetMaxPowers(mPowers);
fit->SetMaxFunctions(1000);
fit->SetMaxStudy(1000);
fit->SetMaxTerms(30);
fit->SetPowerLimit(1);
fit->SetMinAngle(10);
fit->SetMaxAngle(10);
fit->SetMinRelativeError(.01);

double d;
double e;

fit->Print("p");

printf("======================================\n");

int i;
for (i = 0; i < nData ; i++) {

   // Make some data
   makeData(x,d,e);

   // Add the row to the fit object
   fit->AddRow(x,d,e);
}

fit->Print("s");

fit->MakeHistograms();

fit->FindParameterization();

fit->Print("rc");

double *xMax = new double[nVars];
double *xMin = new double[nVars];
for (i = 0; i < nVars; i++) {
   xMax[i] = (*fit->GetMaxVariables())(i);
   xMin[i] = (*fit->GetMinVariables())(i);
}

nData = fit->GetNCoefficients() * 100;
int j;

for (i = 0; i < nData ; i++) {
   // Make some data
   makeData(x,d,e);

   for (j = 0; j < nVars; j++)
      if (x[j] < xMin[j] || x[j] > xMax[j])
   break;

   // If we get through the loop above, all variables are in range
   if (j == nVars)
      // Add the row to the fit object
      fit->AddTestRow(x,d,e);
   else
      i--;
}

if (doFit)
   fit->Fit("M");

fit->Print("fc v");

fit->MakeCode();

output->Write();
output->Close();
delete output;

int compare = CompareResults(fit, doFit);
if (!compare) {
   printf("\nmultidimfit ..............................................  OK\n");
} else {
   printf("\nmultidimfit ..............................................  fails case %d\n",compare);
}

delete fit;
delete [] xMin;
delete [] xMax;
return compare;