#include <getopt.h>
#include <fstream>
#include <sstream>
#include <string>

#include "Examples/JpsiToOmegaPhiGam/JpsiToOmegaPhiGamHist.hh"
#include "Examples/JpsiToOmegaPhiGam/JpsiToOmegaPhiGamEventList.hh"
//#include "Examples/JpsiGamPipPimPi0KK/JpsiGamPipPimPi0KKData.hh"
//#include "Examples/JpsiGamPipPimPi0KK/JpsiGamPipPimPi0KKProdLh.hh"
#include "PwaUtils/KinUtils.hh"
//#include "Examples/JpsiGamPipPimPi0KK/JpsiGamPipPimPi0KKFitParams.hh"

#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TMath.h"
#include "TNtuple.h"
#include "TLorentzVector.h"
#include "ErrLogger/ErrLogger.hh"


JpsiToOmegaPhiGamHist::JpsiToOmegaPhiGamHist(boost::shared_ptr<const EvtDataBaseListNew> theEvtList) :
  _theTFile(0),
  _dalitzDataHist(0),
  _dalitzMcHist(0),
  _dalitzFittedHist(0),
  _PhiPhiMassDataHist(0),
  _PhiPhiMassMcHist(0),
  _PhiPhiMassFittedHist(0),
  _KpKmMassDataHist(0),
  _KpKmMassMcHist(0),
  _KpKmMassFittedHist(0),
  _PipPimPi0MassDataHist(0),
  _PipPimPi0MassMcHist(0),
  _PipPimPi0MassFittedHist(0),
  _costPhi_KpKmDataHist(0), 
  _costPhi_KpKmMcHist(0), 
  _costPhi_KpKmFittedHist(0),
  _phiPhi_KpKmDataHist(0), 
  _phiPhi_KpKmMcHist(0), 
  _phiPhi_KpKmFittedHist(0), 
  _chiDataHist(0), 
  _chiMcHist(0), 
  _chiFittedHist(0),  
  _dataTuple(0),
  _mcTuple(0),
  _massIndepTuple(0),
  _massRange(make_pair(0.,0.) )
{
  if(0==theEvtList){
    Alert <<"JpsiGamPipPimPi0KKEventList* theEvtList is a 0 pointer !!!!" ;  // << endmsg;
    exit(1);
  }

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

  std::vector<EvtDataNew*>::const_iterator it=dataList.begin();
  while(it!=dataList.end())
    {
      const double evtWeight=(*it)->evtWeight;
      plotDalitz(_dalitzDataHist, (*it), evtWeight);
      plotPhiPhi(_PhiPhiMassDataHist, (*it), evtWeight  );
      plotPipPimPi0( _PipPimPi0MassDataHist, (*it), evtWeight );
      plotKpKm( _KpKmMassDataHist, (*it), evtWeight );
      plotCostPhiPip( _costPip_PipPimPi0HeliDataHist,  _phiPip_PipPimPi0HeliDataHist, (*it), evtWeight );
      plotCostPhiKp( _costKp_KpKmHeliDataHist, _phiKp_KpKmHeliDataHist, (*it), evtWeight );
      plotCostGam( _costGamCmDataHist, (*it), evtWeight );
      plotCostPhi_PhiPhiHeli(_costPhi_KpKmDataHist, _phiPhi_KpKmDataHist, (*it)->FourVecsDec[enumJpsiGamX4V::V4_KpKm_HeliPipPimPi0KpKm], evtWeight);  
       plotChi(_chiDataHist, (*it), evtWeight );      
      fillTuple(_dataTuple, (*it), evtWeight);
      
      ++it;
    }

 const std::vector<EvtDataNew*> mcList=theEvtList->getMcVecs();  
  it=mcList.begin();
  while(it!=mcList.end())
    { 
      plotDalitz(_dalitzMcHist, (*it), 1.);
      plotPhiPhi(_PhiPhiMassMcHist, (*it), 1.  );
      plotPipPimPi0( _PipPimPi0MassMcHist, (*it), 1. );
      plotKpKm( _KpKmMassMcHist, (*it), 1. );
      plotCostPhiPip( _costPip_PipPimPi0HeliMcHist, _phiPip_PipPimPi0HeliMcHist, (*it), 1. );
      plotCostPhiKp( _costKp_KpKmHeliMcHist, _phiKp_KpKmHeliMcHist,(*it), 1. ); 
      plotCostGam( _costGamCmMcHist,(*it), 1. );
      plotCostPhi_PhiPhiHeli(_costPhi_KpKmMcHist, _phiPhi_KpKmMcHist, (*it)->FourVecsDec[enumJpsiGamX4V::V4_KpKm_HeliPipPimPi0KpKm], 1.);  
      plotChi(_chiMcHist, (*it), 1. );   
      
      // double evtWeight= _theJpsiGamPipPimPi0KKLh->calcEvtIntensity((*it), _fitParam);
      double evtWeight=1.;
      plotDalitz(_dalitzFittedHist, (*it),evtWeight );
      plotPhiPhi(_PhiPhiMassFittedHist, (*it), evtWeight  );
      plotPipPimPi0( _PipPimPi0MassFittedHist, (*it), evtWeight );
      plotKpKm( _KpKmMassFittedHist, (*it), evtWeight );
      plotCostPhiPip( _costPip_PipPimPi0HeliFittedHist, _phiPip_PipPimPi0HeliFittedHist,(*it), evtWeight );
      plotCostPhiKp( _costKp_KpKmHeliFittedHist, _phiKp_KpKmHeliFittedHist,(*it), evtWeight );
      plotCostGam( _costGamCmFittedHist,(*it), evtWeight );
      plotCostPhi_PhiPhiHeli(_costPhi_KpKmFittedHist, _phiPhi_KpKmFittedHist, (*it)->FourVecsDec[enumJpsiGamX4V::V4_KpKm_HeliPipPimPi0KpKm], evtWeight); 
      plotChi(_chiFittedHist, (*it), evtWeight  );
        
      fillTuple(_mcTuple, (*it), evtWeight);
      ++it;
    }

  double integralDataWoWeight=(double) theEvtList->getDataVecs().size();
  Info <<"No of data events without weight " << integralDataWoWeight ;  // << endmsg;

  double integralDataWWeight = theEvtList->NoOfWeightedDataEvts();
  Info <<"No of data events with weight " << integralDataWWeight ;  // << endmsg;
  
  double integralMC = theEvtList->NoOfWeightedMcEvts();
  Info <<"No of MC events " << integralMC ;  // << endmsg; 

  Info <<"scaling factor  " << integralDataWWeight/integralMC ;  // << endmsg; 

  double scaleFactor = integralDataWWeight/integralMC;
  
  _dalitzFittedHist->Scale(scaleFactor);
  _PhiPhiMassFittedHist->Scale(scaleFactor);
  _PipPimPi0MassFittedHist->Scale(scaleFactor);
  _KpKmMassFittedHist->Scale(scaleFactor);
  _costPip_PipPimPi0HeliFittedHist->Scale(scaleFactor);
  _phiPip_PipPimPi0HeliFittedHist->Scale(scaleFactor);
  _costKp_KpKmHeliFittedHist->Scale(scaleFactor);
  _phiKp_KpKmHeliFittedHist->Scale(scaleFactor);
  _costGamCmFittedHist->Scale(scaleFactor);
  _costPhi_KpKmFittedHist->Scale(scaleFactor);
  _phiPhi_KpKmFittedHist->Scale(scaleFactor);
  _chiFittedHist->Scale(scaleFactor);


}

