42 {
43 double qqmin, qqmax;
44 double cos1min, cos1max, cos2min, cos2max, cos3min, cos3max;
45 double dsigm1, dsigm2, sigma1, sigma2, sigma, dsigm, Ar[14], Ar_r[14];
46 int nm, i, s_seed[105];
47 long int nges, k, j;
48 char outfile[20];
49
50
51 ifstream seeds( "seed.dat" );
52 if ( seeds.is_open() )
53 {
54 int ii = 0;
55 while ( !seeds.eof() ) seeds >> s_seed[ii++];
56 }
57 else cerr << "Cannot open file seed.dat for reading" << endl;
58
60
61
62
63
64 INPUT( nges, nm, outfile );
65
66
67
68
69
70 cout <<
"----------------------------------------------------" <<
FLAGS.pion << endl;
71 if (
FLAGS.pion == 0 ) cout <<
" PHOKHARA 6.0 : e^+ e^- -> mu^+ mu^- gamma" << endl;
72 else if (
FLAGS.pion == 1 ) cout <<
" PHOKHARA 6.0: e^+ e^- -> pi^+ pi^- gamma" << endl;
73 else if (
FLAGS.pion == 2 )
74 cout << " PHOKHARA 6.0: e^+ e^- -> pi^+ pi^- 2pi^0 gamma" << endl;
75 else if (
FLAGS.pion == 3 ) cout <<
" PHOKHARA 6.0: e^+ e^- -> 2pi^+ 2pi^- gamma" << endl;
76 else if (
FLAGS.pion == 4 ) cout <<
" PHOKHARA 6.0: e^+ e^- -> p pbar gamma" << endl;
77 else if (
FLAGS.pion == 5 ) cout <<
" PHOKHARA 6.0: e^+ e^- -> n nbar gamma" << endl;
78 else if (
FLAGS.pion == 6 ) cout <<
" PHOKHARA 6.0: e^+ e^- -> K^+ K^- gamma" << endl;
79 else if (
FLAGS.pion == 7 ) cout <<
" PHOKHARA 6.0: e^+ e^- -> K_0 K_0bar gamma" << endl;
80 else if (
FLAGS.pion == 8 )
81 cout << " PHOKHARA 6.0: e^+ e^- -> pi^+ pi^- pi^0 gamma" << endl;
82 else if (
FLAGS.pion == 9 )
83 {
84 cout << "PHOKHARA 6.0 : e^+ e^- ->" << endl;
85 cout << " Lambda (-> pi^- p) Lambda bar (-> pi^+ pbar) gamma" << endl;
86 }
87 else cout << " PHOKHARA 6.0: not yet implemented" << endl;
88
89
90 cout << "----------------------------------------------------" << endl;
91 printf(
" %s %f %s\n",
"cms total energy = ", sqrt(
CTES.Sp ),
92 " GeV " );
93 if (
FLAGS.tagged == 0 )
94 {
95 if ( (
CUTS.gmin / 2.0 /
CTES.ebeam ) < 0.0098 )
96 {
97 cerr << " minimal missing energy set to small" << endl;
98 return 0;
99 }
100 printf(
" %s %f %s\n",
"minimal tagged photon energy = ",
CUTS.gmin,
" GeV " );
101 printf(
" %s %f,%f\n",
"angular cuts on tagged photon = ",
CUTS.phot1cut,
103 }
104 else
105 {
106 if ( (
CUTS.gmin / 2.0 /
CTES.ebeam ) < 0.0098 )
107 {
108 cerr << " minimal missing energy set to small" << endl;
109 return 0;
110 }
111 printf(
" %s %f %s\n",
"minimal missing energy = ",
CUTS.gmin,
" GeV " );
112 printf(
" %s %f,%f\n",
"angular cuts on missing momentum = ",
CUTS.phot1cut,
114 }
115
116
117 if (
FLAGS.pion == 0 )
118 printf(
" %s %f,%f\n",
"angular cuts on muons = ",
CUTS.pi1cut,
120 else if (
FLAGS.pion == 4 )
121 printf(
" %s %f,%f\n",
"angular cuts on protons = ",
CUTS.pi1cut,
123 else if (
FLAGS.pion == 5 )
124 printf(
" %s %f,%f\n",
"angular cuts on neutrons = ",
CUTS.pi1cut,
126 else if ( (
FLAGS.pion == 6 ) || (
FLAGS.pion == 7 ) )
127 printf(
" %s %f,%f\n",
"angular cuts on kaons = ",
CUTS.pi1cut,
129 else if (
FLAGS.pion == 9 )
130 printf(
" %s %f,%f\n",
"angular cuts on pions and protons = ",
CUTS.pi1cut,
132 else
133 printf(
" %s %f,%f\n",
"angular cuts on pions = ",
CUTS.pi1cut,
135
136 if (
FLAGS.tagged == 0 )
137 {
138 if (
FLAGS.pion == 0 )
139 printf(
" %s %f %s\n",
"min. muons-tagged photon inv.mass^2 = ",
CUTS.q2min,
140 " GeV^2" );
141 else if (
FLAGS.pion == 4 )
142 printf(
" %s %f %s\n",
"min. protons-tagged photon inv.mass^2 = ",
CUTS.q2min,
143 " GeV^2" );
144 else if (
FLAGS.pion == 5 )
145 printf(
" %s %f %s\n",
"min. neutrons-tagged photon inv.mass^2 = ",
CUTS.q2min,
146 " GeV^2" );
147 else if ( (
FLAGS.pion == 6 ) || (
FLAGS.pion == 7 ) )
148 printf(
" %s %f %s\n",
"min. kaons-tagged photon inv.mass^2 = ",
CUTS.q2min,
149 " GeV^2" );
150 else if (
FLAGS.pion == 9 )
151 printf(
" %s %f %s\n",
" min. lambdas-tagged photon inv.mass^2 = ",
CUTS.q2min,
152 " GeV^2" );
153 else
154 printf(
" %s %f %s\n",
"min. pions-tagged photon inv.mass^2 = ",
CUTS.q2min,
155 " GeV^2" );
156 }
157
158
159
161
162
163 if (
FLAGS.tagged == 0 )
164 {
167 }
168 else
169 {
170 cos1min = -1.0;
171 cos1max = 1.0;
173 }
174 cos2min = -1.0;
175 cos2max = 1.0;
176 cos3min = -1.0;
177 cos3max = 1.0;
178 if (
FLAGS.pion == 0 )
180 else if (
FLAGS.pion == 1 ) qqmin = 4.0 *
CTES.mpi *
CTES.mpi;
181 else if (
FLAGS.pion == 2 )
183 else if (
FLAGS.pion == 3 ) qqmin = 16.0 *
CTES.mpi *
CTES.mpi;
185 else if (
FLAGS.pion == 5 ) qqmin = 4.0 *
CTES.mnt *
CTES.mnt;
186 else if (
FLAGS.pion == 6 ) qqmin = 4.0 *
CTES.mKp *
CTES.mKp;
187 else if (
FLAGS.pion == 7 ) qqmin = 4.0 *
CTES.mKn *
CTES.mKn;
188 else if (
FLAGS.pion == 8 )
190 else if (
FLAGS.pion == 9 ) qqmin = 4.0 *
CTES.mlamb *
CTES.mlamb;
191
192
193
195 if (
CUTS.q2_max_c < qqmax ) qqmax =
CUTS.q2_max_c;
196
197
198 if ( (
CUTS.q2_min_c > qqmin ) &&
201 qqmin =
CUTS.q2_min_c;
202 else
203 {
204 cerr << "------------------------------" << endl;
205 cerr << " Q^2_min TOO SMALL" << endl;
206 cerr << " Q^2_min CHANGE BY PHOKHARA = " << qqmin << " GeV^2" << endl;
207 cerr << "------------------------------" << endl;
208 }
209
210 if ( qqmax <= qqmin )
211 {
212 cerr << " Q^2_max to small " << endl;
213 cerr << " Q^2_max = " << qqmax << endl;
214 cerr << " Q^2_min = " << qqmin << endl;
215 return 0;
216 }
217
218
219 if (
FLAGS.pion == 0 )
220 {
221 printf( " %s %f %s\n", "minimal muon-pair invariant mass^2 = ", qqmin, " GeV^2" );
222 printf( " %s %f %s\n", "maximal muon-pair invariant mass^2 = ", qqmax, " GeV^2" );
223 }
224 else if (
FLAGS.pion == 1 )
225 {
226 printf( " %s %f %s\n", "minimal pion-pair invariant mass^2 = ", qqmin, " GeV^2" );
227 printf( " %s %f %s\n", "maximal pion-pair invariant mass^2 = ", qqmax, " GeV^2" );
228 }
229 else if (
FLAGS.pion == 4 )
230 {
231 printf( " %s %f %s\n", "minimal proton-pair invariant mass^2 = ", qqmin, " GeV^2" );
232 printf( " %s %f %s\n", "maximal proton-pair invariant mass^2 = ", qqmax, " GeV^2" );
233 }
234 else if (
FLAGS.pion == 5 )
235 {
236 printf( " %s %f %s\n", "minimal neutron-pair invariant mass^2 = ", qqmin, " GeV^2" );
237 printf( " %s %f %s\n", "maximal neutron-pair invariant mass^2 = ", qqmax, " GeV^2" );
238 }
239 else if ( (
FLAGS.pion == 6 ) || (
FLAGS.pion == 7 ) )
240 {
241 printf( " %s %f %s\n", "minimal kaon-pair invariant mass^2 = ", qqmin, " GeV^2" );
242 printf( " %s %f %s\n", "maximal kaon-pair invariant mass^2 = ", qqmax, " GeV^2" );
243 }
244 else if (
FLAGS.pion == 8 )
245 {
246 printf( " %s %f %s\n", "minimal three-pion invariant mass^2 = ", qqmin, " GeV^2" );
247 printf( " %s %f %s\n", "maximal three-pion invariant mass^2 = ", qqmax, " GeV^2" );
248 }
249 else if (
FLAGS.pion == 9 )
250 {
251 printf( " %s %f %s\n", "minimal lambda-pair invariant mass^2 = ", qqmin, " GeV^2" );
252 printf( " %s %f %s\n", "maximal lambda-pair invariant mass^2 = ", qqmax, " GeV^2" );
253 }
254 else
255 {
256 printf( " %s %f %s\n", "minimal four-pion invariant mass^2 = ", qqmin, " GeV^2" );
257 printf( " %s %f %s\n", "maximal four-pion invariant mass^2 = ", qqmax, " GeV^2" );
258 }
259
260 if (
FLAGS.nlo == 0 )
261 {
262 cout << "Born" << endl;
263 if (
FLAGS.fsrnlo != 0 )
264 {
265 cerr << "WRONG FSRNLO flag - only fsrnlo=0 allowed for Born" << endl;
266 return 0;
267 }
268 }
269
270 if ( (
FLAGS.pion == 9 ) && (
FLAGS.nlo != 0 ) )
271 {
272 cerr << "WRONG NLO flag - only Born allowed for Lambdas" << endl;
273 cerr << "If you feel that you need NLO" << endl;
274 cerr << "please contact the authors" << endl;
275 return 0;
276 }
277
278 if (
FLAGS.nlo == 1 )
279 printf(
" %s %f\n",
"NLO: soft photon cutoff w = ",
CUTS.w );
280 if ( (
FLAGS.pion <= 1 ) || (
FLAGS.pion == 6 ) )
281 {
282
283 if ( !( (
FLAGS.fsr == 1 ) || (
FLAGS.fsr == 2 ) || (
FLAGS.fsrnlo == 0 ) ||
285 (
FLAGS.fsrnlo == 0 ) ) )
286 {
287 cerr << "WRONG combination of FSR, FSRNLO flags" << endl;
288 return 0;
289 }
290
291
292 if (
FLAGS.fsr == 0 ) cout <<
"ISR only" << endl;
293 else if (
FLAGS.fsr == 1 ) cout <<
"ISR+FSR" << endl;
294 else if (
FLAGS.fsr == 2 )
295 {
296 if (
FLAGS.nlo == 0 ) cout <<
"ISR+INT+FSR" << endl;
297 else
298 {
299 cerr << "WRONG FSR flag: interference is included only for nlo=0" << endl;
300 return 0;
301 }
302 }
303 else
304 {
305 cerr <<
"WRONG FSR flag" <<
FLAGS.fsr << endl;
306 return 0;
307 }
308
309 if (
FLAGS.fsrnlo == 1 ) cout <<
"IFSNLO included" << endl;
310 }
311 else
312 {
313 if ( (
FLAGS.fsr == 0 ) && (
FLAGS.fsrnlo == 0 ) ) cout <<
"ISR only" << endl;
314 else
315 {
316 cerr << "FSR is implemented only for pi+pi-, mu+mu- and K+K- modes" << endl;
317 return 0;
318 }
319 }
320
321
322 if (
FLAGS.ivac == 0 ) { cout <<
"Vacuum polarization is NOT included" << endl; }
323 else if (
FLAGS.ivac == 1 ) { cout <<
"Vacuum polarization is included" << endl; }
324 else
325 {
326 cout << "WRONG vacuum polarization switch" << endl;
327 return 0;
328 }
329
330 if (
FLAGS.pion == 1 )
331 {
332 if (
FLAGS.FF_pion == 0 ) cout <<
"Kuhn-Santamaria PionFormFactor" << endl;
333 else if (
FLAGS.FF_pion == 1 ) cout <<
"Gounaris-Sakurai PionFormFactor" << endl;
334 else
335 {
336 cout << "WRONG PionFormFactor switch" << endl;
337 return 0;
338 }
339
340 if (
FLAGS.fsr != 0 )
341 {
342 if (
FLAGS.f0_model == 0 ) cout <<
"f0+f0(600): K+K- model" << endl;
343 else if (
FLAGS.f0_model == 1 ) cout <<
"f0+f0(600): \"no structure\" model" << endl;
344 else if (
FLAGS.f0_model == 2 ) cout <<
"NO f0+f0(600)" << endl;
345 else if (
FLAGS.f0_model == 3 )
346 cout << "only f0, KLOE: Cesare Bini-private communication" << endl;
347 else
348 {
349 cout << "WRONG f0+f0(600) switch" << endl;
350 return 0;
351 }
352 }
353 }
354
355
356
357 k = nm;
358 for ( i = 0; i < 2; i++ )
359 {
363 }
365
366 for ( i = 1; i <= 2; i++ )
367 {
372
373
374
375 for ( j = 1; j <= k; j++ )
376 {
378 Ar[1] = Ar_r[0];
379
381 {
383 GEN_1PH( i, qqmin, qqmax, cos1min, cos1max, cos3min, cos3max );
384 }
385 else
386 {
388 GEN_2PH( i, qqmin, cos1min, cos1max, cos2min, cos2max, cos3min, cos3max );
389 }
390 }
391
392
393
394 k = nges;
395 if ( i == 1 )
396 {
400
401 if ( (
FLAGS.pion == 1 ) && (
FLAGS.fsrnlo == 1 ) )
403 if ( (
FLAGS.pion == 0 ) && (
FLAGS.fsrnlo == 1 ) )
405
406 if ( (
FLAGS.pion == 0 ) && (
FLAGS.fsr == 1 ) && (
FLAGS.fsrnlo == 0 ) )
408
409 if ( (
FLAGS.pion == 2 ) || (
FLAGS.pion == 3 ) )
410 {
413 }
414
415 if (
FLAGS.pion == 8 )
416 {
419 }
420 }
421 }
422
423
424
425 if (
FLAGS.pion == 9 )
429
430
431
433
434
435 if (
FLAGS.nlo == 0 )
436 {
442 }
443 else
444 {
447 dsigm1 =
452 dsigm2 =
457
458 sigma = sigma1 + sigma2;
459 dsigm = sqrt( dsigm1 * dsigm1 + dsigm2 * dsigm2 );
460 }
461
462
463 cout << "----------------------------------------------------" << endl;
464 cout << int(
MAXIMA.tr[0] +
MAXIMA.tr[1] ) <<
" total events accepted of " << endl;
465 cout << int( nges ) << " total events generated" << endl;
466 cout << int(
MAXIMA.tr[0] ) <<
" one photon events accepted of " << endl;
467 cout << int(
MAXIMA.count[0] ) <<
" events generated" << endl;
468 cout << int(
MAXIMA.tr[1] ) <<
" two photon events accepted of " << endl;
469 cout << int(
MAXIMA.count[1] ) <<
" events generated" << endl;
470 cout << endl;
471 if (
FLAGS.nlo != 0 ) cout <<
"sigma1(nbarn) = " << sigma1 <<
" +- " << dsigm1 << endl;
472 if (
FLAGS.nlo != 0 ) cout <<
"sigma2(nbarn) = " << sigma2 <<
" +- " << dsigm2 << endl;
473 cout << "sigma (nbarn) = " << sigma << " +- " << dsigm << endl;
474 cout << endl;
475 cout <<
"maximum1 = " <<
MAXIMA.gross[0] <<
" minimum1 = " <<
MAXIMA.klein[0] << endl;
476 cout <<
"Mmax1 = " <<
MAXIMA.Mmax[0] << endl;
477 cout <<
"maximum2 = " <<
MAXIMA.gross[1] <<
" minimum2 = " <<
MAXIMA.klein[1] << endl;
478 cout <<
"Mmax2 = " <<
MAXIMA.Mmax[1] << endl;
479 cout << "----------------------------------------------------" << endl;
480
481
482
483
484
485
486
487
488
489
490
491
492 return 0;
493}
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)