#include <getopt.h>
#include <fstream>
#include <sstream>
#include <string>
#include "Examples/JpsiGamKsKlKK/JpsiGamKsKlKKHist.hh"
#include "Examples/JpsiGamKsKlKK/JpsiGamKsKlKKEventList.hh"
#include "Examples/JpsiGamKsKlKK/AbsJpsiGamKsKlKKLh.hh"
#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TMath.h"
#include "TNtuple.h"
#include "ErrLogger/ErrLogger.hh"

JpsiGamKsKlKKHist::JpsiGamKsKlKKHist(boost::shared_ptr<const JpsiGamKsKlKKEventList> theEvtList) :
  _theTFile(0),
  _PhiPhiMassDataHist(0),
  _PhiPhiMassMcHist(0),
  _PhiPhiMassFittedHist(0),
  _KpKmMassDataHist(0),
  _KpKmMassMcHist(0),
  _KpKmMassFittedHist(0),
  _KsKlMassDataHist(0),
  _KsKlMassMcHist(0),
  _KsKlMassFittedHist(0),
  
  _dataTuple(0),
  _mcTuple(0)
{
  if(0==theEvtList){
    Alert <<"JpsiGamKsKlKKEventList* theEvtList is a 0 pointer !!!!" ;  // << endmsg;
    exit(1);
  }

  initRootStuff();

  const std::vector<JpsiGamKsKlKKEvtData*> dataList=theEvtList->getDataVecs();

  std::vector<JpsiGamKsKlKKEvtData*>::const_iterator it=dataList.begin();
  while(it!=dataList.end())
    {
      plotDalitz(_dalitzDataHist, (*it), 1.);
      plotPhiPhi(_PhiPhiMassDataHist, (*it), 1.  );
      plotKsKl( _KsKlMassDataHist, (*it), 1. );
      plotKpKm( _KpKmMassDataHist, (*it), 1. );
      
      fillTuple(_dataTuple, (*it), 1.);

      ++it;
    }

  const std::vector<JpsiGamKsKlKKEvtData*> mcList=theEvtList->getMcVecs();
  it=mcList.begin();
  while(it!=mcList.end())
    { 
      plotDalitz(_dalitzMcHist, (*it), 1.);
      plotDalitz(_dalitzMcHist, (*it), 1.);
      plotPhiPhi(_PhiPhiMassMcHist, (*it), 1.  );
      plotKsKl( _KsKlMassMcHist, (*it), 1. );
      plotKpKm( _KpKmMassMcHist, (*it), 1. );
     fillTuple(_mcTuple, (*it), 1.);

      ++it;
    }
}

