//std
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <cmath>
#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>

/*
  CP discovery potential (Test delta = 0, and delta = pi)
  Xi2 as the function of the true delta and the true theta_13 to test if the 
  experiment may be fit by CP conserving oscillations 
  the true values are used to build the referenced experiments
  use chi2Delta to marginalized over all parameters except delta
  
  NO degeneracy included
  29/6/05 JEC clean the code
 */

//FIXME: transfert the #defined variables to Data Cards
#define TH12 ( 0.5 * asin(sqrt(0.82)) )
#define TH23 ( M_PI/4.0 )

#define DELTA_0   0.0
#define DELTA_PI  M_PI 

#define DMQ21 8.2e-5
#define DMQ31 2.5e-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("SPLCPDiscovery.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","CP violation discovery",column,coltype);
  
 
  //Init Globes
  glbInit(argv[0]);
  glbInitExperiment("../data/SPL.glb",
		    &glb_experiment_list[0],
		    &glb_num_of_exps);

  // input Prior errors: 
  // 5% on Solar parameters, the atmospheric parameters are 
  // better determined by disappearence channels here 30% to drive the
  // minimizer
  glb_params input_errors = glbAllocParams();
  glbDefineParams(input_errors, 
		  TH12*0.05, 0., TH23*0.30, 0., 
		  DMQ21*0.05, fabs(DMQ31)*0.30);
  glbSetDensityParams(input_errors, 0.05, GLB_ALL);
  glbSetInputErrors(input_errors);

  
  
  //change ND=31 [0,PI/2] to 121 [0,2PI] scan 
#define ND  121   
#define NTH 41
#define LOG10MIN -4.0
#define LOG10MAX -1.0

  std::cout << "Start to loop...." << std::endl;


  double res,res0,resPi;
  double fit_theta12, fit_theta13,fit_theta23,fit_deltaCP,fit_dmSol,fit_dmAtm;
  double delta, log10s2, sq2th, th;


  glb_params true_values = glbAllocParams();
  glb_params start_values = glbAllocParams();
  glb_params test_values = glbAllocParams();
    
  for(int k = 0; k < NTH; k++) {      
          
    log10s2 = LOG10MIN + (LOG10MAX-LOG10MIN) * double(k) / (NTH-1);
    sq2th = pow(10.,log10s2);
    th = 0.5 * asin(sqrt(sq2th));
    
    for(int j = 0; j < ND; j++) {
      //scan de 0 a +Pi pour delta_CP
      delta =  2. * M_PI * double(j) / (ND-1);
      
      // true values (reference for future Chi2 computations) 
      glbDefineParams(true_values, TH12, th, TH23, delta, DMQ21, DMQ31);
      
      glbSetOscillationParameters(true_values);  
      glbSetRates();
      
      // set Starting Values
      glbDefineParams(start_values, TH12, th, TH23, delta, DMQ21, DMQ31);
      glbSetStartingValues(start_values);
      
      // init test values with current TRUE Theta13 but Fixed DELTA=0
      glbDefineParams(test_values, TH12, th, TH23, DELTA_0, DMQ21, DMQ31);

      //check whether the data can be fitted with the CP conserving values
      res0 = glbChiDelta(test_values,NULL,GLB_ALL);

      // init test values with current TRUE Theta13 but Fixed DELTA=PI
      glbDefineParams(test_values, TH12, th, TH23, DELTA_PI, DMQ21, DMQ31);

      //check whether the data can be fitted with the CP conserving values
      resPi = glbChiDelta(test_values,NULL,GLB_ALL);


      //store the true_values theta13,delta_CP, the other variables are there
      //just for compatibility for the Tuple definition
      fit_theta12 = glbGetOscParams(true_values, GLB_THETA_12);
      fit_theta13 = glbGetOscParams(true_values, GLB_THETA_13);
      fit_theta23 = glbGetOscParams(true_values, GLB_THETA_23);
      fit_deltaCP = glbGetOscParams(true_values, GLB_DELTA_CP);
      fit_dmSol   = glbGetOscParams(true_values, GLB_DM_SOL);
      fit_dmAtm   = glbGetOscParams(true_values, GLB_DM_ATM);

      //Get the Minimum Chi2 between the 2 hypothesis Delta=0 Delta=PI
      res = res0;
      if (resPi < res0) res = resPi;

//       std::cout << "Theta_12(rad) " <<  fit_theta12 << " "
// 		<< "Theta_13(rad) " <<  fit_theta13 << " "
//        		<< "Theta_23(deg) " <<  fit_theta23*180./M_PI << " "
//        		<< "Delta_CP(deg) " <<  fit_deltaCP*180./M_PI << " "
//        		<< "DM_21    " <<  fit_dmSol   << " "
//        		<< "DM_31    " <<  fit_dmAtm   << " "
//        		<< "Chi2 (fin) " << res
//        		<< std::endl;       

      //Store in the 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();
      
      
    }//Loop on Theta13
  }//Loop on DeltaCP
  
  //free GLoBES arrays
  glbFreeParams(true_values);
  glbFreeParams(test_values);
  glbFreeParams(start_values);


  //save Tuple
  fTree->commit();

  //Clean Tuple management
  delete fTree;
  delete aida;

  return 0;
}