#include <getopt.h>
#include <fstream>
#include <sstream>
#include <string>
#include "Examples/MATpbarpToOmegaPi/OmegaPiHist.hh"
#include "Examples/MATpbarpToOmegaPi/OmegaPiEventList.hh"
#include "Examples/MATpbarpToOmegaPi/OmegaPiLh.hh"
#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TMath.h"
#include "ErrLogger/ErrLogger.hh"

OmegaPiHist::OmegaPiHist(boost::shared_ptr<const OmegaPiEventList> theEvtList) :
  _theTFile(0),
  _cosOmegaHeliDataHist(0),
  _cosOmegaHeliMcHist(0),
  _cosOmegaHeliFittedHist(0),
  _cosOmegaAccCorHist(0),
  _cosPi0FromOmegaDataHeli(0),
  _cosPi0FromOmegaMcHeli(0),
  _cosPi0FromOmegaFittedHeli(0),
  _cosPi0FromOmegaAccCorHeli(0),
  _treimanYangDataHist(0),
  _treimanYangMcHist(0),
  _treimanYangFittedHist(0)
{
  if(0==theEvtList){
    Alert <<"OmegaPiEventList* theEvtList is a 0 pointer !!!!" << endmsg;
    exit(1);
  }

  _jmax=theEvtList->jMax();
  _pbarmom=theEvtList->pbarMom();

  initRootStuff();

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

  std::vector<OmPiEvtData*>::const_iterator it=dataList.begin();
  while(it!=dataList.end())
    {
      plotCosOmegaHeli(_cosOmegaHeliDataHist, (*it), 1.);
      plotCosPi0FromOmegaHeli(_cosPi0FromOmegaDataHeli, (*it), 1.);
      plotTreimanYang(_treimanYangDataHist, (*it), 1.);
      ++it;
    }

  const std::vector<OmPiEvtData*> mcList=theEvtList->getMcVecs();
  it=mcList.begin();
  while(it!=mcList.end())
    {
      plotCosOmegaHeli(_cosOmegaHeliMcHist, (*it), 1.);
      plotCosPi0FromOmegaHeli(_cosPi0FromOmegaMcHeli, (*it), 1.);
      plotTreimanYang(_treimanYangMcHist, (*it), 1.);
      ++it;
    }
  _cosOmegaHeliDataHist->Sumw2();
  _cosOmegaHeliMcHist->Sumw2();
  _cosOmegaAccCorHist=(TH1F*) _cosOmegaHeliDataHist->Clone();
  _cosOmegaAccCorHist->Divide(_cosOmegaHeliMcHist);
  _cosOmegaAccCorHist->SetName("_cosOmegaAccCorHist");
  _cosOmegaAccCorHist->SetTitle("cos(#Theta) #omega heli acc cor");

  _cosPi0FromOmegaDataHeli->Sumw2();
  _cosPi0FromOmegaMcHeli->Sumw2();
  _cosPi0FromOmegaAccCorHeli=(TH1F*) _cosPi0FromOmegaDataHeli->Clone();
  _cosPi0FromOmegaAccCorHeli->Divide(_cosPi0FromOmegaMcHeli);
  _cosPi0FromOmegaAccCorHeli->SetName("_cosPi0FromOmegaAccCorHeli");
  _cosPi0FromOmegaAccCorHeli->SetTitle("cos(#Theta) #pi^0 heli (#omega decay) acc cor");
}

