#include <Riostream>
#include <TROOT>
#include <TApplication>
#include <TCanvas>
#include <TH1>
#include <TSystem>
#include <TBrowser>
#include <TFile>
#include <TRandom>
#include <TMultiDimFit>
#include <TVectorD>
#include <TChain>
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>

void multifit3() {


//   ////open new file////
//   TFile *output = new TFile("output.root","recreate");


  ////prepare datasample////
  TChain chainTP1("T");
  chainTP1.Add("wholedataTP1E4.root");

  const Double_t pi = 3.14159265; //define the number PI
  const Double_t tp1 = 0.914;

  Double_t en, sz, xc, yc, ximp, yimp, el, id;

  Double_t dist; //distance from center of field of view
  Double_t logsz; //logarithm of size to the base 10
  Double_t logen; //logarith of energy*1000 to the base 10
  //Double_t imp; //impact distance
  Double_t elrad; //elevation in radians
  Double_t sinel; //sinus of elevation
  Double_t logsinel;
  Double_t x[3];

  chainTP1.SetBranchAddress("id", &id);
  chainTP1.SetBranchAddress("en", &en);
  chainTP1.SetBranchAddress("sz", &sz);
  chainTP1.SetBranchAddress("xc", &xc);
  chainTP1.SetBranchAddress("yc", &yc);
  chainTP1.SetBranchAddress("ximp", &ximp);
  chainTP1.SetBranchAddress("yimp", &yimp);
  chainTP1.SetBranchAddress("el", &el);


  ////create multidim fitter////
  TMultiDimFit *multifit = new TMultiDimFit(3,TMultiDimFit::kMonomials,"");

  Int_t mPowers[] = {5,2,1};
  multifit->SetMaxPowers(mPowers);
  multifit->SetMaxFunctions(500);
  multifit->SetMaxStudy(500);
  multifit->SetMaxTerms(15);
  multifit->SetPowerLimit(3);
  //multifit->SetMinAngle(45);
  //multifit->SetMaxAngle(5);
  multifit->SetMinRelativeError(0.1);

  Double_t d;

  ////fitting////
  Int_t nEntriesTP1 = (Int_t)chainTP1.GetEntries();
  for (Int_t i=0; i<nEntriesTP1; i++) {

    chainTP1.GetEntry(i);

    dist = sqrt(xc*xc + yc*yc);
    logsz = log10(sz/tp1);
    elrad = (el/180.)*pi;
    sinel = sin(elrad);
    logsinel = log10(sinel);
    logen = log10(1000*en);

    x[0] = dist;
    x[1] = logsz;
    x[2] = logsinel;

    d = logen;

    multifit->AddRow(x,d);

  }

  multifit->Print("s");
  multifit->MakeHistograms();
  multifit->FindParameterization("");
  multifit->Print("rc");

  //TBrowser browser = new TBrowser("browser","browsing multifit");
  //multifit->Browse(browser);

  multifit->MakeCode();

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

  //delete output;
  //delete multifit;

}