115 {
116 m_pion=1;
117
118
119
120
121 if(m_pion<0 || m_pion>9){std::cout<<"mode index for phokhar 0~9, but you give "<<m_pion<<std::endl;abort();}
124
125 m_tagged = 0;
126 m_nm = 50000 ;
127 m_nlo = 1;
128 m_w = 0.0001;
129 m_fsr = 1;
130 m_fsrnlo = 1 ;
131 m_NarrowRes = 0 ;
132 m_FF_Kaon = 1 ;
133 m_ivac = 0;
134 m_FF_Pion = 0 ;
135 m_f0_model = 0 ;
136 m_q2min = 0.0;
137 m_q2_min_c = 0.0447 ;
138 m_q2_max_c = m_E*m_E;
139 m_gmin = 0.001;
140 m_phot1cut = 0.0;
141 m_phot2cut = 180.0;
142 m_pi1cut = 0.0 ;
143 m_pi2cut = 180.0;
144
145 if(!(m_pion==0 || m_pion==1 || m_pion==6)){m_fsr = 0; m_fsrnlo = 0 ;}
146 if( m_pion==9 ){m_nlo = 0 ;}
147
148 m_initSeed = 123456789;
156 flags_.FF_pion = m_FF_Pion;
157 flags_.f0_model = m_f0_model;
158 flags_.FF_kaon = m_FF_Kaon;
159 flags_.narr_res = m_NarrowRes;
160
161 ctes_.Sp = m_E*m_E; ;
162
164 cuts_.q2min = m_q2min;
165 cuts_.q2_min_c = m_q2_min_c;
166 cuts_.q2_max_c = m_q2_max_c;
168 cuts_.phot1cut = m_phot1cut;
169 cuts_.phot2cut = m_phot2cut;
170 cuts_.pi1cut = m_pi1cut;
171 cuts_.pi2cut = m_pi2cut;
172
174
175
176 cout << "-------------------------------------------------------------" << endl;
178 cout << " PHOKHARA 7.0 : e^+ e^- -> mu^+ mu^- gamma" << endl;
179 else if (
flags_.pion == 1)
180 cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- gamma" << endl;
181 else if (
flags_.pion == 2)
182 cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- 2pi^0 gamma" << endl;
183 else if (
flags_.pion == 3)
184 cout << " PHOKHARA 7.0: e^+ e^- -> 2pi^+ 2pi^- gamma" << endl;
185 else if (
flags_.pion == 4)
186 cout << " PHOKHARA 7.0: e^+ e^- -> p pbar gamma" << endl;
187 else if (
flags_.pion == 5)
188 cout << " PHOKHARA 7.0: e^+ e^- -> n nbar gamma" << endl;
189 else if (
flags_.pion == 6)
190 cout << " PHOKHARA 7.0: e^+ e^- -> K^+ K^- gamma" << endl;
191 else if (
flags_.pion == 7)
192 cout << " PHOKHARA 7.0: e^+ e^- -> K_0 K_0bar gamma" << endl;
193 else if (
flags_.pion == 8)
194 cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- pi^0 gamma" << endl;
195 else if (
flags_.pion == 9) {
196 cout << "PHOKHARA 7.0 : e^+ e^- ->" << endl;
197 cout << " Lambda (-> pi^- p) Lambda bar (-> pi^+ pbar) gamma" << endl;
198 } else
199 cout << " PHOKHARA 7.0: not yet implemented" << endl;
200
201
202 cout << "--------------------------------------------------------------" << endl;
203 printf(
" %s %f %s\n",
"cms total energy = ",sqrt(
ctes_.Sp),
" GeV ");
204
205
206 if(
cuts_.gmin<0.001){
207 cerr << " minimal missing energy set too small" << endl;
208 abort();
209 }
210 printf(
" %s %f %s\n",
"minimal tagged photon energy = ",
cuts_.gmin,
" GeV ");
211 printf(
" %s %f,%f\n",
"angular cuts on tagged photon = ",
cuts_.phot1cut,
cuts_.phot2cut);
212
213
215 printf(
" %s %f,%f\n",
"angular cuts on muons = ",
cuts_.pi1cut,
cuts_.pi2cut);
216 else if (
flags_.pion == 4)
217 printf(
" %s %f,%f\n",
"angular cuts on protons = ",
cuts_.pi1cut,
cuts_.pi2cut);
218 else if (
flags_.pion == 5)
219 printf(
" %s %f,%f\n",
"angular cuts on neutrons = ",
cuts_.pi1cut,
cuts_.pi2cut);
221 printf(
" %s %f,%f\n",
"angular cuts on kaons = ",
cuts_.pi1cut,
cuts_.pi2cut);
222 else if (
flags_.pion == 9)
223 printf(
" %s %f,%f\n",
"angular cuts on pions and protons = ",
cuts_.pi1cut,
cuts_.pi2cut);
224 else
225 printf(
" %s %f,%f\n",
"angular cuts on pions = ",
cuts_.pi1cut,
cuts_.pi2cut);
227 printf(
" %s %f %s\n",
"min. muons-tagged photon inv.mass^2 = ",
cuts_.q2min,
" GeV^2");
228 else if (
flags_.pion == 4)
229 printf(
" %s %f %s\n",
"min. protons-tagged photon inv.mass^2 = ",
cuts_.q2min,
" GeV^2" );
230 else if (
flags_.pion == 5)
231 printf(
" %s %f %s\n",
"min. neutrons-tagged photon inv.mass^2 = ",
cuts_.q2min,
" GeV^2" );
233 printf(
" %s %f %s\n",
"min. kaons-tagged photon inv.mass^2 = ",
cuts_.q2min,
" GeV^2" );
234 else if (
flags_.pion == 9)
235 printf(
" %s %f %s\n",
" min. lambdas-tagged photon inv.mass^2 = ",
cuts_.q2min,
" GeV^2" );
236 else
237 printf(
" %s %f %s\n",
"min. pions-tagged photon inv.mass^2 = ",
cuts_.q2min,
" GeV^2" );
240 cos2min = -1.0;
241 cos2max = 1.0;
242 cos3min = -1.0;
243 cos3max = 1.0;
246 else if (
flags_.pion == 1)
248 else if (
flags_.pion == 2)
250 else if (
flags_.pion == 3)
252 else if (
flags_.pion == 4)
254 else if (
flags_.pion == 5)
256 else if (
flags_.pion == 6)
258 else if (
flags_.pion == 7)
260 else if (
flags_.pion == 8)
262 else if (
flags_.pion == 9)
265 if (
cuts_.q2_max_c < qqmax)
266 qqmax=
cuts_.q2_max_c;
267
268
270 qqmin =
cuts_.q2_min_c;
271 else {
272 cerr << "------------------------------" << endl;
273 cerr << " Q^2_min TOO SMALL" << endl;
274 cerr << " Q^2_min CHANGED BY PHOKHARA = " << qqmin << " GeV^2" << endl;
275 cerr << "------------------------------" << endl;
276 }
277
278 if(qqmax <= qqmin){
279 cerr << " Q^2_max to small " << endl;
280 cerr << " Q^2_max = " << qqmax << endl;
281 cerr << " Q^2_min = " << qqmin << endl;
282 abort();
283 }
284
285
287 printf(" %s %f %s\n", "minimal muon-pair invariant mass^2 = ", qqmin," GeV^2");
288 printf(" %s %f %s\n", "maximal muon-pair invariant mass^2 = ", qqmax," GeV^2");
289 }
else if (
flags_.pion == 1) {
290 printf(" %s %f %s\n", "minimal pion-pair invariant mass^2 = ", qqmin," GeV^2");
291 printf(" %s %f %s\n", "maximal pion-pair invariant mass^2 = ", qqmax," GeV^2");
292 }
else if (
flags_.pion == 4) {
293 printf(" %s %f %s\n", "minimal proton-pair invariant mass^2 = ", qqmin," GeV^2");
294 printf(" %s %f %s\n", "maximal proton-pair invariant mass^2 = ", qqmax," GeV^2");
295 }
else if (
flags_.pion == 5) {
296 printf(" %s %f %s\n", "minimal neutron-pair invariant mass^2 = ", qqmin," GeV^2");
297 printf(" %s %f %s\n", "maximal neutron-pair invariant mass^2 = ", qqmax," GeV^2");
299 printf(" %s %f %s\n", "minimal kaon-pair invariant mass^2 = ", qqmin," GeV^2");
300 printf(" %s %f %s\n", "maximal kaon-pair invariant mass^2 = ", qqmax," GeV^2");
301 }
else if(
flags_.pion == 8){
302 printf(" %s %f %s\n", "minimal three-pion invariant mass^2 = ", qqmin," GeV^2");
303 printf(" %s %f %s\n", "maximal three-pion invariant mass^2 = ", qqmax," GeV^2");
304 }
else if(
flags_.pion == 9){
305 printf(" %s %f %s\n", "minimal lambda-pair invariant mass^2 = ", qqmin," GeV^2");
306 printf(" %s %f %s\n", "maximal lambda-pair invariant mass^2 = ", qqmax," GeV^2");
307 } else {
308 printf(" %s %f %s\n", "minimal four-pion invariant mass^2 = ", qqmin," GeV^2" );
309 printf(" %s %f %s\n", "maximal four-pion invariant mass^2 = ", qqmax," GeV^2");
310 }
311
313 cout << "Born" << endl;
315 cerr << "WRONG FSRNLO flag - only fsrnlo=0 allowed for Born" << endl;
316 abort();
317 }
318 }
319
321 cerr << "WRONG NLO flag - only Born allowed for Lambdas"<< endl;
322 cerr << "If you feel that you need NLO"<< endl;
323 cerr << "please contact the authors"<< endl;
324 abort();
325 }
326
327 if (
flags_.nlo == 1) printf(
" %s %f\n",
"NLO: soft photon cutoff w = ",
cuts_.w);
329 {
330
334 cerr << "WRONG combination of FSR, FSRNLO flags" <<endl;
335 abort();
336 }
337
339 cout << "ISR only" << endl;
341 cout << "ISR+FSR" << endl;
342 else if (
flags_.fsr == 2) {
344 cout << "ISR+INT+FSR" << endl;
345 else {
346 cerr << "WRONG FSR flag: interference is included only for nlo=0" << endl;
347 abort();
348 }
349 }
350 else {
351 cerr <<
"WRONG FSR flag" <<
flags_.fsr << endl;
352 abort();
353 }
355 cout << "IFSNLO included" << endl;
356 }
357 else
358 {
360 cout << "ISR only" << endl;
361 else {
362 cerr << "FSR is implemented only for pi+pi-, mu+mu- and K+K- modes" << endl;
363 abort();
364 }
365 }
366
368 cout << "Vacuum polarization is NOT included" << endl;}
369 else if(
flags_.ivac == 1){
370 cout << "Vacuum polarization by Fred Jegerlehner (http://www-com.physik.hu-berlin.de/fjeger/alphaQEDn.uu)" << endl;}
371 else if(
flags_.ivac == 2){
372 cout << "Vacuum polarization (VP_HLMNT_v1_3nonr) by Daisuke Nomura and Thomas Teubner" << endl;}
373 else {
374 cout << "WRONG vacuum polarization switch" << endl;
375 abort();
376 }
377
378
381 cout << "Kuhn-Santamaria PionFormFactor" << endl;
382 else if(
flags_.FF_pion == 1)
383 cout << "Gounaris-Sakurai PionFormFactor old" << endl;
384 else if(
flags_.FF_pion == 2)
385 cout << "Gounaris-Sakurai PionFormFactor new" << endl;
386 else {
387 cout << "WRONG PionFormFactor switch" << endl;
388 abort();
389 }
392 cout << "f0+f0(600): K+K- model" << endl;
393 else if(
flags_.f0_model == 1)
394 cout << "f0+f0(600): \"no structure\" model" << endl;
395 else if(
flags_.f0_model == 2)
396 cout << "NO f0+f0(600)" << endl;
397 else if(
flags_.f0_model == 3)
398 cout << "only f0, KLOE: Cesare Bini-private communication" << endl;
399 else {
400 cout << "WRONG f0+f0(600) switch" << endl;
401 abort();
402 }
403 }
404 }
405
408 cout << "constrained KaonFormFactor" << endl;
409 else if(
flags_.FF_kaon == 1) {
410 cout << "unconstrained KaonFormFactor" << endl;}
411 else if(
flags_.FF_kaon == 2) {
412 cout << "Kuhn-Khodjamirian-Bruch KaonFormFactor" << endl;}
413 else{
414 cout << "WRONG KaonFormFactor switch" << endl;
415 abort();
416 }
417 }
418
421 cout << "THE NARROW RESONANCE J/PSI INCLUDED" << endl;}
422 else if(
flags_.narr_res == 2){
423 cout << "THE NARROW RESONANCE PSI(2S) INCLUDED" << endl;}
424 else if(
flags_.narr_res != 0){
425 cout << "wrong flag narr_res: only 0, 1 or 2 allowed" << endl;
426 abort();
427 }
428 }
429
430 cout << "--------------------------------------------------------------" << endl;
431
432
433
434 for(int i = 0; i<2; i++)
435 {
439 }
440
443
448
449
450
451 for(int j = 1; j <= m_nm; j++)
452 {
454 Ar[1] = Ar_r[0];
457 GEN_1PH(1,qqmin,qqmax,cos1min,cos1max,cos3min,cos3max);
458 }
459 else {
461 GEN_2PH(1,qqmin,cos1min,cos1max,cos2min,cos2max,cos3min,cos3max);
462 }
463 }
464
465
466
475
478
482 }
483
487 }
488
489
494}
#define RLXDINIT(LUXURY, SEED)