source: GLBFrejus/HEAD/src/discoCP.cxx @ 150

//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
  
  
  29/6/05 JEC clean the code
  7/7/05  JEC introduce the Mass hierachy and Octant degeneracy from Th.Schwetz
  8/7/05  JEC introduce data card management
  17/5/06 JEC add  TEST_OCTANT to avoid pb if theta23=pi/4. 
 */


#define DELTA_0   0.0
#define DELTA_PI  M_PI 
#define TEST_OCTANT false


int WRONG_TH23, WRONG_HIER; //condition on octant/hierarchy
float TH12, TH23,DMQ21,DMQ31; //theta12,theta23,Dm^2_21,Dm^2_31

int ND; //number of bins to scan delta parameter [0,2pi]
int NTH; //number of bins to scan theta13 parameter
float LOG10MIN; //minimal Log10(sin^2(2theta13))
float LOG10MAX; //minimal Log10(sin^2(2theta13))

//---------------------------------------------------------------------
bool test_th23(glb_params fit_values) {
  if(TEST_OCTANT) {
    double fit_th23 = glbGetOscParams(fit_values, GLB_THETA_23);
    
    // testing if fit_th23 is at the same side as the true value
    
    return  ( 
             (WRONG_TH23 == 0 && 
              (TH23 > M_PI/4. && fit_th23 > M_PI/4.) ||
              (TH23 < M_PI/4. && fit_th23 < M_PI/4.))
             ||
             (WRONG_TH23 == 1 && 
              (TH23 > M_PI/4. && fit_th23 < M_PI/4.) ||
              (TH23 < M_PI/4. && fit_th23 > M_PI/4.))
             );
  }
  return true;
}//same_th23

//---------------------------------------------------------------------
bool test_hier(glb_params fit_values) {
  double fit_dmq = glbGetOscParams(fit_values, GLB_DM_ATM);
  
  // testing if fit_dmq has the same sign as the true value
  return ( 
          (WRONG_HIER == 0 && fit_dmq * DMQ31 > 0) || 
          (WRONG_HIER == 1 && fit_dmq * DMQ31 < 0) 
          );
}//same_hier

//---------------------------------------------------------------------

bool test(glb_params fit_values) {    
  // testing Octant and Hierarchy
  return test_th23(fit_values) && test_hier(fit_values);
}//test

//---------------------------------------------------------------------

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