OmegaPiHist::OmegaPiHist(boost::shared_ptr<OmegaPiLh> omegaPiLh, OmegaPiData::fitParamVal& fitParam) :
  _theTFile(0),
  _cosOmegaHeliDataHist(0),
  _cosOmegaHeliMcHist(0),
  _cosOmegaHeliFittedHist(0),
  _cosOmegaAccCorHist(0),
  _cosPi0FromOmegaDataHeli(0),
  _cosPi0FromOmegaMcHeli(0),
  _cosPi0FromOmegaFittedHeli(0),
  _cosPi0FromOmegaAccCorHeli(0),
  _treimanYangDataHist(0),
  _treimanYangMcHist(0),
  _treimanYangFittedHist(0)
{
  if(0==omegaPiLh){
    Alert <<"OmegaPiLh* omegaPiLh is a 0 pointer !!!!" << endmsg;
    exit(1);
  }

  boost::shared_ptr<const OmegaPiEventList> theEvtList=omegaPiLh->getEventList();  
  if(0==theEvtList){
    Alert <<"OmegaPiEventList* theEvtList is a 0 pointer !!!!" << endmsg;
    exit(1);
  }

  _jmax=theEvtList->jMax();
  _pbarmom=theEvtList->pbarMom();

  initRootStuff();

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

  std::vector<OmPiEvtData*>::iterator it=dataList.begin();
  while(it!=dataList.end())
    {
      plotCosOmegaHeli(_cosOmegaHeliDataHist, (*it), 1.);
      plotCosPi0FromOmegaHeli(_cosPi0FromOmegaDataHeli, (*it), 1.);
      plotTreimanYang(_treimanYangDataHist, (*it), 1.);
      ++it;
    }

  std::vector<OmPiEvtData*> mcList=theEvtList->getMcVecs();
  it=mcList.begin();
  while(it!=mcList.end())
    {
      plotCosOmegaHeli(_cosOmegaHeliMcHist, (*it), 1.);
      plotCosPi0FromOmegaHeli(_cosPi0FromOmegaMcHeli, (*it), 1.);
      plotTreimanYang(_treimanYangMcHist, (*it), 1.);
      double evtWeight= omegaPiLh->calcEvtIntensity((*it), fitParam);
      plotCosOmegaHeli(_cosOmegaHeliFittedHist, (*it), evtWeight);
      plotCosPi0FromOmegaHeli(_cosPi0FromOmegaFittedHeli, (*it), evtWeight);
      plotTreimanYang(_treimanYangFittedHist, (*it), evtWeight);
      ++it;
    }
  _cosOmegaHeliDataHist->Sumw2();
  _cosOmegaHeliMcHist->Sumw2();
  _cosOmegaAccCorHist=(TH1F*) _cosOmegaHeliDataHist->Clone();
  _cosOmegaAccCorHist->Divide(_cosOmegaHeliMcHist);
  _cosOmegaAccCorHist->SetName("_cosOmegaAccCorHist");
  _cosOmegaAccCorHist->SetTitle("cos(#Theta) #omega heli acc cor");

  _cosPi0FromOmegaDataHeli->Sumw2();
  _cosPi0FromOmegaMcHeli->Sumw2();
  _cosPi0FromOmegaAccCorHeli=(TH1F*) _cosPi0FromOmegaDataHeli->Clone();
  _cosPi0FromOmegaAccCorHeli->Divide(_cosPi0FromOmegaMcHeli);
  _cosPi0FromOmegaAccCorHeli->SetName("_cosPi0FromOmegaAccCorHeli");
  _cosPi0FromOmegaAccCorHeli->SetTitle("cos(#Theta) #pi^0 heli (#omega decay) acc cor");


  double integralData=_cosOmegaHeliDataHist->Integral();
  DebugMsg << "integralData= " << integralData << endmsg;

  double integralFitted=_cosOmegaHeliFittedHist->Integral();  
  DebugMsg << "integralFittedHists= " << integralFitted << endmsg;

  _cosOmegaHeliFittedHist->Scale(integralData/integralFitted);
  _cosPi0FromOmegaFittedHeli->Scale(integralData/integralFitted);
  _treimanYangFittedHist->Scale(integralData/integralFitted);
}

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

void OmegaPiHist::initRootStuff()
{
  std::stringstream jmaxStrStr;
  std::stringstream pbarMomStrStr;

  jmaxStrStr << _jmax;
  pbarMomStrStr << _pbarmom;

  std::string rootFileName="./OmegaPi0Fit_jmax"+jmaxStrStr.str()+"_mom"+pbarMomStrStr.str()+".root";

  _theTFile=new TFile(rootFileName.c_str(),"recreate");
  _cosOmegaHeliDataHist= new TH1F("_cosOmegaHeliDataHist","cos(#Theta) #omega heli data",101, -1.0, 1.0);
  _cosOmegaHeliMcHist= new TH1F("_cosOmegaHeliMcHist","cos(#Theta) #omega heli MC",101, -1.0, 1.0);
  _cosOmegaHeliFittedHist= new TH1F("_cosOmegaHeliFittedHist","cos(#Theta) #omega heli fit result",101, -1.0, 1.0);
  _cosPi0FromOmegaDataHeli= new TH1F("_cosPi0FromOmegaDataHeli","cos(#Theta) #pi^{0} heli (#omega decay) data",101, -1.0, 1.0);
  _cosPi0FromOmegaMcHeli= new TH1F("_cosPi0FromOmegaMcHeli","cos(#Theta) #pi^0 heli (#omega decay) MC",101, -1.0, 1.0);
  _cosPi0FromOmegaFittedHeli= new TH1F("_cosPi0FromOmegaFittedHeli","cos(#Theta) #pi^0 heli (#omega decay) fit result",101, -1.0, 1.0);
  _treimanYangDataHist= new TH1F("_treimanYangDataHist","Treiman Yang angle data",101, -TMath::Pi(), TMath::Pi());
  _treimanYangMcHist= new TH1F("_treimanYangMcHist","Treiman Yang angle MC",101, -TMath::Pi(), TMath::Pi());
  _treimanYangFittedHist= new TH1F("_treimanYangFittedHist","Treiman Yang angle fit result",101, -TMath::Pi(), TMath::Pi());
}

void OmegaPiHist::plotCosOmegaHeli(TH1F* theHisto, const OmPiEvtData* theEvtData, double weight)
{
  Vector4<float> omegaHeli4V=theEvtData->omegaHeliCm4Vec;
  theHisto->Fill(omegaHeli4V.CosTheta(), weight);  
}

void OmegaPiHist::plotCosPi0FromOmegaHeli(TH1F* theHisto, const OmPiEvtData* theEvtData, double weight)
{
  Vector4<float> piFromOmegaHeli4V=theEvtData->pi0HeliOmega4Vec;
  theHisto->Fill(piFromOmegaHeli4V.CosTheta(), weight);  
}

void OmegaPiHist::plotTreimanYang(TH1F* theHisto, const OmPiEvtData* theEvtData, double weight)
{
  Vector4<float> piFromOmegaHeli4V=theEvtData->pi0HeliOmega4Vec;
  theHisto->Fill(piFromOmegaHeli4V.Phi(), weight);  
}