source: GLBFrejus/HEAD/src/measTH13DELTA.cxx@ 516

Last change on this file since 516 was 150, checked in by campagne, 19 years ago

add new appli
  • Property svn:executable set to *
File size: 9.6 KB
RevLine 
[150]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// #define TH12 ( asin(sqrt(0.3)) )
24// #define TH23 ( asin(sqrt(0.6)) )
25// #define TH13 ( 0.5 * asin(sqrt(0.03)) )
26// #define DELTA (-0.85 * M_PI )
27// #define DMQ21 8.1e-5
28// #define DMQ31 2.2e-3
29// #define ND 31
30// #define NTH 41
31// #define SQ2THMAX 0.08
32
33
34int WRONG_TH23, WRONG_HIER; //condition on octant/hierarchy
35float TH12,TH13,TH23,DELTA,DMQ21,DMQ31; //theta12,theta13,theta23,deltaCP,Dm^2_21,Dm^2_31
36
37int ND; //number of bins to scan delta parameter
38float DELTAMIN, DELTAMAX; //minimum and maximum values of delta to scan (fraction de Pi radian)
39int NTH13; //number of bins to scan theta13 parameter
40float S22TH13MIN, S22TH13MAX; //minimum and maximum values of sin^2(2theta13) to scan
41
42
43bool test_th23(glb_params fit_values) {
44 double fit_th23 = glbGetOscParams(fit_values, GLB_THETA_23);
45
46 // testing if fit_th23 is at the same side as the true value
47
48 if (fabs(TH23 - M_PI/4) < 1e-6) return true;
49
50 return (
51 (WRONG_TH23 == 0 &&
52 (TH23 > M_PI/4. && fit_th23 > M_PI/4.) ||
53 (TH23 < M_PI/4. && fit_th23 < M_PI/4.))
54 ||
55 (WRONG_TH23 == 1 &&
56 (TH23 > M_PI/4. && fit_th23 < M_PI/4.) ||
57 (TH23 < M_PI/4. && fit_th23 > M_PI/4.))
58 );
59}//same_th23
60
61//---------------------------------------------------------------------
62bool test_hier(glb_params fit_values) {
63 double fit_dmq = glbGetOscParams(fit_values, GLB_DM_ATM);
64
65 // testing if fit_dmq has the same sign as the true value
66 return (
67 (WRONG_HIER == 0 && fit_dmq * DMQ31 > 0) ||
68 (WRONG_HIER == 1 && fit_dmq * DMQ31 < 0)
69 );
70}//same_hier
71
72//---------------------------------------------------------------------
73
74bool test(glb_params fit_values) {
75 return test_th23(fit_values) && test_hier(fit_values);
76}
77
78
79/**************************************/
80/* main */
81/**************************************/
82
83int main(int argc, char *argv[])
84{
85// if(argc != 3){
86// std::cerr << "usage:" << argv[0] << "WRONG_TH23 WRONG_HIER" << std::endl;
87// std::cerr<< "WRONG_* = '1' for wrong * and '0' for right *" << std::endl;
88// return 0;
89// }
90
91// if(sscanf(argv[1], "%d", &WRONG_TH23) != 1 ||
92// (WRONG_TH23 != 1 && WRONG_TH23 != 0) ||
93// sscanf(argv[2], "%d", &WRONG_HIER) != 1 ||
94// (WRONG_HIER != 1 && WRONG_HIER != 0)){
95// std::cerr<< "cannot read parameters" << std::endl;
96// return 0;
97// }
98
99 //Process data card file
100 std::string ROOTDIR = getenv("GLBFREJUSROOT");
101 std::string DATACARD;
102 DATACARD = ROOTDIR + "/run/measTh13Delta.data";
103 std::ifstream fileDataCard;
104 std::cout << "Open DataCard file : " << DATACARD <<std::endl;
105 fileDataCard.open(DATACARD.data());
106 std::string DUMMYSTR; //comment string of the variables
107
108 std::string ROOTFILENAME;
109
110
111 fileDataCard >> DUMMYSTR >> WRONG_TH23;
112 std::cout << DUMMYSTR << " " << WRONG_TH23 << std::endl;
113 fileDataCard >> DUMMYSTR >> WRONG_HIER;
114 std::cout << DUMMYSTR << " " << WRONG_HIER << std::endl;
115 fileDataCard >> DUMMYSTR >> TH12;
116 std::cout << DUMMYSTR << " " << TH12 << std::endl;
117 fileDataCard >> DUMMYSTR >> TH13;
118 std::cout << DUMMYSTR << " " << TH13 << std::endl;
119 fileDataCard >> DUMMYSTR >> TH23;
120 std::cout << DUMMYSTR << " " << TH23 << std::endl;
121 fileDataCard >> DUMMYSTR >> DELTA;
122 std::cout << DUMMYSTR << " " << DELTA << std::endl;
123 fileDataCard >> DUMMYSTR >> DMQ21;
124 std::cout << DUMMYSTR << " " << DMQ21 << std::endl;
125 fileDataCard >> DUMMYSTR >> DMQ31;
126 std::cout << DUMMYSTR << " " << DMQ31 << std::endl;
127 fileDataCard >> DUMMYSTR >> ROOTFILENAME;
128
129 //add to ROOT file name the conditions on the Octant/Hierarchy
130 ROOTFILENAME += "-h";
131 char hier[40];
132 sprintf(hier,"%d",WRONG_HIER);
133 ROOTFILENAME += hier;
134 char octo[40];
135 sprintf(octo,"%d",WRONG_TH23);
136 ROOTFILENAME += "-o";
137 ROOTFILENAME += octo;
138 ROOTFILENAME += ".root";
139 std::cout << DUMMYSTR << " " <<ROOTFILENAME << std::endl;
140
141 fileDataCard >> DUMMYSTR >> ND;
142 std::cout << DUMMYSTR << " " << ND << std::endl;
143 fileDataCard >> DUMMYSTR >> NTH13;
144 std::cout << DUMMYSTR << " " << NTH13 << std::endl;
145 fileDataCard >> DUMMYSTR >> DELTAMIN;
146 std::cout << DUMMYSTR << " " << DELTAMIN << std::endl;
147 fileDataCard >> DUMMYSTR >> DELTAMAX;
148 std::cout << DUMMYSTR << " " << DELTAMAX << std::endl;
149 fileDataCard >> DUMMYSTR >> S22TH13MIN;
150 std::cout << DUMMYSTR << " " << S22TH13MIN << std::endl;
151 fileDataCard >> DUMMYSTR >> S22TH13MAX;
152 std::cout << DUMMYSTR << " " << S22TH13MAX << std::endl;
153
154 //init AIDA for Histo/Tuple
155 AIDA::IAnalysisFactory* aida = AIDA_createAnalysisFactory();
156 if(!aida) {
157 std::cerr << " AIDA not found." << std::endl;
158 return 0;
159 }
160
161 //ROOT tree :
162 AIDA::ITreeFactory* treeFactory = aida->createTreeFactory();
163 std::string opts = "export=root";
164 AIDA::ITree* fTree =
165 treeFactory->create(ROOTFILENAME,"root",false,true,opts);
166 delete treeFactory;
167
168 //Booking Tuple
169 AIDA::ITupleFactory* tf = aida->createTupleFactory(*fTree);
170
171 int NumberOfColumn = 7;
172 std::vector<std::string> column(NumberOfColumn);
173 const char* c_column[] = {"theta12", "theta13", "theta23",
174 "deltaCP", "dm21", "dm31",
175 "chi2"};
176 std::vector<std::string> coltype(NumberOfColumn);
177 const char* c_coltype[] = {"double","double","double",
178 "double","double","double",
179 "double"};
180 for (int icol = 0; icol<NumberOfColumn; ++icol) {
181 column[icol] = c_column[icol];
182 coltype[icol] = c_coltype[icol];
183 }
184
185 AIDA::ITuple* myTuple =
186 tf->create("SplGlb","Measure SPL Hier=T, Oct=T",column,coltype);
187
188
189 //Init Globes
190 glbInit(argv[0]);
191 glbInitExperiment("../data/SPL.glb",
192 &glb_experiment_list[0],
193 &glb_num_of_exps);
194
195 // true values (reference for future Chi2 computations)
196 glb_params true_values = glbAllocParams();
197 glbDefineParams(true_values, TH12, TH13, TH23, DELTA, DMQ21, DMQ31);
198 std::cout << "The Initial True_values" << std::endl;
199 glbPrintParams(stdout,true_values);
200
201 glbSetOscillationParameters(true_values);
202 glbSetRates();
203
204 // set Starting Values
205 glb_params start_values = glbAllocParams();
206 glbDefineParams(start_values, TH12, TH13, TH23, DELTA, DMQ21, DMQ31);
207
208 if(WRONG_HIER == 1)
209 glbSetOscParams(start_values, -DMQ31+DMQ21, GLB_DM_ATM);
210
211 if(WRONG_TH23 == 1)
212 glbSetOscParams(start_values, 0.5*M_PI - TH23, GLB_THETA_23);
213
214 glbSetStartingValues(start_values);
215 std::cout << "The Starting values 1st trail" << std::endl;
216 glbPrintParams(stdout,start_values);
217 glbFreeParams(start_values);
218
219
220 // init test values
221 glb_params test_values = glbAllocParams();
222 glbDefineParams(test_values, TH12, TH13, TH23, DELTA, DMQ21, DMQ31);
223
224 if(WRONG_HIER == 1)
225 glbSetOscParams(test_values, -DMQ31+DMQ21, GLB_DM_ATM);
226
227 if(WRONG_TH23 == 1)
228 glbSetOscParams(test_values, 0.5*M_PI - TH23, GLB_THETA_23);
229 std::cout << "The Test_values" << std::endl;
230 glbPrintParams(stdout,test_values);
231
232 glb_params fit_values = glbAllocParams(); //final minimization
233
234 // Set theta23 error very low to constraint the minimizer in the right octant
235 // input errors
236 glb_params input_errors = glbAllocParams();
237 glbDefineParams(input_errors,
238 TH12*0.05, 0., TH23*0.30, 0.,
239 DMQ21*0.05, fabs(DMQ31)*0.30);
240 glbSetDensityParams(input_errors, 0.05, GLB_ALL);
241 glbSetInputErrors(input_errors);
242
243
244 std::cout << "Start to loop...." << std::endl;
245
246 double res;
247 double fit_theta12, fit_theta13,fit_theta23,fit_deltaCP,fit_dmSol,fit_dmAtm;
248 double delta, sq2th, th;
249
250 for(int j = 0; j < ND; j++) {
251 delta = DELTAMIN + (DELTAMAX-DELTAMIN) * double(j) / (ND-1);
252 delta *= M_PI; //rescale in radian
253
254 glbSetOscParams(test_values, delta, GLB_DELTA_CP);
255
256 for(int k = 0; k < NTH13; k++) {
257 sq2th = S22TH13MIN + (S22TH13MAX - S22TH13MIN) * double(k) / (NTH13-1);
258 th = 0.5 * asin(sqrt(sq2th));
259
260 std::cout << "delta j:" << j << " theta13 k:" << k << std::endl;
261
262 glbSetOscParams(test_values, th, GLB_THETA_13);
263
264
265 //Find minimum in the projection Theta13 Delta_CP
266 res = glbChiThetaDelta(test_values,fit_values,0);
267
268 //Check if the Parameter Vector respect degeneracy Octant/Hierarchy
269 fit_theta12 = glbGetOscParams(fit_values, GLB_THETA_12);
270 fit_theta13 = glbGetOscParams(fit_values, GLB_THETA_13);
271 fit_theta23 = glbGetOscParams(fit_values, GLB_THETA_23);
272 fit_deltaCP = glbGetOscParams(fit_values, GLB_DELTA_CP);
273 fit_dmSol = glbGetOscParams(fit_values, GLB_DM_SOL);
274 fit_dmAtm = glbGetOscParams(fit_values, GLB_DM_ATM);
275
276 if (!test(fit_values)) res += 1e6;
277
278
279
280// std::cout << "Theta_12 " << fit_theta12 << " "
281// << "Theta_13 " << fit_theta13 << " "
282// << "Theta_23 " << fit_theta23*180./M_PI << " "
283// << "Delta_CP " << fit_deltaCP << " "
284// << "DM_21 " << fit_dmSol << " "
285// << "DM_31 " << fit_dmAtm << " "
286// << "Chi2 (mini) " << res
287// << std::endl;
288
289
290 myTuple->fill(0,fit_theta12);
291 myTuple->fill(1,fit_theta13);
292 myTuple->fill(2,fit_theta23);
293 myTuple->fill(3,fit_deltaCP);
294 myTuple->fill(4,fit_dmSol);
295 myTuple->fill(5,fit_dmAtm);
296 myTuple->fill(6,res);
297 myTuple->addRow();
298
299
300
301 }//Loop on Theta13
302 }//Loop on DeltaCP
303
304 glbFreeParams(true_values);
305 glbFreeParams(test_values);
306 glbFreeParams(fit_values);
307
308
309 //save Tuple
310 fTree->commit();
311 delete fTree;
312
313 delete aida;
314 return 0;
315}
316
317
318
319
320
Note: See TracBrowser for help on using the repository browser.