//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 TH13 (0.5*asin(sqrt(0.1))) 
#define TH23 ( M_PI/4.0 )

#define DELTA ( M_PI/2 )

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

#define ENERGY 0.300   //GeV
#define NSTEP  1000    
#define MAXBASELINE 1000.0 //km 


double max(double x) {
// #define MINVALUE 1e-10
//   return (x > MINVALUE ? x : MINVALUE);
  return x;
}


/**************************************/
/*           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 = 25;
  std::vector<std::string> column(NumberOfColumn);
  const char* c_column[] = { "Pee","Pemu","Petau",
			     "Pmue","Pmumu","Pmutau",
			     "Pbee","Pbemu","Pbetau",
			     "Pbmue","Pbmumu","Pbmutau",
			     "PMee","PMemu","PMetau",
			     "PMmue","PMmumu","PMmutau",
			     "PMbee","PMbemu","PMbetau",
			     "PMbmue","PMbmumu","PMbmutau",
			     "LovE" };
  std::vector<std::string> coltype(NumberOfColumn);
  const char* c_coltype[] = {"double","double","double",
			     "double","double","double",
			     "double","double","double",
			     "double","double","double",
			     "double","double","double",
			     "double","double","double",
			     "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("ProbaOscill","dump",column,coltype);
  
 
  //Init Globes: Attention SPL.glb may be a symbolic link
  glbInit(argv[0]);
  glbSetVerbosityLevel(2);
  glbInitExperiment("../data/Bidon.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);
  glbSetOscillationParameters(true_values);  
  glbSetRates();

  std::cout << "Profile  of exp 0: " << glbGetProfileTypeInExperiment(0) << std::endl;


  double baseLine;
  //Vacuum
  //Proba with neutrinos
  double Pee;
  double Pemu;
  double Petau;
  double Pmue;
  double Pmumu;
  double Pmutau;

  //Proba with anti-neutrinos
  double Pbee;
  double Pbemu;
  double Pbetau;
  double Pbmue;
  double Pbmumu;
  double Pbmutau;

  //Matter
  //Proba with neutrinos
  double PMee;
  double PMemu;
  double PMetau;
  double PMmue;
  double PMmumu;
  double PMmutau;

  //PMroba with anti-neutrinos
  double PMbee;
  double PMbemu;
  double PMbetau;
  double PMbmue;
  double PMbmumu;
  double PMbmutau;

  //Proba
  for (int iLen=1; iLen<NSTEP; ++iLen) {
    baseLine = iLen/double(NSTEP) * MAXBASELINE;
    
    //Vaccum
    Pee     = glbVacuumProbability(1,1,+1,ENERGY,baseLine);
    Pemu    = glbVacuumProbability(1,2,+1,ENERGY,baseLine);
    Petau   = glbVacuumProbability(1,3,+1,ENERGY,baseLine);
    Pmue    = glbVacuumProbability(2,1,+1,ENERGY,baseLine);
    Pmumu   = glbVacuumProbability(2,2,+1,ENERGY,baseLine);
    Pmutau  = glbVacuumProbability(2,3,+1,ENERGY,baseLine);

    Pbee    = glbVacuumProbability(1,1,-1,ENERGY,baseLine);
    Pbemu   = glbVacuumProbability(1,2,-1,ENERGY,baseLine);
    Pbetau  = glbVacuumProbability(1,3,-1,ENERGY,baseLine);
    Pbmue   = glbVacuumProbability(2,1,-1,ENERGY,baseLine);
    Pbmumu  = glbVacuumProbability(2,2,-1,ENERGY,baseLine);
    Pbmutau = glbVacuumProbability(2,3,-1,ENERGY,baseLine);
    
    //Matter
    glbSetBaselineInExperiment(0,baseLine); 

    PMee     = glbProfileProbability(0,1,1,+1,ENERGY);
    PMemu    = glbProfileProbability(0,1,2,+1,ENERGY);
    PMetau   = glbProfileProbability(0,1,3,+1,ENERGY);
    PMmue    = glbProfileProbability(0,2,1,+1,ENERGY);
    PMmumu   = glbProfileProbability(0,2,2,+1,ENERGY);
    PMmutau  = glbProfileProbability(0,2,3,+1,ENERGY);

    PMbee    = glbProfileProbability(0,1,1,-1,ENERGY);
    PMbemu   = glbProfileProbability(0,1,2,-1,ENERGY);
    PMbetau  = glbProfileProbability(0,1,3,-1,ENERGY);
    PMbmue   = glbProfileProbability(0,2,1,-1,ENERGY);
    PMbmumu  = glbProfileProbability(0,2,2,-1,ENERGY);
    PMbmutau = glbProfileProbability(0,2,3,-1,ENERGY);
        
    //save into Tuple
    myTuple->fill(0,max(Pee));
    myTuple->fill(1,max(Pemu));
    myTuple->fill(2,max(Petau));
    myTuple->fill(3,max(Pmue));
    myTuple->fill(4,max(Pmumu));
    myTuple->fill(5,max(Pmutau));

    myTuple->fill(6,max(Pbee));
    myTuple->fill(7,max(Pbemu));
    myTuple->fill(8,max(Pbetau));
    myTuple->fill(9,max(Pbmue));
    myTuple->fill(10,max(Pbmumu));
    myTuple->fill(11,max(Pbmutau));

    myTuple->fill(12,max(PMee));
    myTuple->fill(13,max(PMemu));
    myTuple->fill(14,max(PMetau));
    myTuple->fill(15,max(PMmue));
    myTuple->fill(16,max(PMmumu));
    myTuple->fill(17,max(PMmutau));

    myTuple->fill(18,max(PMbee));
    myTuple->fill(19,max(PMbemu));
    myTuple->fill(20,max(PMbetau));
    myTuple->fill(21,max(PMbmue));
    myTuple->fill(22,max(PMbmumu));
    myTuple->fill(23,max(PMbmutau));

    myTuple->fill(24,baseLine/ENERGY);
    
    myTuple->addRow();
  }
  //Clear GLobES
  glbFreeParams(true_values);

  //save Tuple
  fTree->commit();

  //Clear OS
  delete fTree;
  delete aida;
  return 0;
}