JpsiGamKsKlKKHist::JpsiGamKsKlKKHist(boost::shared_ptr<AbsJpsiGamKsKlKKLh> theJpsiGamKsKlKKLh, paramGamKsKlKK& fitParam) :
  _theTFile(0),
  _dalitzDataHist(0),
  _dalitzMcHist(0),
  _dalitzFittedHist(0),
  _PhiPhiMassDataHist(0),
  _PhiPhiMassMcHist(0),
  _PhiPhiMassFittedHist(0),
  _KpKmMassDataHist(0),
  _KpKmMassMcHist(0),
  _KpKmMassFittedHist(0),
  _KsKlMassDataHist(0),
  _KsKlMassMcHist(0),
  _KsKlMassFittedHist(0),
  
  _dataTuple(0),
  _mcTuple(0)
{
  if(0==theJpsiGamKsKlKKLh){
    Alert <<"JpsiGamKsKlKKLh* theJpsiGamKsKlKKLh is a 0 pointer !!!!" ;  // << endmsg;
    exit(1);
  }

  initRootStuff();

  boost::shared_ptr<const JpsiGamKsKlKKEventList> theEvtList=theJpsiGamKsKlKKLh->getEventList();
  const std::vector<JpsiGamKsKlKKEvtData*> dataList=theEvtList->getDataVecs();

  std::vector<JpsiGamKsKlKKEvtData*>::const_iterator it=dataList.begin();
  while(it!=dataList.end())
    {
      plotDalitz(_dalitzDataHist, (*it), 1.);
      plotPhiPhi(_PhiPhiMassDataHist, (*it), 1.  );
      plotKsKl( _KsKlMassDataHist, (*it), 1. );
      plotKpKm( _KpKmMassDataHist, (*it), 1. );
      
      fillTuple(_dataTuple, (*it), 1.);
      
      ++it;
    }
  
  const std::vector<JpsiGamKsKlKKEvtData*> mcList=theEvtList->getMcVecs();
  it=mcList.begin();
  while(it!=mcList.end())
    { 
      plotDalitz(_dalitzMcHist, (*it), 1.);
      plotPhiPhi(_PhiPhiMassMcHist, (*it), 1.  );
      plotKsKl( _KsKlMassMcHist, (*it), 1. );
      plotKpKm( _KpKmMassMcHist, (*it), 1. );
      
      double evtWeight= theJpsiGamKsKlKKLh->calcEvtIntensity((*it), fitParam);
      plotDalitz(_dalitzFittedHist, (*it),evtWeight );
      plotPhiPhi(_PhiPhiMassFittedHist, (*it), evtWeight  );
      plotKsKl( _KsKlMassFittedHist, (*it), evtWeight );
      plotKpKm( _KpKmMassFittedHist, (*it), evtWeight );
      
      fillTuple(_mcTuple, (*it), evtWeight);
      ++it;
    }
  
  
  double integralData=(double) theEvtList->getDataVecs().size();
  Info <<"No of fit data events  " << integralData ;  // << endmsg;
  
  double integralFitted=(double) theEvtList->getMcVecs().size();
  Info <<"No of fit events " << integralFitted ;  // << endmsg; 

  Info <<"scaling factor  " << integralData/integralFitted ;  // << endmsg; 

  double scaleFactor = integralData/integralFitted;
  
  _dalitzMcHist->Scale(scaleFactor);
  _PhiPhiMassFittedHist->Scale(scaleFactor);
  _KsKlMassFittedHist->Scale(scaleFactor);
  _KpKmMassFittedHist->Scale(scaleFactor);
  
}

JpsiGamKsKlKKHist::~JpsiGamKsKlKKHist()
{
  _theTFile->Write();
  _theTFile->Close();
}

void JpsiGamKsKlKKHist::initRootStuff()
{
  std::string rootFileName="./JpsiGamKsKlKK.root";

  _theTFile=new TFile(rootFileName.c_str(),"recreate");
  int xbins = 50;
  double xmin=0.8;
  double xmax =10.;
  int ybins=xbins;
  double ymin=xmin;
  double ymax=xmax;
  
  _dalitzDataHist= new TH2F("_dalitzDataHist","Dpl K+K- K+#pi^{0} data",xbins, xmin, xmax, ybins, ymin, ymax );
  _dalitzMcHist= new TH2F("_dalitzMcHist","Dpl K+K- K+#pi^{0} MC",xbins, xmin, xmax, ybins, ymin, ymax);
  _dalitzFittedHist= new TH2F("_dalitzFittedHist","Dpl K+K- K+#pi^{0} fit",xbins, xmin, xmax, ybins, ymin, ymax );
  
  int nbins=50;
   xmin=2;
   xmax=3.2;
  _PhiPhiMassDataHist = new TH1F("_PhiPhiMassDataHist", "_PhiPhiMassDataHist", nbins, xmin, xmax);
  _PhiPhiMassMcHist = new TH1F("_PhiPhiMassMcHist", "_PhiPhiMassMcHist", nbins, xmin, xmax);
  _PhiPhiMassFittedHist = new TH1F("_PhiPhiMassFittedHist", "_PhiPhiMassFittedHist", nbins, xmin, xmax);
  
  xmin=0.8;
  xmax=1.2;
  _KsKlMassDataHist = new TH1F("_KsKlMassDataHist", "_KsKlMassDataHist", nbins, xmin, xmax);
  _KsKlMassMcHist = new TH1F("_KsKlMassMcHist", "_KsKlMassMcHist", nbins, xmin, xmax);
  _KsKlMassFittedHist = new TH1F("_KsKlMassFittedHist", "_KsKlMassFittedHist", nbins, xmin, xmax);
  
  _KpKmMassDataHist = new TH1F("_KpKmMassDataHist", "_KpKmMassDataHist", nbins, xmin, xmax);
  _KpKmMassMcHist = new TH1F("_KpKmMassMcHist", "_KpKmMassMcHist", nbins, xmin, xmax);
  _KpKmMassFittedHist = new TH1F("_KpKmMassFittedHist", "_KpKmMassFittedHist", nbins, xmin, xmax);
  
  
   std::string tupleVariables = "mKsKlKpKm:mKsKl:mKpKm:KsKlKpKmCostTheta:gamCosTheta:KsKlCosTheta:KpKmCosTheta:KsCosTheta:KpCosTheta:weight";
   
  
   _dataTuple=new TNtuple("_dataTuple", "data ntuple", tupleVariables.data());
   _mcTuple=new TNtuple("_mcTuple", "mc ntuple", tupleVariables.data());
}

