41#include "GeneratorObject/McGenEvent.h"
42#include "BesKernel/IBesRndmGenSvc.h"
43#include "CLHEP/Random/RandomEngine.h"
44#include "cfortran/cfortran.h"
53using std::resetiosflags;
54using std::setiosflags;
57int EvtPhokhara_ppbar::nevtgen=0;
58int EvtPhokhara_ppbar::nphokharadecays=0;
60int EvtPhokhara_ppbar::ntable=0;
62int EvtPhokhara_ppbar::ncommand=0;
63int EvtPhokhara_ppbar::lcommand=0;
64std::string* EvtPhokhara_ppbar::commands=0;
65int EvtPhokhara_ppbar::nevt=0;
72 if (nphokharadecays==0) {
78 for(i=0;i<nphokharadecays;i++){
79 if (phokharadecays[i]==
this){
80 phokharadecays[i]=phokharadecays[nphokharadecays-1];
82 if (nphokharadecays==0) {
90 report(
ERROR,
"EvtGen") <<
"Error in destroying Phokhara model!"<<endl;
97 model_name=
"PHOKHARA_ppbar";
121 if(m_pion<0 || m_pion>9){std::cout<<
"mode index for phokhar 0~9, but you give "<<m_pion<<std::endl;abort();}
137 m_q2_min_c = 0.0447 ;
138 m_q2_max_c = m_E*m_E;
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 ;}
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;
161 ctes_.Sp = m_E*m_E; ;
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;
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;
199 cout <<
" PHOKHARA 7.0: not yet implemented" << endl;
202 cout <<
"--------------------------------------------------------------" << endl;
203 printf(
" %s %f %s\n",
"cms total energy = ",sqrt(
ctes_.Sp),
" GeV ");
206 if(
cuts_.gmin<0.001){
207 cerr <<
" minimal missing energy set too small" << endl;
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);
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);
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" );
237 printf(
" %s %f %s\n",
"min. pions-tagged photon inv.mass^2 = ",
cuts_.q2min,
" GeV^2" );
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;
270 qqmin =
cuts_.q2_min_c;
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;
279 cerr <<
" Q^2_max to small " << endl;
280 cerr <<
" Q^2_max = " << qqmax << endl;
281 cerr <<
" Q^2_min = " << qqmin << endl;
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");
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");
313 cout <<
"Born" << endl;
315 cerr <<
"WRONG FSRNLO flag - only fsrnlo=0 allowed for Born" << endl;
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;
327 if (
flags_.nlo == 1) printf(
" %s %f\n",
"NLO: soft photon cutoff w = ",
cuts_.w);
334 cerr <<
"WRONG combination of FSR, FSRNLO flags" <<endl;
339 cout <<
"ISR only" << endl;
341 cout <<
"ISR+FSR" << endl;
342 else if (
flags_.fsr == 2) {
344 cout <<
"ISR+INT+FSR" << endl;
346 cerr <<
"WRONG FSR flag: interference is included only for nlo=0" << endl;
351 cerr <<
"WRONG FSR flag" <<
flags_.fsr << endl;
355 cout <<
"IFSNLO included" << endl;
360 cout <<
"ISR only" << endl;
362 cerr <<
"FSR is implemented only for pi+pi-, mu+mu- and K+K- modes" << endl;
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;}
374 cout <<
"WRONG vacuum polarization switch" << endl;
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;
387 cout <<
"WRONG PionFormFactor switch" << endl;
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;
400 cout <<
"WRONG f0+f0(600) switch" << endl;
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;}
414 cout <<
"WRONG KaonFormFactor switch" << endl;
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;
430 cout <<
"--------------------------------------------------------------" << endl;
434 for(
int i = 0; i<2; i++)
451 for(
int j = 1; j <= m_nm; j++)
457 GEN_1PH(1,qqmin,qqmax,cos1min,cos1max,cos3min,cos3max);
461 GEN_2PH(1,qqmin,cos1min,cos1max,cos2min,cos2max,cos3min,cos3max);
501 std::string locvp=getenv(
"BESEVTGENROOT");
502 system(
"cat $BESEVTGENROOT/share/phokhara.param>phokhara.param");
507 report(
ERROR,
"EvtGen") <<
"EvtPhokhara finds that you are decaying the"<<endl
508 <<
" aliased particle "
510 <<
" with the Phokhara model"<<endl
511 <<
" this does not work, please modify decay table."
513 report(
ERROR,
"EvtGen") <<
"Will terminate execution!"<<endl;
525 return std::string(
"PhokharaPar");
531 if (ncommand==lcommand){
533 lcommand=10+2*lcommand;
535 std::string* newcommands=
new std::string[lcommand];
539 for(i=0;i<ncommand;i++){
540 newcommands[i]=commands[i];
545 commands=newcommands;
549 commands[ncommand]=cmd;
559 if(p->
getId()!=myvpho) {std::cout<<
"Parent particle is required to be vpho for Phokhara model"<<std::endl;abort();}
562 std::cout<<
"PHOKHARA : ppbar mode "<<std::endl;
568 tr_old[0] = (int)
maxima_.tr[0];
569 tr_old[1] = (
int)
maxima_.tr[1];
571 while( ntrials < 10000)
579 GEN_1PH(2,qqmin,qqmax,cos1min,cos1max,cos3min,cos3max);
583 GEN_2PH(2,qqmin,cos1min,cos1max,cos2min,cos2max,cos3min,cos3max);
586 if( ((
int)
maxima_.tr[0]+(
int)
maxima_.tr[1]) > (tr_old[0]+tr_old[1]) )
593 std::cout <<
"FATAL: Could not satisfy cuts after " << ntrials <<
"trials. Terminate." <<std::endl;
599 EvtId evtnumstable[100];
736 if(
ctes_.momenta[0][3] != 0 )
745 if(more) {std::cout<<
"Existence of mode "<<channel<<
" in exclusive decay list has the same final state as this one"<<std::endl;abort(); }
752 for(
int i=0;i<numstable;i++){
756 if ( ndaugFound == 0 ) {
757 report(
ERROR,
"EvtGen") <<
"Phokhara has failed to do a decay ";
771 if (nphokharadecays==ntable){
775 for(i=0;i<ntable;i++){
776 newphokharadecays[i]=phokharadecays[i];
779 delete [] phokharadecays;
780 phokharadecays=newphokharadecays;
783 phokharadecays[nphokharadecays++]=jsdecay;
809 if(m_pion<0 || m_pion>9){std::cout<<
"mode index for phokhar 0~9, but you give "<<m_pion<<std::endl;abort();}
825 m_q2_min_c = 0.0447 ;
826 m_q2_max_c = m_E*m_E;
833 if(!(m_pion==0 || m_pion==1 || m_pion==6)){m_fsr = 0; m_fsrnlo = 0 ;}
834 if( m_pion==9 ){m_nlo = 0 ;}
842 flags_.FF_pion = m_FF_Pion;
843 flags_.f0_model = m_f0_model;
844 flags_.FF_kaon = m_FF_Kaon;
845 flags_.narr_res = m_NarrowRes;
847 ctes_.Sp = m_E*m_E; ;
850 cuts_.q2min = m_q2min;
851 cuts_.q2_min_c = m_q2_min_c;
852 cuts_.q2_max_c = m_q2_max_c;
854 cuts_.phot1cut = m_phot1cut;
855 cuts_.phot2cut = m_phot2cut;
856 cuts_.pi1cut = m_pi1cut;
857 cuts_.pi2cut = m_pi2cut;
869 else if (
flags_.pion == 1)
871 else if (
flags_.pion == 2)
873 else if (
flags_.pion == 3)
875 else if (
flags_.pion == 4)
877 else if (
flags_.pion == 5)
879 else if (
flags_.pion == 6)
881 else if (
flags_.pion == 7)
883 else if (
flags_.pion == 8)
885 else if (
flags_.pion == 9)
888 if (
cuts_.q2_max_c < qqmax)
889 qqmax=
cuts_.q2_max_c;
893 qqmin =
cuts_.q2_min_c;
900 for(
int i = 0; i<2; i++)
#define RLXDINIT(LUXURY, SEED)
#define GEN_1PH(I, QQMIN, QQMAX, COS1MIN, COS1MAX, COS3MIN, COS3MAX)
#define GEN_2PH(I, QQMIN, COS1MIN, COS1MAX, COS2MIN, COS2MAX, COS3MIN, COS3MAX)
ostream & report(Severity severity, const char *facility)
double cos(const BesAngle a)
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
static int inChannelList(EvtId parent, int ndaug, EvtId *daugs)
static int getStdHep(EvtId id)
static EvtId evtIdFromStdHep(int stdhep)
static std::string name(EvtId i)
static EvtId getId(const std::string &name)
void makeDaughters(int ndaug, EvtId *id)
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
EvtParticle * getDaug(int i)
void command(std::string cmd)
void decay(EvtParticle *p)
void init_evt(EvtParticle *p)
void getName(std::string &name)
void PhokharaInit(int dummy)
virtual ~EvtPhokhara_ppbar()
std::string commandName()
void init_mode(EvtParticle *p)
void set(int i, double d)