JpsiToOmegaPhiGamHist::JpsiToOmegaPhiGamHist(boost::shared_ptr<AbsLhNew> theLh, fitParamsNew& theFitParms):
  _theTFile(0),
  _dalitzDataHist(0),
  _dalitzMcHist(0),
  _dalitzFittedHist(0),
  _PhiPhiMassDataHist(0),
  _PhiPhiMassMcHist(0),
  _PhiPhiMassFittedHist(0),
  _KpKmMassDataHist(0),
  _KpKmMassMcHist(0),
  _KpKmMassFittedHist(0),
  _PipPimPi0MassDataHist(0),
  _PipPimPi0MassMcHist(0),
  _PipPimPi0MassFittedHist(0),
  _costPhi_KpKmDataHist(0), 
  _costPhi_KpKmMcHist(0), 
  _costPhi_KpKmFittedHist(0),
  _phiPhi_KpKmDataHist(0), 
  _phiPhi_KpKmMcHist(0), 
  _phiPhi_KpKmFittedHist(0), 
  _chiDataHist(0), 
  _chiMcHist(0), 
  _chiFittedHist(0),  
  _dataTuple(0),
  _mcTuple(0),
  _massIndepTuple(0),
  _massRange(make_pair(0.,0.) )
{ 
  if(0==theLh){
    Alert <<"AbsLh* is a 0 pointer !!!!" ;  // << endmsg;
    exit(1);
  }

   initRootStuff();
   boost::shared_ptr<const EvtDataBaseListNew> theEvtList=theLh->getEventList();
   const std::vector<EvtDataNew*> dataList=theEvtList->getDataVecs();  

   std::vector<EvtDataNew*>::const_iterator it=dataList.begin();
  while(it!=dataList.end())
    {
      const double evtWeight=(*it)->evtWeight;
      plotDalitz(_dalitzDataHist, (*it), evtWeight);
      plotPhiPhi(_PhiPhiMassDataHist, (*it), evtWeight  );
      plotPipPimPi0( _PipPimPi0MassDataHist, (*it), evtWeight );
      plotKpKm( _KpKmMassDataHist, (*it), evtWeight );
      plotCostPhiPip( _costPip_PipPimPi0HeliDataHist,  _phiPip_PipPimPi0HeliDataHist, (*it), evtWeight );
      plotCostPhiKp( _costKp_KpKmHeliDataHist, _phiKp_KpKmHeliDataHist, (*it), evtWeight );
      plotCostGam( _costGamCmDataHist, (*it), evtWeight );
      plotCostPhi_PhiPhiHeli(_costPhi_KpKmDataHist, _phiPhi_KpKmDataHist, (*it)->FourVecsDec[enumJpsiGamX4V::V4_KpKm_HeliPipPimPi0KpKm], evtWeight);  
       plotChi(_chiDataHist, (*it), evtWeight );      
      fillTuple(_dataTuple, (*it), evtWeight);
      
      ++it;
    }

 const std::vector<EvtDataNew*> mcList=theEvtList->getMcVecs();  
  it=mcList.begin();
  while(it!=mcList.end())
    { 
      plotDalitz(_dalitzMcHist, (*it), 1.);
      plotPhiPhi(_PhiPhiMassMcHist, (*it), 1.  );
      plotPipPimPi0( _PipPimPi0MassMcHist, (*it), 1. );
      plotKpKm( _KpKmMassMcHist, (*it), 1. );
      plotCostPhiPip( _costPip_PipPimPi0HeliMcHist, _phiPip_PipPimPi0HeliMcHist, (*it), 1. );
      plotCostPhiKp( _costKp_KpKmHeliMcHist, _phiKp_KpKmHeliMcHist,(*it), 1. ); 
      plotCostGam( _costGamCmMcHist,(*it), 1. );
      plotCostPhi_PhiPhiHeli(_costPhi_KpKmMcHist, _phiPhi_KpKmMcHist, (*it)->FourVecsDec[enumJpsiGamX4V::V4_KpKm_HeliPipPimPi0KpKm], 1.);  
      plotChi(_chiMcHist, (*it), 1. );   
      
      // double evtWeight= _theJpsiGamPipPimPi0KKLh->calcEvtIntensity((*it), _fitParam);
      double evtWeight= theLh->calcEvtIntensity( (*it), theFitParms );
      plotDalitz(_dalitzFittedHist, (*it),evtWeight );
      plotPhiPhi(_PhiPhiMassFittedHist, (*it), evtWeight  );
      plotPipPimPi0( _PipPimPi0MassFittedHist, (*it), evtWeight );
      plotKpKm( _KpKmMassFittedHist, (*it), evtWeight );
      plotCostPhiPip( _costPip_PipPimPi0HeliFittedHist, _phiPip_PipPimPi0HeliFittedHist,(*it), evtWeight );
      plotCostPhiKp( _costKp_KpKmHeliFittedHist, _phiKp_KpKmHeliFittedHist,(*it), evtWeight );
      plotCostGam( _costGamCmFittedHist,(*it), evtWeight );
      plotCostPhi_PhiPhiHeli(_costPhi_KpKmFittedHist, _phiPhi_KpKmFittedHist, (*it)->FourVecsDec[enumJpsiGamX4V::V4_KpKm_HeliPipPimPi0KpKm], evtWeight); 
      plotChi(_chiFittedHist, (*it), evtWeight  );
        
      fillTuple(_mcTuple, (*it), evtWeight);
      ++it;
    }

  double integralDataWoWeight=(double) theEvtList->getDataVecs().size();
  Info <<"No of data events without weight " << integralDataWoWeight ;  // << endmsg;

  double integralDataWWeight=(double) theEvtList->NoOfWeightedDataEvts();
  Info <<"No of data events with weight " << integralDataWWeight ;  // << endmsg;


  double integralMC=(double) theEvtList->NoOfWeightedMcEvts();
  Info <<"No of MC events " << integralMC ;  // << endmsg; 

  Info <<"scaling factor  " << integralDataWWeight/integralMC ;  // << endmsg;

  double scaleFactor = integralDataWWeight/integralMC;
  
  _dalitzFittedHist->Scale(scaleFactor);
  _PhiPhiMassFittedHist->Scale(scaleFactor);
  _PipPimPi0MassFittedHist->Scale(scaleFactor);
  _KpKmMassFittedHist->Scale(scaleFactor);
  _costPip_PipPimPi0HeliFittedHist->Scale(scaleFactor);
  _phiPip_PipPimPi0HeliFittedHist->Scale(scaleFactor);
  _costKp_KpKmHeliFittedHist->Scale(scaleFactor);
  _phiKp_KpKmHeliFittedHist->Scale(scaleFactor);
  _costGamCmFittedHist->Scale(scaleFactor);
  _costPhi_KpKmFittedHist->Scale(scaleFactor);
  _phiPhi_KpKmFittedHist->Scale(scaleFactor);
  _chiFittedHist->Scale(scaleFactor);

  double integralFitted=(double) _costGamCmFittedHist->Integral();
  Info <<"No of fit events " << integralFitted ;  // << endmsg; 
}


