//	"bw_send01.C"	08/13/2008

//  auto time stamp: (year/month/day time)
//	Time-stamp: <08/08/13 15:56:13 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);

TH1F* piplus;
TH1F* piminus;

Int_t verbose = 0;

const Double_t kPionMass   = .1395679;
const Double_t kMass[6] = {0.93827231,0.493646,0.1395679,0.1395679,0.493646,0.93827231};
Double_t radius;

TH1F *pip010;
TH1F *pim010;

TF2* flowfunc;

void pionpt()
{

	Double_t xAxispi[6] = { 0.38, 0.428, 0.4765, 0.5255, 0.5745, 0.6235 };

	pip010 = new TH1F("pip010", "pip010", 5, xAxispi);	//	pip010 = new TH1F("pip010", "pip010", 11, xAxispi);

    pip010->SetBinContent( 1, 89.98);	//  pip010->SetBinContent( 7, 89.98);
    pip010->SetBinContent( 2, 71.1);	//  pip010->SetBinContent( 8, 71.1);
    pip010->SetBinContent( 3, 56.1);	//  pip010->SetBinContent( 9, 56.1);
    pip010->SetBinContent( 4, 44.3);	//  pip010->SetBinContent(10, 44.3);
    pip010->SetBinContent( 5, 35.8);	//  pip010->SetBinContent(11, 35.8);

    pip010->SetBinError( 1, 0.92);
    pip010->SetBinError( 2, 1.4);
    pip010->SetBinError( 3, 1.6);
    pip010->SetBinError( 4, 1.3);
    pip010->SetBinError( 5, 1.4);

	pim010 = new TH1F("pim010", "pim010", 5, xAxispi);	//	pim010 = new TH1F("pim010", "pim010", 11, xAxispi);

    pim010->SetBinContent( 1, 90.68);	//  pim010->SetBinContent( 7, 90.68);
    pim010->SetBinContent( 2, 72.0);	//  pim010->SetBinContent( 8, 72.0);
    pim010->SetBinContent( 3, 56.7);	//  pim010->SetBinContent( 9, 56.7);
    pim010->SetBinContent( 4, 44.8);	//  pim010->SetBinContent(10, 44.8);
    pim010->SetBinContent( 5, 35.7);	//  pim010->SetBinContent(11, 35.7);

    pim010->SetBinError( 1, 0.93);
    pim010->SetBinError( 2, 1.4);
    pim010->SetBinError( 3, 1.7);
    pim010->SetBinError( 4, 1.3);
    pim010->SetBinError( 5, 1.4);

}

Double_t PipFit(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[2];


//+++++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(2);

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


Double_t PimFit(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[3];


//+++++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(3);

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

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.,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;
}


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

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

  Float_t sta_err = 1.0;
  Float_t err0;
  Float_t err1;

  radius = rad;

//  TH2F * multi = new TH2F("multi","multi",40,0,4.0,6,0,6);	// original
  TH2F * multi = new TH2F("multi", "multi", 20, 0, 2.0, 6, 0, 6);
 
	TCanvas *c1 = new TCanvas( "c1", "c1", 700, 500);
	c1->Clear();

TH1F* piplus  = 0;
TH1F* piminus = 0;

	if(rapidity == 0)
	{

	piplus = pip010;	//new TH1F(*pip010);
	piminus = pim010;	//new TH1F(*pim010);

	}




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

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

		}
	}




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

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

		}
	}


    
  TFile* chi2depend = new TFile(FileName,"RECREATE");
  flowfunc = new TF2("flowfit", func2D, 0, 2.5, 0, 6, 9);
  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 " << flowfunc->GetParameter(1) << endl;

  flowfunc->Print();

	//	(Re-scale of y-axis by normalization parameters)

	flowfunc->SetParameter(2, flowfunc->GetParameter(2) * 1.0);		//	pi+
	flowfunc->SetParameter(3, flowfunc->GetParameter(3) * 0.1);		//	pi-

  flowfunc->Print();

//	Double_t xAxispi[12] = { 0.1045, 0.1475, 0.1915, 0.2375, 0.2845, 0.332, 0.38, 0.428, 0.4765, 0.5255, 0.5745, 0.6235 };
	Double_t xAxispi[6] = { 0.38, 0.428, 0.4765, 0.5255, 0.5745, 0.6235 };


//  TH1F* pipproj = new TH1F("pipproj", "#pi^{+} projection", 11, xAxispi);
  TH1F* pipproj = new TH1F("pipproj", "#pi^{+} projection", 5, xAxispi);

//  TH1F* pimproj = new TH1F("pimproj", "#pi^{-} projection", 11, xAxispi);
  TH1F* pimproj = new TH1F("pimproj", "#pi^{-} projection", 5, xAxispi);



//  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 = 0; i <= multi->GetNbinsX(); i++)
  {
    pipproj->SetBinContent(i,multi->GetBinContent(i,3));
    pipproj->SetBinError(i,multi->GetBinError(i,3));

    pimproj->SetBinContent(i,multi->GetBinContent(i,4));
    pimproj->SetBinError(i,multi->GetBinError(i,4));

  }
 
  //Calculating the projected result of the simultaneous fit.
//  TF1* pipfit = new TF1("pipfit", PipFit, 0.0, 0.6, 7);	// original is TF1("~", ~, 0.0, 1.0, 0);

  TF1* pipfit = new TF1("pipfit", PipFit, 0.0, 1.0, 0);
  pipfit->SetLineWidth(1);
  pipproj->GetListOfFunctions()->Add(pipfit);

  TF1* pimfit = new TF1("pimfit", PimFit, 0.0, 1.0, 0);
  pimfit->SetLineWidth(1);
  pimproj->GetListOfFunctions()->Add(pimfit);

 
  pipproj->GetXaxis()->SetRangeUser(0, 1.0);
  pimproj->GetXaxis()->SetRangeUser(0, 1.0);


//	Re-scale of y-axis

	pipproj->TH1::Scale(1.0);
	pimproj->TH1::Scale(0.1);

  c2->cd();
  c2->SetFrameFillColor(0);
  gPad->SetLogy();
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);
     
//  pipproj->SetMinimum(0.000000001);
  pipproj->SetMarkerColor(2);
  pipproj->SetMarkerStyle(24);
  pipproj->SetMarkerSize(0.5);
//  pimproj->SetMinimum(0.001);
  pimproj->SetMinimum(0.00001);
  pimproj->SetMaximum(1000.0);
  pimproj->SetMarkerColor(2);
  pimproj->SetMarkerStyle(20);
  pimproj->SetMarkerSize(0.5);
  pimproj->SetTitle("Projection of multifit");
//  pipproj->SetTitle("Projection of multifit");



  pimproj->Draw();
  pipproj->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);
     
  pipproj->Write();

  pimproj->Write();

  multi->Write();

  chi2depend->Write();


}