Index: GLBFrejus/HEAD/data/SPL.glb
===================================================================
--- GLBFrejus/HEAD/data/SPL.glb	(revision 4)
+++ GLBFrejus/HEAD/data/SPL.glb	(revision 4)
@@ -0,0 +1,1 @@
+link SPL11Oct05.glb
Index: GLBFrejus/HEAD/run/TDSyst.C
===================================================================
--- GLBFrejus/HEAD/run/TDSyst.C	(revision 3)
+++ 	(revision )
@@ -1,60 +1,0 @@
-#define TDSyst_cxx
-#include "TDSyst.h"
-#include <TH2D.h>
-#include <TStyle.h>
-#include <TCanvas.h>
-#include <TGraph.h>
-
-#include <iostream>
-
-void TDSyst::Loop()
-{
-//   In a ROOT session, you can do:
-//      Root > .L TDSyst.C
-//      Root > TDSyst t
-//      Root > t.GetEntry(12); // Fill t data members with entry number 12
-//      Root > t.Show();       // Show values of entry 12
-//      Root > t.Show(16);     // Read and show values of entry 16
-//      Root > t.Loop();       // Loop on all entries
-//
-
-//     This is the loop skeleton where:
-//    jentry is the global entry number in the chain
-//    ientry is the entry number in the current Tree
-//  Note that the argument to GetEntry must be:
-//    jentry for TChain::GetEntry
-//    ientry for TTree::GetEntry and TBranch::GetEntry
-//
-//       To read only selected branches, Insert statements like:
-// METHOD1:
-//    fChain->SetBranchStatus("*",0);  // disable all branches
-//    fChain->SetBranchStatus("branchname",1);  // activate branchname
-// METHOD2: replace line
-//    fChain->GetEntry(jentry);       //read all branches
-//by  b_branchname->GetEntry(ientry); //read only this branch
-   if (fChain == 0) return;
-
-   const Double_t PI = TMath::Pi(); 
-
-
-   Long64_t nentries = fChain->GetEntriesFast();
-
-   Double_t sin2Cur = -999.;
-   Double_t deltaCur = -999.;
-
-   Int_t nbytes = 0, nb = 0;
-   for (Long64_t jentry=0; jentry<nentries;jentry++) {
-      Long64_t ientry = LoadTree(jentry);
-      if (ientry < 0) break;
-      nb = fChain->GetEntry(jentry);   nbytes += nb;
-      // if (Cut(ientry) < 0) continue;
-      sin2Cur = TMath::Power(TMath::Sin(2.0 * theta13),2.0);
-      deltaCur= deltaCP/PI*180.0;
-
-//       std::cout << deltaCur << " " << TMath::Log10(sin2Cur)
-// 		<< ": chi2 " << chi2
-// 		<< std::endl;
-      histo->Fill(deltaCur,TMath::Log10(sin2Cur),chi2);
-   }
-}//Loop
-
Index: GLBFrejus/HEAD/run/discover3.C
===================================================================
--- GLBFrejus/HEAD/run/discover3.C	(revision 3)
+++ 	(revision )
@@ -1,83 +1,0 @@
-{
-
-    Int_t nLogSin22q13 = 41;
-   Double_t logSin2Min = -4;
-   Double_t logSin2Max = -1;
-   Double_t logSin2BinWidth = 
-     (logSin2Max - logSin2Min)/double(nLogSin22q13-1.);
-   Double_t logSin2EdgeMax = logSin2Max + logSin2BinWidth/2.;
-   Double_t logSin2EdgeMin = logSin2Min - logSin2BinWidth/2.;
-
-   Int_t nDeltaCP  = 61;
-   Double_t deltaMin = 0.; //degree
-   Double_t deltaMax =  180.;
-   Double_t deltaBinWidth = (deltaMax - deltaMin)/double(nDeltaCP-1.);
-   Double_t deltaEdgeMax = deltaMax + deltaBinWidth/2.;
-   Double_t deltaEdgeMin = deltaMin - deltaBinWidth/2.;
-
-
-
- 
-
-  TCanvas* c1 = new TCanvas("c1","Contour",100,100,600,600);
-  c1->SetGridx();
-  c1->SetGridy();
-  c1->SetLogx();
-
-  TH2D* hdum = new TH2D("hdum","#delta_{CP} 3#sigma discovery",
-			nLogSin22q13,TMath::Power(10.,logSin2Min),TMath::Power(10.,logSin2Max),
-			nDeltaCP,deltaMin,deltaMax
-			);
-  hdum->SetYTitle("#delta_{CP} (deg)");
-  hdum->SetXTitle("sin^{2}(2#theta_{13})");
-  hdum->Draw();
-
-  TPad *overlay = new TPad("overlay","",
-			   gStyle->GetPadLeftMargin(),
-			   gStyle->GetPadBottomMargin(),
-			   1-gStyle->GetPadRightMargin(),
-			   1-gStyle->GetPadTopMargin());
-  overlay->SetFillStyle(4000);
-  overlay->Draw();
-  overlay->Range(logSin2Min,deltaMin,logSin2Max,deltaMax);
-
-
-  TFile f1("./OptiFlux/SPLDeltaCorr3.root");
-  TTree* tree1 = (TTree*)f1.Get("SplGlb");
-  TH2D* histo1 = new TH2D("histo1"," Theta vs Delta",
-			  nLogSin22q13,logSin2EdgeMin,logSin2EdgeMax,
-			  nDeltaCP,deltaEdgeMin,deltaEdgeMax
-			  );
-
-  tree1->Project(histo1->GetName(),"deltaCP/pi*180.:TMath::Log10(pow(sin(2*theta13),2))","chi2");
-
-  histo1->SetContour(1);
-  histo1->SetContourLevel(0,9); //Chi2 for 1 dof and 3sigma
-  histo1->SetLineColor(2); //red
-  overlay->cd();
-  histo1->Draw("CONT3 SAME");
-
-
-
-  TFile f2("./SPLDeltaCorrOldFlux3.root");
-  TTree* tree2 = (TTree*)f2.Get("SplGlb");
-  TH2D* histo2 = new TH2D("histo2"," Theta vs Delta",
-			  nLogSin22q13,logSin2EdgeMin,logSin2EdgeMax,
-			  nDeltaCP,deltaEdgeMin,deltaEdgeMax
-			  );
-  tree2->Project(histo2->GetName(),"deltaCP/pi*180.:TMath::Log10(pow(sin(2*theta13),2))","chi2");
-
-  histo2->SetContour(1);
-  histo2->SetContourLevel(0,9); //Chi2 for 1 dof and 3sigma
-  histo2->SetLineColor(4); //blue
-  overlay->cd();
-  histo2->Draw("CONT3 SAME");
-
-
-
-
-  overlay->Update();
-
-
-
-}
Index: GLBFrejus/HEAD/run/dmtheta.C
===================================================================
--- GLBFrejus/HEAD/run/dmtheta.C	(revision 3)
+++ 	(revision )
@@ -1,70 +1,0 @@
-{
-
-  TCanvas* c1 = new TCanvas("c1","Surf 1",100,100,600,600);
-
-  TFile* file1 = new TFile("./OptiFlux/SPLDmThetaCorr.root");
-  TTree* tree1 = (TTree*)file1->Get("SplGlb");
-  tree1->Draw("log10(dm31):log10(pow(sin(2.0*theta13),2.0))>>h1(41,-3.3325,-1.9675,41,-3.3325,-1.967)","chi2","SURF");
-
-  //  tree1->Draw("log10(dm31):log10(pow(sin(2.0*theta13),2.0))>>h1(21,-3.3065,-1.935,21,-3.3065,-1.935)","chi2","SURF");
-
-
-
-
-  TCanvas* c3 = new TCanvas("c3","Contour",100,100,600,600);
-  c3->SetGridx();
-  c3->SetGridy();
-  c3->SetLogx();
-  c3->SetLogy();
-
-
-  TH2D* hdum = new TH2D("hdum","Theta vs Delta",10,5e-4,1e-2,10,5e-4,1e-2);
-  hdum->SetXTitle("sin^{2}(2#theta_{13})");
-  hdum->SetYTitle("#Delta m^{2}_{31} (eV^{2})");
-  hdum->Draw();
-
-  TPad *overlay = new TPad("overlay","",
-			   gStyle->GetPadLeftMargin(),
-			   gStyle->GetPadBottomMargin(),
-			   1-gStyle->GetPadRightMargin(),
-			   1-gStyle->GetPadTopMargin());
-  overlay->SetFillStyle(4000);
-  overlay->Draw();
-  overlay->Range(-3.3,-3.3,-2.0,-2.0);
-
-
-  TH2F *h1 = (TH2F*)gDirectory->Get("h1");
-  if (h1) {
-    h1->SetContour(1);
-    h1->SetContourLevel(0,4.605); //90%CL 2dof
-    //  h1->SetContourLevel(0,11.83); //3sigma 2dof 
-    h1->SetLineColor(2); //red
-    overlay->cd();
-    h1->Draw("CONT3 SAME");
-  }
-
-  TCanvas* c2 = new TCanvas("c2","Surf 2",100,100,600,600);
-
-  TFile* file2 = new TFile("./SPLDmThetaCorr.root");
-  TTree* tree2 = (TTree*)file2->Get("SplGlb");
-  tree2->Draw("log10(dm31):log10(pow(sin(2.0*theta13),2.0))>>h2(41,-3.3325,-1.9675,41,-3.3325,-1.967)","chi2","SURF");
-
-//   //  tree2->Draw("log10(dm31):log10(pow(sin(2.0*theta13),2.0))>>h2(21,-3.3065,-1.935,21,-3.3065,-1.935)","chi2","SURF");
-
-  TH2F *h2 = (TH2F*)gDirectory->Get("h2");
-  if (h2) {
-    h2->SetContour(1);
-    h2->SetContourLevel(0,4.605); //90%CL 2dof
-    //  h2->SetContourLevel(0,11.83); //3sigma 2dof 
-    h2->SetLineColor(4); //blue
-    overlay->cd();
-    h2->Draw("CONT3 SAME");
-  }
-//   TH2D* h1bis = h1->Clone("h1bis");
-//   h1bis->SetContour(1);
-//   h1bis->SetContourLevel(0,4.605); //90%CL 2dof
-//   h1bis->SetLineColor(4); //blue
-//   h1bis->Draw("CONT3 SAME");
-
-
-}
Index: GLBFrejus/HEAD/run/makeComp.C
===================================================================
--- GLBFrejus/HEAD/run/makeComp.C	(revision 3)
+++ 	(revision )
@@ -1,92 +1,0 @@
-{
-  gROOT->LoadMacro("Compare.C++");
- 
-
-  
-
-
-  TCanvas* c1 = new TCanvas("c1","Contour",100,100,600,600);
-  c1->SetGridx();
-  c1->SetGridy();
-  
-
-  //Original
-  TFile* file1 = new TFile("./Mauro/Original/SPLHierTOctTv2.root");
-  TTree* tree1 = (TTree*)file1->Get("SplGlb");
-  Compare* comp1 = new Compare(tree1,"h1");
-  comp1->Loop();
-  TH2D* h1 = comp1->GetHisto();
-  //  h1->SetTitle("Hier^{tr} Oct^{tr}");
-  h1->SetXTitle("sin^{2}(2#theta_{13})");
-  h1->SetYTitle("#delta_{CP}/#pi");
-  h1->SetContour(1);
-  h1->SetContourLevel(0,4.605); //Chi2 for 2 dof and 90% CL
-  h1->SetLineColor(2); //red
-  h1->Draw("CONT3");
-  
-//   TH2D h1bis; 
-//   h1->Copy(h1bis);
-//   h1bis.SetContour(1);
-//   h1bis.SetContourLevel(0,6.0); //Chi2 for 2 dof and 95% CL
-//   h1bis.SetLineColor(2); //red
-//   h1bis.SetLineStyle(2); //dashed
-//   h1bis.Draw("CONT3 SAME");
-
-
-
-  //Second to be compared
-  TFile* file2 = new TFile("./SPLHierTOctFv2.root");
-  TTree* tree2 = (TTree*)file2->Get("SplGlb");
-  Compare* comp2 = new Compare(tree2,"h2");
-  comp2->Loop();
-  TH2D* h2 = comp2->GetHisto();
-  h2->SetContour(1);
-  h2->SetContourLevel(0,4.605);
-  h2->SetLineStyle(1); //full
-  h2->SetLineColor(4); //blue
-  h2->Draw("CONT3 SAME");
-
-//   TH2D h2bis; 
-//   h2->Copy(h2bis);
-//   h2bis.SetContour(1);
-//   h2bis.SetContourLevel(0,6.0); //Chi2 for 2 dof and 95% CL
-//   h2bis.SetLineColor(4); //blue
-//   h2bis.SetLineStyle(2); //dashed
-//   h2bis.Draw("CONT3 SAME");
-
-
-
-//   //Third to be compared
-//   TFile* file3 = new TFile("./Mauro/Original/SPLHierFOctT.root");
-//   TTree* tree3 = (TTree*)file3->Get("SplGlb");
-//   Compare* comp3 = new Compare(tree3,"h3");
-//   comp3->Loop();
-//   TH2D* h3 = comp3->GetHisto();
-//   h3->SetContour(1);
-//   h3->SetContourLevel(0,4.605);
-//   h3->SetLineColor(3); //green
-//   h3->Draw("CONT3 SAME");
-
-//   //Fourth to be compared
-//   TFile* file4 = new TFile("./Mauro/Original/SPLHierFOctF.root");
-//   TTree* tree4 = (TTree*)file4->Get("SplGlb");
-//   Compare* comp4 = new Compare(tree4,"h4");
-//   comp4->Loop();
-//   TH2D* h4 = comp4->GetHisto();
-//   h4->SetContour(1);
-//   h4->SetContourLevel(0,4.605);
-//   h4->SetLineColor(5); //yellow
-//   h4->Draw("CONT3 SAME");
-
-
-
-  //Draw True sin2(2theta) delta/pi (Hierarchy = T && Octant = T)
-  double x[1] = {0.03};
-  double y[1] = {-0.85};
-  TGraph* true_value = new TGraph(1,x,y);
-  true_value->Draw("*");
-
-
-  c1->Update();
-  
-}
Index: GLBFrejus/HEAD/run/measTh13Delta.data
===================================================================
--- GLBFrejus/HEAD/run/measTh13Delta.data	(revision 3)
+++ 	(revision )
@@ -1,15 +1,0 @@
-WRONG-OCTANT-OPTION     0
-WRONG-HIERAR-OPTION     0
-THETA12-RADIAN          0.594816801
-THETA13-RADIAN          0.087041505
-THETA23-RADIAN          0.785398163
-DELTA-RADIAN            -2.670353756
-DMQ21-EV2               7.92e-5
-DMQ31-EV2               2.4e-3
-ROOT-OUTPUT-FILE        SPLTh13DeltaMeas
-NBIN-DELTA              21
-NBIN-THETA13            21
-MIN-DELTA   		-1.
-MAX-DELTA		0.5
-MIN-SIN22TH13   	0
-MAX-SIN22TH13		0.08
Index: GLBFrejus/HEAD/run/newComp.C
===================================================================
--- GLBFrejus/HEAD/run/newComp.C	(revision 3)
+++ 	(revision )
@@ -1,47 +1,0 @@
-{
-  TCanvas* c1 = new TCanvas("c1","Contour",100,100,600,600);
-  c1->SetGridx();
-  c1->SetGridy();
-  
-  //Reference plot
-  TFile* file1 = new TFile("./Mauro/Original/SPLHierTOctTv2.root");
-  TTree* tree1 = (TTree*)file1->Get("SplGlb");
-  tree1->Draw("deltaCP:sin(2*theta13)^2>>h1(41,-0.001,0.081,31,-1.025,0.525)","chi2","SURF");
-  TH2F *h1 = (TH2F*)gDirectory->Get("h1");
-  h1->SetXTitle("sin^{2}(2#theta_{13})");
-  h1->SetYTitle("#delta_{CP}/#pi");
-  h1->SetContour(1);
-  h1->SetContourLevel(0,4.605);
-  h1->SetLineColor(2); //red
-  h1->Draw("CONT3");
-
-//   TH2D h1bis; 
-//   h1->Copy(h1bis);
-//   h1bis.SetContour(1);
-//   h1bis.SetContourLevel(0,6.0); //Chi2 for 2 dof and 95% CL
-//   h1bis.SetLineColor(2); //red
-//   h1bis.SetLineStyle(2); //dashed
-//   h1bis.Draw("CONT3 SAME");
-
-
-  //Second to be compared
-  TFile* file2 = new TFile("./SPLHierTOctFv2.root");
-  TTree* tree2 = (TTree*)file2->Get("SplGlb");
-  tree2->Draw("deltaCP:sin(2*theta13)^2>>h2(41,-0.001,0.081,31,-1.025,0.525)","chi2","SURF");
-  TH2F *h2 = (TH2F*)gDirectory->Get("h2");
-  h2->SetContour(1);
-  h2->SetContourLevel(0,4.605);
-  h2->SetLineStyle(1); //full
-  h2->SetLineColor(4); //blue
-  h2->Draw("CONT3 SAME");
-
-//   TH2D h2bis; 
-//   h2->Copy(h2bis);
-//   h2bis.SetContour(1);
-//   h2bis.SetContourLevel(0,6.0); //Chi2 for 2 dof and 95% CL
-//   h2bis.SetLineColor(4); //blue
-//   h2bis.SetLineStyle(2); //dashed
-//   h2bis.Draw("CONT3 SAME");
-
-
-}
Index: GLBFrejus/HEAD/run/rate5p.C
===================================================================
--- GLBFrejus/HEAD/run/rate5p.C	(revision 3)
+++ 	(revision )
@@ -1,49 +1,0 @@
-{
- 
-  TCanvas* c1 = new TCanvas("c1","SPL Flux",100,100,600,600);
-  c1->SetLogy();
-
-  TLegend* leg = new TLegend(0.20,0.20,0.45,0.55);
-  leg->SetTextFont(gStyle->GetTextFont());
-  leg->SetTextSize(0.04);
-
-
-
-  TH1D* hdum = new TH1D("hdum","Rates 5yrs (+)",5,0.08,1.08);
-  hdum->SetMaximum(1e3);
-  hdum->SetMinimum(1e-3);
-  hdum->SetXTitle("E_{#nu} (GeV)");
-  hdum->SetYTitle("Number of Events");
-  hdum->Draw();
-  
-  TTree T1("t1","SPL Rate for 5yrs +");
-  T1.ReadFile("./rate5+.dat","ene/F:numuCC/F:nueB/F:nubare/F:pi0/F:nue1d/F:nue3d/F");
-
-  T1.Draw("ene>>hnueB(5,0.08,1.08)","nueB","SAME");
-
-  T1.SetLineColor(6);
-  T1.Draw("ene>>hnumuCC(5,0.08,1.08)","numuCC","SAME");
-  
-  T1.SetLineColor(7);
-  T1.Draw("ene>>hnubare(5,0.08,1.08)","nubare","SAME");
-  
-  T1.SetLineColor(3);
-  T1.Draw("ene>>hpi0(5,0.08,1.08)","pi0","SAME");
-
-  T1.SetLineColor(4);
-  T1.Draw("ene>>hnue3d(5,0.08,1.08)","nue3d","SAME");
-
-  T1.SetLineColor(4);
-  T1.SetLineStyle(2);
-  T1.Draw("ene>>hnue1d(5,0.08,1.08)","nue1d","SAME");
-  
-
-  leg->AddEntry(hnue3d,"#nu_{e} (#theta_{13} = 3^{o})","l");
-  leg->AddEntry(hnue1d,"#nu_{e} (#theta_{13} = 1^{o})","l");
-  leg->AddEntry(hnueB,"#nu_{e} Beam","l");
-  leg->AddEntry(hpi0,"#pi^{0} (NC)","l");
-  leg->AddEntry(hnumuCC,"#mu missId ","l");
-  leg->AddEntry(hnubare,"#bar{#nu_{e}}","l");
-  leg->Draw();
-
-}
Index: GLBFrejus/HEAD/src/degener.cxx
===================================================================
--- GLBFrejus/HEAD/src/degener.cxx	(revision 3)
+++ 	(revision )
@@ -1,291 +1,0 @@
-//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>
-
-
-
-#define TH12 ( asin(sqrt(0.3)) )
-#define TH23 ( asin(sqrt(0.6)) )
-#define TH13 ( 0.5 * asin(sqrt(0.03)) )
-
-#define DELTA (-0.85 * M_PI )
-
-#define DMQ21 8.1e-5
-#define DMQ31 2.2e-3
-
-
-int WRONG_TH23, WRONG_HIER;
-
-
-bool test_th23(glb_params fit_values) {
-  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.))
-	   );
-}//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) {    
-  return test_th23(fit_values) && test_hier(fit_values);
-}
-
-
-/**************************************/
-/*           main                     */
-/**************************************/
-
-int main(int argc, char *argv[])
-{
-  if(argc != 3){
-    std::cerr << "usage:" << argv[0] << "WRONG_TH23 WRONG_HIER" << std::endl;
-    std::cerr<< "WRONG_* = '1' for wrong * and '0' for right *" << std::endl;
-    return 0;
-  }
-  
-  if(sscanf(argv[1], "%d", &WRONG_TH23) != 1 ||
-     (WRONG_TH23 != 1 && WRONG_TH23 != 0) ||
-     sscanf(argv[2], "%d", &WRONG_HIER) != 1 ||
-     (WRONG_HIER != 1 && WRONG_HIER != 0)){
-    std::cerr<< "cannot read parameters" << std::endl;
-    return 0;
-  }
- 
-  //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("degener.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","Theta23 [0,90deg]:Dm31>0 scan",column,coltype);
-  
- 
-  //Init Globes
-  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 (1st trial)
-  glb_params start_values = glbAllocParams();
-  glbDefineParams(start_values, TH12, TH13, 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);
-  std::cout << "The Starting values 1st trail" << std::endl;
-  glbPrintParams(stdout,start_values);
-  glbFreeParams(start_values);
-
-  
-  // init test values
-  glb_params test_values = glbAllocParams();
-  glbDefineParams(test_values, TH12, TH13, TH23, DELTA, 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);
-  std::cout << "The Test_values" << std::endl;
-  glbPrintParams(stdout,test_values);
-
-  glb_params fit_values = glbAllocParams(); //final minimization
-
-  // Set theta23 error very low to constraint the minimizer in the right octant
-  // input errors
-  glb_params input_errors = glbAllocParams();
-  glbDefineParams(input_errors, 
-		  TH12*0.1, 0., TH23*0.05, 0., 
-		  DMQ21*0.05, fabs(DMQ31)*0.05);
-  glbSetDensityParams(input_errors, 0.05, GLB_ALL);
-  glbSetInputErrors(input_errors);
-
-  //double err_fit_theta12 = glbGetOscParams(input_errors, GLB_THETA_12);
-  //double err_fit_theta13 = glbGetOscParams(input_errors, GLB_THETA_13);
-  double err_fit_theta23 = glbGetOscParams(input_errors, GLB_THETA_23);
-  //double err_fit_deltaCP = glbGetOscParams(input_errors, GLB_DELTA_CP);
-  //double err_fit_dmSol   = glbGetOscParams(input_errors, GLB_DM_SOL);
-  double err_fit_dmAtm   = glbGetOscParams(input_errors, GLB_DM_ATM);
-  std::cout << "Theta23 error " << err_fit_theta23/TH23 
-	    << " and Dm31 error " << err_fit_dmAtm/fabs(DMQ31) 
-	    << std::endl;
-  
-  static const int NDM31 =26;  
-  static const int NTH23 =91;
-  static const double DM31MIN = 1e-3;
-  static const double DM31MAX = 5e-3;
-  static const double THETA23MIN =  0.0; //degree
-  static const double THETA23MAX = 90.0;    
-
-  std::cout << "Start to loop...." << std::endl;
-
-  double res;
-  double fit_theta12, fit_theta13,fit_theta23,fit_deltaCP,fit_dmSol,fit_dmAtm;
-  double dm31cur, theta23cur;
-
-  
-  //set a DELTA_CP test value
-  glbSetOscParams(test_values,-0.12*M_PI,GLB_DELTA_CP);
-  glbSetOscParams(test_values,0.5*asin(sqrt(0.025)),GLB_THETA_13);
-
-
-  //Scan in theta23 and Dm31 to find minimum Chi2
-  for(int j = 0; j < NDM31; j++)
-  {
-    dm31cur = DM31MIN 
-      + (DM31MAX-DM31MIN) * double(j)/double(NDM31-1);
-
-    glbSetOscParams(test_values, dm31cur, GLB_DM_ATM);
-
-    for(int k = 0; k < NTH23; k++)
-    {      
-      theta23cur = THETA23MIN
-	+ (THETA23MAX - THETA23MIN) * double(k)/double(NTH23-1);
-      
-      std::cout << "delta j:" << j << " theta13 k:" << k << std::endl;
-      
-      glbSetOscParams(test_values,theta23cur*M_PI/180.0, GLB_THETA_23);
-      
-
-      //Find Chi2 associated
-      res = glbChiSys(test_values,0,GLB_ALL);
-      
-
-      fit_theta12 = glbGetOscParams(test_values, GLB_THETA_12);
-      fit_theta13 = glbGetOscParams(test_values, GLB_THETA_13);
-      fit_theta23 = glbGetOscParams(test_values, GLB_THETA_23);
-      fit_deltaCP = glbGetOscParams(test_values, GLB_DELTA_CP);
-      fit_dmSol   = glbGetOscParams(test_values, GLB_DM_SOL);
-      fit_dmAtm   = glbGetOscParams(test_values, GLB_DM_ATM);
-
-//       if (!test(fit_values)) {
-// 	std::cout << "Problem at (j,k): " << j << " , " << k << std::endl;
-// 	std::cout << "Theta_12 " <<  fit_theta12 << " "
-// 		  << "Theta_13 " <<  fit_theta13 << " "
-// 		  << "Theta_23 " <<  fit_theta23*180./M_PI << " "
-// 		  << "Delta_CP " <<  fit_deltaCP << " "
-// 		  << "DM_21    " <<  fit_dmSol   << " "
-// 		  << "DM_31    " <<  fit_dmAtm   << " "
-// 		  << "Chi2 (mini) " << res
-// 		  << std::endl;
-// 	continue;
-//       }
-
-      
-//       std::cout << "Theta_12 " <<  fit_theta12 << " "
-//  		<< "Theta_13 " <<  fit_theta13 << " "
-// 		<< "Theta_23 " <<  fit_theta23*180./M_PI << " "
-// 		<< "Delta_CP " <<  fit_deltaCP << " "
-// 		<< "DM_21    " <<  fit_dmSol   << " "
-// 		<< "DM_31    " <<  fit_dmAtm   << " "
-// 		<< "Chi2 (mini) " << res
-// 		<< std::endl;       
-
-
-      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
-  
-  glbFreeParams(true_values);
-  glbFreeParams(test_values);
-  glbFreeParams(fit_values);
-
-
-  //save Tuple
-  fTree->commit();
-  delete fTree;
-  
-  delete aida;
-  return 0;
-}
-    
-
-
-
-
Index: GLBFrejus/HEAD/src/discoTheta13.cxx
===================================================================
--- GLBFrejus/HEAD/src/discoTheta13.cxx	(revision 3)
+++ GLBFrejus/HEAD/src/discoTheta13.cxx	(revision 4)
@@ -24,8 +24,8 @@
   experiment may be fit by Theta13 = 0 
   the true values are used to build the referenced experiments
-  use chi2Delta to marginalized over all parameters except delta
-  
+  use chi2Theta to marginalized over all parameters except Theta
   
   11/7/05 JEC creation
+  31/10/05 JEC correct the comment chi2Theta
  */
 
Index: GLBFrejus/HEAD/src/discovery3.cxx
===================================================================
--- GLBFrejus/HEAD/src/discovery3.cxx	(revision 3)
+++ 	(revision )
@@ -1,209 +1,0 @@
-//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;
-}
-    
-
-
-
-
Index: GLBFrejus/HEAD/src/examine.cxx
===================================================================
--- GLBFrejus/HEAD/src/examine.cxx	(revision 3)
+++ 	(revision )
@@ -1,427 +1,0 @@
-//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 )
-//theta13 used by Thomas hep-ph/0501037
-//#define TH13 ( 0.5 * asin(sqrt(0.05)) )
-//theta13 = 1 degree, 3 degree
-#define TH13 ( 1.0/180.0 * M_PI) 
-
-
-#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[])
-{
-//   if(argc != 3){
-//     std::cerr << "usage:" << argv[0] << "WRONG_TH23 WRONG_HIER [debug]" << std::endl;
-//     std::cerr<< "WRONG_* = '1' for wrong * and '0' for right *" << std::endl;
-//     return 0;
-//   }
-  
-//   if(sscanf(argv[1], "%d", &WRONG_TH23) != 1 ||
-//      (WRONG_TH23 != 1 && WRONG_TH23 != 0) ||
-//      sscanf(argv[2], "%d", &WRONG_HIER) != 1 ||
-//      (WRONG_HIER != 1 && WRONG_HIER != 0)){
-//     std::cerr<< "cannot read parameters" << std::endl;
-//     return 0;
-//   }
-
-//   int debug = 0;
-//   if(4 == argc) {
-//     if(sscanf(argv[3],"%d",&debug) !=1 ) {
-//       std::cerr<< "cannot read debug" << std::endl;
-//       return 0;
-//     }
-//   }
- 
-  //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("bidon.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","dump",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);
-
-  // 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);
-
-  
-
-  //Fluxes
-  std::cout << "Flux Name(0):<" << glbValueToName(0,"flux",0) << ">" << std::endl; 
-  int splFluxPosFocus = glbNameToValue(0,"flux","#SPLplus");
-  std::cout << "Flux Id(#SPLplus):<"<<splFluxPosFocus << ">" << std::endl;
-  int splFluxNegFocus = glbNameToValue(0,"flux","#SPLminus");
-
-
-  std::cout << "Fluxes at 130km for 300MeV neutrinos" << std::endl;
-  double energy = 0.300;
-  double baseline = 1.0; //the Flux are already defined at 130km...
-  std::cout << "Flux(ve)  (+foc): " << glbFlux(0,splFluxPosFocus,energy,baseline,1,1)  << std::endl;
-  std::cout << "Flux(vm)  (+foc): " << glbFlux(0,splFluxPosFocus,energy,baseline,2,1)  << std::endl;
-  std::cout << "Flux(ave) (+foc): " << glbFlux(0,splFluxPosFocus,energy,baseline,1,-1) << std::endl;
-  std::cout << "Flux(avm) (+foc): " << glbFlux(0,splFluxPosFocus,energy,baseline,2,-1) << std::endl;
-
-  std::cout << "Flux(ve)  (-foc): " << glbFlux(0,splFluxNegFocus,energy,baseline,1,1)  << std::endl;
-  std::cout << "Flux(vm)  (-foc): " << glbFlux(0,splFluxNegFocus,energy,baseline,2,1)  << std::endl;
-  std::cout << "Flux(ave) (-foc): " << glbFlux(0,splFluxNegFocus,energy,baseline,1,-1) << std::endl;
-  std::cout << "Flux(avm) (-foc): " << glbFlux(0,splFluxNegFocus,energy,baseline,2,-1) << std::endl;
-
-  //Proba
-  double ene[5] = {0.18,0.38,0.58,0.78,0.98};
-  std::cout << "Energy, mu->e, mu->mu, mu->tau" << std::endl;
-  for (int iene=0; iene<5; ++iene) {
-    std::cout << ene[iene] << " " 
-	      << glbVacuumProbability(2,1,+1,ene[iene],130.0) << " "
-	      << glbVacuumProbability(2,2,+1,ene[iene],130.0) << " "
-	      << glbVacuumProbability(2,3,+1,ene[iene],130.0) << " "
-	      << std::endl;
-  }
-
-
- //biodn
-  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;
-
-
-  //Rules
-
-
-  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;
-  
-
-
-
-  //NU_E_app
-  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)
-  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
-  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)
-  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;
-  }
-
-  /*
-    fit_theta12 = glbGetOscParams(fit_values, GLB_THETA_12);
-    fit_theta13 = glbGetOscParams(fit_values, GLB_THETA_13);
-    fit_theta23 = glbGetOscParams(fit_values, GLB_THETA_23);
-    fit_deltaCP = glbGetOscParams(fit_values, GLB_DELTA_CP);
-    fit_dmSol   = glbGetOscParams(fit_values, GLB_DM_SOL);
-    fit_dmAtm   = glbGetOscParams(fit_values, GLB_DM_ATM);
-  */
-      
-
-  //store result in 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();
-  
-      
-
-  
-  glbFreeParams(true_values);
-  glbFreeParams(start_values);
-
-  //save Tuple
-  fTree->commit();
-  delete fTree;
-  
-  delete aida;
-  return 0;
-}
-    
-
-
-
-
Index: GLBFrejus/HEAD/src/measTH13DELTA.cxx
===================================================================
--- GLBFrejus/HEAD/src/measTH13DELTA.cxx	(revision 3)
+++ 	(revision )
@@ -1,320 +1,0 @@
-//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>
-
-
-
-// #define TH12 ( asin(sqrt(0.3)) )
-// #define TH23 ( asin(sqrt(0.6)) )
-// #define TH13 ( 0.5 * asin(sqrt(0.03)) )
-// #define DELTA (-0.85 * M_PI )
-// #define DMQ21 8.1e-5
-// #define DMQ31 2.2e-3
-// #define ND  31  
-// #define NTH 41
-// #define SQ2THMAX 0.08
-
-
-int WRONG_TH23, WRONG_HIER; //condition on octant/hierarchy
-float TH12,TH13,TH23,DELTA,DMQ21,DMQ31; //theta12,theta13,theta23,deltaCP,Dm^2_21,Dm^2_31
-
-int ND; //number of bins to scan delta parameter
-float DELTAMIN, DELTAMAX; //minimum and maximum values of delta to scan (fraction de Pi radian)
-int NTH13; //number of bins to scan theta13 parameter
-float S22TH13MIN, S22TH13MAX; //minimum and maximum values of sin^2(2theta13) to scan
-
-
-bool test_th23(glb_params fit_values) {
-  double fit_th23 = glbGetOscParams(fit_values, GLB_THETA_23);
-  
-  // testing if fit_th23 is at the same side as the true value
-  
-  if (fabs(TH23 - M_PI/4) < 1e-6) return true;
-  
-  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.))
-	   );
-}//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) {    
-  return test_th23(fit_values) && test_hier(fit_values);
-}
-
-
-/**************************************/
-/*           main                     */
-/**************************************/
-
-int main(int argc, char *argv[])
-{
-//   if(argc != 3){
-//     std::cerr << "usage:" << argv[0] << "WRONG_TH23 WRONG_HIER" << std::endl;
-//     std::cerr<< "WRONG_* = '1' for wrong * and '0' for right *" << std::endl;
-//     return 0;
-//   }
-  
-//   if(sscanf(argv[1], "%d", &WRONG_TH23) != 1 ||
-//      (WRONG_TH23 != 1 && WRONG_TH23 != 0) ||
-//      sscanf(argv[2], "%d", &WRONG_HIER) != 1 ||
-//      (WRONG_HIER != 1 && WRONG_HIER != 0)){
-//     std::cerr<< "cannot read parameters" << std::endl;
-//     return 0;
-//   }
- 
-  //Process data card file
-  std::string ROOTDIR = getenv("GLBFREJUSROOT");
-  std::string DATACARD;
-  DATACARD = ROOTDIR + "/run/measTh13Delta.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 >> TH13;  
-  std::cout << DUMMYSTR << " " << TH13 << std::endl; 
-  fileDataCard >> DUMMYSTR >> TH23;  
-  std::cout << DUMMYSTR << " " << TH23 << std::endl; 
-  fileDataCard >> DUMMYSTR >> DELTA;  
-  std::cout << DUMMYSTR << " " << DELTA << 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 >> NTH13;  
-  std::cout << DUMMYSTR << " " << NTH13  << std::endl; 
-  fileDataCard >> DUMMYSTR >> DELTAMIN;  
-  std::cout << DUMMYSTR << " " << DELTAMIN << std::endl; 
-  fileDataCard >> DUMMYSTR >> DELTAMAX;  
-  std::cout << DUMMYSTR << " " << DELTAMAX << std::endl; 
-  fileDataCard >> DUMMYSTR >> S22TH13MIN;  
-  std::cout << DUMMYSTR << " " << S22TH13MIN << std::endl; 
-  fileDataCard >> DUMMYSTR >> S22TH13MAX;  
-  std::cout << DUMMYSTR << " " << S22TH13MAX << 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","Measure SPL Hier=T, Oct=T",column,coltype);
-  
- 
-  //Init Globes
-  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);
-    
-  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);
-  std::cout << "The Starting values 1st trail" << std::endl;
-  glbPrintParams(stdout,start_values);
-  glbFreeParams(start_values);
-
-  
-  // init test values
-  glb_params test_values = glbAllocParams();
-  glbDefineParams(test_values, TH12, TH13, TH23, DELTA, 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);
-  std::cout << "The Test_values" << std::endl;
-  glbPrintParams(stdout,test_values);
-
-  glb_params fit_values = glbAllocParams(); //final minimization
-
-  // Set theta23 error very low to constraint the minimizer in the right octant
-  // input errors
-  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);
-
-
-  std::cout << "Start to loop...." << std::endl;
-
-  double res;
-  double fit_theta12, fit_theta13,fit_theta23,fit_deltaCP,fit_dmSol,fit_dmAtm;
-  double delta, sq2th, th;
-
-  for(int j = 0; j < ND; j++) {
-    delta = DELTAMIN + (DELTAMAX-DELTAMIN)  * double(j) / (ND-1);
-    delta *= M_PI; //rescale in radian
-
-    glbSetOscParams(test_values, delta, GLB_DELTA_CP);
-
-    for(int k = 0; k < NTH13; k++) {      
-      sq2th = S22TH13MIN + (S22TH13MAX - S22TH13MIN) * double(k) / (NTH13-1);
-      th = 0.5 * asin(sqrt(sq2th));
-      
-      std::cout << "delta j:" << j << " theta13 k:" << k << std::endl;
-      
-      glbSetOscParams(test_values, th, GLB_THETA_13);
-      
-
-      //Find minimum in the projection Theta13 Delta_CP
-      res = glbChiThetaDelta(test_values,fit_values,0);
-      
-      //Check if the Parameter Vector respect degeneracy Octant/Hierarchy
-      fit_theta12 = glbGetOscParams(fit_values, GLB_THETA_12);
-      fit_theta13 = glbGetOscParams(fit_values, GLB_THETA_13);
-      fit_theta23 = glbGetOscParams(fit_values, GLB_THETA_23);
-      fit_deltaCP = glbGetOscParams(fit_values, GLB_DELTA_CP);
-      fit_dmSol   = glbGetOscParams(fit_values, GLB_DM_SOL);
-      fit_dmAtm   = glbGetOscParams(fit_values, GLB_DM_ATM);
-
-      if (!test(fit_values)) res += 1e6;
-
-
-      
-//       std::cout << "Theta_12 " <<  fit_theta12 << " "
-//  		<< "Theta_13 " <<  fit_theta13 << " "
-// 		<< "Theta_23 " <<  fit_theta23*180./M_PI << " "
-// 		<< "Delta_CP " <<  fit_deltaCP << " "
-// 		<< "DM_21    " <<  fit_dmSol   << " "
-// 		<< "DM_31    " <<  fit_dmAtm   << " "
-// 		<< "Chi2 (mini) " << res
-// 		<< std::endl;       
-
-
-      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
-  
-  glbFreeParams(true_values);
-  glbFreeParams(test_values);
-  glbFreeParams(fit_values);
-
-
-  //save Tuple
-  fTree->commit();
-  delete fTree;
-  
-  delete aida;
-  return 0;
-}
-    
-
-
-
-
Index: GLBFrejus/HEAD/src/proba.cxx
===================================================================
--- GLBFrejus/HEAD/src/proba.cxx	(revision 3)
+++ 	(revision )
@@ -1,124 +1,0 @@
-//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;
-}
-    
-
-
-
-