// JpsiToOmegaPhiGamHist::JpsiToOmegaPhiGamHist(JpsiGamPipPimPi0KKProdLh* theJpsiGamPipPimPi0KKLh, fitParams& fitParam, FitParamErrorMatrix* theErrorMatrix) :
//   _theTFile(0),
//   _dalitzDataHist(0),
//   _dalitzMcHist(0),
//   _dalitzFittedHist(0),
//   _PhiPhiMassDataHist(0),
//   _PhiPhiMassMcHist(0),
//   _PhiPhiMassFittedHist(0),
//   _KpKmMassDataHist(0),
//   _KpKmMassMcHist(0),
//   _KpKmMassFittedHist(0),
//   _PipPimPi0MassDataHist(0),
//   _PipPimPi0MassMcHist(0),
//   _PipPimPi0MassFittedHist(0),
//   _costPhi_KpKmDataHist(0), 
//   _costPhi_KpKmMcHist(0), 
//   _costPhi_KpKmFittedHist(0),
//   _phiPhi_KpKmDataHist(0), 
//   _phiPhi_KpKmMcHist(0), 
//   _phiPhi_KpKmFittedHist(0), 
//   _chiDataHist(0), 
//   _chiMcHist(0), 
//   _chiFittedHist(0),  
//   _dataTuple(0),
//   _mcTuple(0),
//   _massIndepTuple(0),
//   _massRange(make_pair(0.,0.) )
// {
//   if(0==theJpsiGamPipPimPi0KKLh){
//     Alert <<"JpsiGamPipPimPi0KKLh* theJpsiGamPipPimPi0KKLh is a 0 pointer !!!!" ;  // << endmsg;
//     exit(1);
//   }
  
//   initRootStuff();
//   _theJpsiGamPipPimPi0KKLh = theJpsiGamPipPimPi0KKLh;
//   _fitParam = fitParam;
//   std::vector<double> data;
//   _errMatrix = theErrorMatrix;
  
// }
  

// void JpsiToOmegaPhiGamHist::fill(){

