#include <iostream>
#include <vector>
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TTree.h"
#include "TGraphErrors.h"
#include "TCanvas.h"
void fitl(TH1F * hist, TF1 * total);
void masspt(TTree *Fitter){
  Int_t step=20;
  Int_t pos=60;
  Int_t i;
  Int_t n=8;
  Double_t ptmin[8]={1.5,  1.5,2,3,4,5,7,10};
  Double_t ptmax[8]={1000.,2.0,3,4,5,7,10,15};
  char name[80];
  char cutstring[80];
  char drawstring[80];
  char canv[80];
  Double_t mean[20];
  Double_t sigma[20];
  Double_t bkg[20];
  Double_t meane[20];
  Double_t sigmae[20];
  Double_t bkge[20];
  Double_t pt[20];
  Double_t pte[20];

  TCanvas *myCanvas;


  TH1F*jmpt=0;


  TF1 * total;

  TCanvas *tc33 = new TCanvas("MassvsPt2","Mass Vs Pt2",400, pos, 800, 600);
  Fitter->Draw("FMass:PtJpsi","PtJpsi< 20");
  pos+=step;
  TCanvas *tc = new TCanvas("MassvsPt","Mass Vs Pt",400, pos, 800, 600);
  pos+=step;
  TH2F* mvp = new TH2F("2dhist","Mass vs Pt",50, 0.,21.,50, 2.6,3.6); 
  Fitter->Draw("FMass:PtJpsi>> 2dhist","PtJpsi< 20");
  mvp->Draw("surf");
  TCanvas *tc2 = new TCanvas("MassvsPt2","Mass Vs Pt:box",400, pos, 800, 600);
  pos+=step;
  mvp->Draw("box");

  for(i=0;i<n;i++){
    sprintf(name,"%s%d","jmpt",i);
    sprintf(cutstring,"%s%lf%s%lf","PtJpsi>",ptmin[i],"&&PtJpsi<",ptmax[i]);
    jmpt=new TH1F(name,cutstring,100,2.9,3.3);
    sprintf(drawstring,"%s%s","FMass>> ",name);

    if(!(i%4)){
      sprintf(canv,"%s%d","PtRange",i/4);
      TCanvas *t = new TCanvas(canv,canv,400, pos, 800, 600);
      myCanvas=t;
      myCanvas->Divide(2,2);

      pos+=step;
    }
    myCanvas->cd(i%4+1);
    Fitter->Draw(drawstring,cutstring);
    total = new TF1("total","gaus(0)+pol0(3)");
    fitl(jmpt, total);
    Double_t par[4];
    total->GetParameters(&par[0]);
    mean[i]=par[1];
    meane[i]=total->GetParError(1);
    sigma[i]=par[2];
    sigmae[i]=total->GetParError(2);
    bkg[i]=par[3];
    bkge[i]=total->GetParError(3);
    if(i!=0){
      pt[i]=ptmin[i]+(ptmax[i]-ptmin[i])/2;
      pte[i]=ptmax[i]-ptmin[i];
    }else{
      pt[i]=20;
      pte[i]=1;
    }

  } 
  
  std::cout << "Mass"<< std::endl;
  for (i=0;i<n;i++)
    std::cout << pt[i] << " +- " << pte[i] 
	      << " :  " <<  mean[i] << " +- " << meane[i] << std::endl;
  std::cout << "Width"<< std::endl;
  for (i=0;i<n;i++)
    std::cout << pt[i] << " +- " << pte[i] << " :  " 
	      <<  sigma[i] << " +- " << sigmae[i] << std::endl;
  std::cout << "Background"<< std::endl;
  for (i=0;i<n;i++)
    std::cout << pt[i] << " +- " << pte[i] << " :  "
	      <<  bkg[i] << " +- " << bkge[i] << std::endl;
  sprintf(canv,"%s","Pt Dependencies");
  TCanvas *t = new TCanvas(canv,canv,400, pos, 800, 600);
  pos+=step;
  myCanvas=t;
  myCanvas->Divide(2,2);
  Int_t j=1;
  // Mass
  myCanvas->cd(j++);
  TGraphErrors * gMass = new TGraphErrors(i,pt,mean,pte,meane);
  gMass->SetTitle("JPsi Mass vs Pt");
  gMass->Draw("AP");
  // sigma
  myCanvas->cd(j++);
  TGraphErrors * gSig = new TGraphErrors(i,pt,sigma,pte,sigmae);
  gSig->SetTitle("JPsi Width vs Pt");
  gSig->Draw("AP");
  // Mass
  myCanvas->cd(j++);
  TGraphErrors * gBkg = new TGraphErrors(i,pt,bkg,pte,bkge);
  gBkg->SetTitle("JPsi Background vs Pt");
  gBkg->Draw("AP");
  float mass0=3.09687;
  for (i=0;i<n;i++)
    std::cout << pt[i] << " +- " << pte[i] << " :  " 
	      << "Sig: " <<  mass0-3*sigma[i] << " to  " << mass0+3*sigma[i]
	      << " Wind: " <<  mass0-7*sigma[i] << " and " << mass0+7*sigma[i]
	      << std::endl;

}
void fitl(TH1F * hist, TF1* total){
  TF1 *g1 = new TF1("g1", "gaus");
  TF1 *l1 = new TF1("l1", "pol0");
  hist->Fit("g1"); 
  hist->Fit("l1");  
  Double_t par[4];
  g1->GetParameters(&par[0]);
  l1->GetParameters(&par[3]);
  total->SetParameters(par);
  hist->Fit("total");     
}