//	"bw_send03.C"	08/14/2008

//  auto time stamp: (year/month/day time)
//	Time-stamp: <08/08/14 17:25:24 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 "TLegend.h"
#include "TString.h"
#include "TFrame.h"
#include "TFile.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 PipFit(Double_t* xp, Double_t* par);
Double_t PimFit(Double_t* xp, Double_t* par);


TH1F* piplus;
TH1F* piminus;

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 integral = 2.0;
  return integral; 
}

Double_t PimFit(Double_t* xp, Double_t* par)
{
  
  Double_t integral = 3.0;
  return integral; 
}


Double_t func2D(Double_t* xp, Double_t* pp) 
{
  // xp  array of dependent variables  
  // pp  array of parameters

  Double_t integral = 1.0;
  return integral;
}


Double_t int_func(Double_t* xx, Double_t* param)
{
	Double_t out = 10.0;
	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", 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);

	}

  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");


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

	//	(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-

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

  TH1F* pipproj = new TH1F("pipproj", "#pi^{+} projection", 5, 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, 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");
     
//  pipproj->Write();
//  pimproj->Write();
//  multi->Write();

}