#include "qvalue_fitter.h"
#include "RooTrace.h" 

using namespace std;
using namespace RooFit ;

Qvalue_Fitter::Qvalue_Fitter()
{
  TString name = "signalfraction";
  signalfraction = new RooRealVar(name,name,0.36,0,1); //sung had fixed signal fraction to 0.36 for his method.
  for(int i = 0 ; i < 2; i++){
    Qvalue[i] = -9.;
    Qvalue_error[i] = -9.;
    seateventvalue[i] = -9.;
    fitrange[i] = -9.;
  }
  Varnames = new vector <TString>  ;
  vartype = new vector <char>  ;
  varibles = new vector<RooRealVar* >  ;
  Nvars = 1;NumNN = -1;
  Varnames->push_back(name);
  vartype->push_back('0');
  //vartype->push_back('1');
  varibles->push_back(signalfraction);
}

Qvalue_Fitter::~Qvalue_Fitter()
{
  
  varibles->clear();
  Varnames->clear();
  vartype->clear();
  
}

void Qvalue_Fitter::SetBackgroundVariable(RooRealVar *var)
{
  if(var!=NULL){
    
    TString name = var->GetName();
    Varnames->push_back(name);
    vartype->push_back('1');     // vartype = 1 -> background variable
    varibles->push_back(var);
    Nvars++;
  }
  
}


void Qvalue_Fitter::SetBackgroundPdf(RooAbsPdf *var)
{
  backgroundpdf = var;    
}

void Qvalue_Fitter::SetSignalVariable(RooRealVar *var)
{
  TString name = var->GetName();
  Varnames->push_back(name);
  vartype->push_back('2');    // vartype = 2 -> signal variable
  varibles->push_back(var);
  Nvars++;
  
}


void Qvalue_Fitter::SetSignalPdf(RooAbsPdf *var)
{
  signalpdf = var;    
  
}

void Qvalue_Fitter::SetSeatEventValue(float value1, float value2 )
{
  seateventvalue[0] = value1;
  seateventvalue[1] = value2; 
  
}


float Qvalue_Fitter::GetQvalue(int index,float *error)
{
  if(error != NULL) *error =  Qvalue_error[index];
  return Qvalue[index];
}

void Qvalue_Fitter::SetDataSet(RooRealVar *DataVar, RooDataSet *DataSet, int NumofNN)
{
  DataVariable = DataVar;
  DataSet_p = DataSet;
  NumNN = NumofNN;
}


void Qvalue_Fitter::SetFitRange(float minrange, float maxrange)
{
  fitrange[0] = minrange;
  fitrange[1] = maxrange;
  
 }

//void Qvalue_Fitter::Fit(float *chi2 ,TCanvas *C1, float *SigFrac, int *NDF, int hist_Binning, int *SigFlag, double *signalFitmean, double *signalFitgamma, double *signalFitsigma, double *Max_total)

//void Qvalue_Fitter::Fit(float *chi2 ,TCanvas *C1, float *SigFrac, int *NDF, int hist_Binning, int *SigFlag, double *signalFitmean, double *signalFitsigma, double *Max_total, double *c1, double *c2) 
//void Qvalue_Fitter::Fit(float *chi2 ,TCanvas *C1, float *SigFrac, int *NDF, int hist_Binning, int *SigFlag, double *signalFitmean, double *signalFitsigma, double *Max_total, double *arguspar)