//   boost::shared_ptr<const EvtDataBaseList> theEvtList=_theJpsiGamPipPimPi0KKLh->getEventList();
//   const std::vector<EvtData*> dataList=theEvtList->getDataVecs();

//   std::vector<EvtDataNew*>::const_iterator it=dataList.begin();
//   while(it!=dataList.end())
//     {
//       plotDalitz(_dalitzDataHist, (*it), 1.);
//       plotPhiPhi(_PhiPhiMassDataHist, (*it), 1.  );
//       plotPipPimPi0( _PipPimPi0MassDataHist, (*it), 1. );
//       plotKpKm( _KpKmMassDataHist, (*it), 1. );
//       plotCostPhiPip( _costPip_PipPimPi0HeliDataHist,  _phiPip_PipPimPi0HeliDataHist, (*it), 1. );
//       plotCostPhiKp( _costKp_KpKmHeliDataHist, _phiKp_KpKmHeliDataHist, (*it), 1. );
//       plotCostGam( _costGamCmDataHist, (*it), 1. );
//       plotCostPhi_PhiPhiHeli(_costPhi_KpKmDataHist, _phiPhi_KpKmDataHist, (*it)->FourVecsDec[enumJpsiGamX4V::V4_KpKm_HeliPipPimPi0KpKm], 1.);  
//        plotChi(_chiDataHist, (*it), 1. );      
//       fillTuple(_dataTuple, (*it), 1.);
      
//       ++it;
//     }
  
//   const std::vector<EvtDataNew*> mcList=theEvtList->getMcVecs();
//   it=mcList.begin();
//   while(it!=mcList.end())
//     { 
//       plotDalitz(_dalitzMcHist, (*it), 1.);
//       plotPhiPhi(_PhiPhiMassMcHist, (*it), 1.  );
//       plotPipPimPi0( _PipPimPi0MassMcHist, (*it), 1. );
//       plotKpKm( _KpKmMassMcHist, (*it), 1. );
//       plotCostPhiPip( _costPip_PipPimPi0HeliMcHist, _phiPip_PipPimPi0HeliMcHist, (*it), 1. );
//       plotCostPhiKp( _costKp_KpKmHeliMcHist, _phiKp_KpKmHeliMcHist,(*it), 1. ); 
//       plotCostGam( _costGamCmMcHist,(*it), 1. );
//       plotCostPhi_PhiPhiHeli(_costPhi_KpKmMcHist, _phiPhi_KpKmMcHist, (*it)->FourVecsDec[enumJpsiGamX4V::V4_KpKm_HeliPipPimPi0KpKm], 1.);  
//       plotChi(_chiMcHist, (*it), 1. );   
      
//       double evtWeight= _theJpsiGamPipPimPi0KKLh->calcEvtIntensity((*it), _fitParam);
//       plotDalitz(_dalitzFittedHist, (*it),evtWeight );
//       plotPhiPhi(_PhiPhiMassFittedHist, (*it), evtWeight  );
//       plotPipPimPi0( _PipPimPi0MassFittedHist, (*it), evtWeight );
//       plotKpKm( _KpKmMassFittedHist, (*it), evtWeight );
//       plotCostPhiPip( _costPip_PipPimPi0HeliFittedHist, _phiPip_PipPimPi0HeliFittedHist,(*it), evtWeight );
//       plotCostPhiKp( _costKp_KpKmHeliFittedHist, _phiKp_KpKmHeliFittedHist,(*it), evtWeight );
//       plotCostGam( _costGamCmFittedHist,(*it), evtWeight );
//       plotCostPhi_PhiPhiHeli(_costPhi_KpKmFittedHist, _phiPhi_KpKmFittedHist, (*it)->FourVecsDec[enumJpsiGamX4V::V4_KpKm_HeliPipPimPi0KpKm], evtWeight); 
//       plotChi(_chiFittedHist, (*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;
  
//   _dalitzFittedHist->Scale(scaleFactor);
//   _PhiPhiMassFittedHist->Scale(scaleFactor);
//   _PipPimPi0MassFittedHist->Scale(scaleFactor);
//   _KpKmMassFittedHist->Scale(scaleFactor);
//   _costPip_PipPimPi0HeliFittedHist->Scale(scaleFactor);
//   _phiPip_PipPimPi0HeliFittedHist->Scale(scaleFactor);
//   _costKp_KpKmHeliFittedHist->Scale(scaleFactor);
//   _phiKp_KpKmHeliFittedHist->Scale(scaleFactor);
//   _costGamCmFittedHist->Scale(scaleFactor);
//   _costPhi_KpKmFittedHist->Scale(scaleFactor);
//   _phiPhi_KpKmFittedHist->Scale(scaleFactor);
//   _chiFittedHist->Scale(scaleFactor);  

  
//   double iEta(0.), iEtaErr(0.), iF0(0.), iF0Err(0.), iEta2(0.), iEta2Err(0.), iF1(0.), iF1Err(0.), iF2(0.), iF2Err(0.);
//   double etaReal(0.), etaImg(0.);
  
//   it=mcList.begin();
//   while(it!=mcList.end()){
    
//     std::pair<double, double> intensityEvent = make_pair(0.,0.);
//     _theJpsiGamPipPimPi0KKLh->calcComponentIntensity((*it), _fitParam, *_errMatrix, "etacHyp", intensityEvent);
//     iEta+= intensityEvent.first*scaleFactor;
//     iEtaErr+= intensityEvent.second*scaleFactor;
    
