source: GLBFrejus/HEAD/src/dm-theta2.cxx@ 4

Last change on this file since 4 was 2, checked in by campagne, 20 years ago

first import

File size: 5.3 KB
Line 
1//std
2#include <stdio.h>
3#include <stdlib.h>
4#include <math.h>
5#include <string.h>
6#include <cmath>
7#include <fstream>
8#include <iostream>
9
10//Globes
11#include <globes/globes.h>
12
13// AIDA :
14#include <AIDA/IAnalysisFactory.h>
15#include <AIDA/ITreeFactory.h>
16#include <AIDA/IHistogramFactory.h>
17#include <AIDA/ITree.h>
18#include <AIDA/ITupleFactory.h>
19#include <AIDA/ITuple.h>
20
21
22/*
23 Discovery plot dm31 and theta_13
24 Theta_13=0 and Delta=0 for the true value of the reference experiment
25 use glbChiNP to marginalized other all parameters except of
26 the current dm31 and theta (delta=0 is also fixed)*/
27
28#define TH12 ( 0.5 * asin(sqrt(0.82)) )
29#define TH23 ( M_PI/4.0 )
30
31#define TH13 0.0
32#define DELTA 0.0
33
34#define DMQ21 8.2e-5
35#define DMQ31 2.5e-3
36
37
38/**************************************/
39/* main */
40/**************************************/
41
42int main(int argc, char *argv[])
43{
44 //init AIDA for Histo/Tuple
45 AIDA::IAnalysisFactory* aida = AIDA_createAnalysisFactory();
46 if(!aida) {
47 std::cerr << " AIDA not found." << std::endl;
48 return 0;
49 }
50
51 //ROOT tree :
52 AIDA::ITreeFactory* treeFactory = aida->createTreeFactory();
53 std::string opts = "export=root";
54 AIDA::ITree* fTree =
55 treeFactory->create("SPLDmThetaCorr.root","root",false,true,opts);
56 delete treeFactory;
57
58 //Booking Tuple
59 AIDA::ITupleFactory* tf = aida->createTupleFactory(*fTree);
60
61 int NumberOfColumn = 7;
62 std::vector<std::string> column(NumberOfColumn);
63 const char* c_column[] = {"theta12", "theta13", "theta23",
64 "deltaCP", "dm21", "dm31",
65 "chi2"};
66 std::vector<std::string> coltype(NumberOfColumn);
67 const char* c_coltype[] = {"double","double","double",
68 "double","double","double",
69 "double"};
70 for (int icol = 0; icol<NumberOfColumn; ++icol) {
71 column[icol] = c_column[icol];
72 coltype[icol] = c_coltype[icol];
73 }
74
75 AIDA::ITuple* myTuple =
76 tf->create("SplGlb","dm31 - theta13",column,coltype);
77
78
79 //Init Globes
80 glbInit(argv[0]);
81 glbInitExperiment("../data/SPL.glb",
82 &glb_experiment_list[0],
83 &glb_num_of_exps);
84
85 // input errors
86 glb_params input_errors = glbAllocParams();
87 glbDefineParams(input_errors,
88 TH12*0.05, 0., TH23*0.30, 0.,
89 DMQ21*0.05, fabs(DMQ31)*0.30);
90 glbSetDensityParams(input_errors, 0.05, GLB_ALL);
91 glbSetInputErrors(input_errors);
92
93#define NDM 41
94#define LOG10DMMIN -3.3
95#define LOG10DMMAX -2.0
96
97#define NTH 41
98#define LOG10THMIN -3.3
99#define LOG10THMAX -2.0
100
101
102 glb_params true_values = glbAllocParams();
103 glb_params start_values = glbAllocParams();
104 glb_params test_values = glbAllocParams();
105
106 //Projection. Theta13-Dm31 (et deltaCP fixed also)
107 glb_projection proj = glbAllocProjection();
108 glbDefineProjection(proj, GLB_FREE, GLB_FIXED,GLB_FREE,
109 GLB_FIXED, GLB_FREE, GLB_FIXED);
110 glbSetProjection(proj);
111
112 std::cout << "Start to loop...." << std::endl;
113
114
115 double res;
116 double fit_theta12, fit_theta13,fit_theta23,fit_deltaCP,fit_dmSol,fit_dmAtm;
117 double log10s2, sq2th, th;
118 double log10dm, dm;
119
120 for(int j = 0; j < NDM; j++)
121 {
122 //scan de DM31
123 log10dm = LOG10DMMIN + (LOG10DMMAX - LOG10DMMIN) * double(j) / (NDM-1);
124 dm = pow(10.,log10dm);
125
126
127 for(int k = 0; k < NTH; k++)
128 {
129 //scan de THETA13
130 log10s2 = LOG10THMIN + (LOG10THMAX-LOG10THMIN) * double(k) / (NTH-1);
131 sq2th = pow(10.,log10s2);
132 th = 0.5 * asin(sqrt(sq2th));
133
134 std::cout << "dm31 j:" << j << " theta13 k:" << k << std::endl;
135 glbDefineParams(true_values, TH12, th, TH23, DELTA, DMQ21, dm);
136 glbSetOscillationParameters(true_values);
137 glbSetRates();
138
139
140 // set Starting Values
141 glbDefineParams(start_values, TH12, th, TH23, DELTA, DMQ21, dm);
142 glbSetStartingValues(start_values);
143
144 // init test values with Theta13=0
145 glbDefineParams(test_values, TH12, TH13, TH23, DELTA, DMQ21, dm);
146 // check whether the data can be fitted with Theta13=0 (delta=0 fixed)
147 res = glbChiNP(test_values,NULL,GLB_ALL);
148
149 //store the true_values theta13,dm31, the other variables are there
150 //just for compatibility for the Tuple definition
151 fit_theta12 = glbGetOscParams(true_values, GLB_THETA_12);
152 fit_theta13 = glbGetOscParams(true_values, GLB_THETA_13);
153 fit_theta23 = glbGetOscParams(true_values, GLB_THETA_23);
154 fit_deltaCP = glbGetOscParams(true_values, GLB_DELTA_CP);
155 fit_dmSol = glbGetOscParams(true_values, GLB_DM_SOL);
156 fit_dmAtm = glbGetOscParams(true_values, GLB_DM_ATM);
157
158
159// std::cout << "Theta_12(rad) " << fit_theta12 << " "
160// << "Theta_13(rad) " << fit_theta13 << " "
161// << "Theta_23(deg) " << fit_theta23*180./M_PI << " "
162// << "Delta_CP(deg) " << fit_deltaCP*180./M_PI << " "
163// << "DM_21 " << fit_dmSol << " "
164// << "DM_31 " << fit_dmAtm << " "
165// << "Chi2 (fin) " << res
166// << std::endl;
167
168
169 //Store in the Tuple
170 myTuple->fill(0,fit_theta12);
171 myTuple->fill(1,fit_theta13);
172 myTuple->fill(2,fit_theta23);
173 myTuple->fill(3,fit_deltaCP);
174 myTuple->fill(4,fit_dmSol);
175 myTuple->fill(5,fit_dmAtm);
176 myTuple->fill(6,res);
177 myTuple->addRow();
178
179
180
181 }//Loop on Theta13
182 }//Loop on Dm31
183
184 //free GLoBES arrays
185 glbFreeParams(true_values);
186 glbFreeParams(test_values);
187 glbFreeParams(start_values);
188
189 //save Tuple
190 fTree->commit();
191 delete fTree;
192
193 delete aida;
194 return 0;
195}
196
197
198
199
200
Note: See TracBrowser for help on using the repository browser.