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

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

add new appli
  • Property svn:executable set to *
File size: 15.1 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 <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.