//     _theJpsiGamPipPimPi0KKLh->calcComponentIntensity((*it), _fitParam, *_errMatrix, "f02020Hyp",intensityEvent);
//     iF0 += intensityEvent.first*scaleFactor;
//     iF0Err += intensityEvent.second*scaleFactor;
    
//     _theJpsiGamPipPimPi0KKLh->calcComponentIntensity((*it), _fitParam, *_errMatrix, "f1Hyp", intensityEvent);
//     iF1 += intensityEvent.first*scaleFactor;
//     iF1Err += intensityEvent.second*scaleFactor;
    
//      _theJpsiGamPipPimPi0KKLh->calcComponentIntensity((*it), _fitParam, *_errMatrix, "f22300Hyp",intensityEvent);
//     iF2 += intensityEvent.first*scaleFactor;
//     iF2Err += intensityEvent.second*scaleFactor;
    
//      _theJpsiGamPipPimPi0KKLh->calcComponentIntensity((*it), _fitParam, *_errMatrix, "eta21870Hyp",intensityEvent);
//     iEta2+= intensityEvent.first*scaleFactor;
//     iEta2Err += intensityEvent.second*scaleFactor;
    


    
//     it++;
//   }
//   double meanMassRange = _massRange.first + 0.5*(_massRange.second-_massRange.first);
  
  
//   Info << "Events for component eta : " << iEta << " +/- " << iEtaErr ;
//   Info << "Events for component f0:   " << iF0 << " +/- " << iF0Err ;
//   Info << "Events for component f1:   " << iF1 << " +/- " << iF1Err ;
//   Info << "Events for component f2:   " << iF2 << " +/- " << iF2Err ;
//   Info << "Events for component eta2: " << iEta2 << " +/- " << iEta2Err ;
  
//   std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess > PsiToEtacGamPhi=_fitParam.Phis[paramEnumJpsiGamPipPimPi0KK::PsiToEtacGamma];
//   std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess > PsiToEtacGamMag=_fitParam.Mags[paramEnumJpsiGamPipPimPi0KK::PsiToEtacGamma];
//   std::map< boost::shared_ptr<const JPCLS>, double, pawian::Collection::SharedPtrLess >::iterator phiIter = PsiToEtacGamPhi.begin();
//   for(; phiIter!=PsiToEtacGamPhi.end(); phiIter++){
//     double thePhase = phiIter->second;
//     double theMag = PsiToEtacGamMag[phiIter->first];
    
//     etaImg = theMag*cos(thePhase) ;
//     etaReal = theMag*sin(thePhase);
//   }
  
//   _massIndepTuple->Fill(meanMassRange, iEta, iEtaErr, iF0, iF0Err, iF1, iF1Err, iF2, iF2Err, iEta2, iEta2Err, etaImg, etaReal );
  
// }

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

