// Plot the TF3 Integral function used for the fit of 3 GeV data
// Fitting one particle species.
// root -x 'BWModified_Proton1.C("proton_pt_data_3Gev.dat")'  

//#include "TF23.h"

// class projectiong a TF3(x,y,z)  in a TF2 (x,z)

struct TF3Projector { 
   TF3Projector(TF3 * f, double y) : fFunc(f), fX(), fP()  { fX[0] = y;}
   double operator() (const double *x, const double *p) { 
      fX[1] = x[0];
      fX[2] = x[1]; 
      fP[0] = p[0];
      fP[1] = p[1];
      fP[2] = p[2];
      fP[3] = p[3];
      fP[4] = p[4];
      return  (*fFunc)( fX,fP );
   }
   TF3 * fFunc; 
   double fX[3];
   double fP[5];
};

double Gaussian(double y, double sigma,double mean){

 double f1 = pow((y-mean),2);
 double f2 = 2*pow(sigma,2);
 double f3 = -f1/f2;
 
 return TMath::Exp(f3);

}

Double_t fRy(Double_t *x,Double_t *par){

  //par[0] = Tkin
  //par[1] = betaS
  //par[2] = n
  //par[3] = sigma
  //par[4] = Normalization
  //x[0] = pT
  //x[1] = y
  //x[2] = r
  const double R = 10.;
  const double m = 0.938;
  const double ycm = 1.045;
  double beta = par[1] * pow(x[2]/R, par[2]);
  double rho = TMath::ATanH(beta);
  double mT = sqrt(x[0]*x[0] + m*m);

  return par[4] * (mT * x[2] * TMath::BesselI0(x[0] * TMath::SinH(rho) / par[0]) * TMath::CosH(x[1]) * TMath::Exp(-mT * TMath::CosH(x[1]) * TMath::CosH(rho) / par[0]) * Gaussian(x[1],par[3],ycm));

}

double getYield(Double_t *x,Double_t *par){

   TF3 * f3 = new TF3("fRy",fRy);
   /*
   TF3Projector * p = new TF3Projector(f3,par[0]);
   TF2 *fpT = new TF2("fpT",p,0.0,10.0,-2.04,2.04,4);//fpT should have only 4 pars. after integration.
   fpT->SetParameter(0,par[1]);
   fpT->SetParameter(1,par[2]);
   fpT->SetParameter(2,par[3]);
   fpT->SetParameter(3,par[4]);
   return  fpT->Integral(0.,10.,-2.04,2.04);
   */
   TF3Projector * p = new TF3Projector(f3,x[0]);
   TF2 *fpT = new TF2("fpT",p,-2.04,2.04,0,10,5);//fpT should have only 5 pars. after integration.
   fpT->SetParameter(0,par[0]);
   fpT->SetParameter(1,par[1]);
   fpT->SetParameter(2,par[2]);
   fpT->SetParameter(3,par[3]);
   fpT->SetParameter(4,par[4]);   
   double f = fpT->Integral(-2.04,2.04,0,10);
   delete fpT;
   delete p;
   delete f3;
   return f;

}


void BWModified_Proton1(Char_t *inFile){

ifstream mFile;

mFile.open(inFile);

  if (mFile.is_open() == 0)
    {
      cout << "  file not found! " << endl;
      return 0;
    } 

Char_t buf[256];


const unsigned int MaxNBins = 100;

Double_t x[MaxNBins];
Double_t y[MaxNBins];
Double_t ey[MaxNBins];

Int_t Bins;
Int_t j=0;

while (!mFile.eof()){
         
         mFile >> x[j] >> y[j] >> ey[j] ; 
         cout << x[j] << " " <<y[j] << " " << ey[j] << endl;
         j++;
  }
  cout << " Number of Bins for Particle"<< " is " << j << endl;
  Bins = j;

TGraphErrors *grBW = new TGraphErrors(Bins,x,y,0,ey);

TCanvas *c1 = new TCanvas("c1","K^{+} Fitting",10,10,1200,1800);
c1->GetFrame()->SetFillColor(21);
c1->GetFrame()->SetBorderSize(12); 

//TH2F *fhist = new TH2F("fhist","",50,0.15,0.7,200,0.5e1,1e2);
TH2F *fhist = new TH2F("fhist","",50,0.6,1.3,200,4,25);
gStyle->SetOptStat(0);
  fhist->SetTitle("Blast Wave Fit for p;p_{t} (GeV/c);#frac{1}{2#pi N_{evt}} #frac{1}{p_{T}} #frac{d^{2}N}{dp_{T}dy} (GeV/c)^{-2};");
  fhist->GetXaxis()->CenterTitle();
  fhist->GetYaxis()->CenterTitle();

  c1->SetLogy();

fhist->Draw("");

Double_t par0[9] = {0.09,0.1,0.1,0.1,0.1, 0.1,0.1,0.1, 0.1}; // Temperature
Double_t par1[9] = {0.8,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5}; // Beta_s
Double_t par2[9] = {1.0,1, 1,1.0,1,1,1,1,1}; // n
Double_t par3[9] = {1e5,1e4,1e5,1e6,1e5,1e5,1e5,1e5,1e5};  // Normalization constant


TVirtualFitter::SetDefaultFitter("Minuit");
TVirtualFitter::SetMaxIterations(500000); 

TF1 *fns = NULL;          //initialize
/*
fns = new TF1("fns",getYield,.1,.5,4);
    fns->SetParNames("T","beta_s" , "n" , "norm");
    fns->SetParameters(par0[0],par1[0],par2[0],par3[0]);
*/
fns = new TF1("fns",getYield,.1,2,5);
    fns->SetParNames("T","beta_s" , "n" , "sigma", "norm");
    fns->SetParameters(0.073,0.63,0.5,0.5,1.1e7);
    
    fns->SetLineColor(kBlue);
    fns->SetParLimits(2,0,2);
    fns->FixParameter(3,0.5);
    
    grBW->Fit("fns","S");
    grBW->Draw("P");
    grBW->SetMarkerColor(kRed);
    grBW->SetMarkerStyle(21);  
}