source: GLBFrejus/HEAD/src/examine.cxx@ 151

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

add new appli

  • Property svn:executable set to *
File size: 15.1 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 <values.h>
8#include <fstream>
9#include <iostream>
10
11//Globes
12#include <globes/globes.h>
13
14// AIDA :
15#include <AIDA/IAnalysisFactory.h>
16#include <AIDA/ITreeFactory.h>
17#include <AIDA/IHistogramFactory.h>
18#include <AIDA/ITree.h>
19#include <AIDA/ITupleFactory.h>
20#include <AIDA/ITuple.h>
21
22
23
24#define TH12 ( asin(sqrt(0.3)) )
25#define TH23 ( M_PI/4.0 )
26//theta13 used by Thomas hep-ph/0501037
27//#define TH13 ( 0.5 * asin(sqrt(0.05)) )
28//theta13 = 1 degree, 3 degree
29#define TH13 ( 3.0/180.0 * M_PI)
30
31
32#define DELTA ( M_PI/2 )
33//#define DELTA ( 0.0 )
34
35#define DMQ21 7.9e-5
36#define DMQ31 2.4e-3
37
38
39
40/**************************************/
41/* main */
42/**************************************/
43
44int main(int argc, char *argv[])
45{
46// if(argc != 3){
47// std::cerr << "usage:" << argv[0] << "WRONG_TH23 WRONG_HIER [debug]" << std::endl;
48// std::cerr<< "WRONG_* = '1' for wrong * and '0' for right *" << std::endl;
49// return 0;
50// }
51
52// if(sscanf(argv[1], "%d", &WRONG_TH23) != 1 ||
53// (WRONG_TH23 != 1 && WRONG_TH23 != 0) ||
54// sscanf(argv[2], "%d", &WRONG_HIER) != 1 ||
55// (WRONG_HIER != 1 && WRONG_HIER != 0)){
56// std::cerr<< "cannot read parameters" << std::endl;
57// return 0;
58// }
59
60// int debug = 0;
61// if(4 == argc) {
62// if(sscanf(argv[3],"%d",&debug) !=1 ) {
63// std::cerr<< "cannot read debug" << std::endl;
64// return 0;
65// }
66// }
67
68 //init AIDA for Histo/Tuple
69 AIDA::IAnalysisFactory* aida = AIDA_createAnalysisFactory();
70 if(!aida) {
71 std::cerr << " AIDA not found." << std::endl;
72 return 0;
73 }
74
75 //ROOT tree :
76 AIDA::ITreeFactory* treeFactory = aida->createTreeFactory();
77 std::string opts = "export=root";
78 AIDA::ITree* fTree =
79 treeFactory->create("bidon.root","root",false,true,opts);
80 delete treeFactory;
81
82 //Booking Tuple
83 AIDA::ITupleFactory* tf = aida->createTupleFactory(*fTree);
84
85 int NumberOfColumn = 7;
86 std::vector<std::string> column(NumberOfColumn);
87 const char* c_column[] = {"theta12", "theta13", "theta23",
88 "deltaCP", "dm21", "dm31",
89 "chi2"};
90 std::vector<std::string> coltype(NumberOfColumn);
91 const char* c_coltype[] = {"double","double","double",
92 "double","double","double",
93 "double"};
94 for (int icol = 0; icol<NumberOfColumn; ++icol) {
95 column[icol] = c_column[icol];
96 coltype[icol] = c_coltype[icol];
97 }
98
99 AIDA::ITuple* myTuple =
100 tf->create("SplGlb","dump",column,coltype);
101
102
103 //Init Globes: Attention SPL.glb may be a symbolic link
104 glbInit(argv[0]);
105 glbInitExperiment("../data/SPL.glb",
106 &glb_experiment_list[0],
107 &glb_num_of_exps);
108
109 // true values (reference for future Chi2 computations)
110 glb_params true_values = glbAllocParams();
111 glbDefineParams(true_values, TH12, TH13, TH23, DELTA, DMQ21, DMQ31);
112 std::cout << "The Initial True_values" << std::endl;
113 glbPrintParams(stdout,true_values);
114
115 glbSetOscillationParameters(true_values);
116 glbSetRates();
117
118 // set Starting Values
119 glb_params start_values = glbAllocParams();
120 glbDefineParams(start_values, TH12, TH13, TH23, DELTA, DMQ21, DMQ31);
121 glbSetStartingValues(start_values);
122 std::cout << "The Starting values" << std::endl;
123 glbPrintParams(stdout,start_values);
124
125
126
127 //Fluxes
128 std::cout << "Flux Name(0):<" << glbValueToName(0,"flux",0) << ">" << std::endl;
129 int splFluxPosFocus = glbNameToValue(0,"flux","#SPLplus");
130 std::cout << "Flux Id(#SPLplus):<"<<splFluxPosFocus << ">" << std::endl;
131 int splFluxNegFocus = glbNameToValue(0,"flux","#SPLminus");
132
133
134 std::cout << "Fluxes at 130km for 300MeV neutrinos" << std::endl;
135 double energy = 0.300;
136 double baseline = 1.0; //the Flux are already defined at 130km...
137 std::cout << "Flux(ve) (+foc): " << glbFlux(0,splFluxPosFocus,energy,baseline,1,1) << std::endl;
138 std::cout << "Flux(vm) (+foc): " << glbFlux(0,splFluxPosFocus,energy,baseline,2,1) << std::endl;
139 std::cout << "Flux(ave) (+foc): " << glbFlux(0,splFluxPosFocus,energy,baseline,1,-1) << std::endl;
140 std::cout << "Flux(avm) (+foc): " << glbFlux(0,splFluxPosFocus,energy,baseline,2,-1) << std::endl;
141
142 std::cout << "Flux(ve) (-foc): " << glbFlux(0,splFluxNegFocus,energy,baseline,1,1) << std::endl;
143 std::cout << "Flux(vm) (-foc): " << glbFlux(0,splFluxNegFocus,energy,baseline,2,1) << std::endl;
144 std::cout << "Flux(ave) (-foc): " << glbFlux(0,splFluxNegFocus,energy,baseline,1,-1) << std::endl;
145 std::cout << "Flux(avm) (-foc): " << glbFlux(0,splFluxNegFocus,energy,baseline,2,-1) << std::endl;
146
147 //Proba
148 double ene[5] = {0.18,0.38,0.58,0.78,0.98};
149 std::cout << "Energy, mu->e, mu->mu, mu->tau" << std::endl;
150 for (int iene=0; iene<5; ++iene) {
151 std::cout << ene[iene] << " "
152 << glbVacuumProbability(2,1,+1,ene[iene],130.0) << " "
153 << glbVacuumProbability(2,2,+1,ene[iene],130.0) << " "
154 << glbVacuumProbability(2,3,+1,ene[iene],130.0) << " "
155 << std::endl;
156 }
157
158
159 //biodn
160 int splChannelNumuNonOsci = glbNameToValue(0,"channel","#numu_non_osci");
161 std::cout << "Rate of Channel(#numuNonOsci) NO_OSC: GLB_PRE, W_EFF, W_BG " << std::endl;
162 glbShowChannelRates(stdout,0,splChannelNumuNonOsci,GLB_PRE,GLB_W_EFF,GLB_W_BG);
163 std::cout << std::endl;
164 std::cout << "Rate of Channel(#numuNonOsci) NO_OSC: GLB_POST, W_EFF, W_BG " << std::endl;
165 glbShowChannelRates(stdout,0,splChannelNumuNonOsci,GLB_POST,GLB_W_EFF,GLB_W_BG);
166 std::cout << std::endl;
167
168
169 //Rules
170
171
172 int splChannelNCplus = glbNameToValue(0,"channel","#NC_bckg_plus");
173 std::cout << "Rate of Channel(#NC_bckg_plus) NO_OSC: GLB_PRE, W_EFF, W_BG " << std::endl;
174 glbShowChannelRates(stdout,0,splChannelNCplus,GLB_PRE,GLB_W_EFF,GLB_W_BG);
175 std::cout << std::endl;
176 std::cout << "Rate of Channel(#NC_bckg_plus) NO_OSC: GLB_POST, W_EFF, W_BG " << std::endl;
177 glbShowChannelRates(stdout,0,splChannelNCplus,GLB_POST,GLB_W_EFF,GLB_W_BG);
178 std::cout << std::endl;
179
180
181
182
183 //NU_E_app
184 int splRuleNueApp = glbNameToValue(0,"rule","#NU_E_app");
185 int nItemRuleNUEappSignal = glbGetLengthOfRule(0,splRuleNueApp,GLB_SIG);
186 int nItemRuleNUEappBkg = glbGetLengthOfRule(0,splRuleNueApp,GLB_BG);
187
188 std::cout << "Length Of Rule(#NU_E_app) Signal: "<< nItemRuleNUEappSignal
189 << " Bkg: " << nItemRuleNUEappBkg
190 << std::endl;
191
192 std::cout << "Normalisation Rule(#NU_E_app) Signal: "
193 << glbGetNormalizationInRule(0,nItemRuleNUEappSignal,GLB_SIG)
194 << " Bkg: "
195 << glbGetNormalizationInRule(0,nItemRuleNUEappSignal,GLB_BG)
196 << std::endl;
197
198 std::cout << "Rate of Rules(#NU_E_app) Signal: W_EFF, W_BG, W_COEFF" << std::endl;
199 glbShowRuleRates(stdout,0,splRuleNueApp,GLB_ALL,GLB_W_EFF,GLB_W_BG,GLB_W_COEFF,GLB_SIG);
200 std::cout << std::endl;
201
202
203 std::cout << "Channels of Rule(#NU_E_app) Signal" << std::endl;
204 for (int item=0; item<nItemRuleNUEappSignal; ++item) {
205 std::cout << " Channel<"
206 << glbValueToName(0,"channel",
207 glbGetChannelInRule(0,splRuleNueApp,item,GLB_SIG))
208 << "> with coeff. "
209 << glbGetCoefficientInRule(0,splRuleNueApp,item,GLB_SIG)
210 << std::endl;
211 std::cout << " Rate: WO_EFF, WO_BG, WO_COEFF: "
212 << glbTotalRuleRate(0,splRuleNueApp,item,GLB_WO_EFF,GLB_WO_BG,GLB_WO_COEFF,GLB_SIG)
213 << std::endl;
214 }
215
216 std::cout << "Rate of Rules(#NU_E_app) BG: W_EFF, W_BG, W_COEFF"
217 << std::endl;
218 glbShowRuleRates(stdout,0,splRuleNueApp,
219 GLB_ALL,GLB_W_EFF,GLB_W_BG,GLB_W_COEFF,GLB_BG);
220
221 std::cout << "Channels of Rule(#NU_E_app) Bkg" << std::endl;
222 for (int item=0; item<nItemRuleNUEappBkg; ++item) {
223 std::cout << " Channel<"
224 << glbValueToName(0,"channel",
225 glbGetChannelInRule(0,splRuleNueApp,item,GLB_BG))
226 << "> with coeff. "
227 << glbGetCoefficientInRule(0,splRuleNueApp,item,GLB_BG)
228 << std::endl;
229 std::cout << " Rate: WO_EFF, WO_BG, WO_COEFF: "
230 << glbTotalRuleRate(0,splRuleNueApp,item,GLB_WO_EFF,GLB_WO_BG,GLB_WO_COEFF,GLB_BG)
231 << std::endl;
232 }
233
234
235 //NUMU_DISA: FIXME: problem of the PRE-SMEARING efficiency effects (JEC 6/6/05)
236 int splRuleNumudisa = glbNameToValue(0,"rule","#NUMU_DISA");
237 int nItemRuleNUMUDisaSignal = glbGetLengthOfRule(0,splRuleNumudisa,GLB_SIG);
238 int nItemRuleNUMUDisaBkg = glbGetLengthOfRule(0,splRuleNumudisa,GLB_BG);
239
240 std::cout << "Length Of Rule(#NUMU_disa) Signal: "<< nItemRuleNUMUDisaSignal
241 << " Bkg: " << nItemRuleNUMUDisaBkg
242 << std::endl;
243
244 std::cout << "Normalisation Rule(#NUMU_disa) Signal: "
245 << glbGetNormalizationInRule(0,nItemRuleNUMUDisaSignal,GLB_SIG)
246 << " Bkg: "
247 << glbGetNormalizationInRule(0,nItemRuleNUMUDisaSignal,GLB_BG)
248 << std::endl;
249
250 std::cout << "Rate of Rules(#NUMU_disa) Signal: W_EFF, W_BG, W_COEFF" << std::endl;
251 glbShowRuleRates(stdout,0,splRuleNumudisa,GLB_ALL,GLB_W_EFF,GLB_W_BG,GLB_W_COEFF,GLB_SIG);
252 std::cout << std::endl;
253
254
255 std::cout << "Channels of Rule(#NUMU_disa) Signal" << std::endl;
256 for (int item=0; item<nItemRuleNUMUDisaSignal; ++item) {
257 std::cout << " Channel<"
258 << glbValueToName(0,"channel",
259 glbGetChannelInRule(0,splRuleNumudisa,item,GLB_SIG))
260 << "> with coeff. "
261 << glbGetCoefficientInRule(0,splRuleNumudisa,item,GLB_SIG)
262 << std::endl;
263 std::cout << " Rate: WO_EFF, WO_BG, WO_COEFF: "
264 << glbTotalRuleRate(0,splRuleNumudisa,item,GLB_WO_EFF,GLB_WO_BG,GLB_WO_COEFF,GLB_SIG)
265 << std::endl;
266 }
267
268 std::cout << "Rate of Rules(#NUMU_disa) BG: W_EFF, W_BG, W_COEFF"
269 << std::endl;
270 glbShowRuleRates(stdout,0,splRuleNumudisa,
271 GLB_ALL,GLB_W_EFF,GLB_W_BG,GLB_W_COEFF,GLB_BG);
272
273 std::cout << "Channels of Rule(#NUMU_disa) Bkg" << std::endl;
274 for (int item=0; item<nItemRuleNUMUDisaBkg; ++item) {
275 std::cout << " Channel<"
276 << glbValueToName(0,"channel",
277 glbGetChannelInRule(0,splRuleNumudisa,item,GLB_BG))
278 << "> with coeff. "
279 << glbGetCoefficientInRule(0,splRuleNumudisa,item,GLB_BG)
280 << std::endl;
281 std::cout << " Rate: WO_EFF, WO_BG, WO_COEFF: "
282 << glbTotalRuleRate(0,splRuleNumudisa,item,GLB_WO_EFF,GLB_WO_BG,GLB_WO_COEFF,GLB_BG)
283 << std::endl;
284 }
285
286
287
288
289 //NU_E_BAR_app
290 int splRuleAntiNueApp = glbNameToValue(0,"rule","#NU_E_BAR_app");
291 int nItemRuleAntiNUEappSignal = glbGetLengthOfRule(0,splRuleAntiNueApp,GLB_SIG);
292 int nItemRuleAntiNUEappBkg = glbGetLengthOfRule(0,splRuleAntiNueApp,GLB_BG);
293 std::cout << "Length Of Rule(#NU_E_BAR_app) Signal: "
294 << nItemRuleAntiNUEappSignal
295 << " Bkg: " << nItemRuleAntiNUEappBkg
296 << std::endl;
297
298 std::cout << "Normalisation Rule(#NU_E_BAR_app) Signal: "
299 << glbGetNormalizationInRule(0,nItemRuleAntiNUEappSignal,GLB_SIG)
300 << " Bkg: "
301 << glbGetNormalizationInRule(0,nItemRuleAntiNUEappSignal,GLB_BG)
302 << std::endl;
303
304 std::cout << "Rate of Rules(#NU_E_BARR_app) Signal: W_EFF, W_BG, W_COEFF" << std::endl;
305 glbShowRuleRates(stdout,0,splRuleAntiNueApp,GLB_ALL,GLB_W_EFF,GLB_W_BG,GLB_W_COEFF,GLB_SIG);
306 std::cout << std::endl;
307
308 std::cout << "Channels of Rule(#NU_E_BAR_app) Signal" << std::endl;
309 for (int item=0; item<nItemRuleAntiNUEappSignal; ++item) {
310 std::cout << " Channel<"
311 << glbValueToName(0,"channel",
312 glbGetChannelInRule(0,splRuleAntiNueApp,item,GLB_SIG))
313 << "> with coeff. "
314 << glbGetCoefficientInRule(0,splRuleAntiNueApp,item,GLB_SIG)
315 << std::endl;
316 std::cout << " Rate: WO_EFF, WO_BG, WO_COEFF: "
317 << glbTotalRuleRate(0,splRuleAntiNueApp,item,GLB_WO_EFF,GLB_WO_BG,GLB_WO_COEFF,GLB_SIG)
318 << std::endl;
319 }
320
321 std::cout << "Rate of Rules(#NU_E_BARR_app) BG: W_EFF, W_BG, W_COEFF"
322 << std::endl;
323 glbShowRuleRates(stdout,0,splRuleAntiNueApp,
324 GLB_ALL,GLB_W_EFF,GLB_W_BG,GLB_W_COEFF,GLB_BG);
325
326 std::cout << "Channels of Rule(#NU_E_BAR_app) Bkg" << std::endl;
327 for (int item=0; item<nItemRuleAntiNUEappBkg; ++item) {
328 std::cout << " Channel<"
329 << glbValueToName(0,"channel",
330 glbGetChannelInRule(0,splRuleAntiNueApp,item,GLB_BG))
331 << "> with coeff. "
332 << glbGetCoefficientInRule(0,splRuleAntiNueApp,item,GLB_BG)
333 << std::endl;
334 std::cout << " Rate: WO_EFF, WO_BG, WO_COEFF: "
335 << glbTotalRuleRate(0,splRuleAntiNueApp,item,GLB_WO_EFF,GLB_WO_BG,GLB_WO_COEFF,GLB_BG)
336 << std::endl;
337 }
338
339 //NUMUBAR_DISA: FIXME: problem of the PRE-SMEARING efficiency effects (JEC 6/6/05)
340 int splRuleNumuBardisa = glbNameToValue(0,"rule","#NUMUBAR_DISA");
341 int nItemRuleNUMUBARDisaSignal = glbGetLengthOfRule(0,splRuleNumuBardisa,GLB_SIG);
342 int nItemRuleNUMUBARDisaBkg = glbGetLengthOfRule(0,splRuleNumuBardisa,GLB_BG);
343
344 std::cout << "Length Of Rule(#NUMUBAR_disa) Signal: "<< nItemRuleNUMUBARDisaSignal
345 << " Bkg: " << nItemRuleNUMUBARDisaBkg
346 << std::endl;
347
348 std::cout << "Normalisation Rule(#NUMUBAR_disa) Signal: "
349 << glbGetNormalizationInRule(0,nItemRuleNUMUBARDisaSignal,GLB_SIG)
350 << " Bkg: "
351 << glbGetNormalizationInRule(0,nItemRuleNUMUBARDisaSignal,GLB_BG)
352 << std::endl;
353
354 std::cout << "Rate of Rules(#NUMUBAR_disa) Signal: W_EFF, W_BG, W_COEFF" << std::endl;
355 glbShowRuleRates(stdout,0,splRuleNumuBardisa,GLB_ALL,GLB_W_EFF,GLB_W_BG,GLB_W_COEFF,GLB_SIG);
356 std::cout << std::endl;
357
358
359 std::cout << "Channels of Rule(#NUMUBAR_disa) Signal" << std::endl;
360 for (int item=0; item<nItemRuleNUMUBARDisaSignal; ++item) {
361 std::cout << " Channel<"
362 << glbValueToName(0,"channel",
363 glbGetChannelInRule(0,splRuleNumuBardisa,item,GLB_SIG))
364 << "> with coeff. "
365 << glbGetCoefficientInRule(0,splRuleNumuBardisa,item,GLB_SIG)
366 << std::endl;
367 std::cout << " Rate: WO_EFF, WO_BG, WO_COEFF: "
368 << glbTotalRuleRate(0,splRuleNumuBardisa,item,GLB_WO_EFF,GLB_WO_BG,GLB_WO_COEFF,GLB_SIG)
369 << std::endl;
370 }
371
372 std::cout << "Rate of Rules(#NUMUBAR_disa) BG: W_EFF, W_BG, W_COEFF"
373 << std::endl;
374 glbShowRuleRates(stdout,0,splRuleNumuBardisa,
375 GLB_ALL,GLB_W_EFF,GLB_W_BG,GLB_W_COEFF,GLB_BG);
376
377 std::cout << "Channels of Rule(#NUMUBAR_disa) Bkg" << std::endl;
378 for (int item=0; item<nItemRuleNUMUBARDisaBkg; ++item) {
379 std::cout << " Channel<"
380 << glbValueToName(0,"channel",
381 glbGetChannelInRule(0,splRuleNumuBardisa,item,GLB_BG))
382 << "> with coeff. "
383 << glbGetCoefficientInRule(0,splRuleNumuBardisa,item,GLB_BG)
384 << std::endl;
385 std::cout << " Rate: WO_EFF, WO_BG, WO_COEFF: "
386 << glbTotalRuleRate(0,splRuleNumuBardisa,item,GLB_WO_EFF,GLB_WO_BG,GLB_WO_COEFF,GLB_BG)
387 << std::endl;
388 }
389
390 /*
391 fit_theta12 = glbGetOscParams(fit_values, GLB_THETA_12);
392 fit_theta13 = glbGetOscParams(fit_values, GLB_THETA_13);
393 fit_theta23 = glbGetOscParams(fit_values, GLB_THETA_23);
394 fit_deltaCP = glbGetOscParams(fit_values, GLB_DELTA_CP);
395 fit_dmSol = glbGetOscParams(fit_values, GLB_DM_SOL);
396 fit_dmAtm = glbGetOscParams(fit_values, GLB_DM_ATM);
397 */
398
399
400 //store result in Tuple
401 // myTuple->fill(0,fit_theta12);
402// myTuple->fill(1,fit_theta13);
403// myTuple->fill(2,fit_theta23);
404// myTuple->fill(3,fit_deltaCP);
405// myTuple->fill(4,fit_dmSol);
406// myTuple->fill(5,fit_dmAtm);
407// myTuple->fill(6,res);
408// myTuple->addRow();
409
410
411
412
413 glbFreeParams(true_values);
414 glbFreeParams(start_values);
415
416 //save Tuple
417 fTree->commit();
418 delete fTree;
419
420 delete aida;
421 return 0;
422}
423
424
425
426
427
Note: See TracBrowser for help on using the repository browser.