void JpsiToOmegaPhiGamHist::initRootStuff()
{
  std::string rootFileName="./JpsiGamPipPimPi0KK.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=400;
   xmin=2;
   xmax=4.;
  _PhiPhiMassDataHist = new TH1F("_PhiPhiMassDataHist", "_PhiPhiMassDataHist", nbins, xmin, xmax);
  _PhiPhiMassMcHist = new TH1F("_PhiPhiMassMcHist", "_PhiPhiMassMcHist", nbins, xmin, xmax);
  _PhiPhiMassFittedHist = new TH1F("_PhiPhiMassFittedHist", "_PhiPhiMassFittedHist", nbins, xmin, xmax);

  nbins=50;  
  xmin=0.8;
  xmax=1.2;
  _PipPimPi0MassDataHist = new TH1F("_PipPimPi0MassDataHist", "_PipPimPi0MassDataHist", nbins, xmin, xmax);
  _PipPimPi0MassMcHist = new TH1F("_PipPimPi0MassMcHist", "_PipPimPi0MassMcHist", nbins, xmin, xmax);
  _PipPimPi0MassFittedHist = new TH1F("_PipPimPi0MassFittedHist", "_PipPimPi0MassFittedHist", 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);
  
  _costPip_PipPimPi0HeliDataHist= new TH1F("_costPip_PipPimPi0HeliDataHist", "cos(#Theta_{Pip}) PipPimPi0Heli data", 100, -1., 1.); 
  _costPip_PipPimPi0HeliMcHist= new TH1F("_costPip_PipPimPi0HeliMcHist", "cos(#Theta_{Pip}) PipPimPi0Heli Mc", 100, -1., 1.); 
  _costPip_PipPimPi0HeliFittedHist= new TH1F("_costPip_PipPimPi0HeliFittedHist", "cos(#Theta_{Pip}) PipPimPi0Heli Fit", 100, -1, 1);
  _phiPip_PipPimPi0HeliDataHist= new TH1F("_phiPip_PipPimPi0HeliDataHist", "#Phi_{Pip} PipPimPi0Heli data", 100, -TMath::Pi(), TMath::Pi()); 
  _phiPip_PipPimPi0HeliMcHist= new TH1F("_phiPip_PipPimPi0HeliMcHist", "#Phi_{Pip} PipPimPi0Heli Mc", 100, -TMath::Pi(), TMath::Pi()); 
  _phiPip_PipPimPi0HeliFittedHist= new TH1F("_phiPip_PipPimPi0HeliFittedHist", "#Phi_{Pip} PipPimPi0Heli Fit", 100, -TMath::Pi(), TMath::Pi());

  _costKp_KpKmHeliDataHist= new TH1F("_costKp_KpKmHeliDataHist", "cos(#Theta_{K+}) PipPimPi0Heli data", 100, -1., 1.); 
  _costKp_KpKmHeliMcHist= new TH1F("_costKp_KpKmHeliMcHist", "cos(#Theta_{K+}) PipPimPi0Heli Mc", 100, -1., 1.); 
  _costKp_KpKmHeliFittedHist= new TH1F("_costKp_KpKmHeliFittedHist", "cos(#Theta_{K+}) PipPimPi0Heli Fit", 100, -1, 1);
  _phiKp_KpKmHeliDataHist= new TH1F("_phiKp_KpKmHeliDataHist", "#Phi_{K+} PipPimPi0Heli data", 100, -TMath::Pi(), TMath::Pi()); 
  _phiKp_KpKmHeliMcHist= new TH1F("_phiKp_KpKmHeliMcHist", "#Phi_{K+} PipPimPi0Heli Mc", 100, -TMath::Pi(), TMath::Pi()); 
  _phiKp_KpKmHeliFittedHist= new TH1F("_phiKp_KpKmHeliFittedHist", "#Phi_{K+} PipPimPi0Heli Fit", 100, -TMath::Pi(), TMath::Pi());


  _costGamCmDataHist= new TH1F("_costGamCmDataHist", "cos(#Theta_{#gamma}) CM data", 100, -1., 1.); 
  _costGamCmMcHist= new TH1F("_costGamCmMcHist", "cos(#Theta_{#gamma}) CM Mc", 100, -1., 1.); 
  _costGamCmFittedHist= new TH1F("_costGamCmFittedHist", "cos(#Theta_{#gamma}) CM Fit", 100, -1, 1);

  _costPhi_KpKmDataHist= new TH1F("_costPhi_KpKmDataHist", "cos(#Theta_{#phi}) K+ K- data", 100, -1., 1.);  
  _costPhi_KpKmMcHist= new TH1F("_costPhi_KpKmMcHist", "cos(#Theta_{#phi}) K+ K- Mc", 100, -1., 1.);   
  _costPhi_KpKmFittedHist= new TH1F("_costPhi_KpKmFittedHist", "cos(#Theta_{#phi}) K+ K- Fit", 100, -1., 1.);   

  _phiPhi_KpKmDataHist= new TH1F("_phiPhi_KpKmDataHist", "cos(#Phi_{#phi}) K+ K- data", 100, -TMath::Pi(), TMath::Pi());  
  _phiPhi_KpKmMcHist= new TH1F("_phiPhi_KpKmMcHist", "cos(#Phi_{#phi}) K+ K- Mc", 100, -TMath::Pi(), TMath::Pi());  
  _phiPhi_KpKmFittedHist= new TH1F("_phiPhi_KpKmFittedHist", "cos(#Phi_{#phi}) K+ K- Fit", 100, -TMath::Pi(), TMath::Pi());

  _chiDataHist= new TH1F("_chiDataHist", "#chi data", 90, 0., 90.);
  _chiMcHist= new TH1F("_chiMcHist", "#chi Mc", 90, 0., 90.); 
  _chiFittedHist= new TH1F("_chiFittedHist", "#chi Fit", 90, 0., 90.);  

   std::string tupleVariables = "mPipPimPi0KpKm:mPipPimPi0:mKpKm:PipPimPi0KpKmCostTheta:gamCosTheta:PipPimPi0CosTheta:KpKmCosTheta:PipCosTheta:KpCosTheta:decPlaneChi:testHeli:weight";
   
  
   _dataTuple=new TNtuple("_dataTuple", "data ntuple", tupleVariables.data());
   _mcTuple=new TNtuple("_mcTuple", "mc ntuple", tupleVariables.data());

   _massIndepTuple = new TNtuple("_massIndepTuple","results from mass indep. fits", ("mass:eta:etaErr:f0:f0Err:f1:f1Err:f2:f2Err:eta2:eta2Err:etaReal:etaImg")   );
}

void JpsiToOmegaPhiGamHist::plotDalitz(TH2F* theHisto,  EvtDataNew* theData, double weight)
{
  Vector4<double>& V4_PipPimPi0_HeliPsi = theData->FourVecsDec[enumJpsiGamX4V::V4_PipPimPi0_HeliPsi]  ;
  Vector4<double>& V4_KpKm_HeliPsi = theData->FourVecsDec[enumJpsiGamX4V::V4_KpKm_HeliPsi] ;
  Vector4<double>& V4_gamma_HeliPsi = theData->FourVecsDec[enumJpsiGamX4V::V4_gamma_HeliPsi] ;
  
  double gphi1 = (V4_gamma_HeliPsi + V4_PipPimPi0_HeliPsi ).M();
  double gphi2 = (V4_gamma_HeliPsi + V4_KpKm_HeliPsi ).M();
  theHisto->Fill( gphi1*gphi1, gphi2*gphi2   ,weight);  
}

void JpsiToOmegaPhiGamHist::plotPhiPhi(TH1F* theHisto, EvtDataNew* theData, double weight){
  Vector4<double>& v4 = theData->FourVecsDec[enumJpsiGamX4V::V4_PipPimPi0KpKm_HeliPsi];
  theHisto->Fill( v4.M(), weight );
}

void JpsiToOmegaPhiGamHist::plotPipPimPi0(TH1F* theHisto, EvtDataNew* theData, double weight){
  Vector4<double>& v4 = theData->FourVecsDec[enumJpsiGamX4V::V4_PipPimPi0_HeliPsi]; 
  theHisto->Fill( v4.M(), weight ); 
}
void JpsiToOmegaPhiGamHist::plotKpKm(TH1F* theHisto, EvtDataNew* theData, double weight){
  Vector4<double>& v4 = theData->FourVecsDec[enumJpsiGamX4V::V4_KpKm_HeliPipPimPi0KpKm]; 
  theHisto->Fill( v4.M(), weight );   
}

