//std
//#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <cmath>
#include <values.h>
#include <fstream>
#include <iostream>

//Globes
#include <globes/globes.h>

// AIDA :
#include <AIDA/IAnalysisFactory.h>
#include <AIDA/ITreeFactory.h>
#include <AIDA/IHistogramFactory.h>
#include <AIDA/ITree.h>
#include <AIDA/ITupleFactory.h>
#include <AIDA/ITuple.h>



#define TH12 ( 0.5 * asin(sqrt(0.82)) )
#define TH23 ( M_PI/4.0 )
//theta13 used by Thomas hep-ph/0501037
//#define TH13 ( 0.5 * asin(sqrt(0.05)) )
//theta13 = 1 degree, 3 degree
#define TH13 ( 1.0/180.0 * M_PI) 


#define DELTA ( M_PI/2 )
//#define DELTA ( 0.0 )

#define DMQ21 8.1e-5
#define DMQ31 2.2e-3



/**************************************/
/*           main                     */
/**************************************/

int main(int argc, char *argv[])
{
//   if(argc != 3){
//     std::cerr << "usage:" << argv[0] << "WRONG_TH23 WRONG_HIER [debug]" << std::endl;
//     std::cerr<< "WRONG_* = '1' for wrong * and '0' for right *" << std::endl;
//     return 0;
//   }
  
//   if(sscanf(argv[1], "%d", &WRONG_TH23) != 1 ||
//      (WRONG_TH23 != 1 && WRONG_TH23 != 0) ||
//      sscanf(argv[2], "%d", &WRONG_HIER) != 1 ||
//      (WRONG_HIER != 1 && WRONG_HIER != 0)){
//     std::cerr<< "cannot read parameters" << std::endl;
//     return 0;
//   }

//   int debug = 0;
//   if(4 == argc) {
//     if(sscanf(argv[3],"%d",&debug) !=1 ) {
//       std::cerr<< "cannot read debug" << std::endl;
//       return 0;
//     }
//   }
 
  //init AIDA for Histo/Tuple
  AIDA::IAnalysisFactory* aida = AIDA_createAnalysisFactory();
  if(!aida) {
    std::cerr  << " AIDA not found." << std::endl;
    return 0;
  }

  //ROOT tree :
  AIDA::ITreeFactory* treeFactory = aida->createTreeFactory();
  std::string opts = "export=root";
  AIDA::ITree* fTree = 
    treeFactory->create("bidon.root","root",false,true,opts);
  delete treeFactory;

 //Booking Tuple
  AIDA::ITupleFactory* tf = aida->createTupleFactory(*fTree);

  int NumberOfColumn = 7;
  std::vector<std::string> column(NumberOfColumn);
  const char* c_column[] = {"theta12", "theta13", "theta23", 
			    "deltaCP", "dm21", "dm31", 
			    "chi2"};
  std::vector<std::string> coltype(NumberOfColumn);
  const char* c_coltype[] = {"double","double","double",
			     "double","double","double",
			     "double"};
  for (int icol = 0; icol<NumberOfColumn; ++icol) {
    column[icol]  = c_column[icol];
    coltype[icol] = c_coltype[icol];
  }

  AIDA::ITuple* myTuple = 
    tf->create("SplGlb","dump",column,coltype);
  
 
  //Init Globes: Attention SPL.glb may be a symbolic link
  glbInit(argv[0]);
  glbInitExperiment("../data/SPL.glb", 
		    &glb_experiment_list[0], 
		    &glb_num_of_exps);

  // true values (reference for future Chi2 computations) 
  glb_params true_values = glbAllocParams();
  glbDefineParams(true_values, TH12, TH13, TH23, DELTA, DMQ21, DMQ31);
  std::cout << "The Initial True_values" << std::endl;
  glbPrintParams(stdout,true_values);
  
  glbSetOscillationParameters(true_values);  
  glbSetRates();
  
  // set Starting Values
  glb_params start_values = glbAllocParams();
  glbDefineParams(start_values, TH12, TH13, TH23, DELTA, DMQ21, DMQ31);
  glbSetStartingValues(start_values);
  std::cout << "The Starting values" << std::endl;
  glbPrintParams(stdout,start_values);

  

  //Fluxes
  std::cout << "Flux Name(0):<" << glbValueToName(0,"flux",0) << ">" << std::endl; 
  int splFluxPosFocus = glbNameToValue(0,"flux","#SPLplus");
  std::cout << "Flux Id(#SPLplus):<"<<splFluxPosFocus << ">" << std::endl;
  int splFluxNegFocus = glbNameToValue(0,"flux","#SPLminus");


  std::cout << "Fluxes at 130km for 300MeV neutrinos" << std::endl;
  double energy = 0.300;
  double baseline = 1.0; //the Flux are already defined at 130km...
  std::cout << "Flux(ve)  (+foc): " << glbFlux(0,splFluxPosFocus,energy,baseline,1,1)  << std::endl;
  std::cout << "Flux(vm)  (+foc): " << glbFlux(0,splFluxPosFocus,energy,baseline,2,1)  << std::endl;
  std::cout << "Flux(ave) (+foc): " << glbFlux(0,splFluxPosFocus,energy,baseline,1,-1) << std::endl;
  std::cout << "Flux(avm) (+foc): " << glbFlux(0,splFluxPosFocus,energy,baseline,2,-1) << std::endl;

  std::cout << "Flux(ve)  (-foc): " << glbFlux(0,splFluxNegFocus,energy,baseline,1,1)  << std::endl;
  std::cout << "Flux(vm)  (-foc): " << glbFlux(0,splFluxNegFocus,energy,baseline,2,1)  << std::endl;
  std::cout << "Flux(ave) (-foc): " << glbFlux(0,splFluxNegFocus,energy,baseline,1,-1) << std::endl;
  std::cout << "Flux(avm) (-foc): " << glbFlux(0,splFluxNegFocus,energy,baseline,2,-1) << std::endl;

  //Proba
  double ene[5] = {0.18,0.38,0.58,0.78,0.98};
  std::cout << "Energy, mu->e, mu->mu, mu->tau" << std::endl;
  for (int iene=0; iene<5; ++iene) {
    std::cout << ene[iene] << " " 
	      << glbVacuumProbability(2,1,+1,ene[iene],130.0) << " "
	      << glbVacuumProbability(2,2,+1,ene[iene],130.0) << " "
	      << glbVacuumProbability(2,3,+1,ene[iene],130.0) << " "
	      << std::endl;
  }


 //biodn
  int splChannelNumuNonOsci = glbNameToValue(0,"channel","#numu_non_osci");
  std::cout << "Rate of Channel(#numuNonOsci) NO_OSC: GLB_PRE, W_EFF, W_BG " << std::endl; 
  glbShowChannelRates(stdout,0,splChannelNumuNonOsci,GLB_PRE,GLB_W_EFF,GLB_W_BG); 
  std::cout << std::endl;
  std::cout << "Rate of Channel(#numuNonOsci) NO_OSC: GLB_POST, W_EFF, W_BG " << std::endl; 
  glbShowChannelRates(stdout,0,splChannelNumuNonOsci,GLB_POST,GLB_W_EFF,GLB_W_BG); 
  std::cout << std::endl;


  //Rules


  int splChannelNCplus = glbNameToValue(0,"channel","#NC_bckg_plus");
  std::cout << "Rate of Channel(#NC_bckg_plus) NO_OSC: GLB_PRE, W_EFF, W_BG " << std::endl; 
  glbShowChannelRates(stdout,0,splChannelNCplus,GLB_PRE,GLB_W_EFF,GLB_W_BG); 
  std::cout << std::endl;
  std::cout << "Rate of Channel(#NC_bckg_plus) NO_OSC: GLB_POST, W_EFF, W_BG " << std::endl; 
  glbShowChannelRates(stdout,0,splChannelNCplus,GLB_POST,GLB_W_EFF,GLB_W_BG); 
  std::cout << std::endl;
  



  //NU_E_app
  int splRuleNueApp     = glbNameToValue(0,"rule","#NU_E_app");
  int nItemRuleNUEappSignal = glbGetLengthOfRule(0,splRuleNueApp,GLB_SIG);
  int nItemRuleNUEappBkg    = glbGetLengthOfRule(0,splRuleNueApp,GLB_BG);
  
  std::cout << "Length Of Rule(#NU_E_app) Signal: "<< nItemRuleNUEappSignal
	    << " Bkg: "                            << nItemRuleNUEappBkg
	    << std::endl;

  std::cout << "Normalisation Rule(#NU_E_app) Signal: " 
	    << glbGetNormalizationInRule(0,nItemRuleNUEappSignal,GLB_SIG)
	    << " Bkg: " 
	    << glbGetNormalizationInRule(0,nItemRuleNUEappSignal,GLB_BG)
	    << std::endl;

  std::cout << "Rate of Rules(#NU_E_app) Signal: W_EFF, W_BG, W_COEFF" << std::endl;
  glbShowRuleRates(stdout,0,splRuleNueApp,GLB_ALL,GLB_W_EFF,GLB_W_BG,GLB_W_COEFF,GLB_SIG); 
  std::cout << std::endl;


  std::cout << "Channels of Rule(#NU_E_app) Signal" << std::endl;
  for (int item=0; item<nItemRuleNUEappSignal; ++item) {
    std::cout << "  Channel<" 
	      << glbValueToName(0,"channel",
				glbGetChannelInRule(0,splRuleNueApp,item,GLB_SIG))
	      << "> with coeff. " 
	      << glbGetCoefficientInRule(0,splRuleNueApp,item,GLB_SIG)
	      << std::endl;
    std::cout << "  Rate: WO_EFF, WO_BG, WO_COEFF: "
	      << glbTotalRuleRate(0,splRuleNueApp,item,GLB_WO_EFF,GLB_WO_BG,GLB_WO_COEFF,GLB_SIG)
	      << std::endl;
  }

  std::cout << "Rate of Rules(#NU_E_app) BG: W_EFF, W_BG, W_COEFF" 
	    << std::endl;
  glbShowRuleRates(stdout,0,splRuleNueApp,
		   GLB_ALL,GLB_W_EFF,GLB_W_BG,GLB_W_COEFF,GLB_BG); 

  std::cout << "Channels of Rule(#NU_E_app) Bkg" << std::endl;
  for (int item=0; item<nItemRuleNUEappBkg; ++item) {
    std::cout << "  Channel<" 
	      << glbValueToName(0,"channel",
				glbGetChannelInRule(0,splRuleNueApp,item,GLB_BG))
	      << "> with coeff. " 
	      << glbGetCoefficientInRule(0,splRuleNueApp,item,GLB_BG)
	      << std::endl;
    std::cout << "  Rate: WO_EFF, WO_BG, WO_COEFF: "
	      << glbTotalRuleRate(0,splRuleNueApp,item,GLB_WO_EFF,GLB_WO_BG,GLB_WO_COEFF,GLB_BG)
	      << std::endl;
  }


  //NUMU_DISA: FIXME: problem of the PRE-SMEARING efficiency effects (JEC 6/6/05)
  int splRuleNumudisa         = glbNameToValue(0,"rule","#NUMU_DISA");
  int nItemRuleNUMUDisaSignal = glbGetLengthOfRule(0,splRuleNumudisa,GLB_SIG);
  int nItemRuleNUMUDisaBkg    = glbGetLengthOfRule(0,splRuleNumudisa,GLB_BG);
  
  std::cout << "Length Of Rule(#NUMU_disa) Signal: "<< nItemRuleNUMUDisaSignal
	    << " Bkg: "                             << nItemRuleNUMUDisaBkg
	    << std::endl;

  std::cout << "Normalisation Rule(#NUMU_disa) Signal: " 
	    << glbGetNormalizationInRule(0,nItemRuleNUMUDisaSignal,GLB_SIG)
	    << " Bkg: " 
	    << glbGetNormalizationInRule(0,nItemRuleNUMUDisaSignal,GLB_BG)
	    << std::endl;

  std::cout << "Rate of Rules(#NUMU_disa) Signal: W_EFF, W_BG, W_COEFF" << std::endl;
  glbShowRuleRates(stdout,0,splRuleNumudisa,GLB_ALL,GLB_W_EFF,GLB_W_BG,GLB_W_COEFF,GLB_SIG); 
  std::cout << std::endl;


  std::cout << "Channels of Rule(#NUMU_disa) Signal" << std::endl;
  for (int item=0; item<nItemRuleNUMUDisaSignal; ++item) {
    std::cout << "  Channel<" 
	      << glbValueToName(0,"channel",
				glbGetChannelInRule(0,splRuleNumudisa,item,GLB_SIG))
	      << "> with coeff. " 
	      << glbGetCoefficientInRule(0,splRuleNumudisa,item,GLB_SIG)
	      << std::endl;
    std::cout << "  Rate: WO_EFF, WO_BG, WO_COEFF: "
	      << glbTotalRuleRate(0,splRuleNumudisa,item,GLB_WO_EFF,GLB_WO_BG,GLB_WO_COEFF,GLB_SIG)
	      << std::endl;
  }

  std::cout << "Rate of Rules(#NUMU_disa) BG: W_EFF, W_BG, W_COEFF" 
	    << std::endl;
  glbShowRuleRates(stdout,0,splRuleNumudisa,
		   GLB_ALL,GLB_W_EFF,GLB_W_BG,GLB_W_COEFF,GLB_BG); 

  std::cout << "Channels of Rule(#NUMU_disa) Bkg" << std::endl;
  for (int item=0; item<nItemRuleNUMUDisaBkg; ++item) {
    std::cout << "  Channel<" 
	      << glbValueToName(0,"channel",
				glbGetChannelInRule(0,splRuleNumudisa,item,GLB_BG))
	      << "> with coeff. " 
	      << glbGetCoefficientInRule(0,splRuleNumudisa,item,GLB_BG)
	      << std::endl;
    std::cout << "  Rate: WO_EFF, WO_BG, WO_COEFF: "
	      << glbTotalRuleRate(0,splRuleNumudisa,item,GLB_WO_EFF,GLB_WO_BG,GLB_WO_COEFF,GLB_BG)
	      << std::endl;
  }


  

  //NU_E_BAR_app
  int splRuleAntiNueApp = glbNameToValue(0,"rule","#NU_E_BAR_app");
  int nItemRuleAntiNUEappSignal = glbGetLengthOfRule(0,splRuleAntiNueApp,GLB_SIG);
  int nItemRuleAntiNUEappBkg    = glbGetLengthOfRule(0,splRuleAntiNueApp,GLB_BG);
  std::cout << "Length Of Rule(#NU_E_BAR_app) Signal: "
	    << nItemRuleAntiNUEappSignal
	    << " Bkg: " <<  nItemRuleAntiNUEappBkg
	    << std::endl;

  std::cout << "Normalisation Rule(#NU_E_BAR_app) Signal: " 
	    << glbGetNormalizationInRule(0,nItemRuleAntiNUEappSignal,GLB_SIG)
	    << " Bkg: " 
	    << glbGetNormalizationInRule(0,nItemRuleAntiNUEappSignal,GLB_BG)
	    << std::endl;

  std::cout << "Rate of Rules(#NU_E_BARR_app) Signal: W_EFF, W_BG, W_COEFF" << std::endl;
  glbShowRuleRates(stdout,0,splRuleAntiNueApp,GLB_ALL,GLB_W_EFF,GLB_W_BG,GLB_W_COEFF,GLB_SIG); 
  std::cout << std::endl;

  std::cout << "Channels of Rule(#NU_E_BAR_app) Signal" << std::endl;
  for (int item=0; item<nItemRuleAntiNUEappSignal; ++item) {
    std::cout << "  Channel<" 
	      << glbValueToName(0,"channel",
				glbGetChannelInRule(0,splRuleAntiNueApp,item,GLB_SIG))
	      << "> with coeff. " 
	      << glbGetCoefficientInRule(0,splRuleAntiNueApp,item,GLB_SIG)
	      << std::endl;
    std::cout << "  Rate: WO_EFF, WO_BG, WO_COEFF: "
	      << glbTotalRuleRate(0,splRuleAntiNueApp,item,GLB_WO_EFF,GLB_WO_BG,GLB_WO_COEFF,GLB_SIG)
	      << std::endl;
  }

  std::cout << "Rate of Rules(#NU_E_BARR_app) BG: W_EFF, W_BG, W_COEFF" 
	    << std::endl;
  glbShowRuleRates(stdout,0,splRuleAntiNueApp,
		   GLB_ALL,GLB_W_EFF,GLB_W_BG,GLB_W_COEFF,GLB_BG); 

  std::cout << "Channels of Rule(#NU_E_BAR_app) Bkg" << std::endl;
  for (int item=0; item<nItemRuleAntiNUEappBkg; ++item) {
    std::cout << "  Channel<" 
	      << glbValueToName(0,"channel",
				glbGetChannelInRule(0,splRuleAntiNueApp,item,GLB_BG))
	      << "> with coeff. " 
	      << glbGetCoefficientInRule(0,splRuleAntiNueApp,item,GLB_BG)
	      << std::endl;
    std::cout << "  Rate: WO_EFF, WO_BG, WO_COEFF: "
	      << glbTotalRuleRate(0,splRuleAntiNueApp,item,GLB_WO_EFF,GLB_WO_BG,GLB_WO_COEFF,GLB_BG)
	      << std::endl;
 }

  //NUMUBAR_DISA: FIXME: problem of the PRE-SMEARING efficiency effects (JEC 6/6/05)
  int splRuleNumuBardisa         = glbNameToValue(0,"rule","#NUMUBAR_DISA");
  int nItemRuleNUMUBARDisaSignal = glbGetLengthOfRule(0,splRuleNumuBardisa,GLB_SIG);
  int nItemRuleNUMUBARDisaBkg    = glbGetLengthOfRule(0,splRuleNumuBardisa,GLB_BG);
  
  std::cout << "Length Of Rule(#NUMUBAR_disa) Signal: "<< nItemRuleNUMUBARDisaSignal
	    << " Bkg: "                             << nItemRuleNUMUBARDisaBkg
	    << std::endl;

  std::cout << "Normalisation Rule(#NUMUBAR_disa) Signal: " 
	    << glbGetNormalizationInRule(0,nItemRuleNUMUBARDisaSignal,GLB_SIG)
	    << " Bkg: " 
	    << glbGetNormalizationInRule(0,nItemRuleNUMUBARDisaSignal,GLB_BG)
	    << std::endl;

  std::cout << "Rate of Rules(#NUMUBAR_disa) Signal: W_EFF, W_BG, W_COEFF" << std::endl;
  glbShowRuleRates(stdout,0,splRuleNumuBardisa,GLB_ALL,GLB_W_EFF,GLB_W_BG,GLB_W_COEFF,GLB_SIG); 
  std::cout << std::endl;


  std::cout << "Channels of Rule(#NUMUBAR_disa) Signal" << std::endl;
  for (int item=0; item<nItemRuleNUMUBARDisaSignal; ++item) {
    std::cout << "  Channel<" 
	      << glbValueToName(0,"channel",
				glbGetChannelInRule(0,splRuleNumuBardisa,item,GLB_SIG))
	      << "> with coeff. " 
	      << glbGetCoefficientInRule(0,splRuleNumuBardisa,item,GLB_SIG)
	      << std::endl;
    std::cout << "  Rate: WO_EFF, WO_BG, WO_COEFF: "
	      << glbTotalRuleRate(0,splRuleNumuBardisa,item,GLB_WO_EFF,GLB_WO_BG,GLB_WO_COEFF,GLB_SIG)
	      << std::endl;
  }

  std::cout << "Rate of Rules(#NUMUBAR_disa) BG: W_EFF, W_BG, W_COEFF" 
	    << std::endl;
  glbShowRuleRates(stdout,0,splRuleNumuBardisa,
		   GLB_ALL,GLB_W_EFF,GLB_W_BG,GLB_W_COEFF,GLB_BG); 

  std::cout << "Channels of Rule(#NUMUBAR_disa) Bkg" << std::endl;
  for (int item=0; item<nItemRuleNUMUBARDisaBkg; ++item) {
    std::cout << "  Channel<" 
	      << glbValueToName(0,"channel",
				glbGetChannelInRule(0,splRuleNumuBardisa,item,GLB_BG))
	      << "> with coeff. " 
	      << glbGetCoefficientInRule(0,splRuleNumuBardisa,item,GLB_BG)
	      << std::endl;
    std::cout << "  Rate: WO_EFF, WO_BG, WO_COEFF: "
	      << glbTotalRuleRate(0,splRuleNumuBardisa,item,GLB_WO_EFF,GLB_WO_BG,GLB_WO_COEFF,GLB_BG)
	      << std::endl;
  }

  /*
    fit_theta12 = glbGetOscParams(fit_values, GLB_THETA_12);
    fit_theta13 = glbGetOscParams(fit_values, GLB_THETA_13);
    fit_theta23 = glbGetOscParams(fit_values, GLB_THETA_23);
    fit_deltaCP = glbGetOscParams(fit_values, GLB_DELTA_CP);
    fit_dmSol   = glbGetOscParams(fit_values, GLB_DM_SOL);
    fit_dmAtm   = glbGetOscParams(fit_values, GLB_DM_ATM);
  */
      

  //store result in Tuple
  // myTuple->fill(0,fit_theta12);
//   myTuple->fill(1,fit_theta13);
//   myTuple->fill(2,fit_theta23);
//   myTuple->fill(3,fit_deltaCP);
//   myTuple->fill(4,fit_dmSol);
//   myTuple->fill(5,fit_dmAtm);
//   myTuple->fill(6,res);
//   myTuple->addRow();
  
      

  
  glbFreeParams(true_values);
  glbFreeParams(start_values);

  //save Tuple
  fTree->commit();
  delete fTree;
  
  delete aida;
  return 0;
}