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

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] = y
  //x[1] = r
  const double R = 10.;
  const double m = 0.938;
  const double ycm = 1.045;
  double pT = par[5];
  double beta = par[1] * pow(x[1]/R, par[2]);
  double rho = TMath::ATanH(beta);
  double mT = sqrt(pT*pT + m*m);

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

}

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

   TF2 *fpT = new TF2("fpT",fRy,-2.04,2.04,0,10,6); // getYield has 5 pars.
   fpT->SetParameters(par[0],par[1],par[2],par[3],par[4],x[0]);
   double f = fpT->Integral(-2.04,2.04,0,10);
   delete fpT;
   return f;

}


void BWModified_Proton2(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 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;
  Int_t 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.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("");

//TVirtualFitter::SetDefaultFitter("Minuit2");//default fitter is "Minuit", both give same result
//TVirtualFitter::SetMaxIterations(500000); 

TF1 *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");
    //fns->Draw("same");
    grBW->SetMarkerColor(kRed);
    grBW->SetMarkerStyle(21);
    
}