int main(int argc, char *argv[]) {

  //Process data card file
  std::string ROOTDIR = getenv("GLBFREJUSROOT");
  std::string DATACARD;
  DATACARD = ROOTDIR + "/run/discoCP.data"; 
  std::ifstream fileDataCard;
  std::cout << "Open DataCard file : " << DATACARD <<std::endl;
  fileDataCard.open(DATACARD.data());
  std::string DUMMYSTR; //comment string of the variables

  std::string ROOTFILENAME;


  fileDataCard >> DUMMYSTR >> WRONG_TH23; 
  std::cout << DUMMYSTR << " " << WRONG_TH23 << std::endl; 
  fileDataCard >> DUMMYSTR >> WRONG_HIER;
  std::cout << DUMMYSTR << " " << WRONG_HIER << std::endl; 
  fileDataCard >> DUMMYSTR >> TH12;  
  std::cout << DUMMYSTR << " " << TH12 << std::endl; 
  fileDataCard >> DUMMYSTR >> TH23;  
  std::cout << DUMMYSTR << " " << TH23 << std::endl; 
  fileDataCard >> DUMMYSTR >> DMQ21;  
  std::cout << DUMMYSTR << " " << DMQ21 << std::endl; 
  fileDataCard >> DUMMYSTR >> DMQ31;  
  std::cout << DUMMYSTR << " " << DMQ31 << std::endl; 
  fileDataCard >> DUMMYSTR >> ROOTFILENAME;  

  //add to ROOT file name the conditions on the Octant/Hierarchy
  ROOTFILENAME += "-h";
  char hier[40];
  sprintf(hier,"%d",WRONG_HIER);
  ROOTFILENAME += hier;
  char octo[40];
  sprintf(octo,"%d",WRONG_TH23);
  ROOTFILENAME += "-o";
  ROOTFILENAME += octo;
  ROOTFILENAME += ".root";
  std::cout << DUMMYSTR <<ROOTFILENAME  << std::endl; 

  fileDataCard >> DUMMYSTR >> ND;  
  std::cout << DUMMYSTR << " "  << ND  << std::endl; 
  fileDataCard >> DUMMYSTR >> NTH;  
  std::cout << DUMMYSTR << " " << NTH  << std::endl; 
  fileDataCard >> DUMMYSTR >> LOG10MIN;  
  std::cout << DUMMYSTR << " " << LOG10MIN << std::endl; 
  fileDataCard >> DUMMYSTR >> LOG10MAX;  
  std::cout << DUMMYSTR << " "  << LOG10MAX << std::endl; 



  
  //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(ROOTFILENAME,"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();
  //JEC 17-05-06 use 10% except for Dm2_21 4%
//   glbDefineParams(input_errors, 
//                TH12*0.05, 0., TH23*0.30, 0., 
//                DMQ21*0.05, fabs(DMQ31)*0.30);
  glbDefineParams(input_errors, 
                  TH12*0.10, 0., TH23*0.10, 0., 
                  DMQ21*0.04, fabs(DMQ31)*0.10);
  glbSetDensityParams(input_errors, 0.05, GLB_ALL);
  glbSetInputErrors(input_errors);

  
  

  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();
  glb_params fit_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 +2Pi pour delta_CP
      delta =  2.0 * M_PI * double(j) / (ND-1);

      std::cout << "Consider (k,j) = (" << k << "," << j << ")" << std::endl;
      
      // 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);
      if(WRONG_HIER == 1)
        glbSetOscParams(start_values, -DMQ31+DMQ21, GLB_DM_ATM);
      if(WRONG_TH23 == 1)
        glbSetOscParams(start_values, 0.5*M_PI - TH23, GLB_THETA_23);

      glbSetStartingValues(start_values);
      
      // init test values with current TRUE Theta13 but Fixed DELTA=0
      glbDefineParams(test_values, TH12, th, TH23, DELTA_0, DMQ21, DMQ31);
      if(WRONG_HIER == 1)
        glbSetOscParams(test_values, -DMQ31+DMQ21, GLB_DM_ATM);
      
      if(WRONG_TH23 == 1)
        glbSetOscParams(test_values, 0.5*M_PI - TH23, GLB_THETA_23);
      

      //check whether the data can be fitted with the CP conserving values
      res0 = glbChiDelta(test_values,fit_values,GLB_ALL);
      //Add a big factor in case if fit_values does not respect the Octant &
      //Mass Hierarchy of the true values
      if (!test(fit_values)) res0 += 1e6; 

      // init test values with current TRUE Theta13 but Fixed DELTA=PI
      glbDefineParams(test_values, TH12, th, TH23, DELTA_PI, DMQ21, DMQ31);
      if(WRONG_HIER == 1)
        glbSetOscParams(test_values, -DMQ31+DMQ21, GLB_DM_ATM);
      
      if(WRONG_TH23 == 1)
        glbSetOscParams(test_values, 0.5*M_PI - TH23, GLB_THETA_23);

      //check whether the data can be fitted with the CP conserving values
      resPi = glbChiDelta(test_values,fit_values,GLB_ALL);
      //Add a big factor in case if fit_values does not respect the Octant &
      //Mass Hierarchy of the true values
      if (!test(fit_values)) resPi += 1e6; 

      //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);
  glbFreeParams(fit_values);


  //save Tuple
  fTree->commit();

  //Clean Tuple management
  delete fTree;
  delete aida;

  return 0;
}