void JpsiToOmegaPhiGamHist::plotCostPhiPip(TH1F* theCostHisto,  TH1F* thePhiHisto, EvtDataNew* theData, double weight){
  Vector4<double>& v4 = theData->FourVecsDec[enumJpsiGamX4V::V4_Pip_HeliPipPimPi0]; 
  theCostHisto->Fill( v4.CosTheta(), weight );
  thePhiHisto->Fill( v4.Phi(), weight );
}

void JpsiToOmegaPhiGamHist::plotCostPhiKp(TH1F* theCostHisto,  TH1F* thePhiHisto, EvtDataNew* theData, double weight){
  Vector4<double>& v4 = theData->FourVecsDec[enumJpsiGamX4V::V4_Kp_HeliKpKm]; 
  theCostHisto->Fill( v4.CosTheta(), weight );
  thePhiHisto->Fill( v4.Phi(), weight );
}

void JpsiToOmegaPhiGamHist::plotCostGam(TH1F* theCostHisto, EvtDataNew* theData, double weight){
  Vector4<double>& v4 = theData->FourVecsDec[enumJpsiGamX4V::V4_gamma_HeliPsi]; 
  theCostHisto->Fill( v4.CosTheta(), weight );
}

void JpsiToOmegaPhiGamHist::plotCostPhi_PhiPhiHeli(TH1F* theCostHisto, TH1F* thePhiHisto, const Vector4<double>& the4Vec, double weight){
  theCostHisto->Fill( the4Vec.CosTheta(), weight);
  thePhiHisto->Fill( the4Vec.Phi(), weight);
}

void JpsiToOmegaPhiGamHist::plotChi(TH1F* theChiHisto, EvtDataNew* theData, double weight){
//   Vector4<double>& V4_PipPimPi0KpKm_HeliPsi = theData->FourVecsDec[enumJpsiGamX4V::V4_PipPimPi0KpKm_HeliPsi]  ;
//   Vector4<double>& V4_Pip_HeliPsi= theData->FourVecsDec[enumJpsiGamX4V::V4_Pip_HeliPsi] ;
//   Vector4<double>& V4_Kl_HeliPsi= theData->FourVecsDec[enumJpsiGamX4V::V4_Kl_HeliPsi] ;
//   Vector4<double>& V4_Kp_HeliPsi= theData->FourVecsDec[enumJpsiGamX4V::V4_Kp_HeliPsi] ;
//   Vector4<double>& V4_Km_HeliPsi= theData->FourVecsDec[enumJpsiGamX4V::V4_Km_HeliPsi] ;

//   double thePhiPhiDecayPlaneAngle = decayAngleChi( V4_PipPimPi0KpKm_HeliPsi, V4_Kp_HeliPsi, V4_Km_HeliPsi, V4_Pip_HeliPsi, V4_Kl_HeliPsi   );


  Vector4<double>& V4_normKpKmDecHeliPipPimPi0KpKm = theData->FourVecsDec[enumJpsiGamX4V::V4_normKpKmDecHeliPipPimPi0KpKm];
  Vector4<double>& V4_normPipPimPi0DecHeliPipPimPi0KpKm = theData->FourVecsDec[enumJpsiGamX4V::V4_normPipPimPi0DecHeliPipPimPi0KpKm];

  double cosChi=(V4_normKpKmDecHeliPipPimPi0KpKm.Px()*V4_normPipPimPi0DecHeliPipPimPi0KpKm.Px()
		 +V4_normKpKmDecHeliPipPimPi0KpKm.Py()*V4_normPipPimPi0DecHeliPipPimPi0KpKm.Py()
		 +V4_normKpKmDecHeliPipPimPi0KpKm.Pz()*V4_normPipPimPi0DecHeliPipPimPi0KpKm.Pz())
    / (V4_normKpKmDecHeliPipPimPi0KpKm.P()*V4_normPipPimPi0DecHeliPipPimPi0KpKm.P());

  double chi=acos(fabs(cosChi));
  theChiHisto->Fill(chi*180./TMath::Pi(),weight);   
}

