source: GLBFrejus/HEAD/src/proba.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: 5.9 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 TH13 (0.5*asin(sqrt(0.1)))
26#define TH23 ( M_PI/4.0 )
27
28#define DELTA ( M_PI/2 )
29
30#define DMQ21 7.9e-5
31#define DMQ31 2.4e-3
32
33#define ENERGY 0.300 //GeV
34#define NSTEP 1000
35#define MAXBASELINE 1000.0 //km
36
37
38double max(double x) {
39// #define MINVALUE 1e-10
40// return (x > MINVALUE ? x : MINVALUE);
41 return x;
42}
43
44
45/**************************************/
46/* main */
47/**************************************/
48
49int main(int argc, char *argv[])
50{
51 //init AIDA for Histo/Tuple
52 AIDA::IAnalysisFactory* aida = AIDA_createAnalysisFactory();
53 if(!aida) {
54 std::cerr << " AIDA not found." << std::endl;
55 return 0;
56 }
57
58 //ROOT tree :
59 AIDA::ITreeFactory* treeFactory = aida->createTreeFactory();
60 std::string opts = "export=root";
61 AIDA::ITree* fTree =
62 treeFactory->create("Proba.root","root",false,true,opts);
63 delete treeFactory;
64
65 //Booking Tuple
66 AIDA::ITupleFactory* tf = aida->createTupleFactory(*fTree);
67
68 int NumberOfColumn = 25;
69 std::vector<std::string> column(NumberOfColumn);
70 const char* c_column[] = { "Pee","Pemu","Petau",
71 "Pmue","Pmumu","Pmutau",
72 "Pbee","Pbemu","Pbetau",
73 "Pbmue","Pbmumu","Pbmutau",
74 "PMee","PMemu","PMetau",
75 "PMmue","PMmumu","PMmutau",
76 "PMbee","PMbemu","PMbetau",
77 "PMbmue","PMbmumu","PMbmutau",
78 "LovE" };
79 std::vector<std::string> coltype(NumberOfColumn);
80 const char* c_coltype[] = {"double","double","double",
81 "double","double","double",
82 "double","double","double",
83 "double","double","double",
84 "double","double","double",
85 "double","double","double",
86 "double","double","double",
87 "double","double","double",
88 "double" };
89 for (int icol = 0; icol<NumberOfColumn; ++icol) {
90 column[icol] = c_column[icol];
91 coltype[icol] = c_coltype[icol];
92 }
93
94 AIDA::ITuple* myTuple = tf->create("ProbaOscill","dump",column,coltype);
95
96
97 //Init Globes: Attention SPL.glb may be a symbolic link
98 glbInit(argv[0]);
99 glbSetVerbosityLevel(2);
100 glbInitExperiment("../data/Bidon.glb", &glb_experiment_list[0], &glb_num_of_exps);
101 // true values (reference for future Chi2 computations)
102 glb_params true_values = glbAllocParams();
103 glbDefineParams(true_values, TH12, TH13, TH23, DELTA, DMQ21, -DMQ31);
104 glbSetOscillationParameters(true_values);
105 glbSetRates();
106
107 std::cout << "Profile of exp 0: " << glbGetProfileTypeInExperiment(0) << std::endl;
108
109
110 double baseLine;
111 //Vacuum
112 //Proba with neutrinos
113 double Pee;
114 double Pemu;
115 double Petau;
116 double Pmue;
117 double Pmumu;
118 double Pmutau;
119
120 //Proba with anti-neutrinos
121 double Pbee;
122 double Pbemu;
123 double Pbetau;
124 double Pbmue;
125 double Pbmumu;
126 double Pbmutau;
127
128 //Matter
129 //Proba with neutrinos
130 double PMee;
131 double PMemu;
132 double PMetau;
133 double PMmue;
134 double PMmumu;
135 double PMmutau;
136
137 //PMroba with anti-neutrinos
138 double PMbee;
139 double PMbemu;
140 double PMbetau;
141 double PMbmue;
142 double PMbmumu;
143 double PMbmutau;
144
145 //Proba
146 for (int iLen=1; iLen<NSTEP; ++iLen) {
147 baseLine = iLen/double(NSTEP) * MAXBASELINE;
148
149 //Vaccum
150 Pee = glbVacuumProbability(1,1,+1,ENERGY,baseLine);
151 Pemu = glbVacuumProbability(1,2,+1,ENERGY,baseLine);
152 Petau = glbVacuumProbability(1,3,+1,ENERGY,baseLine);
153 Pmue = glbVacuumProbability(2,1,+1,ENERGY,baseLine);
154 Pmumu = glbVacuumProbability(2,2,+1,ENERGY,baseLine);
155 Pmutau = glbVacuumProbability(2,3,+1,ENERGY,baseLine);
156
157 Pbee = glbVacuumProbability(1,1,-1,ENERGY,baseLine);
158 Pbemu = glbVacuumProbability(1,2,-1,ENERGY,baseLine);
159 Pbetau = glbVacuumProbability(1,3,-1,ENERGY,baseLine);
160 Pbmue = glbVacuumProbability(2,1,-1,ENERGY,baseLine);
161 Pbmumu = glbVacuumProbability(2,2,-1,ENERGY,baseLine);
162 Pbmutau = glbVacuumProbability(2,3,-1,ENERGY,baseLine);
163
164 //Matter
165 glbSetBaselineInExperiment(0,baseLine);
166
167 PMee = glbProfileProbability(0,1,1,+1,ENERGY);
168 PMemu = glbProfileProbability(0,1,2,+1,ENERGY);
169 PMetau = glbProfileProbability(0,1,3,+1,ENERGY);
170 PMmue = glbProfileProbability(0,2,1,+1,ENERGY);
171 PMmumu = glbProfileProbability(0,2,2,+1,ENERGY);
172 PMmutau = glbProfileProbability(0,2,3,+1,ENERGY);
173
174 PMbee = glbProfileProbability(0,1,1,-1,ENERGY);
175 PMbemu = glbProfileProbability(0,1,2,-1,ENERGY);
176 PMbetau = glbProfileProbability(0,1,3,-1,ENERGY);
177 PMbmue = glbProfileProbability(0,2,1,-1,ENERGY);
178 PMbmumu = glbProfileProbability(0,2,2,-1,ENERGY);
179 PMbmutau = glbProfileProbability(0,2,3,-1,ENERGY);
180
181 //save into Tuple
182 myTuple->fill(0,max(Pee));
183 myTuple->fill(1,max(Pemu));
184 myTuple->fill(2,max(Petau));
185 myTuple->fill(3,max(Pmue));
186 myTuple->fill(4,max(Pmumu));
187 myTuple->fill(5,max(Pmutau));
188
189 myTuple->fill(6,max(Pbee));
190 myTuple->fill(7,max(Pbemu));
191 myTuple->fill(8,max(Pbetau));
192 myTuple->fill(9,max(Pbmue));
193 myTuple->fill(10,max(Pbmumu));
194 myTuple->fill(11,max(Pbmutau));
195
196 myTuple->fill(12,max(PMee));
197 myTuple->fill(13,max(PMemu));
198 myTuple->fill(14,max(PMetau));
199 myTuple->fill(15,max(PMmue));
200 myTuple->fill(16,max(PMmumu));
201 myTuple->fill(17,max(PMmutau));
202
203 myTuple->fill(18,max(PMbee));
204 myTuple->fill(19,max(PMbemu));
205 myTuple->fill(20,max(PMbetau));
206 myTuple->fill(21,max(PMbmue));
207 myTuple->fill(22,max(PMbmumu));
208 myTuple->fill(23,max(PMbmutau));
209
210 myTuple->fill(24,baseLine/ENERGY);
211
212 myTuple->addRow();
213 }
214 //Clear GLobES
215 glbFreeParams(true_values);
216
217 //save Tuple
218 fTree->commit();
219
220 //Clear OS
221 delete fTree;
222 delete aida;
223 return 0;
224}
225
226
227
228
229
Note: See TracBrowser for help on using the repository browser.