void Qvalue_Fitter::Fit(float *chi2 ,TCanvas *C1, float *SigFrac, int *NDF, int hist_Binning, int OmegaFlag, double *par1, double *par2, double *par3, float *Qvalue_error_total, int flag)//here C0 is the first coefficient of Chebychev polynomial
{
  float b,s;
  float sf,spdf,bpdf;
  //  double b_max,s_max;
  //double bRatio_max[2];
  if(seateventvalue[0]< 0|| fitrange[0] < 0 || fitrange[1]  < 0  ){
    cerr << " seateventvalue or fitrange is not set. Exiting." << endl;
    exit(1);    
  }
  
  //   if(FitResult != NULL) FitResult->Delete();

  /////////////////////////////////////////////////////////////////////////////////////////
  // totalpdf = signal + background 
  // [Q factor = signal/(signal + background)] at the pion mass of the seed event
  // [signal = signalpdf * weight which we are calling signal fraction] at the pion mass
  //for the seed event
  // [background = backgroundpdf * (1-signal fraction)] at the pion mass for the seed event
  //////////////////////////////////////////////////////////////////////////////////////////

  Int_t numfloatingpars; 
  RooArgList pdfList(*signalpdf ,*backgroundpdf );
  RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); //to prevent messages about DEBUG, INFO, PROGRESS 
  //and WARNING from getting printed in the logfile (to make the logfile smaller and the running time faster)


  //cout<<"  totalpdf = signalpdf*signalfraction + backgroundpdf*(1-signalfraction) "<<endl<<endl;
  RooAbsPdf* totalpdf;
  totalpdf = new RooAddPdf("totalpdf","totalpdf",pdfList,*signalfraction);

  //Tracing memory

  //RooTrace::active(1);
  //RooTrace::verbose(true);
  //cerr<<endl<<endl<<"this dump is right before the fitto statement in the qfitter code."<<endl;
  //RooTrace::dump();  
  //RooTrace::mark();
 RooFitResult *FitResult = totalpdf->fitTo(*DataSet_p,Range(fitrange[0],fitrange[1]), PrintLevel(-1),Verbose(false),Save());
  /*
  RooAbsReal *nll = totalpdf.createNLL(*DataSet_p,Range(fitrange[0],fitrange[1]), PrintLevel(-1),Verbose(false),Save());
  RooMinimizer m(*nll);
  m.migrad();
  m.hesse();
  delete nll;
  */
 // cerr<<"This dump is right after the fitto statement in qfitter code"<<endl;
 //RooTrace::dump();

  RooArgSet* value =  totalpdf->getObservables(*DataVariable);
  Qvalue[0] = 0.0;   
    
    if(seateventvalue[0] > fitrange[0] && seateventvalue[0] < fitrange[1]){
      
      value->setRealValue(DataVariable->GetName(),seateventvalue[0]);
       sf = signalfraction->getVal();
      *SigFrac = sf;
      spdf = signalpdf->getVal(value);
      bpdf = backgroundpdf->getVal(value);
      b = (1.- signalfraction->getVal()) * backgroundpdf->getVal(value);
      s =  signalfraction->getVal() * signalpdf->getVal(value);

      Qvalue[0] = s /(s+b);
      if(s < 1e-4) Qvalue[0] = s;
      if(TMath::IsNaN(Qvalue[0])) Qvalue[0] = -0.5;

      //if(event>0 && event<20){
      
      //cout<<"  value : ";
      //value->Print("1");
      
      //cout<<"    Qvalue [s/(s+b)]               : "<<Qvalue[0]<<endl;
      
    }
 
    // if(event>0 && event<20){
    //cout<<" ___________________________________________________ "<<endl<<endl;
    //FitResult->Print();
    //cout<<endl<<"___________________________________________________ "<<endl<<endl;
    //}

  //Area under the distributions plotted wrt the datavariable   

  //RooAbsReal* intbckgrd = backgroundpdf->createIntegral(*DataVariable);
  //RooAbsReal* intsignal = signalpdf->createIntegral(*DataVariable);
  RooAbsReal* inttotal = totalpdf->createIntegral(*DataVariable);
  
  // // // // // // // // // // // // // // // // // // // // // // // //   d
  // // Setting Qvalue error
  // // // // // // // // // // // // // // // // // // // // // // // // // // // 


  // double corrmatrix[Nvars][Nvars]
  double sigfrac, dQ_dCs, dQ_dB,dQ_dS, derivatives[Nvars], errorlocal;
  // RooDerivative roodervar[Nvars];
  double chi2_local;
  TString name;
   
    Qvalue_error[0] = 0.;

    if(seateventvalue[0] > fitrange[0] && seateventvalue[0] < fitrange[1]){

      value->setRealValue(DataVariable->GetName(),seateventvalue[0]);
      sigfrac = signalfraction->getVal();
      
      //Changing the definition of s and b here (basically this time s and b don't have the signal fraction included in their def.)

      s=signalpdf->getVal(value);//This is the value of signalfunction [normalized to total no. of events in the mass range i.e. inttotal->getVal()] at the seedmass         
                                 //When we call the fitTo function, the normalization of totalpdf is set to the total # events in the mass plot. 
                                 //To know more about normalization of pdfs and the usage of getVal(): http://roofit.sourceforge.net/docs/classref/examples/intro9.cc.html

      b =  backgroundpdf->getVal(value);//This is the value of background function [normalized to total no. of events in the mass range i.e. inttotal->getVal()] at the seedmass

      dQ_dCs = b * Qvalue[0] * Qvalue[0] / ( sigfrac * sigfrac * s  );//will be used to calculate dQ/d(sigfrac)                                                                   
      dQ_dS = (1. - sigfrac)* b * Qvalue[0] * Qvalue[0]/ ( sigfrac * s * s );//will be used to calculate dQ/d(sigma)and dQ/d(Gaus_mean)                                            
      dQ_dB = -1. * (1. - sigfrac) * Qvalue[0] * Qvalue[0] / (sigfrac * s );//will be used to calculate dQ/d(c1 or c2 or c3)       

      //****************************************************************************************************************************************
      //Verifying the values of s and b by scaling them properly and then comparing them with the values observed in corresponding mass histograms -
     
      //cerr<<"**********This is the qvalue code -"<<endl;
      // cerr<<"The height of the signal function at the seed mass as seen in the plot = "<<(signalfraction->getVal())*(signalpdf->getVal(value))*(inttotal->getVal())*(450/hist_Binning)<<endl;
      //450/hist_binning = bin_width (in my test plot). Change the no. 450 to appropriate no. for different plots.

      //cerr<<"The height of the back function  at the seed mass as seen in the plot = "<<(1-signalfraction->getVal())*(backgroundpdf->getVal(value))*(inttotal->getVal())*(450/hist_Binning)<<endl;
      //cerr<<"The height of the totalpdf at the seed mass as seen in the plot = "<<(totalpdf->getVal(value))*(inttotal->getVal())*(450/hist_Binning)<<endl;
      //****************************************************************************************************************************************
      /*
      cerr<<"total # of events in the mass range = "<<inttotal->getVal()<<endl;

      cerr<<"sig_fraction = "<<signalfraction->getVal()<<endl;
      cerr<<"normalized value of signalpdf at the seedevent's mass (normalized to mass variable) = "<<signalpdf->getVal(value)<<endl;//pdf normalised wrt mass variable (i.e. total # events in the mass plot)

      cerr<<"total no. of events under signal function = intsignal->getVal()"<<intsignal->getVal()<<endl;
      cerr<<"total no. of events in the mass range = intttotal->getVal() ="<<inttotal->getVal()<<endl;
      cerr<<"normalized value of backgrdpdf  at the seedevent's mass (normalized to mass variable)= "<<backgroundpdf->getVal(value)<<endl;
      cerr<<"total no. of events under background function = intbckgrd->getVal() = "<<intbckgrd->getVal()<<endl;
      cerr<<"total no. of events in the mass range = intttotal->getVal() ="<<inttotal->getVal()<<endl;
      */
     
      // floatParsFinal :  Return list of floating parameters after fit                                                     
      /*

      RooArgList variablelistforcovarinacematrix(FitResult->floatParsFinal());
      int varcounter = 0;
      */

      int varmap = 0;
    
      //////////////////////////////////////////////
      //identifying signal and background parameters 
      //////////////////////////////////////////////

      //cerr<<"*******identifying signal and background*********"<<endl;
      //cerr<<"# of fit parameters = "<<FitResult->floatParsFinal().getSize()<<endl<<endl;           

      for(int Vs = 0; Vs < FitResult->floatParsFinal().getSize(); Vs++){
	    name = "variable";
	    name += Vs;
	    varmap = -9;

	    //In this loop we are finding out whether the parameter i is measured or not

	    for(int Vsl = 0; Vsl < Nvars; Vsl++){
	      if( *(FitResult->floatParsFinal().at(Vs)) == *(varibles->at(Vsl))) {
		     varmap = Vsl;
	             break;
	        }
	    }

	   if(varmap == -9){
	         cerr << "Variable not found." << endl;
	         exit(1);
	    }
	
	   if(( (varibles->at(varmap))->isConstant())) exit(1);
	 
	   //identifying signal and background parameters

	   if(vartype->at(varmap) == '1') {

	     if(OmegaFlag == 0){
 	       if(varmap==1) *par3 = varibles->at(varmap)->getVal();//C1 variable of Chebychev
	     }

	     RooDerivative roodervar(name, name,*backgroundpdf,*(varibles->at(varmap)) );

	     derivatives[Vs] = (roodervar.getVal()/inttotal->getVal()) * dQ_dB;
             //(roodervar.getVal()/inttotal->getVal()) : value of the derivative of the signalpdf [normalized to total # of events in the mass range] w.r.t Vs.
             //getVal() returns the value of the function after automatically substituting the values of the parameters that the qvalue fit returns and the seedevnt's mass.
	     //The usage of getVal(): http://roofit.sourceforge.net/docs/classref/examples/intro9.cc.html

	     //cerr<<"(FitResult->floatParsFinal().at(Vs))->GetName() = "<<(FitResult->floatParsFinal().at(Vs))->GetName()<<endl;
	     //cerr<<"(varibles->at(varmap))->GetName()="<<(varibles->at(varmap))->GetName()<<endl;
	     //cerr <<"    "<< (varibles->at(varmap))->GetName() << " is a background variable " <<endl
	     } 

	   else if (vartype->at(varmap) == '2') {
	       /*
	        if (varmap==1){ 
		  *par1 = varibles->at(varmap)->getVal(); //signalFitMean (fot both omega as well as 2pion analysis)
		}

		if (varmap==2){
                  *par2 = varibles->at(varmap)->getVal();//signalFitSigma for 2pion analysis
                                                         //signalFitGamma for omega analysis
		}

		if(OmegaFlag == 1){
		  if (varmap==3)*par3 = varibles->at(varmap)->getVal();//signalFitSigma for omega analysis
		}
	       */ 
	     RooDerivative roodervar(name, name,*signalpdf,*(varibles->at(varmap)) );
	     derivatives[Vs] = (roodervar.getVal()/inttotal->getVal()) * dQ_dS;      
	       
	   }

	  else if (vartype->at(varmap) == '0'){
	    //cerr<<"vartype = "<<vartype->at(varmap)<<endl<<"Vs = "<<Vs<<"  (signal fraction)"<<endl;
	        derivatives[Vs] = dQ_dCs;
	   }
	   //varcounter++;
      }   

      numfloatingpars= FitResult->floatParsFinal().getSize();
      TMatrixDSym covarmatrix( (FitResult->covarianceMatrix()));

      /*
          cout<<" ___________________________________________________ "<<endl<<endl;
          cout<<" covarmatrix : "<<endl<<endl;
          cout<<" [0]["<<variablelistforcovarinacematrix[0].GetName()<<"]   [1]["<<variablelistforcovarinacematrix[1].GetName()<<"]   [2]["<<variablelistforcovarinacematrix[2].GetName()<<"]   [3]["<<variablelistforcovarinacematrix[3].GetName()<<"]  [4]["<<variablelistforcovarinacematrix[4].GetName()<<"]   [5]["<<variablelistforcovarinacematrix[5].GetName()<<"]"<<endl;
          cout<<" (0,0)("<<covarmatrix(0,0)<<")  (0,1)("<<covarmatrix(0,1)<<")  (0,2)("<<covarmatrix(0,2)<<")  (0,3)("<<covarmatrix(0,3)<<")  (0,4)("<<covarmatrix(0,4)<<")  (0,5)("<<covarmatrix(0,5)<<")"<<endl;
          cout<<" (1,0)("<<covarmatrix(1,0)<<")  (1,1)("<<covarmatrix(1,1)<<")  (1,2)("<<covarmatrix(1,2)<<")  (1,3)("<<covarmatrix(1,3)<<")  (1,4)("<<covarmatrix(1,4)<<")  (1,5)("<<covarmatrix(1,5)<<")"<<endl;
          cout<<"."<<endl<<endl<<"."<<endl<<endl;
          cout<<" ___________________________________________________ "<<endl<<endl;
      */


   
      // for(int v1 = 0 ; v1 < varcounter; v1++) {
      //  for(int v2 = 0 ; v2 < varcounter; v2++) {
      for(int v1 = 0 ; v1 < FitResult->floatParsFinal().getSize(); v1++) {                                                                               
         for(int v2 = 0 ; v2 < FitResult->floatParsFinal().getSize(); v2++) {   
	         errorlocal = derivatives[v1] *covarmatrix(v1,v2) * derivatives[v2];
	         Qvalue_error[0] += errorlocal;

	   }
      }

      Qvalue_error[0] = TMath::Sqrt(Qvalue_error[0]);
      if(TMath::IsNaN(Qvalue_error[0])) Qvalue_error[0] = -0.5;
      *Qvalue_error_total = *Qvalue_error_total + Qvalue_error[0]; 
    }
    

  int RooPlot_Binning = hist_Binning;
  int NDF_local;
  RooPlot* backgroundFrame = DataVariable->frame(); //define a frame with DataVariable on the X axis
  DataSet_p->plotOn(backgroundFrame,Name("data"),Binning(RooPlot_Binning)); 
  totalpdf->plotOn(backgroundFrame,LineColor(kBlue),Name("totalpdf")); //plot DataSet_p and totalpdf on this frame

  // Method RooPlot::chiSquare() calculated the reduced chi^2
  chi2_local = backgroundFrame->chiSquare("totalpdf","data",numfloatingpars);
  //NDF_local = RooPlot_Binning + numfloatingpars -1; // no. of degrees of freedom 
  NDF_local = RooPlot_Binning - numfloatingpars; // no. of degrees of freedom 
  *NDF = NDF_local;
  if(chi2 != NULL) *chi2 = chi2_local;
  
  int Yaxis_Max = backgroundFrame->GetMaximum();
  if(C1 != NULL && flag>0){
    totalpdf->plotOn(backgroundFrame,Name("totalpdf"),Components(*backgroundpdf),LineStyle(kDashed));
    totalpdf->plotOn(backgroundFrame,Components(*signalpdf),LineColor(kRed));
    totalpdf->paramOn(backgroundFrame,Layout(0.55));
    stringstream plotinfo1;
    stringstream plotinfo2;
    stringstream plotChi2;
    stringstream plotinfo3;
    stringstream plotinfo4;
    stringstream plotinfo5;
    stringstream plotinfo6;
    stringstream plotinfo7;
    plotinfo1.str("");
    plotinfo2.str("");
    plotChi2.str("");
    plotinfo3.str("");
    plotinfo4.str("");
    plotinfo5.str("");
    plotinfo6.str("");
    plotinfo7.str("");
    plotinfo1 << " QV = " << setprecision(3) << Qvalue[0] << " E = " << setprecision(3) << Qvalue_error[0] ;
    plotinfo2 << " seed = " << setprecision(4) <<seateventvalue[0] ;
    plotChi2  << " chi2 = " << setprecision(3) << chi2_local;
    //plotinfo3 << " total_max = " << setprecision(4) << *Max_total;
    plotinfo4 << " sig-pdf = " << setprecision(4) << spdf ;
    plotinfo5 << " back-pdf = " << setprecision(4) << bpdf ;
    plotinfo6 << " background = " << setprecision(4) << b ;
    plotinfo7 << " signal = " << setprecision(4) << s ;

    if(OmegaFlag == 1 ){
      TText *textbox1 = new TText(675,Yaxis_Max*0.90,plotinfo1.str().c_str());//675 since the X axis starts at 650                                               
      TText *textbox2 = new TText(675,Yaxis_Max*0.85,plotinfo2.str().c_str());
      TText *textbox3 = new TText(675,Yaxis_Max*0.80,plotChi2.str().c_str());

      textbox1->SetTextSize(0.03);
      textbox2->SetTextSize(0.03);
      textbox3->SetTextSize(0.03);

      C1->cd();
      backgroundFrame->Draw();
      textbox1->Draw();
      textbox2->Draw();
      textbox3->Draw();
    }
    else if(OmegaFlag == 0){
      TText *textbox_1 = new TText(10,Yaxis_Max*0.90,plotinfo1.str().c_str());
      TText *textbox_2 = new TText(10,Yaxis_Max*0.85,plotinfo2.str().c_str());
      TText *textbox_3 = new TText(10,Yaxis_Max*0.80,plotChi2.str().c_str());

      textbox_1->SetTextSize(0.03);
      textbox_2->SetTextSize(0.03);
      textbox_3->SetTextSize(0.03);

      C1->cd();
      backgroundFrame->Draw();
      textbox_1->Draw();
      textbox_2->Draw();
      textbox_3->Draw();
    }
  }

else {
    backgroundFrame->Delete();
 }
  //RooTrace::mark();
  value->Delete();
  totalpdf->Delete();  
  FitResult->Delete();
  //intbckgrd->Delete();
  //intsignal->Delete();
  inttotal->Delete();
  //RooTrace::mark();
  //  RooTrace::active(1);                                                                                               
  //RooTrace::verbose(true);  
  //  cerr<<endl<<endl<<endl<<"this dump is right after deleting all pointers, and is at the end of the qfitter code"<<endl<<endl;
  //RooTrace::dump();
  //cerr<<"end of qfitter code"<<endl<<endl;  
}
// end of file