//	"bw_fit_simple02_sent.C"	08/26/2008

//  auto time stamp: (year/month/day time)
//	Time-stamp: <08/08/26 16:56:39 tsuchiya>

#include <iostream>
#include <fstream>
#include "TH1F.h"
#include "TH2F.h"
#include "TF1.h"
#include "TF2.h"
#include "TMath.h"
#include "TMinuit.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TROOT.h"
#include "TAxis.h"
#include "TGraph.h"
#include "TLatex.h"
#include "TText.h"
#include "TGraphErrors.h"
#include "TLegend.h"
#include "TString.h"
#include "TFrame.h"
#include "TFile.h"
#include "TAttLine.h"

using namespace std;

#ifdef __CINT__
#endif

Double_t int_func(Double_t* xx, Double_t* param);
Double_t func2D(Double_t* xp, Double_t* pp);

Double_t ProFit(Double_t* xp, Double_t* par);

Double_t radius;

TH1F* proton;
TH1F *proton010;

TF2* flowfunc;

Int_t verbose = 0;

const Double_t kProtonMass = 0.93827231;

//                              0         1         2         3       4        5
const Double_t kMass[6] = {0.93827231,0.493646,0.1395679,0.1395679,0.493646,0.93827231};

void ppbarpt()
{

//  proton
//	y = 0 (-0.10~0.10), centrality 0-5% Au+Au @ STAR data

	Double_t xAxispro[17] = { 0.081, 0.103, 0.1255, 0.1495, 0.1755, 0.2035, 0.233, 0.2635, 0.295, 0.328,
							0.362, 0.397, 0.4335, 0.4705, 0.508, 0.5465, 0.5855 };

	proton010 = new TH1F("proton010", "proton010", 16, xAxispro);

//	real data
    proton010->SetBinContent( 1, 6.42);
    proton010->SetBinContent( 2, 6.45);
    proton010->SetBinContent( 3, 6.18);
    proton010->SetBinContent( 4, 5.96);
    proton010->SetBinContent( 5, 5.75);
    proton010->SetBinContent( 6, 5.44);
    proton010->SetBinContent( 7, 5.07);
    proton010->SetBinContent( 8, 4.89);
    proton010->SetBinContent( 9, 4.67);
    proton010->SetBinContent(10, 4.436);
    proton010->SetBinContent(11, 4.16);
    proton010->SetBinContent(12, 3.95);
    proton010->SetBinContent(13, 3.68);
    proton010->SetBinContent(14, 3.36);
    proton010->SetBinContent(15, 3.24);
    proton010->SetBinContent(16, 2.89);

    proton010->SetBinError( 1, 0.31);
    proton010->SetBinError( 2, 0.25);
    proton010->SetBinError( 3, 0.20);
    proton010->SetBinError( 4, 0.16);
    proton010->SetBinError( 5, 0.13);
    proton010->SetBinError( 6, 0.13);
    proton010->SetBinError( 7, 0.12);
    proton010->SetBinError( 8, 0.11);
    proton010->SetBinError( 9, 0.10);
    proton010->SetBinError(10, 0.095);
    proton010->SetBinError(11, 0.13);
    proton010->SetBinError(12, 0.12);
    proton010->SetBinError(13, 0.12);
    proton010->SetBinError(14, 0.11);
    proton010->SetBinError(15, 0.14);
    proton010->SetBinError(16, 0.14);

}


//===========================================================

Double_t ProFit(Double_t* xp, Double_t* par)
{

    Double_t beta  = flowfunc->GetParameter(6); 
	Double_t temp  = flowfunc->GetParameter(7);
//  Double_t rmin = 0.0;
	Double_t alpha = flowfunc->GetParameter(8);
//  Double_t alpha = 0.5; 
	Double_t integral;
	Double_t mass = kMass[0];

//+++++07/18/2008 changed+++++
	Double_t mt = xp[0] + mass;//	xp = mt - mass, so mt = xp + mass
//+++++end changed+++++

//  Double_t mt   = sqrt(xp[0]*xp[0]+mass*mass); 
	Double_t con = flowfunc->GetParameter(0);

	TF1* func_pro = new TF1( "func_pro", int_func, 0.0, radius+2.0, 7);

//                         0  1    2    3     4     5    6  
	Double_t params[7] = {mt,con,mass,beta,radius,temp,alpha};
	func_pro->SetParameters(params); 
	integral = func_pro->Integral(0.0, radius);
	if(integral <= 0 )integral = 0.00000001;
	return integral; 

}


//______________________________________________________________________
//The simultaneous fit function. 
//

