#include <algorithm>
#include <string>
#include "TFile.h"
#include "RooRealVar.h"

using namespace RooFit;

void fitmZstarSignal(string fileName="../dataFiles/ZZ4muAnalysisTree_H160_wSelections.root",double mH=160){

  gSystem->Load("libRooFit");  
  gROOT->ProcessLine(".L HiggsCSandWidth.cc+");
  gROOT->ProcessLine(".L ParametersZZ4l.cc+");
  gROOT->ProcessLine(".L ../PDFs/HZZ4LRooPdfs.cc+");
  gROOT->ProcessLine(".L ../PDFs/RooMZstar.cc+");

  //================ Parameters for Signal PDF ===============
            //=====================================
  HiggsCSandWidth *myCSW = new HiggsCSandWidth();
  double sqrts = 7.; 
  double g_sqrts = sqrts;
  double widthHVal =  myCSW->HiggsWidth(0,mH);
  double high_M=1000;
  double windowVal = max( widthHVal, 1. );
  double lowside = 100.;
  double g_channel=1;
  if (mH >= 275){ lowside = 180.; }
  else { lowside = 100.; }
  double low_M = max( (mH - 20.*windowVal), lowside) ;
  double high_M = min( (mH + 15.*windowVal), 900.) ;
  char* name="";

  //================== measurables

  RooRealVar Z2Mass("Z2Mass","m_{Z*}",12,80);
  RooRealVar mZStar("mZStar","m'Z*",12,80);
  RooRealVar ZZMass("ZZMass","m_{ZZ}",low_M,high_M);
  ZZMass.setBins(100000,"fft") ;

  //================== P_1(m'Z*,mZZ) ==================
	//=============================

  RooRealVar mZ("mZ","m_{Z}",60,120);
  RooRealVar Gamma("Gamma","#Gamma",0,100);
  
  mZ.setVal(91.2);   mZ.setConstant(kTRUE);
  Gamma.setVal(2.5); Gamma.setConstant(kTRUE);

  RooMZstar idealmZshape("idealmZshape","idealmZshape",mZStar,ZZMass,mZ,Gamma);

  //================= G(m'Z*) ========================
 	//=============================

  RooRealVar a1("a1","a1",1,-1000,1000);
  a1.setConstant(kTRUE);
  RooRealVar a2("a2","a2",0,-1000,1000);
  a2.setConstant(kTRUE);

  RooPolynomial accept("accept","accept",mZStar,RooArgList(a1,a2));

  //================= P_1*G ===========================
	//===============================
	
  RooProdPdf mZshape("mZshape","mZshape",idealmZshape,accept);		

  //================= P_2(mZZ) ========================
 	//=============================
  sprintf( name, "CMS_zz4l_mean_sig_%i", g_channel );
  RooRealVar CMS_zz4l_mean_sig(name,"CMS_zz4l_mean_sig",0.,-10.,10.); 
  CMS_zz4l_mean_sig.setConstant(kTRUE);
  sprintf( name, "CMS_zz4l_sigma_sig_%i", g_channel );
  RooRealVar CMS_zz4l_sigma_sig(name,"CMS_zz4l_sigma_sig",3.,0.,30.);
  CMS_zz4l_sigma_sig.setConstant(kTRUE);
  sprintf( name, "CMS_zz4l_alpha_%i", g_channel );
  RooRealVar CMS_zz4l_alpha(name,"CMS_zz4l_alpha",1.,-10.,10.);
  CMS_zz4l_alpha.setConstant(kTRUE);
  sprintf( name, "CMS_zz4l_n_%i", g_channel );
  RooRealVar CMS_zz4l_n(name,"CMS_zz4l_n",2.,-10.,10.);
  CMS_zz4l_n.setConstant(kTRUE);
  RooCBShape signalCB("signalCB","signalCB",ZZMass,CMS_zz4l_mean_sig,CMS_zz4l_sigma_sig,CMS_zz4l_alpha,CMS_zz4l_n);
  
  sprintf( name, "CMS_zz4l_mean_BW_%i", g_channel );
  RooRealVar CMS_zz4l_mean_BW(name,"CMS_zz4l_mean_BW",mH,low_M,high_M);
  sprintf( name, "CMS_zz4l_gamma_BW_%i", g_channel );
  RooRealVar CMS_zz4l_gamma_BW(name,"CMS_zz4l_gamma_BW",widthHVal);
  //RooBreitWigner signalBW("signalBW", "signalBW",ZZMass,CMS_zz4l_mean_BW,CMS_zz4l_gamma_BW);
  //RooRelBW1 signalBW("signalBW", "signalBW",ZZMass,CMS_zz4l_mean_BW,CMS_zz4l_gamma_BW);
  sprintf( name, "CMS_zz4l_widthScale_%i", g_channel );
  RooRealVar CMS_zz4l_widthScale(name,"CMS_zz4Z2l_widthScale",1.);
  RooRelBWUFParam signalBW("signalBW", "signalBW",ZZMass,CMS_zz4l_mean_BW,CMS_zz4l_widthScale);
  
  CMS_zz4l_widthScale.setConstant(kTRUE);  //--> leave this floating as a nuisance parameter for the theoretical line shape
  
  //Construct BW (x) CB
  RooFFTConvPdf sig_ggH("sig_ggH","BW (X) CB",ZZMass,signalBW,signalCB, 2);
  // Buffer fraction for cyclical behavior
  
  CMS_zz4l_mean_BW.setConstant(kTRUE);
  CMS_zz4l_gamma_BW.setConstant(kTRUE);
  
  double n_CB_d = 0.;
  double alpha_CB_d = 0.;
  double mean_CB_d = 0.;
  double sigma_CB_d = 0.;
  double mean_BW_d = mH;
  ParametersZZ4l *parameters = new ParametersZZ4l( g_channel, mH, 180);//should be 1000 for above threshold, I think
  parameters->get_massDependence_BWxCB( mH, n_CB_d, alpha_CB_d, mean_CB_d, sigma_CB_d);
  
  CMS_zz4l_mean_sig.setVal( mean_CB_d );
  cout << "mean_CB_d: " << mean_CB_d << endl;
  CMS_zz4l_sigma_sig.setVal( sigma_CB_d );
  cout << "sigma_CB_d: " << sigma_CB_d << endl;
  CMS_zz4l_alpha.setVal( alpha_CB_d );
  cout << "alpha_CB_d: " << alpha_CB_d << endl;
  CMS_zz4l_n.setVal( n_CB_d );
  cout << "n_CB_d: " << n_CB_d << endl;
  CMS_zz4l_mean_BW.setVal( mean_BW_d );
  cout << "mean_BW_d: " << mean_BW_d << endl;
  CMS_zz4l_gamma_BW.setVal( widthHVal );
  cout << "widthHVal: " << widthHVal << endl;

  //============== P(m'Z*,mZZ) ======================
	//============================

  RooProdPdf totalPDF("totalPDF","totalPDF",sig_ggH,mZshape);

  //============== load data ========================
	//============================
  TFile *file = new TFile(fileName.c_str());
  char* cutString="";
  sprintf(cutString,"ZZMass<%i && ZZMass>%i && mZStar>12 && mZStar<80",high_M,low_M);
  RooDataSet data("data","data",(TTree*)file->Get("outTree"),RooArgSet(ZZMass,mZStar),cutString);

  //totalPDF.fitTo(data);

  RooPlot *plot=mZStar.frame(12,80,40);
  RooPlot *plot2=ZZMass.frame(low_M,high_M,40);
  data.plotOn(plot);
  data.plotOn(plot2);
  //sig_ggH.plotOn(plot2,LineColor(2));
  //totalPDF.plotOn(plot2);
  totalPDF.createPlotProjection(RooArgSet(mZStar),RooArgSet(ZZMass))->plotOn(plot);
  TCanvas*c = new TCanvas("c","c",800,400);
  c->Divide(2,1);
  c->cd(1);
  plot->Draw();
  c->cd(2);  
  plot2->Draw();	

}