35 {
36 double qqmin,qqmax;
37 double cos1min,cos1max,cos2min,cos2max,cos3min,cos3max;
38 double dsigm1,dsigm2,sigma1,
sigma2,sigma,dsigm,Ar[14],Ar_r[14];
39 int nm,i,s_seed[105];
40 long int nges,k,j;
41 char outfile[20];
42
43
45 if( seeds.is_open() )
46 {
47 int ii=0;
48 while( ! seeds.eof() )
49 seeds >> s_seed[ii++];
50 }
51 else
52 cerr << "Cannot open file seed.dat for reading" << endl;
53
55
56
57
58
59 INPUT(nges, nm, outfile);
60
61
62
63
64
65 cout <<
"----------------------------------------------------" <<
FLAGS.pion<< endl;
67 cout << " PHOKHARA 6.0 : e^+ e^- -> mu^+ mu^- gamma" << endl;
68 else if (
FLAGS.pion == 1)
69 cout << " PHOKHARA 6.0: e^+ e^- -> pi^+ pi^- gamma" << endl;
70 else if (
FLAGS.pion == 2)
71 cout << " PHOKHARA 6.0: e^+ e^- -> pi^+ pi^- 2pi^0 gamma" << endl;
72 else if (
FLAGS.pion == 3)
73 cout << " PHOKHARA 6.0: e^+ e^- -> 2pi^+ 2pi^- gamma" << endl;
74 else if (
FLAGS.pion == 4)
75 cout << " PHOKHARA 6.0: e^+ e^- -> p pbar gamma" << endl;
76 else if (
FLAGS.pion == 5)
77 cout << " PHOKHARA 6.0: e^+ e^- -> n nbar gamma" << endl;
78 else if (
FLAGS.pion == 6)
79 cout << " PHOKHARA 6.0: e^+ e^- -> K^+ K^- gamma" << endl;
80 else if (
FLAGS.pion == 7)
81 cout << " PHOKHARA 6.0: e^+ e^- -> K_0 K_0bar gamma" << endl;
82 else if (
FLAGS.pion == 8)
83 cout << " PHOKHARA 6.0: e^+ e^- -> pi^+ pi^- pi^0 gamma" << endl;
84 else if (
FLAGS.pion == 9) {
85 cout << "PHOKHARA 6.0 : e^+ e^- ->" << endl;
86 cout << " Lambda (-> pi^- p) Lambda bar (-> pi^+ pbar) gamma" << endl;
87 } else
88 cout << " PHOKHARA 6.0: not yet implemented" << endl;
89
90
91 cout << "----------------------------------------------------" << endl;
92 printf(
" %s %f %s\n",
"cms total energy = ",sqrt(
CTES.Sp),
" GeV ");
93 if (
FLAGS.tagged == 0) {
94 if((
CUTS.gmin/2.0/
CTES.ebeam) < 0.0098){
95 cerr << " minimal missing energy set to small" << endl;
96 return 0;
97 }
98 printf(
" %s %f %s\n",
"minimal tagged photon energy = ",
CUTS.gmin,
" GeV ");
99 printf(
" %s %f,%f\n",
"angular cuts on tagged photon = ",
CUTS.phot1cut,
CUTS.phot2cut);
100 }
101 else {
102 if((
CUTS.gmin/2.0/
CTES.ebeam) < 0.0098){
103 cerr << " minimal missing energy set to small" << endl;
104 return 0;
105 }
106 printf(
" %s %f %s\n",
"minimal missing energy = ",
CUTS.gmin,
" GeV ");
107 printf(
" %s %f,%f\n",
"angular cuts on missing momentum = ",
CUTS.phot1cut,
CUTS.phot2cut);
108 }
109
110
111
113 printf(
" %s %f,%f\n",
"angular cuts on muons = ",
CUTS.pi1cut,
CUTS.pi2cut);
114 else if (
FLAGS.pion == 4)
115 printf(
" %s %f,%f\n",
"angular cuts on protons = ",
CUTS.pi1cut,
CUTS.pi2cut);
116 else if (
FLAGS.pion == 5)
117 printf(
" %s %f,%f\n",
"angular cuts on neutrons = ",
CUTS.pi1cut,
CUTS.pi2cut);
118 else if ((
FLAGS.pion == 6)||(
FLAGS.pion == 7))
119 printf(
" %s %f,%f\n",
"angular cuts on kaons = ",
CUTS.pi1cut,
CUTS.pi2cut);
120 else if (
FLAGS.pion == 9)
121 printf(
" %s %f,%f\n",
"angular cuts on pions and protons = ",
CUTS.pi1cut,
CUTS.pi2cut);
122 else
123 printf(
" %s %f,%f\n",
"angular cuts on pions = ",
CUTS.pi1cut,
CUTS.pi2cut);
124
125 if (
FLAGS.tagged == 0) {
127 printf(
" %s %f %s\n",
"min. muons-tagged photon inv.mass^2 = ",
CUTS.q2min,
" GeV^2");
128 else if (
FLAGS.pion == 4)
129 printf(
" %s %f %s\n",
"min. protons-tagged photon inv.mass^2 = ",
CUTS.q2min,
" GeV^2" );
130 else if (
FLAGS.pion == 5)
131 printf(
" %s %f %s\n",
"min. neutrons-tagged photon inv.mass^2 = ",
CUTS.q2min,
" GeV^2" );
132 else if ((
FLAGS.pion == 6)||(
FLAGS.pion == 7))
133 printf(
" %s %f %s\n",
"min. kaons-tagged photon inv.mass^2 = ",
CUTS.q2min,
" GeV^2" );
134 else if (
FLAGS.pion == 9)
135 printf(
" %s %f %s\n",
" min. lambdas-tagged photon inv.mass^2 = ",
CUTS.q2min,
" GeV^2" );
136 else
137 printf(
" %s %f %s\n",
"min. pions-tagged photon inv.mass^2 = ",
CUTS.q2min,
" GeV^2" );
138 }
139
140
141
143
144
145 if (
FLAGS.tagged == 0) {
148 } else {
149 cos1min = -1.0;
150 cos1max = 1.0;
152 }
153 cos2min = -1.0;
154 cos2max = 1.0;
155 cos3min = -1.0;
156 cos3max = 1.0;
159 else if (
FLAGS.pion == 1)
161 else if (
FLAGS.pion == 2)
163 else if (
FLAGS.pion == 3)
165 else if (
FLAGS.pion == 4)
167 else if (
FLAGS.pion == 5)
169 else if (
FLAGS.pion == 6)
171 else if (
FLAGS.pion == 7)
173 else if (
FLAGS.pion == 8)
175 else if (
FLAGS.pion == 9)
177
178
179
181 if (
CUTS.q2_max_c < qqmax)
183
184
186 qqmin =
CUTS.q2_min_c;
187 else {
188 cerr << "------------------------------" << endl;
189 cerr << " Q^2_min TOO SMALL" << endl;
190 cerr << " Q^2_min CHANGE BY PHOKHARA = " << qqmin << " GeV^2" << endl;
191 cerr << "------------------------------" << endl;
192 }
193
194 if(qqmax <= qqmin){
195 cerr << " Q^2_max to small " << endl;
196 cerr << " Q^2_max = " << qqmax << endl;
197 cerr << " Q^2_min = " << qqmin << endl;
198 return 0;
199 }
200
201
202 if (
FLAGS.pion == 0) {
203 printf(" %s %f %s\n", "minimal muon-pair invariant mass^2 = ", qqmin," GeV^2");
204 printf(" %s %f %s\n", "maximal muon-pair invariant mass^2 = ", qqmax," GeV^2");
205 }
else if (
FLAGS.pion == 1) {
206 printf(" %s %f %s\n", "minimal pion-pair invariant mass^2 = ", qqmin," GeV^2");
207 printf(" %s %f %s\n", "maximal pion-pair invariant mass^2 = ", qqmax," GeV^2");
208 }
else if (
FLAGS.pion == 4) {
209 printf(" %s %f %s\n", "minimal proton-pair invariant mass^2 = ", qqmin," GeV^2");
210 printf(" %s %f %s\n", "maximal proton-pair invariant mass^2 = ", qqmax," GeV^2");
211 }
else if (
FLAGS.pion == 5) {
212 printf(" %s %f %s\n", "minimal neutron-pair invariant mass^2 = ", qqmin," GeV^2");
213 printf(" %s %f %s\n", "maximal neutron-pair invariant mass^2 = ", qqmax," GeV^2");
214 }
else if ((
FLAGS.pion == 6) || (
FLAGS.pion == 7)) {
215 printf(" %s %f %s\n", "minimal kaon-pair invariant mass^2 = ", qqmin," GeV^2");
216 printf(" %s %f %s\n", "maximal kaon-pair invariant mass^2 = ", qqmax," GeV^2");
217 }
else if(
FLAGS.pion == 8){
218 printf(" %s %f %s\n", "minimal three-pion invariant mass^2 = ", qqmin," GeV^2");
219 printf(" %s %f %s\n", "maximal three-pion invariant mass^2 = ", qqmax," GeV^2");
220 }
else if(
FLAGS.pion == 9){
221 printf(" %s %f %s\n", "minimal lambda-pair invariant mass^2 = ", qqmin," GeV^2");
222 printf(" %s %f %s\n", "maximal lambda-pair invariant mass^2 = ", qqmax," GeV^2");
223 } else {
224 printf(" %s %f %s\n", "minimal four-pion invariant mass^2 = ", qqmin," GeV^2" );
225 printf(" %s %f %s\n", "maximal four-pion invariant mass^2 = ", qqmax," GeV^2");
226 }
227
228 if (
FLAGS.nlo == 0) {
229 cout << "Born" << endl;
230 if(
FLAGS.fsrnlo != 0){
231 cerr << "WRONG FSRNLO flag - only fsrnlo=0 allowed for Born" << endl;
232 return 0;
233 }
234 }
235
237 cerr << "WRONG NLO flag - only Born allowed for Lambdas"<< endl;
238 cerr << "If you feel that you need NLO"<< endl;
239 cerr << "please contact the authors"<< endl;
240 return 0;
241 }
242
243 if (
FLAGS.nlo == 1) printf(
" %s %f\n",
"NLO: soft photon cutoff w = ",
CUTS.w);
245 {
246
249 || (
FLAGS.fsr == 0) || (
FLAGS.fsrnlo == 0))) {
250 cerr << "WRONG combination of FSR, FSRNLO flags" <<endl;
251 return 0;
252 }
253
254
256 cout << "ISR only" << endl;
257 else if (
FLAGS.fsr == 1)
258 cout << "ISR+FSR" << endl;
259 else if (
FLAGS.fsr == 2) {
261 cout << "ISR+INT+FSR" << endl;
262 else {
263 cerr << "WRONG FSR flag: interference is included only for nlo=0" << endl;
264 return 0;
265 }
266 }
267 else {
268 cerr <<
"WRONG FSR flag" <<
FLAGS.fsr << endl;
269 return 0;
270 }
271
272 if(
FLAGS.fsrnlo == 1)
273 cout << "IFSNLO included" << endl;
274 }
275 else
276 {
278 cout << "ISR only" << endl;
279 else {
280 cerr << "FSR is implemented only for pi+pi-, mu+mu- and K+K- modes" << endl;
281 return 0;
282 }
283 }
284
285
287 cout << "Vacuum polarization is NOT included" << endl;
288 }
else if(
FLAGS.ivac == 1){
289 cout << "Vacuum polarization is included" << endl;
290 } else {
291 cout << "WRONG vacuum polarization switch" << endl;
292 return 0;
293 }
294
296 if(
FLAGS.FF_pion == 0)
297 cout << "Kuhn-Santamaria PionFormFactor" << endl;
298 else if(
FLAGS.FF_pion == 1)
299 cout << "Gounaris-Sakurai PionFormFactor" << endl;
300 else {
301 cout << "WRONG PionFormFactor switch" << endl;
302 return 0;
303 }
304
306 if(
FLAGS.f0_model == 0)
307 cout << "f0+f0(600): K+K- model" << endl;
308 else if(
FLAGS.f0_model == 1)
309 cout << "f0+f0(600): \"no structure\" model" << endl;
310 else if(
FLAGS.f0_model == 2)
311 cout << "NO f0+f0(600)" << endl;
312 else if(
FLAGS.f0_model == 3)
313 cout << "only f0, KLOE: Cesare Bini-private communication" << endl;
314 else {
315 cout << "WRONG f0+f0(600) switch" << endl;
316 return 0;
317 }
318 }
319 }
320
321
322
323 k = nm;
324 for( i = 0; i<2; i++)
325 {
329 }
332
333 for( i = 1; i<=2; i++) {
338
339
340
341 for(j = 1; j <= k; j++)
342 {
344 Ar[1] = Ar_r[0];
345
348 GEN_1PH(i,qqmin,qqmax,cos1min,cos1max,cos3min,cos3max);
349 }
350 else {
352 GEN_2PH(i,qqmin,cos1min,cos1max,cos2min,cos2max,cos3min,cos3max);
353 }
354 }
355
356
357
358 k = nges;
359 if (i == 1) {
362
367
370
374 }
375
379 }
380 }
381 }
382
383
384
388
389
390
392
393
394 if (
FLAGS.nlo == 0) {
397 } else {
402
404 dsigm = sqrt(dsigm1*dsigm1 + dsigm2*dsigm2);
405 }
406
407
408 cout << "----------------------------------------------------" << endl;
409 cout << int(
MAXIMA.tr[0]+
MAXIMA.tr[1]) <<
" total events accepted of " << endl;
410 cout << int(nges) << " total events generated" << endl;
411 cout << int(
MAXIMA.tr[0]) <<
" one photon events accepted of " << endl;
412 cout << int(
MAXIMA.count[0]) <<
" events generated" << endl;
413 cout << int(
MAXIMA.tr[1]) <<
" two photon events accepted of " << endl;
414 cout << int(
MAXIMA.count[1]) <<
" events generated" << endl;
415 cout << endl;
416 if (
FLAGS.nlo != 0) cout <<
"sigma1(nbarn) = " << sigma1 <<
" +- " << dsigm1 << endl;
417 if (
FLAGS.nlo != 0) cout <<
"sigma2(nbarn) = " <<
sigma2 <<
" +- " << dsigm2 << endl;
418 cout << "sigma (nbarn) = " << sigma << " +- " <<dsigm << endl;
419 cout << endl;
420 cout <<
"maximum1 = " <<
MAXIMA.gross[0] <<
" minimum1 = " <<
MAXIMA.klein[0] << endl;
421 cout <<
"Mmax1 = " <<
MAXIMA.Mmax[0] << endl;
422 cout <<
"maximum2 = " <<
MAXIMA.gross[1] <<
" minimum2 = " <<
MAXIMA.klein[1] << endl;
423 cout <<
"Mmax2 = " <<
MAXIMA.Mmax[1] << endl;
424 cout << "----------------------------------------------------" << endl;
425
426
427
428
429
430
431
432
433
434
435
436
437return 0;
438}
double cos(const BesAngle a)
#define GEN_1PH(I, QQMIN, QQMAX, COS1MIN, COS1MAX, COS3MIN, COS3MAX)
#define INPUT(NGES, NM, OUTFILE)
#define GEN_2PH(I, QQMIN, COS1MIN, COS1MAX, COS2MIN, COS2MAX, COS3MIN, COS3MAX)