Double_t func2D(Double_t* xp, Double_t* pp) 
{

//  xp  array of dependent variables  
//  pp  array of parameters

    Double_t beta  = pp[6]; 
	Double_t temp  = pp[7];
//  Double_t rmin = 0.0;

	Double_t alpha = pp[8];
//  Double_t alpha = 0.5;
	Double_t integral;

	Int_t y = (Int_t) xp[1];

	Double_t mass = kMass[y];
	Double_t mt = xp[0] + mass;
//  Double_t mt   = sqrt(xp[0]*xp[0]+mass*mass); 

	TF1* funcx = new TF1("funcx",int_func,0.0,radius+2.0, 7);

//                       0   1     2    3     4     5    6  
  Double_t params[7] = {mt,pp[y],mass,beta,radius,temp,alpha};
  funcx->SetParameters(params); 

  integral = funcx->Integral(0.0, radius);
  if(integral <= 0.0 ) integral = 0.00000001;
  return integral;
}


Double_t int_func(Double_t* xx, Double_t* param)
{
//	Double_t arg = 0.0;
	Double_t out;

//+++++06/30/2008 changed+++++  
	Double_t rho = TMath::ATanH(param[3]*pow(xx[0]/param[4],param[6]));
  Double_t pt = sqrt(param[0]*param[0]-param[2]*param[2]);	// original
//	Double_t pt = param[0];
  Double_t mt = param[0];	// original
//	Double_t mt = sqrt(param[0]*param[0] + param[2]*param[2]);  
//+++++end changed+++++

	Double_t con = param[1]; 
	Double_t temp = param[5];

	if( (mt * TMath::CosH(rho)/temp) == 0.0 ) return 0.0;	//	avoiding BesselK1(0.0) error message!

	out =   con * xx[0] * mt				//	the original is "out = con * mt", xx[0] is r(radius)
		* TMath::BesselI0(pt * TMath::SinH(rho)/temp) 
		* TMath::BesselK1(mt * TMath::CosH(rho)/temp);
	return out;
}


//______________________________________________________________________
//The method to create simultaneous blast wave fit and corresponding
//minuit curves.