void  JpsiToOmegaPhiGamHist::fillTuple( TNtuple* theTuple, EvtDataNew* theData, double weight)
{
  
  Vector4<double>& V4_PipPimPi0KpKm_HeliPsi = theData->FourVecsDec[enumJpsiGamX4V::V4_PipPimPi0KpKm_HeliPsi]  ;
  Vector4<double>& V4_PipPimPi0_HeliPsi= theData->FourVecsDec[enumJpsiGamX4V::V4_PipPimPi0_HeliPsi] ;
  Vector4<double>& V4_KpKm_HeliPsi= theData->FourVecsDec[enumJpsiGamX4V::V4_KpKm_HeliPsi] ;
  Vector4<double>& V4_gamma_HeliPsi= theData->FourVecsDec[enumJpsiGamX4V::V4_gamma_HeliPsi] ;
  
  Vector4<double>& V4_PipPimPi0_HeliPipPimPi0KpKm= theData->FourVecsDec[enumJpsiGamX4V::V4_PipPimPi0_HeliPipPimPi0KpKm] ;
  Vector4<double>& V4_KpKm_HeliPipPimPi0KpKm= theData->FourVecsDec[enumJpsiGamX4V::V4_KpKm_HeliPipPimPi0KpKm] ;
  Vector4<double>& V4_Pip_HeliPipPimPi0= theData->FourVecsDec[enumJpsiGamX4V::V4_Pip_HeliPipPimPi0] ;
  Vector4<double>& V4_Kp_HeliKpKm= theData->FourVecsDec[enumJpsiGamX4V::V4_Kp_HeliKpKm] ;
  
  Vector4<double>& V4_Pip_HeliPsi= theData->FourVecsDec[enumJpsiGamX4V::V4_Pip_HeliPsi] ;
  Vector4<double>& V4_Kl_HeliPsi= theData->FourVecsDec[enumJpsiGamX4V::V4_Kl_HeliPsi] ;
  Vector4<double>& V4_Kp_HeliPsi= theData->FourVecsDec[enumJpsiGamX4V::V4_Kp_HeliPsi] ;
  Vector4<double>& V4_Km_HeliPsi= theData->FourVecsDec[enumJpsiGamX4V::V4_Km_HeliPsi] ;
  
  
  Vector4<double>& V4_normKpKmDecHeliPipPimPi0KpKm = theData->FourVecsDec[enumJpsiGamX4V::V4_normKpKmDecHeliPipPimPi0KpKm];
  Vector4<double>& V4_normPipPimPi0DecHeliPipPimPi0KpKm = theData->FourVecsDec[enumJpsiGamX4V::V4_normPipPimPi0DecHeliPipPimPi0KpKm];
  double cosChi=(V4_normKpKmDecHeliPipPimPi0KpKm.Px()*V4_normPipPimPi0DecHeliPipPimPi0KpKm.Px()
		 +V4_normKpKmDecHeliPipPimPi0KpKm.Py()*V4_normPipPimPi0DecHeliPipPimPi0KpKm.Py()
		 +V4_normKpKmDecHeliPipPimPi0KpKm.Pz()*V4_normPipPimPi0DecHeliPipPimPi0KpKm.Pz())
    / (V4_normKpKmDecHeliPipPimPi0KpKm.P()*V4_normPipPimPi0DecHeliPipPimPi0KpKm.P());

  double chi=acos(fabs(cosChi));
  double thePhiPhiDecayPlaneAngle = chi/TMath::Pi()*180.;
  
  //double thePhiPhiDecayPlaneAngle = decayAngleChi( V4_PipPimPi0KpKm_HeliPsi, V4_Kp_HeliPsi, V4_Km_HeliPsi, V4_Pip_HeliPsi, V4_Kl_HeliPsi   );
  double testHeli = costDecHeli( V4_PipPimPi0KpKm_HeliPsi+V4_gamma_HeliPsi, V4_Pip_HeliPsi+V4_Kl_HeliPsi+V4_Km_HeliPsi+V4_Kp_HeliPsi, V4_Pip_HeliPsi+V4_Kl_HeliPsi );
  
  theTuple->Fill( V4_PipPimPi0KpKm_HeliPsi.M(),
		  V4_PipPimPi0_HeliPsi.M(),
		  V4_KpKm_HeliPsi.M(),
		  V4_PipPimPi0KpKm_HeliPsi.CosTheta(),
		  V4_gamma_HeliPsi.CosTheta(),
		  V4_PipPimPi0_HeliPipPimPi0KpKm.CosTheta(),
		  V4_KpKm_HeliPipPimPi0KpKm.CosTheta(),
		  V4_Pip_HeliPipPimPi0.CosTheta(),
		  V4_Kp_HeliKpKm.CosTheta(),
		  thePhiPhiDecayPlaneAngle,		  
		  testHeli,
		  weight);
}


double JpsiToOmegaPhiGamHist::decayAngleChi(const Vector4<double>& v4_p,const Vector4<double>& v4_d1,
		const Vector4<double>& v4_d2,const Vector4<double>& v4_h1,
		const Vector4<double>& v4_h2 ) {
  
  TLorentzVector p4_p(  v4_p.Px(), v4_p.Py(), v4_p.Pz(), v4_p.E() );
  TLorentzVector p4_d1p( v4_d1.Px(), v4_d1.Py(), v4_d1.Pz(), v4_d1.E() );
  TLorentzVector p4_d2p( v4_d2.Px(), v4_d2.Py(), v4_d2.Pz(), v4_d2.E() );
  TLorentzVector p4_h1p( v4_h1.Px(), v4_h1.Py(), v4_h1.Pz(), v4_h1.E() );
  TLorentzVector p4_h2p( v4_h2.Px(), v4_h2.Py(), v4_h2.Pz(), v4_h2.E() );
  
  
  // boost all vectors parent restframe
  
  p4_d1p.Boost( -p4_p.BoostVector() );
  p4_d2p.Boost( -p4_p.BoostVector() );
  p4_h1p.Boost( -p4_p.BoostVector() );
  p4_h2p.Boost( -p4_p.BoostVector() );
  
  
  TVector3 d1_perp,d1_prime,h1_perp;
  TVector3 D;
  
  D=(p4_d1p+p4_d2p).Vect();
  
  d1_perp=p4_d1p.Vect()-(D.Dot(p4_d1p.Vect())/D.Dot(D))*D;
  h1_perp=p4_h1p.Vect()-(D.Dot(p4_h1p.Vect())/D.Dot(D))*D;
  
  // orthogonal to both D and d1_perp
  
  d1_prime=D.Cross(d1_perp);
  
  d1_perp= d1_perp* (1./d1_perp.Mag());
  d1_prime= d1_prime * (1./d1_prime.Mag());
  
  double x,y;
 
 x=d1_perp.Dot(h1_perp);
 y=d1_prime.Dot(h1_perp);
 
 double chi=atan2(y,x);
 
 if (chi<0.0) chi+=2.*TMath::Pi();
 
 return chi;
 
}