#include <iostream>
#include <iomanip>
#include <string>
#include <TROOT.h>
#include <TChain.h>
#include <TString.h>
#include <TFile.h>
#include <TH2D.h>
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooConstVar.h"
#include "RooKeysPdf.h"
#include "RooNDKeysPdf.h"
#include "RooProdPdf.h"
#include "RooPlot.h"
#include "../CoreAnalysisFiles/CPtree.h"
#include "analyzer.h"
#include <string.h>
#include <ostream>
#include <time.h>

//RooNDKeysPdf *kest4 ;
RooArgSet *observables ;
RooDataSet *data ;
// RooAbsArg *dZ, *eZ ;

using namespace std ;

TH2D *hToy, *hData ;

RooPlot *kmP, *kdZP, *keZP ;
RooPlot *mPlot, *dZPlot, *eZPlot ;
 
TChain *t ;
using namespace std ;
Int_t AddFilesToChain( char* buf , Int_t Run ) 
{
  char* RootFileName = new char[100];

  Int_t fstRun(1),lstRun(6), step(1) ;

  string sbuf (buf) ;

  if( Run != 0 ) { fstRun = Run ; lstRun = Run; }
  for( Int_t iRun = fstRun ; iRun <= lstRun ; iRun++ )
    {
      Int_t aRun = abs(iRun) ;
      Int_t a = sprintf( RootFileName ,"../%s_run%d%s", buf, aRun,".root" ) ;
	  TFile f( RootFileName,"READ" ) ;
	  if( f.IsZombie() ) 
	    { cout << RootFileName << " does not exist " << endl ;break ;}
	  f.Close() ;
	  cout << " Add file " << RootFileName << endl ;
	  t->Add( RootFileName ) ;
    }
  if( t->GetEntries() == 0 )
    { return 0 ; }
  else
    { return 1 ; }
}
void AnaCP( char* Str , Int_t Run, int nevmax )
{
  gDirectory->DeleteAll() ;
  if( t != 0 ) delete t ;
  t = new TChain( "CPtree" ) ;

  // Add file(s) to chain
  Int_t Ok = AddFilesToChain( Str , Run ) ; 
  if( Ok == 0 ) 
    { 
      printf( " Cannot Open Chain for Set %s \n", Str ) ;
      return ;
    }

  cout << " Load Tree " ;
  int s = t->GetEntries() ;
  cout << " entries " << s  ;
  if( s>nevmax ) s=nevmax ;
  cout << " process " << s << " events " << endl ;

  RooRealVar Mnu2("Mnu2","Squared Missing Mass",-12.5,2.5,"GeV/c^{2}") ;
  RooRealVar ZVDec("ZVDec","Z Decay Vertex",-5,5,"cm") ;
  RooRealVar eVDec("eVDec","eZ Decay Vertex",0.,0.05,"cm") ;
  RooRealVar ZVTag("ZVTag","Z Tag Vertex",-5,5,"cm") ;
  RooRealVar eVTag("eVTag","eZ Tag Vertex",0.,0.05,"cm") ;

  RooRealVar dZ("dZ","delta Z",-0.25,0.25 ) ;
  RooRealVar eZ("eZ","sigma delta Z",0.,0.05 ) ;
  RooFormulaVar dZfunc("dZ","ZVDec-ZVTag",RooArgList(ZVDec,ZVTag)) ;
  RooFormulaVar eZfunc("eZ","sqrt(eVDec*eVDec+eVTag*eVTag)",RooArgList(eVDec,eVTag)) ;

  observables = new RooArgSet( Mnu2,ZVDec,eVDec,ZVTag,eVTag ) ;

  data = new RooDataSet( "data","data", (TTree*) t, *observables ) ;
  data->addColumn( dZfunc ) ;
  data->addColumn( eZfunc ) ;

  int ndata = data->numEntries() ;

  cout << " Number of data events " << ndata << endl ;

  dZPlot = dZ.frame(50) ;
  data->plotOn( dZPlot ) ;
  dZPlot->Draw() ;

  eZPlot = eZ.frame(50) ;
  data->plotOn( eZPlot ) ;
  eZPlot->Draw() ;

  mPlot = Mnu2.frame(50) ;
  data->plotOn( mPlot ) ;
  mPlot->Draw() ;

  RooNDKeysPdf kest4("kest4","kest4",RooArgSet(Mnu2,dZ),*data,"am") ;

  kmP = Mnu2.frame(50) ;
  kest4.plotOn( kmP ) ;
  kmP->Draw() ;

  kdZP = dZ.frame(50) ;
  kest4.plotOn( kdZP ) ;
  kdZP->Draw() ;

  // generate toy events
  hToy = new TH2D("hToy","dZ vs MM2 toy",60,-12.5,2.5, 50,-0.25,0.25) ;
  RooRealVar M2gen("M2gen","Generated Squared Missing Mass",-12.5,2.5,"GeV/c^{2}") ;
  RooRealVar dZgen("dZgen","Generated Delta Z",-0.25,0.25,"cm") ;  
  const RooArgSet toy( M2gen, dZgen ) ;
  for( int i=0; i<1000 ; i++ ) 
    {
      kest4.generate( toy, 1 ) ;
      cout << i << "  " << M2gen.getValV() << "  " << dZgen.getValV() << endl ;
      hToy->Fill( M2gen.getValV(), dZgen.getValV() ) ;
    }

  return ;

}
int main ( int narg, char* args[])
{ 
  cout << " " << narg << endl ;
  if( narg<2 ) 
    { AnaCP("b0b0",1,1000000000 ) ;}
  else if( narg<3 ) 
    { AnaCP(args[1],1,1000000000 ) ;}
  else
    {
      cout << " Analyse File "<< args[1] << "  -- Run " <<  args[2] << " -- Events " << args[3] << endl ;
      int irun(1) ;
      sscanf( args[2],"%d",&irun ) ;
      int nevmax ;
      sscanf( args[3],"%d",&nevmax ) ;
      AnaCP(args[1],irun,nevmax) ;
    }
  return 0 ;
}