44#include "CLHEP/Random/RandomEngine.h"
45#include "cfortran/cfortran.h"
54using std::resetiosflags;
55using std::setiosflags;
58int EvtPhokhara_Lambda::nevtgen=0;
59int EvtPhokhara_Lambda::nphokharadecays=0;
61int EvtPhokhara_Lambda::ntable=0;
63int EvtPhokhara_Lambda::ncommand=0;
64int EvtPhokhara_Lambda::lcommand=0;
65std::string* EvtPhokhara_Lambda::commands=0;
66int EvtPhokhara_Lambda::nevt=0;
73 if (nphokharadecays==0) {
79 for(i=0;i<nphokharadecays;i++){
80 if (phokharadecays[i]==
this){
81 phokharadecays[i]=phokharadecays[nphokharadecays-1];
83 if (nphokharadecays==0) {
91 report(
ERROR,
"EvtGen") <<
"Error in destroying Phokhara model!"<<endl;
98 model_name=
"PHOKHARA_LAMBDA";
122#include "Phokhara_init_mode.txt"
130 std::string locvp=getenv(
"BESEVTGENROOT");
131 system(
"cat $BESEVTGENROOT/share/phokhara_10.0.param>phokhara_10.0.param");
132 system(
"cat $BESEVTGENROOT/share/phokhara_10.0.fferr>phokhara_10.0.fferr");
133 system(
"cat $BESEVTGENROOT/share/phokhara_10.0.ffwarn>phokhara_10.0.ffwarn");
138 report(
ERROR,
"EvtGen") <<
"EvtPhokhara finds that you are decaying the"<<endl
139 <<
" aliased particle "
141 <<
" with the Phokhara model"<<endl
142 <<
" this does not work, please modify decay table."
144 report(
ERROR,
"EvtGen") <<
"Will terminate execution!"<<endl;
156 return std::string(
"PhokharaPar");
162 if (ncommand==lcommand){
164 lcommand=10+2*lcommand;
166 std::string* newcommands=
new std::string[lcommand];
170 for(i=0;i<ncommand;i++){
171 newcommands[i]=commands[i];
176 commands=newcommands;
180 commands[ncommand]=cmd;
191 if(p->
getId()!=myvpho) {std::cout<<
"Parent particle is required to be vpho for Phokhara model"<<std::endl;abort();}
193 std::cout<<
"PHOKHARA : Lambda anti-Lambda mode "<<std::endl;
196 if(debug)cout<<
"init_evt works!"<<endl;
200 std::cout<<
"flags :ph0,nlo,pion,fsr,fsrnlo,ivac,FF_pion,f0_model,FF_kaon,narr_res,FF_pp,chi_sw,chi_pion,FF_Pgg,nlo2 "<<std::endl;
201 std::cout<<
"= "<<
flags_.ph0<<
","<<
flags_.nlo<<
","<<
flags_.pion<<
","<<
flags_.fsr<<
","<<
flags_.fsrnlo<<
","<<
flags_.ivac<<
","<<
flags_.FF_pion<<
","<<
flags_.f0_model<<
","<<
flags_.FF_kaon<<
","<<
flags_.narr_res<<
","<<
flags_.FF_pp<<
","<<
flags_.chi_sw<<
","<<
flags_.be_r<<
","<<
flags_.FF_Pgg<<
","<<
flags_.nlo2<<std::endl;
202 std::cout<<
"ctes: Sp = "<<
ctes_.Sp<<std::endl;
203 std::cout<<
"cuts : w,q2min,q2_min_c,gmin,phot1cut,phot2cut,pi1cut,pi2cut"<<std::endl<<
"= "<<
cuts_.w<<
","<<
cuts_.q2min<<
","<<
cuts_.q2_min_c<<
","<<
cuts_.gmin<<
","<<
cuts_.phot1cut<<
","<<
cuts_.phot2cut<<
","<<
cuts_.pi1cut<<
","<<
cuts_.pi2cut<<std::endl;
211 tr_old[0] = (int)
maxima_.tr[0];
212 tr_old[1] = (int)
maxima_.tr[1];
213 tr_old[2] = (int)
maxima_.tr[2];
215 while( ntrials < 1000000)
227 GEN_1PH(2,qqmin,qqmax,cos1min,cos1max,cos3min,cos3max);
231 GEN_2PH(2,qqmin,cos1min,cos1max,cos2min,cos2max,cos3min,cos3max);
240 std::cout <<
"FATAL: Could not satisfy cuts after " << ntrials <<
"trials. Terminate." <<std::endl;
246 EvtId evtnumstable[100];
379 EvtId LambdaMode[100];
388 if (debug)std::cout<<
"Phokhara_Lambda: anti-Lambda0 p4[numstable] = "<<
p4[numstable]<<std::endl;
393 if (debug) std::cout<<
"Phokhara_Lambda: Lambda0 p4[numstable] = "<<
p4[numstable]<<std::endl;
399 if (debug) std::cout<<
"Phokhara_Lambda: Pi+ p4[numstable] = "<<Lambdap4[LambdaDaus]<<std::endl;
402 tmp.
set(
p4[0].get(0),-
p4[0].get(1),-
p4[0].get(2),-
p4[0].get(3));
403 Lambdap4[LambdaDaus] =
boostTo(Lambdap4[LambdaDaus],tmp);
404 if (debug) std::cout<<
"Phokhara_Lambda:Boosted Pi+ p4[numstable] = "<<Lambdap4[LambdaDaus]<<std::endl;
409 if (debug) std::cout<<
"Phokhara_Lambda: pbar p4[numstable] = "<<Lambdap4[LambdaDaus]<<std::endl;
413 tmp.
set(
p4[0].get(0),-
p4[0].get(1),-
p4[0].get(2),-
p4[0].get(3));
414 Lambdap4[LambdaDaus] =
boostTo(Lambdap4[LambdaDaus],tmp);
415 if (debug) std::cout<<
"Phokhara_Lambda:Boosted pbar p4[numstable] = "<<Lambdap4[LambdaDaus]<<std::endl;
420 if (debug) std::cout<<
"Phokhara_Lambda: pi- p4[numstable] = "<<Lambdap4[LambdaDaus]<<std::endl;
424 tmp.
set(
p4[1].get(0),-
p4[1].get(1),-
p4[1].get(2),-
p4[1].get(3));
425 Lambdap4[LambdaDaus] =
boostTo(Lambdap4[LambdaDaus],tmp);
426 if (debug) std::cout<<
"Phokhara_Lambda:Boosted pi- p4[numstable] = "<<Lambdap4[LambdaDaus]<<std::endl;
430 Lambdap4[LambdaDaus].
set(
ctes_.momenta[0][10],
ctes_.momenta[1][10],
ctes_.momenta[2][10],
ctes_.momenta[3][10]);
431 if (debug) std::cout<<
"Phokhara_Lambda: p p4[numstable] = "<<Lambdap4[LambdaDaus]<<std::endl;
435 tmp.
set(
p4[1].get(0),-
p4[1].get(1),-
p4[1].get(2),-
p4[1].get(3));
436 Lambdap4[LambdaDaus] =
boostTo(Lambdap4[LambdaDaus],tmp);
437 if (debug) std::cout<<
"Phokhara_Lambda:Boosted p p4[numstable] = "<<Lambdap4[LambdaDaus]<<std::endl;
444 if (debug) std::cout<<
"Phokhara_Lambda: first gamma p4[numstable] = "<<
p4[numstable]<<std::endl;
446 if(
ctes_.momenta[0][3] != 0 )
450 if (debug) std::cout<<
"Phokhara_Lambda: second gamma p4[numstable] = "<<
p4[numstable]<<std::endl;
456 if(more) {std::cout<<
"EvtPhokhara_Lambda:Existence of mode "<<channel<<
" in exclusive decay list has the same final state as this one"<<std::endl;abort(); }
463 for(
int i = 0; i< LambdaDaus;i++){
470 for(
int i=0;i<numstable;i++){
474 if ( ndaugFound == 0 ) {
475 report(
ERROR,
"EvtGen") <<
"Phokhara has failed to do a decay ";
479 if (debug) std::cout<<
"EvtPhokhara_Lambda SUMMARY: part p4"<<p->
getP4Lab()<<std::endl;
480 if (debug) std::cout<<
"EvtPhokhara_Lambda SUMMARY: Daug0 p4"<<p->
getDaug(0)->
getP4Lab()<<std::endl;
481 if (debug) std::cout<<
"EvtPhokhara_Lambda SUMMARY: Daug1 p4"<<p->
getDaug(1)->
getP4Lab()<<std::endl;
482 if (debug) std::cout<<
"EvtPhokhara_Lambda SUMMARY: Daug2 p4"<<p->
getDaug(2)->
getP4Lab()<<std::endl;
493 if (nphokharadecays==ntable){
497 for(i=0;i<ntable;i++){
498 newphokharadecays[i]=phokharadecays[i];
501 delete [] phokharadecays;
502 phokharadecays=newphokharadecays;
505 phokharadecays[nphokharadecays++]=jsdecay;
531#include "Phokhara_init_evt.txt"
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
#define GEN_1PH(I, QQMIN, QQMAX, COS1MIN, COS1MAX, COS3MIN, COS3MAX)
#define GEN_0PH(I, QQMIN, SP, COS3MIN, COS3MAX)
#define GEN_2PH(I, QQMIN, COS1MIN, COS1MAX, COS2MIN, COS2MAX, COS3MIN, COS3MAX)
ostream & report(Severity severity, const char *facility)
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
static EvtId _NextLevelId[20]
static int _NextLevelDauNum
EvtParticle * getDaug(int i)
static EvtVector4R _NextLevelP4[20]
void PhokharaInit(int dummy)
std::string commandName()
void getName(std::string &name)
void init_evt(EvtParticle *p)
virtual ~EvtPhokhara_Lambda()
void decay(EvtParticle *p)
void command(std::string cmd)
void init_mode(EvtParticle *p)
void set(int i, double d)
double double double * p4