22#include "EvtGenBase/EvtPatches.hh"
23#include "EvtGenBase/EvtPatches.hh"
24#include "EvtGenBase/EvtParticle.hh"
25#include "EvtGenBase/EvtStringParticle.hh"
26#include "EvtGenBase/EvtDecayTable.hh"
27#include "EvtGenBase/EvtPDL.hh"
28#include "EvtGenModels/EvtLunda.hh"
29#include "EvtGenBase/EvtReport.hh"
31#include "EvtGenBase/EvtId.hh"
43using std::resetiosflags;
44using std::setiosflags;
48int EvtLunda::nlundadecays=0;
50int EvtLunda::ntable=0;
52int EvtLunda::ncommand=0;
53int EvtLunda::lcommand=0;
54std::string* EvtLunda::commands=0;
59 extern void lundar_(
int *,
int *,
float *,
int *,
int *,
int *,
60 double *,
double *,
double *,
double *);
68 extern void lugive0_(
const char *cnfgstr,
int length);
99 if (nlundadecays==0) {
105 for(i=0;i<nlundadecays;i++){
106 if (lundadecays[i]==
this){
107 lundadecays[i]=lundadecays[nlundadecays-1];
109 if (nlundadecays==0) {
117 report(
ERROR,
"EvtGen") <<
"Error in destroying Lunda model!"<<endl;
145 std::string strpdg=getenv(
"BESEVTGENROOT");
146 strpdg +=
"/share/r_pdg.dat";
148 std::cout<<
"mypdg= "<<strpdg<<std::endl;
150 if(
getNArg()>2){std::cout<<
"Parameter can be 1 or 2, You have "<<
getNArg()<<std::endl; ::abort();}
154 report(
ERROR,
"EvtGen") <<
"EvtLunda finds that you are decaying the"<<endl
155 <<
" aliased particle "
157 <<
" with the Lunda model"<<endl
158 <<
" this does not work, please modify decay table."
160 report(
ERROR,
"EvtGen") <<
"Will terminate execution!"<<endl;
172 return std::string(
"LundaPar");
177 if (ncommand==lcommand){
178 lcommand=10+2*lcommand;
179 std::string* newcommands=
new std::string[lcommand];
181 for(i=0;i<ncommand;i++){
182 newcommands[i]=commands[i];
185 commands=newcommands;
187 commands[ncommand]=cmd;
195 static int iniflag=1;
227 EvtId evtnumstable[100],evtnumparton[100];
228 int stableindex[100],partonindex[100];
237 static double px[4000],py[4000],pz[4000],e[4000];
238 if (iniflag==1)
lundar_(&isr,&iniflag,&xmp,&ndaugjs,kf,km,px,py,pz,e);
250 lundar_(&isr,&iniflag,&xmp,&ndaugjs,kf,km,px,py,pz,e);
253 if(mtg)std::cout<<
"checktag_.decaytag="<<
checktag_.decaytag<<std::endl;
267 for(i=0;i<ndaugjs;i++){
268 if (
abs(kf[i])==11 || kf[i]==92 || kf[i]==22)
continue;
271 report(
ERROR,
"EvtGen") <<
"Lunda returned particle:"<<kf[i]<<endl;
272 report(
ERROR,
"EvtGen") <<
"This can not be translated to evt number"<<endl;
273 report(
ERROR,
"EvtGen") <<
"and the decay will be rejected!"<<endl;
274 report(
ERROR,
"EvtGen") <<
"The decay was of particle:"<<ip<<endl;
276 std::cout<<
"xmp= "<<xmp<<std::endl;
281 if (
abs(kf[i])<=6||kf[i]==21){
282 partonindex[numparton]=i;
287 stableindex[numstable]=i;
297 if (px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i]>=e[i]*e[i]){
299 e[i]=sqrt(px[i]*px[i]+py[i]*py[i]+pz[i]*pz[i])+0.0000000001;
303 p4[i].
set(e[i],px[i],py[i],pz[i]);
314 if(
abs(esum-xmp)>0.001 ){
315 std::cout<<
"Unexpected final states from Lunda with total energy "<<esum<<
" unequal to "<<xmp<<std::endl;
317 std::cout<<
"xmp= "<<xmp<<std::endl;
322 }
while( more && (
count<10000) && mtg );
326 report(
INFO,
"EvtGen") <<
"Too many loops in EvtLunda!!!"<<endl;
328 for(i=0;i<numstable;i++){
340 for(i=0;i<numstable;i++){
341 p->
getDaug(i)->
init(evtnumstable[i],p4[stableindex[i]]);
344 if ( ndaugFound == 0 ) {
345 report(
ERROR,
"EvtGen") <<
"Lunda has failed to do a decay ";
363 for(i=0;i<numparton;i++){
364 p4string+=p4[partonindex[i]];
369 for(i=0;i<numstable;i++){
370 if (km[stableindex[i]]==0){
371 type[nprimary++]=evtnumstable[i];
381 for(i=0;i<numparton;i++){
382 p4partons[i]=p4[partonindex[i]];
391 for(i=0;i<numstable;i++){
393 if (km[stableindex[i]]==0){
394 p->
getDaug(nprimary++)->
init(evtnumstable[i],p4[stableindex[i]]);
400 for(i=0;i<numstable;i++){
401 if (km[stableindex[i]]!=0){
402 type[nsecond++]=evtnumstable[i];
409 EvtVector4R p4stringboost(p4string.get(0),-p4string.get(1),
410 -p4string.get(2),-p4string.get(3));
413 for(i=0;i<numstable;i++){
414 if (km[stableindex[i]]!=0){
415 p4[stableindex[i]]=
boostTo(p4[stableindex[i]],p4stringboost);
423 if ( nsecond == 0 ) {
424 report(
ERROR,
"EvtGen") <<
"Jetset has failed to do a decay ";
447 for(i=0;i<ndaug;i++){
468 double beta=acos(p4Psi.
get(3)/p4Psi.
d3mag());
481 if (nlundadecays==ntable){
485 for(i=0;i<ntable;i++){
486 newlundadecays[i]=lundadecays[i];
489 delete [] lundadecays;
490 lundadecays=newlundadecays;
493 lundadecays[nlundadecays++]=jsdecay;
506 for(
int i=0;i<ncommand;i++)
522 for(
int i=0;i< _ndaugs;i++){
528 if(totMass>p->
mass())
goto checkA;
void lundar_(int *, int *, float *, int *, int *, int *, double *, double *, double *, double *)
void lugive0_(const char *cnfgstr, int length)
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
ostream & report(Severity severity, const char *facility)
void lundar_(int *, int *, float *, int *, int *, int *, double *, double *, double *, double *)
void lugive0_(const char *cnfgstr, int length)
DOUBLE_PRECISION count[3]
void setDaughterSpinDensity(int daughter)
static int inChannelList(EvtId parent, int ndaug, EvtId *daugs)
void LundaInit(int dummy)
std::string commandName()
void decay(EvtParticle *p)
void ExclusiveDecay(EvtParticle *p)
void command(std::string cmd)
void getName(std::string &name)
static int getStdHep(EvtId id)
static EvtId evtIdFromStdHep(int stdhep)
static std::string name(EvtId i)
static double getMass(EvtId i)
static EvtId getId(const std::string &name)
void setSpinDensityForwardHelicityBasis(const EvtSpinDensity &rho)
void makeDaughters(int ndaug, EvtId *id)
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
void setDiagonalSpinDensity()
const EvtVector4R & getP4() const
EvtParticle * getDaug(int i)
void deleteDaughters(bool keepChannel=false)
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
void Set(int i, int j, const EvtComplex &rhoij)
void set(int i, double d)