void JpsiGamKsKlKKHist::plotDalitz(TH2F* theHisto, const JpsiGamKsKlKKEvtData* theData, double weight)
{
  Vector4<float> V4_KsKlKpKm_HeliPsi;
  Vector4<float> V4_KsKl_HeliPsi = theData->V4_KsKl_HeliPsi  ;
  Vector4<float> V4_KpKm_HeliPsi = theData->V4_KpKm_HeliPsi ;
  Vector4<float> V4_gamma_HeliPsi = theData->V4_gamma_HeliPsi ;
  
  double gphi1 = (V4_gamma_HeliPsi + V4_KsKl_HeliPsi ).M();
  double gphi2 = (V4_gamma_HeliPsi + V4_KpKm_HeliPsi ).M();
  theHisto->Fill( gphi1*gphi1, gphi2*gphi2   ,weight);  
}

void JpsiGamKsKlKKHist::plotPhiPhi(TH1F* theHisto, const JpsiGamKsKlKKEvtData* theData, double weight){
  Vector4<float> v4 = theData->V4_KsKlKpKm_HeliPsi;
  theHisto->Fill( v4.M() );
}

void JpsiGamKsKlKKHist::plotKsKl(TH1F* theHisto, const JpsiGamKsKlKKEvtData* theData, double weight){
  
}
void JpsiGamKsKlKKHist::plotKpKm(TH1F* theHisto, const JpsiGamKsKlKKEvtData* theData, double weight){
  
}

void  JpsiGamKsKlKKHist::fillTuple( TNtuple* theTuple, const JpsiGamKsKlKKEvtData* theData, double weight)
{
  
  Vector4<float> V4_KsKlKpKm_HeliPsi = theData->V4_KsKlKpKm_HeliPsi  ;
  Vector4<float> V4_KsKl_HeliPsi= theData->V4_KsKl_HeliPsi ;
  Vector4<float> V4_KpKm_HeliPsi= theData->V4_KpKm_HeliPsi ;
  Vector4<float> V4_gamma_HeliPsi= theData->V4_gamma_HeliPsi ;
  
  Vector4<float> V4_KsKl_HeliKsKlKpKm= theData->V4_KsKl_HeliKsKlKpKm ;
  Vector4<float> V4_KpKm_HeliKsKlKpKm= theData->V4_KpKm_HeliKsKlKpKm ;
  Vector4<float> V4_Ks_HeliKsKl= theData->V4_Ks_HeliKsKl ;
  Vector4<float> V4_Kp_HeliKpKm= theData->V4_Kp_HeliKpKm ;
  
  
  theTuple->Fill( V4_KsKlKpKm_HeliPsi.M(),
		  V4_KsKl_HeliPsi.M(),
		  V4_KpKm_HeliPsi.M(),
		  V4_KsKlKpKm_HeliPsi.CosTheta(),
		  V4_gamma_HeliPsi.CosTheta(),
		  V4_KsKl_HeliKsKlKpKm.CosTheta(),
		  V4_KpKm_HeliKsKlKpKm.CosTheta(),
		  V4_Ks_HeliKsKl.CosTheta(),
		  V4_Kp_HeliKpKm.CosTheta(),
		  weight);
}