//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 ( asin(sqrt(0.3)) )
#define TH23 ( M_PI/4.0 )
//theta13 = 1 degree, 3 degree
#define TH13 ( 1.0/180.0 * M_PI) 
//#define TH13 ( 3.0/180.0 * M_PI) 
//#define TH13 ( 0.5 * asin(sqrt(0.05)) )

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

#define DMQ21 7.9e-5
#define DMQ31 2.4e-3



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

int main(int argc, char *argv[])
{
 
  //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);

  
  //Rates of some channels

  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;




  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;
  

  //Rules

  //NU_E_app
  std::cout << "----------------------------------------------" << std::endl;
  std::cout << ">>>>>>>>>>>> NU_E_APP <<<<<<<<<<<<<<<<<<" << std::endl;
  std::cout << "----------------------------------------------" << std::endl;

  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)
  std::cout << "----------------------------------------------" << std::endl;
  std::cout << ">>>>>>>>>>>> NU_MU_DISA <<<<<<<<<<<<<<<<<<" << std::endl;
  std::cout << "----------------------------------------------" << std::endl;

  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
  std::cout << "----------------------------------------------" << std::endl;
  std::cout << ">>>>>>>>>>>> NU_E__BARR_APP <<<<<<<<<<<<<<<<<<" << std::endl;
  std::cout << "----------------------------------------------" << std::endl;

  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)
  std::cout << "----------------------------------------------" << std::endl;
  std::cout << ">>>>>>>>>>>> NU_MU_BAR_DISA <<<<<<<<<<<<<<<<<<" << std::endl;
  std::cout << "----------------------------------------------" << std::endl;

  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;
  }

      
      

  
  glbFreeParams(true_values);
  glbFreeParams(start_values);

  return 0;
}