void flowfit(Int_t rapidity, Double_t rad, Char_t* FileName)
{

  Float_t sys_err = 0.05;			// original sys_err = 0.1;
  //  Float_t sys_err_proton = 0.05;	// original sys_err_proton = 0.1;

  Float_t sta_err = 1.0;
  Float_t err0;
  Float_t err1;
//  Int_t bin;

  radius = rad;

//  TH2F * multi = new TH2F("multi","multi",40,0,4.0,6,0,6);	// original
  TH2F * multi = new TH2F("multi", "multi", 18, 0, 2.0, 6, 0, 6);
  //    multi->GetXaxis()->SetRangeUser(0.0,2.5);

  //Get the different spectras and add cuts to them. 
  //     TFile* hfile  = new TFile(FileName);

//  FindCanvas("c1",600,800);
	TCanvas *c1 = new TCanvas( "c1", "c1", 700, 500);
	c1->Clear();

TH1F* proton  = 0;

	if(rapidity == 0)
	{
	    proton = proton010;	//new TH1F(*proton010);
	}


	cout << "\t<proton>" << endl;
	cout << "\tchecking TH1F * proton" << endl;
    cout << "\tproton = " << proton << endl;

	for (Int_t i=1; i <= proton->GetNbinsX(); i++)
	{
		if(proton->GetBinContent(i)<100000.0 && proton->GetBinContent(i) > 0.01)
		{
			err0 = proton->GetBinError(i)*sta_err;
			err1 = proton->GetBinContent(i)*sys_err;
			multi->SetBinContent(i,1,proton->GetBinContent(i));
			multi->SetBinError(i,1,sqrt(err0*err0 + err1*err1));

			cout << "\t" << i << "\t" << proton->GetBinContent(i) << endl;

		}
	}


  //Creating a file to write down the fits in.
    
  TFile* chi2depend = new TFile(FileName,"RECREATE");
  flowfunc = new TF2("flowfunc", func2D, 0, 2.5, 0, 6, 9);    // already decrear "TF2* flowfunc;" on the top of the code
                                                              // do not write "TF2* flowfunc = new TF2(..." again, error messages...

  flowfunc->SetParNames( "Norm p^{+}", "Norm K^{+}", "Norm#pi^{+}", "Norm#pi^{-}", "Norm k{-}", "Norm #bar{p}", "#beta", "T", "alpha" );
  flowfunc->SetParameters( 32427, 3754, 4337, 4150, 2618, 7355, 0.8, 0.12, 0.5 );

//	The last parameter is alpha, 0.5 or 1.0 ...

  if(rapidity == 0) flowfunc->SetParameters(32427,3754,4337,4150,2618,7355,0.8,0.10,0.8);	// original

  flowfunc->SetParLimits(6,0.5,0.95);	//	parameter 6 is "beta" original

  flowfunc->SetParLimits(7, 0.01, 0.9);	//	parameter 7 is "Temp"

  flowfunc->SetParLimits(8,0.1,1.3);	//	original

  multi->Fit("flowfunc");
  multi->Fit("flowfunc","+");

  cout << "Checking parameters(0) " << flowfunc->GetParameter(0) << endl;
  cout << "Checking parameters(1) " << flowfunc->GetParameter(1) << endl;
  cout << "Checking parameters(2) " << flowfunc->GetParameter(2) << endl;
  cout << "Checking parameters(3) " << flowfunc->GetParameter(3) << endl;
  cout << "Checking parameters(4) " << flowfunc->GetParameter(4) << endl;
  cout << "Checking parameters(5) " << flowfunc->GetParameter(5) << endl;
  cout << "Checking parameters(6) " << flowfunc->GetParameter(6) << endl;
  cout << "Checking parameters(7) " << flowfunc->GetParameter(7) << endl;
  cout << "Checking parameters(8) " << flowfunc->GetParameter(8) << endl;

  flowfunc->Print();

	//	(Re-scale of y-axis by normalization parameters)
	flowfunc->SetParameter(0, flowfunc->GetParameter(0) * 0.0001);	//	proton

  flowfunc->Print();

  //=====08/26/2008=====
  //  multi->Draw(); ok, but no change!
  //=====
 
	Double_t xAxispro[17] = { 0.081, 0.103, 0.1255, 0.1495, 0.1755, 0.2035, 0.233, 0.2635, 0.295, 0.328,
							0.362, 0.397, 0.4335, 0.4705, 0.508, 0.5465, 0.5855 };

  TH1F* proproj = new TH1F("proproj", "pro projection",     16, xAxispro);

//  FindCanvas("c2",500,500);
	TCanvas *c2 = new TCanvas( "c2", "c2", 500, 500);
	c2->Clear();


  //Finding the spectras used in the fit.

  for(Int_t i = 1; i <= 16; i++)
  {
    proproj->SetBinContent(i,multi->GetBinContent(i,1));
	cout << "\t" << i << "\t" << proproj->GetBinContent(i) << endl;

    proproj->SetBinError(i,multi->GetBinError(i,1));
  }


  //Calculating the projected result of the simultaneous fit.

  TF1* profit = new TF1("profit", ProFit, 0.0, 1.0, 9);
  profit->SetLineWidth(1);
  proproj->GetListOfFunctions()->Add(profit);

  //  proproj->Fit(profit);     does not work!
  //  proproj->Fit(profit,"+"); does not work!

//	Re-scale of y-axis

	proproj->TH1::Scale(0.0001);

  c2->cd();
  c2->SetFrameFillColor(0);
  gPad->SetLogy();
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);

//  gPad->DrawFrame(xmin,ymin,xmax,ymax);
  gPad->DrawFrame(0.0, 0.00001, 1.0, 0.001);
  gPad->Update();

  proproj->SetMarkerColor(4);
  proproj->SetMarkerStyle(25);
  proproj->SetMarkerSize(0.5);

  proproj->Draw("same");
  
  Float_t T         = flowfunc->GetParameter(7);
  Float_t T_err     = flowfunc->GetParError(7);
  Float_t chi2      = flowfunc->GetChisquare();
//  Float_t ndof      = flowfunc->GetNDF();
  Int_t ndof      = flowfunc->GetNDF();
  Float_t beta      = flowfunc->GetParameter(6);
  Float_t beta_err  = flowfunc->GetParError(6);
  Float_t prof      = flowfunc->GetParameter(8);
  Float_t prof_err  = flowfunc->GetParError(8);
          
  TLatex l;
  Char_t fittext[50];
  l.SetNDC();
  l.SetTextColor(2);

  sprintf(fittext, "#beta_{s} = %3.3f #pm %2.3f", beta, beta_err);
  l.DrawLatex(0.4, 0.85, fittext);
  sprintf(fittext, "Temp = %3.3f #pm %2.3f", T, T_err);
  l.DrawLatex(0.4, 0.8, fittext);
  sprintf(fittext, "#alpha = %3.3f #pm %2.3f", prof, prof_err);
  l.DrawLatex(0.4, 0.75, fittext);
  sprintf(fittext, "#chi^{2}/NDF = %3.1f / %d", chi2, ndof);
  l.DrawLatex(0.4, 0.70, fittext);

  //=====08/26/2008=====
  profit->Write();
  //=====

  proproj->Write();

  multi->Write();

  chi2depend->Write();

}