//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 )
#define TH13 ( 0.0 ) 

//#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[])
{
  //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("proba.root","root",false,true,opts);
  delete treeFactory;

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

  int NumberOfColumn = 3;
  std::vector<std::string> column(NumberOfColumn);
  const char* c_column[] = {"channel","theta23","proba"};
  std::vector<std::string> coltype(NumberOfColumn);
  const char* c_coltype[] = {"int","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","Proba",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);

  double ene = 0.260; //neutrino energy (GeV)
  double baseline = 130.; //baseline (km)
  double theta23;// current value of Theta23
  double sin2th; // current sin2(theta23) (and not sin2(2theta23))
  double sin2ThMin = 0.3;
  double sin2ThMax = 0.7;
  double stepSin2Th = 0.01;
  int muonType = 2; //nu_mu init
  int channel; //oscillation channel
  double probaV; //vaccum oscillation probability  

  glb_params true_values = glbAllocParams();
  //nu_mu->nu_e: 1; ->nu_mu: 2; ->nu_tau: 3
  for (channel = 1; channel<=3 ; channel++) {
    for (sin2th = sin2ThMin; sin2th<=sin2ThMax; sin2th += stepSin2Th ) {
      theta23 = asin(sqrt(sin2th));
      // true values (reference for future Chi2 computations) 
      glbDefineParams(true_values, TH12, TH13, theta23, DELTA, DMQ21, DMQ31);
      glbSetOscillationParameters(true_values);  
      //      glbSetRates();
      probaV = glbVacuumProbability(muonType,channel,+1,ene,baseline);

      //store in Tuple
      myTuple->fill(0,channel);
      myTuple->fill(1,theta23);
      myTuple->fill(2,probaV);
      myTuple->addRow();

    }//loop on sin2th
  }//loop on channel
  
  glbFreeParams(true_values);
  //  glbFreeParams(start_values);

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