1#include "Phokhara/PhokharaDef.h"
7#include "Phokhara/ranlxd.h"
12#define RLXDRESETF(SEED) CCALLSFSUB1(RLXDRESETF,rlxdresetf,INTV, SEED)
15#define INPUT(NGES,NM,OUTFILE) CCALLSFSUB3(INPUT,input,PLONG,PINT,PSTRING,NGES,NM,OUTFILE)
18#define INITHISTO() CCALLSFSUB0(INITHISTO,inithisto)
21#define ENDHISTO() CCALLSFSUB0(ENDHISTO,endhisto)
24#define WRITEEVENT() CCALLSFSUB0(WRITEEVENT,writeevent)
27#define RANLXDF(AR, VAL) CCALLSFSUB2(RANLXDF,ranlxdf,DOUBLEV, INT, AR, VAL)
30#define GEN_1PH(I,QQMIN,QQMAX,COS1MIN,COS1MAX,COS3MIN,COS3MAX) CCALLSFSUB7(GEN_1PH,gen_1ph_,INT, PDOUBLE, PDOUBLE, PDOUBLE, PDOUBLE, PDOUBLE, PDOUBLE,I,QQMIN,QQMAX,COS1MIN,COS1MAX,COS3MIN,COS3MAX)
33#define GEN_2PH(I,QQMIN,COS1MIN,COS1MAX,COS2MIN,COS2MAX,COS3MIN,COS3MAX) CCALLSFSUB8(GEN_2PH,gen_2ph_,INT, PDOUBLE, PDOUBLE, PDOUBLE, PDOUBLE, PDOUBLE, PDOUBLE, PDOUBLE,I,QQMIN,COS1MIN,COS1MAX,COS2MIN,COS2MAX,COS3MIN,COS3MAX)
37 double cos1min,cos1max,cos2min,cos2max,cos3min,cos3max;
38 double dsigm1,dsigm2,sigma1,sigma2,sigma,dsigm,Ar[14],Ar_r[14];
44 ifstream seeds(
"seed.dat");
48 while( ! seeds.eof() )
49 seeds >> s_seed[ii++];
52 cerr <<
"Cannot open file seed.dat for reading" << endl;
59 INPUT(nges, nm, outfile);
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;
88 cout <<
" PHOKHARA 6.0: not yet implemented" << endl;
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;
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);
102 if((
CUTS.gmin/2.0/
CTES.ebeam) < 0.0098){
103 cerr <<
" minimal missing energy set to small" << endl;
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);
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);
123 printf(
" %s %f,%f\n",
"angular cuts on pions = ",
CUTS.pi1cut,
CUTS.pi2cut);
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" );
137 printf(
" %s %f %s\n",
"min. pions-tagged photon inv.mass^2 = ",
CUTS.q2min,
" GeV^2" );
145 if (
FLAGS.tagged == 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)
181 if (
CUTS.q2_max_c < qqmax)
186 qqmin =
CUTS.q2_min_c;
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;
195 cerr <<
" Q^2_max to small " << endl;
196 cerr <<
" Q^2_max = " << qqmax << endl;
197 cerr <<
" Q^2_min = " << qqmin << endl;
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");
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");
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;
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;
243 if (
FLAGS.nlo == 1) printf(
" %s %f\n",
"NLO: soft photon cutoff w = ",
CUTS.w);
249 || (
FLAGS.fsr == 0) || (
FLAGS.fsrnlo == 0))) {
250 cerr <<
"WRONG combination of FSR, FSRNLO flags" <<endl;
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;
263 cerr <<
"WRONG FSR flag: interference is included only for nlo=0" << endl;
268 cerr <<
"WRONG FSR flag" <<
FLAGS.fsr << endl;
272 if(
FLAGS.fsrnlo == 1)
273 cout <<
"IFSNLO included" << endl;
278 cout <<
"ISR only" << endl;
280 cerr <<
"FSR is implemented only for pi+pi-, mu+mu- and K+K- modes" << endl;
287 cout <<
"Vacuum polarization is NOT included" << endl;
288 }
else if(
FLAGS.ivac == 1){
289 cout <<
"Vacuum polarization is included" << endl;
291 cout <<
"WRONG vacuum polarization switch" << endl;
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;
301 cout <<
"WRONG PionFormFactor switch" << endl;
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;
315 cout <<
"WRONG f0+f0(600) switch" << endl;
324 for( i = 0; i<2; i++)
333 for( i = 1; i<=2; i++) {
341 for(j = 1; j <= k; j++)
348 GEN_1PH(i,qqmin,qqmax,cos1min,cos1max,cos3min,cos3max);
352 GEN_2PH(i,qqmin,cos1min,cos1max,cos2min,cos2max,cos3min,cos3max);
394 if (
FLAGS.nlo == 0) {
403 sigma = sigma1+sigma2;
404 dsigm = sqrt(dsigm1*dsigm1 + dsigm2*dsigm2);
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;
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;
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;
#define PROTOCCALLSFSUB3(UN, LN, T1, T2, T3)
#define PROTOCCALLSFSUB2(UN, LN, T1, T2)
#define PROTOCCALLSFSUB7(UN, LN, T1, T2, T3, T4, T5, T6, T7)
#define PROTOCCALLSFSUB0(UN, LN)
#define PROTOCCALLSFSUB1(UN, LN, T1)
#define PROTOCCALLSFSUB8(UN, LN, T1, T2, T3, T4, T5, T6, T7, T8)
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)