EvtDDalitz Class Reference

#include <EvtDDalitz.hh>

Inheritance diagram for EvtDDalitz:

List of all members.

Public Member Functions
void	checkNArg (int a1, int a2=-1, int a3=-1, int a4=-1)
void	checkNDaug (int d1, int d2=-1)
void	checkQ ()
void	checkSpinDaughter (int d1, EvtSpinType::spintype sp)
void	checkSpinParent (EvtSpinType::spintype sp)
EvtDecayBase *	clone ()
virtual void	command (std::string cmd)
virtual std::string	commandName ()
void	decay (EvtParticle *p)
void	disableCheckQ ()
	EvtDDalitz ()
double	getArg (int j)
double *	getArgs ()
std::string *	getArgsStr ()
std::string	getArgStr (int j)
double	getBranchingFraction ()
EvtId	getDaug (int i)
EvtId *	getDaugs ()
int	getDSum ()
std::string	getModelName ()
void	getName (std::string &name)
int	getNArg ()
int	getNDaug ()
EvtId	getParentId ()
int	getPHOTOS ()
double	getProbMax (double prob)
void	init ()
void	initProbMax ()
void	makeDecay (EvtParticle *p)
virtual bool	matchingDecay (const EvtDecayBase &other) const
void	noProbMax ()
virtual int	nRealDaughters ()
void	printSummary ()
double	resetProbMax (double prob)
void	saveDecayInfo (EvtId ipar, int ndaug, EvtId *daug, int narg, std::vector< std::string > &args, std::string name, double brfr)
void	setPHOTOS ()
void	setProbMax (double prbmx)
void	setSummary ()
void	setVerbose ()
void	setWeight (double weight)
int	summary ()
int	verbose ()
void	vertex (int *i1, const EvtComplex &amp)
void	vertex (int i1, int i2, int i3, const EvtComplex &amp)
void	vertex (int i1, int i2, const EvtComplex &amp)
void	vertex (int i1, const EvtComplex &amp)
void	vertex (const EvtComplex &amp)
virtual	~EvtDDalitz ()
Static Public Member Functions
void	findMass (EvtParticle *p)
void	findMasses (EvtParticle *p, int ndaugs, EvtId daugs[10], double masses[10])
double	findMaxMass (EvtParticle *p)
Protected Member Functions
bool	daugsDecayedByParentModel ()
Protected Attributes
EvtAmp	_amp2
bool	_daugsDecayedByParentModel
Private Member Functions
EvtComplex	amplDtoK0KK (EvtVector4R p4_p, EvtVector4R moms1, EvtVector4R moms2, EvtVector4R moms3)
EvtComplex	amplDtoK0PiPi (EvtVector4R p4_p, EvtVector4R moms1, EvtVector4R moms2, EvtVector4R moms3)
Private Attributes
int	_d1
int	_d2
int	_d3
int	_flag
vector< EvtFlatteParam >	_kkpi_params

---

Constructor & Destructor Documentation

EvtDDalitz::EvtDDalitz (   )  [inline]

{}

EvtDDalitz::~EvtDDalitz (   )  [virtual]

{}

---

Member Function Documentation

EvtComplex EvtDDalitz::amplDtoK0KK (  EvtVector4R  p4_p, EvtVector4R  moms1, EvtVector4R  moms2, EvtVector4R  moms3 )  [private]

{ //phi EvtResonance DK0KKRes1( p4_p, moms2, moms3, 113.75, -40.0, 0.0043, 1.019456, 1 ) ; //a0(980) EvtResonance DK0KKRes2( p4_p, moms2, moms3, 152.25, 69.0, 0.1196 , 0.9847, 0 ) ; //f0(980) EvtResonance DK0KKRes3( p4_p, moms2, moms3, 30.5, -201.0, 0.05, 0.980 , 0 ) ; //a0(980)+ EvtResonance DK0KKRes4( p4_p, moms1, moms2, 85.75, -93.0, 0.1196 , 0.9847, 0 ) ; //a0(980)- EvtResonance DK0KKRes5( p4_p, moms3, moms1, 8. , -53.0 ,0.1196, 0.9847, 0 ) ; return EvtComplex(1.0,0.0) + DK0KKRes1.resAmpl() + DK0KKRes2.resAmpl() + DK0KKRes3.resAmpl() + DK0KKRes4.resAmpl() + DK0KKRes5.resAmpl() ; }

EvtComplex EvtDDalitz::amplDtoK0PiPi (  EvtVector4R  p4_p, EvtVector4R  moms1, EvtVector4R  moms2, EvtVector4R  moms3 )  [private]

{ //K*(892)- EvtResonance2 DK2piRes1(p4_p,moms1,moms2,1.418,-190.0,0.0508,0.89166,1); //K0*(1430) EvtResonance2 DK2piRes2(p4_p,moms1,moms2,1.818,-337.0,0.294 ,1.412 ,0); //K2*(1430) EvtResonance2 DK2piRes3(p4_p,moms1,moms2,0.909, -5.0,0.0985,1.4256 ,2); //K*(1680) EvtResonance2 DK2piRes4(p4_p,moms1,moms2,5.091,-166.0,0.322 ,1.717 ,1); //DCS K*(892) EvtResonance2 DK2piRes5(p4_p,moms1,moms3,0.100, -19.0,0.0508,0.89166,1); //Rho EvtResonance2 DK2piRes6(p4_p,moms3,moms2,0.909,-340.0,0.1502,0.7693,1); //Omega EvtResonance2 DK2piRes7(p4_p,moms3,moms2,.0336,-226.0,0.00844,0.78257,1); //f0(980) EvtResonance2 DK2piRes8(p4_p,moms3,moms2,0.309,-152.0,0.05,0.977,0); //f0(1370) EvtResonance2 DK2piRes9(p4_p,moms3,moms2,1.636,-255.0,0.272,1.31,0); //f2(1270) EvtResonance2 DK2piRes10(p4_p,moms3,moms2,0.636,-32.0,0.1851,1.2754,2); return EvtComplex(1.0,0.0) + DK2piRes1.resAmpl() + DK2piRes2.resAmpl() + DK2piRes3.resAmpl() + DK2piRes4.resAmpl() + DK2piRes5.resAmpl() + DK2piRes6.resAmpl() + DK2piRes7.resAmpl() + DK2piRes8.resAmpl() + DK2piRes9.resAmpl() + DK2piRes10.resAmpl(); }

void EvtDecayBase::checkNArg (  int  a1, int  a2 = -1, int  a3 = -1, int  a4 = -1 )  [inherited]

{ if ( _narg != a1 && _narg != a2 && _narg != a3 && _narg != a4 ) { report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected "<<endl; report(ERROR,"EvtGen") << a1<<endl;; if ( a2>-1) { report(ERROR,"EvtGen") << " or " << a2<<endl; } if ( a3>-1) { report(ERROR,"EvtGen") << " or " << a3<<endl; } if ( a4>-1) { report(ERROR,"EvtGen") << " or " << a4<<endl; } report(ERROR,"EvtGen") << " arguments but found:"<< _narg << endl; printSummary(); report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } }

void EvtDecayBase::checkNDaug (  int  d1, int  d2 = -1 )  [inherited]

{ if ( _ndaug != d1 && _ndaug != d2 ) { report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected "; report(ERROR,"EvtGen") << d1; if ( d2>-1) { report(ERROR,"EvtGen") << " or " << d2; } report(ERROR,"EvtGen") << " daughters but found:"<< _ndaug << endl; printSummary(); report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } }

void EvtDecayBase::checkQ (   )  [inherited]

{ int i; int q=0; int qpar; //If there are no daughters (jetset etc) then we do not //want to do this test. Why? Because sometimes the parent //will have a nonzero charge. if ( _ndaug != 0) { for(i=0; i<_ndaug; i++ ) { q += EvtPDL::chg3(_daug[i]); } qpar = EvtPDL::chg3(_parent); if ( q != qpar ) { report(ERROR,"EvtGen") <<_modelname.c_str()<< " generator expected " << " charge to be conserved, found:"<<endl; report(ERROR,"EvtGen") << "Parent charge of "<<(qpar/3)<<endl; report(ERROR,"EvtGen") << "Sum of daughter charge of "<<(q/3)<<endl; report(ERROR,"EvtGen") << "The parent is "<< EvtPDL::name(_parent).c_str()<<endl; for(i=0; i<_ndaug; i++ ) { report(ERROR,"EvtGen") << "Daughter "<< EvtPDL::name(_daug[i]).c_str()<<endl; } report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } } }

void EvtDecayBase::checkSpinDaughter (  int  d1, EvtSpinType::spintype  sp )  [inherited]

{ EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getDaug(d1)); if ( parenttype != sp ) { report(ERROR,"EvtGen") << _modelname.c_str() << " did not get the correct daughter spin d=" << d1 << endl; printSummary(); report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } }

void EvtDecayBase::checkSpinParent (  EvtSpinType::spintype  sp  )  [inherited]

{ EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId()); if ( parenttype != sp ) { report(ERROR,"EvtGen") << _modelname.c_str() << " did not get the correct parent spin\n"; printSummary(); report(ERROR,"EvtGen") << "Will terminate execution!"<<endl; ::abort(); } }

EvtDecayBase * EvtDDalitz::clone (   )  [virtual]

Implements EvtDecayBase. { return new EvtDDalitz; }

void EvtDecayBase::command (  std::string  cmd  )  [virtual, inherited]

Reimplemented in EvtJetSet, EvtLunda, EvtLundCharm, EvtPythia, and EvtTauola. { report(ERROR,"EvtGen") << "Should never call EvtDecayBase::command"<<endl; ::abort(); }

std::string EvtDecayBase::commandName (   )  [virtual, inherited]

Reimplemented in EvtJetSet, EvtLunda, EvtLundCharm, EvtPythia, and EvtTauola. { return std::string(""); }

bool EvtDecayBase::daugsDecayedByParentModel (   )  [inline, protected, inherited]

{return _daugsDecayedByParentModel;}

void EvtDDalitz::decay (  EvtParticle *  p  )  [virtual]

Implements EvtDecayBase. { //added by Lange Jan4,2000 static EvtId BP = EvtPDL::getId("B+"); static EvtId BM = EvtPDL::getId("B-"); static EvtId B0 = EvtPDL::getId("B0"); static EvtId B0B = EvtPDL::getId("anti-B0"); static EvtId DM=EvtPDL::getId("D-"); static EvtId DP=EvtPDL::getId("D+"); static EvtId D0=EvtPDL::getId("D0"); static EvtId D0B=EvtPDL::getId("anti-D0"); static EvtId KM=EvtPDL::getId("K-"); static EvtId KP=EvtPDL::getId("K+"); static EvtId K0=EvtPDL::getId("K0"); static EvtId KB=EvtPDL::getId("anti-K0"); static EvtId PIM=EvtPDL::getId("pi-"); static EvtId PIP=EvtPDL::getId("pi+"); static EvtId PI0=EvtPDL::getId("pi0"); double oneby2 = 0.707106782; bool isBToDK=false; if ( p -> getParent () ) { std::string name1; EvtDecayTable::getDecayFunc(p->getParent())->getName(name1); EvtId parId = p -> getParent()->getId (); if ( ( BP == parId ) || ( BM == parId ) || ( B0 == parId ) || ( B0B == parId ) ) if (name1 == "BTODDALITZCPK") isBToDK=true; } //same structure for all of these decays p->initializePhaseSpace(getNDaug(),getDaugs()); EvtVector4R moms1 = p->getDaug(_d1)->getP4(); EvtVector4R moms2 = p->getDaug(_d2)->getP4(); EvtVector4R moms3 = p->getDaug(_d3)->getP4(); EvtVector4R p4_p; p4_p.set(p->mass(),0.0,0.0,0.0); EvtComplex amp(1.0,0.0); //now determine which D and which decay //data from Anjos et al, Phys.Rev.D 1993, v.48,num.1,p.56 (E691 resuls) //for D+ -> K- pi+ pi+, and from Adler et al, Phys.Lett. B196 (1987), 107 //(Mark III results) for D+ -> K0bar pi+ pi0. //CLEO results for D0->k-pi+pi0 if ( _flag==1) { //have a D+ -> K- pi+ pi+ decay, or charge conjugate //Anjos etal e691 - Phys Rev D48, 56 (1993) EvtResonance DplusRes11(p4_p,moms1,moms2,0.78,-60.0,0.0498,0.89610,1); EvtResonance DplusRes12(p4_p,moms3,moms1,0.78,-60.0,0.0498,0.89610,1);//K*(892) EvtResonance DplusRes21(p4_p,moms1,moms2,0.53,132.0,0.287,1.429,0); EvtResonance DplusRes22(p4_p,moms3,moms1,0.53,132.0,0.287,1.429,0);//K*(1430) EvtResonance DplusRes31(p4_p,moms1,moms2,0.47,-51.0,0.323,1.714,1); EvtResonance DplusRes32(p4_p,moms3,moms1,0.47,-51.0,0.323,1.714,1);//K*(1680) amp = amp + oneby2*(-DplusRes11.resAmpl()+DplusRes12.resAmpl()) + oneby2*(DplusRes21.resAmpl() + DplusRes22.resAmpl()) + oneby2*(-DplusRes31.resAmpl()+ DplusRes32.resAmpl()); } if ( _flag==2) { //have a D+ -> K0bar pi+ pi0 decay (daughters must be listed in this order!) //adler etal MarkIII - Phys Lett B196, 107 (1987) // Results in this paper: // Kbar rho+ FitFraction = 68+/-8+/-12 Phase 0 // Kbar* pi+ 19+/-6+/-6 43+/-23 // nonres 13+/-7+/-8 250+/-19 // These numbers below seem not to be exactly the same // the phases are equiv to -106=254 and 41 // EvtResonance DplusKpipi0Res1(p4_p,moms2,moms3,1.00,0.00,0.1512,0.7699,1); //rho+ EvtResonance DplusKpipi0Res2(p4_p,moms3,moms1,0.8695,0.7191,0.0498,0.89159,1); //K*0 amp = 0.9522*EvtComplex(cos(-1.8565),sin(-1.8565)) + 1.00*DplusKpipi0Res1.relBrWig(0) + 0.8695*EvtComplex(cos(0.7191),sin(0.7191))*DplusKpipi0Res2.relBrWig(1); } if(_flag==3) { //know it's a D0 -> K0_bar pi+ pi- or charge conjugate // EvtResonance DK2piRes1(p4_p,moms1,moms2,2.31,109.0,0.0498,0.89610,1);//K*(892) // EvtResonance DK2piRes2(p4_p,moms3,moms2,1.59,-123.0,0.1491,0.7683,1);//RHO(770) // amp = amp + DK2piRes1.resAmpl() + DK2piRes2.resAmpl(); // D0 -> K0 pi+ pi- + CC // If it does not come from a B->DK, decay it as D0 or D0bar separatly // if p4_p is D0, moms1 is K0, moms2 is pi-, moms3 is pi+ // if p4_p is D0bar, moms1 is K0, moms2 is pi+, moms3 is pi- if ( isBToDK ) { // Gamma angle in rad. double gamma = EvtDecayTable::getDecayFunc( p->getParent() ) -> getArg( 0 ) ; // Strong phase in rad. double delta = EvtDecayTable::getDecayFunc( p->getParent() ) -> getArg( 1 ) ; // Ratio between B->D0K and B->D0barK double A = EvtDecayTable::getDecayFunc( p->getParent() ) -> getArg( 2 ) ; EvtComplex Factor( fabs( A ) * cos ( delta ) , fabs( A ) * sin ( delta ) ) ; if ( ( p->getParent()->getId() == BP ) || ( p->getParent()->getId() == B0 ) ) { // the ratio D/Dbar Factor = Factor * EvtComplex( cos ( gamma ) , sin ( gamma ) ) ; if ( p->getId() == D0 ) { // the flavor of the particle has no meaning. But we need // it to know which daughter is pi+ or pi- // M( B+ or B0 ) = f(Dbar) + factor * f(D) // f(Dbar) = amplDtoK0PiPi(pD, K0, pi+, pi-) // f(D) = amplDtoK0PiPi(pD, K0, pi-, pi+) // Then ... amp = amplDtoK0PiPi( p4_p , moms1 , moms3 , moms2 ) + Factor * amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) ; } else { amp = amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) + Factor * amplDtoK0PiPi( p4_p , moms1 , moms3 , moms2 ) ; } } else if ( ( p->getParent() -> getId() == BM ) || ( p->getParent() -> getId() == B0B ) ) { Factor = Factor * EvtComplex( cos ( gamma ) , - sin ( gamma ) ) ; // here M( B- or B0bar ) = f(D) + factor * f(Dbar) then ... if ( p->getId() == D0 ) { amp = amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) + Factor * amplDtoK0PiPi( p4_p , moms1 , moms3 , moms2 ) ; } else { amp = amplDtoK0PiPi( p4_p , moms1 , moms3 , moms2 ) + Factor * amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) ; } } } else { amp = amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) ; } } if(_flag==4) { EvtResonance2 DKpipi0Res1(p4_p,moms2,moms3,1.0 ,0.0 ,0.1507,0.770 ,1); //rho EvtResonance2 DKpipi0Res2(p4_p,moms1,moms2,0.39, -0.2 ,0.0505,0.8961,1); //k*0 EvtResonance2 DKpipi0Res3(p4_p,moms1,moms3,0.44, 163.0 ,0.050 ,0.8915,1); //k*- EvtResonance2 DKpipi0Res4(p4_p,moms1,moms3,0.77 ,55.5 ,0.294 ,1.412 ,0); //k01430- EvtResonance2 DKpipi0Res5(p4_p,moms1,moms2,0.85 ,166.0 ,0.294 ,1.412 ,0); //k01430bar EvtResonance2 DKpipi0Res6(p4_p,moms2,moms3,2.5 ,171.0 ,0.240 ,1.700 ,1); //rho1700 EvtResonance2 DKpipi0Res7(p4_p,moms1,moms3,2.5 ,103.0 ,0.322 ,1.717 ,1); //K*1680- double pi180inv = 1.0/EvtConst::radToDegrees; amp = EvtComplex(1.75*cos(31.2*pi180inv),1.75*sin(31.2*pi180inv)) + DKpipi0Res1.resAmpl() + DKpipi0Res2.resAmpl() + DKpipi0Res3.resAmpl() + DKpipi0Res4.resAmpl() + DKpipi0Res5.resAmpl() + DKpipi0Res6.resAmpl() + DKpipi0Res7.resAmpl(); } if(_flag==5) { // D0 -> K0 K+ K- + CC // If it does not come from a B->DK, decay it as D0 or D0bar separatly // if p4_p is D0, moms1 is K0, moms2 is pi-, moms3 is pi+ // if p4_p is D0bar, moms1 is K0, moms2 is pi+, moms3 is pi- if ( isBToDK ){ // Gamma angle in rad. double gamma = EvtDecayTable::getDecayFunc( p->getParent() ) -> getArg( 0 ) ; // Strong phase in rad. double delta = EvtDecayTable::getDecayFunc( p->getParent() ) -> getArg( 1 ) ; // Ratio between B->D0K and B->D0barK double A = EvtDecayTable::getDecayFunc( p->getParent() ) -> getArg( 2 ) ; EvtComplex Factor( fabs( A ) * cos ( delta ) , fabs( A ) * sin ( delta ) ) ; if ( ( p->getParent()->getId() == BP ) || ( p->getParent()->getId() == B0 ) ) { // the ratio D/Dbar Factor = Factor * EvtComplex( cos ( gamma ) , sin ( gamma ) ) ; if ( p->getId() == D0 ) { // the flavor of the particle has no meaning. But we need // it to know which daughter is pi+ or pi- // M( B+ or B0 ) = f(Dbar) + factor * f(D) // f(Dbar) = amplDtoK0PiPi(pD, K0, K+, K-) // f(D) = amplDtoK0PiPi(pD, K0, K-, K+) // Then ... amp = amplDtoK0KK( p4_p , moms1 , moms3 , moms2 ) + Factor * amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) ; } else { amp = amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) + Factor * amplDtoK0KK( p4_p , moms1 , moms3 , moms2 ) ; } } else if ( ( p->getParent() -> getId() == BM ) || ( p->getParent() -> getId() == B0B ) ) { Factor = Factor * EvtComplex( cos ( gamma ) , - sin ( gamma ) ) ; // here M( B- or B0bar ) = f(D) + factor * f(Dbar) then ... if ( p->getId() == D0 ) { amp = amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) + Factor * amplDtoK0KK( p4_p , moms1 , moms3 , moms2 ) ; } else { amp = amplDtoK0KK( p4_p , moms1 , moms3 , moms2 ) + Factor * amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) ; } } } else { amp = amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) ; } } // Ds -> K K pi if(_flag==6) { EvtResonance2 DsKKpiRes1(p4_p, moms3, moms1, 1.0, 0.0, 0.0455, 0.8944, 1, true); // K*(892) EvtResonance2 DsKKpiRes2(p4_p, moms3, moms1, 1.48, 138., 0.290, 1.414, 0); // K*_0(1430) EvtFlatte DsKKpiRes3(p4_p, moms1, moms2, 5.07, 156., 0.965, _kkpi_params); // f_0(980) EvtResonance2 DsKKpiRes4(p4_p, moms1, moms2, 1.15, -10., 0.00426, 1.019455, 1, true); // phi(1020) EvtResonance2 DsKKpiRes5(p4_p, moms1, moms2, 1.28, 53., 0.265, 1.350, 0); // f_0(1370) EvtResonance2 DsKKpiRes6(p4_p, moms1, moms2, 1.19, 87., 0.137, 1.724, 0); // f_0(1710) amp = DsKKpiRes1.resAmpl() + DsKKpiRes2.resAmpl() + DsKKpiRes3.resAmpl() + DsKKpiRes4.resAmpl() + DsKKpiRes5.resAmpl() + DsKKpiRes6.resAmpl(); } vertex(amp); return ; }

void EvtDecayBase::disableCheckQ (   )  [inline, inherited]

{_chkCharge=0;};

void EvtDecayBase::findMass (  EvtParticle *  p  )  [static, inherited]

{ //Need to also check that this mass does not screw //up the parent //This code assumes that for the ith daughter, 0..i-1 //already have a mass double maxOkMass=findMaxMass(p); int count=0; double mass; bool massOk=false; int i; while (!massOk) { count++; if ( count > 10000 ) { report(INFO,"EvtGen") << "Can not find a valid mass for: " << EvtPDL::name(p->getId()).c_str() <<endl; report(INFO,"EvtGen") << "Now printing parent and/or grandparent tree\n"; if ( p->getParent() ) { if ( p->getParent()->getParent() ) { p->getParent()->getParent()->printTree(); report(INFO,"EvtGen") << p->getParent()->getParent()->mass() <<endl; report(INFO,"EvtGen") << p->getParent()->mass() <<endl; } else{ p->getParent()->printTree(); report(INFO,"EvtGen") << p->getParent()->mass() <<endl; } } else p->printTree(); report(INFO,"EvtGen") << "maxokmass=" << maxOkMass << " " << EvtPDL::getMinMass(p->getId()) << " " << EvtPDL::getMaxMass(p->getId())<<endl; if ( p->getNDaug() ) { for (i=0; i<p->getNDaug(); i++) { report(INFO,"EvtGen") << p->getDaug(i)->mass()<<" "; } report(INFO,"EvtGen") << endl; } if ( maxOkMass >= EvtPDL::getMinMass(p->getId()) ) { report(INFO,"EvtGen") << "taking a default value\n"; p->setMass(maxOkMass); return; } assert(0); } mass = EvtPDL::getMass(p->getId()); //Just need to check that this mass is > than //the mass of all daughters double massSum=0.; if ( p->getNDaug() ) { for (i=0; i<p->getNDaug(); i++) { massSum+= p->getDaug(i)->mass(); } } //some special cases are handled with 0 (stable) or 1 (k0->ks/kl) daughters if (p->getNDaug()<2) massOk=true; if ( p->getParent() ) { if ( p->getParent()->getNDaug()==1 ) massOk=true; } if ( !massOk ) { if (massSum < mass) massOk=true; if ( mass> maxOkMass) massOk=false; } } p->setMass(mass); }

void EvtDecayBase::findMasses (  EvtParticle *  p, int  ndaugs, EvtId  daugs[10], double  masses[10] )  [static, inherited]

{ int i; double mass_sum; int count=0; if (!( p->firstornot() )) { for (i = 0; i < ndaugs; i++ ) { masses[i] = p->getDaug(i)->mass(); } //for } //if else { p->setFirstOrNot(); // if only one daughter do it if (ndaugs==1) { masses[0]=p->mass(); return; } //until we get a combo whose masses are less than _parent mass. do { mass_sum = 0.0; for (i = 0; i < ndaugs; i++ ) { masses[i] = EvtPDL::getMass(daugs[i]); mass_sum = mass_sum + masses[i]; } count++; if(count==10000) { report(ERROR,"EvtGen") <<"Decaying particle:"<< EvtPDL::name(p->getId()).c_str()<<" (m="<<p->mass()<<")"<<endl; report(ERROR,"EvtGen") <<"To the following daugthers"<<endl; for (i = 0; i < ndaugs; i++ ) { report(ERROR,"EvtGen") << EvtPDL::name(daugs[i]).c_str() << endl; } report(ERROR,"EvtGen") << "Has been rejected "<<count << " times, will now take minimal masses " << " of daugthers"<<endl; mass_sum=0.; for (i = 0; i < ndaugs; i++ ) { masses[i] = EvtPDL::getMinMass(daugs[i]); mass_sum = mass_sum + masses[i]; } if (mass_sum > p->mass()){ report(ERROR,"EvtGen") << "Parent mass="<<p->mass() << "to light for daugthers."<<endl << "Will throw the event away."<<endl; //dont terminate - start over on the event. EvtStatus::setRejectFlag(); mass_sum=0.; // ::abort(); } } } while ( mass_sum > p->mass()); } //else return; }

double EvtDecayBase::findMaxMass (  EvtParticle *  p  )  [static, inherited]

{ double maxOkMass=EvtPDL::getMaxMass(p->getId()); //protect against vphotons if ( maxOkMass < 0.0000000001 ) return 10000000.; //and against already determined masses if ( p->hasValidP4() ) maxOkMass=p->mass(); EvtParticle *par=p->getParent(); if ( par ) { double maxParMass=findMaxMass(par); int i; double minDaugMass=0.; for(i=0;i<par->getNDaug();i++){ EvtParticle *dau=par->getDaug(i); if ( dau!=p) { // it might already have a mass if ( dau->isInitialized() || dau->hasValidP4() ) minDaugMass+=dau->mass(); else //give it a bit of phase space minDaugMass+=1.000001*EvtPDL::getMinMass(dau->getId()); } } if ( maxOkMass>(maxParMass-minDaugMass)) maxOkMass=maxParMass-minDaugMass; } return maxOkMass; }

double EvtDecayBase::getArg (  int  j  )  [inherited]

{ // Verify string const char* str = _args[j].c_str(); int i = 0; while(str[i]!=0){ if (isalpha(str[i]) && str[i]!='e') { report(INFO,"EvtGen") << "String " << str << " is not a number" << endl; assert(0); } i++; } char** tc=0; return strtod(_args[j].c_str(),tc); }

double * EvtDecayBase::getArgs (   )  [inherited]

{ if ( _argsD ) return _argsD; //The user has asked for a list of doubles - the arguments //better all be doubles... if ( _narg==0 ) return _argsD; _argsD = new double[_narg]; int i; char * tc; for(i=0;i<_narg;i++) { _argsD[i] = strtod(_args[i].c_str(),&tc); } return _argsD; }

std::string* EvtDecayBase::getArgsStr (   )  [inline, inherited]

{return _args;}

std::string EvtDecayBase::getArgStr (  int  j  )  [inline, inherited]

{return _args[j];}

double EvtDecayBase::getBranchingFraction (   )  [inline, inherited]

{return _brfr;}

EvtId EvtDecayBase::getDaug (  int  i  )  [inline, inherited]

{return _daug[i];}

EvtId* EvtDecayBase::getDaugs (   )  [inline, inherited]

{return _daug;}

int EvtDecayBase::getDSum (   )  [inline, inherited]

{return _dsum; }

std::string EvtDecayBase::getModelName (   )  [inline, inherited]

{return _modelname; }

void EvtDDalitz::getName (  std::string &  name  )  [virtual]

Implements EvtDecayBase. { model_name="D_DALITZ"; }

int EvtDecayBase::getNArg (   )  [inline, inherited]

{return _narg;}

int EvtDecayBase::getNDaug (   )  [inline, inherited]

{return _ndaug;}

EvtId EvtDecayBase::getParentId (   )  [inline, inherited]

{return _parent;}

int EvtDecayBase::getPHOTOS (   )  [inline, inherited]

{return _photos;}

double EvtDecayBase::getProbMax (  double  prob  )  [inherited]

{ int i; //diagnostics sum_prob+=prob; if (prob>max_prob) max_prob=prob; if ( defaultprobmax && ntimes_prob<=500 ) { //We are building up probmax with this iteration ntimes_prob += 1; if ( prob > probmax ) { probmax = prob;} if (ntimes_prob==500) { probmax*=1.2; } return 1000000.0*prob; } if ( prob> probmax*1.0001) { report(INFO,"EvtGen") << "prob > probmax:("<<prob<<">"<<probmax<<")"; report(INFO,"") << "("<<_modelname.c_str()<<") "; report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> "; for(i=0;i<_ndaug;i++){ report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " "; } report(INFO,"") << endl; if (defaultprobmax) probmax = prob; } ntimes_prob += 1; return probmax; } //getProbMax

void EvtDDalitz::init (   )  [virtual]

Reimplemented from EvtDecayBase. { // check that there are 0 arguments static EvtId DM=EvtPDL::getId("D-"); static EvtId DP=EvtPDL::getId("D+"); static EvtId D0=EvtPDL::getId("D0"); static EvtId D0B=EvtPDL::getId("anti-D0"); static EvtId DSP=EvtPDL::getId("D_s+"); static EvtId DSM=EvtPDL::getId("D_s-"); static EvtId KM=EvtPDL::getId("K-"); static EvtId KP=EvtPDL::getId("K+"); static EvtId K0=EvtPDL::getId("K0"); static EvtId KB=EvtPDL::getId("anti-K0"); static EvtId KL=EvtPDL::getId("K_L0"); static EvtId KS=EvtPDL::getId("K_S0"); static EvtId PIM=EvtPDL::getId("pi-"); static EvtId PIP=EvtPDL::getId("pi+"); static EvtId PI0=EvtPDL::getId("pi0"); static double MPI = EvtPDL::getMeanMass(PI0); static double MKP = EvtPDL::getMeanMass(KP); checkNArg(0); checkNDaug(3); checkSpinParent(EvtSpinType::SCALAR); checkSpinDaughter(0,EvtSpinType::SCALAR); checkSpinDaughter(1,EvtSpinType::SCALAR); checkSpinDaughter(2,EvtSpinType::SCALAR); EvtId parnum=getParentId(); EvtId d1=getDaug(0); EvtId d2=getDaug(1); EvtId d3=getDaug(2); _flag=0; if ( parnum == D0 ) { //look for either a K- pi+ pi0 or K0bar pi+ pi- if ( d1==KM && d2==PIP && d3==PI0 ) { _flag=4; _d1=0; _d2=1; _d3=2;} if ( d1==KM && d3==PIP && d2==PI0 ) { _flag=4; _d1=0; _d2=2; _d3=1;} if ( d2==KM && d1==PIP && d3==PI0 ) { _flag=4; _d1=1; _d2=0; _d3=2;} if ( d2==KM && d3==PIP && d1==PI0 ) { _flag=4; _d1=1; _d2=2; _d3=0;} if ( d3==KM && d1==PIP && d2==PI0 ) { _flag=4; _d1=2; _d2=0; _d3=1;} if ( d3==KM && d2==PIP && d1==PI0 ) { _flag=4; _d1=2; _d2=1; _d3=0;} if ( d1==KB && d2==PIP && d3==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;} if ( d1==KB && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;} if ( d2==KB && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;} if ( d2==KB && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;} if ( d3==KB && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;} if ( d3==KB && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;} if ( d1==KL && d2==PIP && d3==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;} if ( d1==KL && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;} if ( d2==KL && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;} if ( d2==KL && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;} if ( d3==KL && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;} if ( d3==KL && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;} if ( d1==KS && d2==PIP && d3==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;} if ( d1==KS && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;} if ( d2==KS && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;} if ( d2==KS && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;} if ( d3==KS && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;} if ( d3==KS && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;} if ( d1==KS && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;} if ( d1==KS && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;} if ( d2==KS && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;} if ( d2==KS && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;} if ( d3==KS && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;} if ( d3==KS && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;} if ( d1==KL && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;} if ( d1==KL && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;} if ( d2==KL && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;} if ( d2==KL && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;} if ( d3==KL && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;} if ( d3==KL && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;} if ( d1==K0 && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;} if ( d1==K0 && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;} if ( d2==K0 && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;} if ( d2==K0 && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;} if ( d3==K0 && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;} if ( d3==K0 && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;} } if ( parnum == D0B ) { //look for either a K+ pi- pi0 or K0 pi+ pi- if ( d1==KP && d2==PIM && d3==PI0 ) { _flag=4; _d1=0; _d2=1; _d3=2;} if ( d1==KP && d3==PIM && d2==PI0 ) { _flag=4; _d1=0; _d2=2; _d3=1;} if ( d2==KP && d1==PIM && d3==PI0 ) { _flag=4; _d1=1; _d2=0; _d3=2;} if ( d2==KP && d3==PIM && d1==PI0 ) { _flag=4; _d1=1; _d2=2; _d3=0;} if ( d3==KP && d1==PIM && d2==PI0 ) { _flag=4; _d1=2; _d2=0; _d3=1;} if ( d3==KP && d2==PIM && d1==PI0 ) { _flag=4; _d1=2; _d2=1; _d3=0;} if ( d1==K0 && d2==PIP && d3==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;} if ( d1==K0 && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;} if ( d2==K0 && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;} if ( d2==K0 && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;} if ( d3==K0 && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;} if ( d3==K0 && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;} if ( d1==KL && d2==PIP && d3==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;} if ( d1==KL && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;} if ( d2==KL && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;} if ( d2==KL && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;} if ( d3==KL && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;} if ( d3==KL && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;} if ( d1==KS && d2==PIP && d3==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;} if ( d1==KS && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;} if ( d2==KS && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;} if ( d2==KS && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;} if ( d3==KS && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;} if ( d3==KS && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;} if ( d1==KS && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;} if ( d1==KS && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;} if ( d2==KS && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;} if ( d2==KS && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;} if ( d3==KS && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;} if ( d3==KS && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;} if ( d1==KL && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;} if ( d1==KL && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;} if ( d2==KL && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;} if ( d2==KL && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;} if ( d3==KL && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;} if ( d3==KL && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;} if ( d1==K0 && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;} if ( d1==K0 && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;} if ( d2==K0 && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;} if ( d2==K0 && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;} if ( d3==K0 && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;} if ( d3==K0 && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;} } if ( parnum == DP ) { //look for K- pi+ pi+ if ( d1==KB && d2==PIP && d3==PI0 ) { _flag=2; _d1=0; _d2=1; _d3=2;} if ( d1==KB && d3==PIP && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;} if ( d2==KB && d1==PIP && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;} if ( d2==KB && d3==PIP && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;} if ( d3==KB && d1==PIP && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;} if ( d3==KB && d2==PIP && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;} if ( d1==KL && d2==PIP && d3==PI0 ) { _flag=2; _d1=0; _d2=1; _d3=2;} if ( d1==KL && d3==PIP && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;} if ( d2==KL && d1==PIP && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;} if ( d2==KL && d3==PIP && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;} if ( d3==KL && d1==PIP && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;} if ( d3==KL && d2==PIP && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;} if ( d1==KS && d2==PIP && d3==PI0 ) { _flag=2; _d1=0; _d2=1; _d3=2;} if ( d1==KS && d3==PIP && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;} if ( d2==KS && d1==PIP && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;} if ( d2==KS && d3==PIP && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;} if ( d3==KS && d1==PIP && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;} if ( d3==KS && d2==PIP && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;} if ( d1==KM && d2==PIP && d3==PIP ) { _flag=1; _d1=0; _d2=1; _d3=2;} if ( d2==KM && d1==PIP && d3==PIP ) { _flag=1; _d1=1; _d2=0; _d3=2;} if ( d3==KM && d1==PIP && d2==PIP ) { _flag=1; _d1=2; _d2=0; _d3=1;} } if ( parnum == DM ) { //look for K- pi+ pi+ if ( d1==K0 && d2==PIM && d3==PI0 ) { _flag=2; _d1=0; _d2=1; _d3=2;} if ( d1==K0 && d3==PIM && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;} if ( d2==K0 && d1==PIM && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;} if ( d2==K0 && d3==PIM && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;} if ( d3==K0 && d1==PIM && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;} if ( d3==K0 && d2==PIM && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;} if ( d1==KL && d2==PIM && d3==PI0 ) { _flag=2; _d1=0; _d2=1; _d3=2;} if ( d1==KL && d3==PIM && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;} if ( d2==KL && d1==PIM && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;} if ( d2==KL && d3==PIM && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;} if ( d3==KL && d1==PIM && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;} if ( d3==KL && d2==PIM && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;} if ( d1==KS && d2==PIM && d3==PI0 ) { _flag=2; _d1=0; _d2=1; _d3=2;} if ( d1==KS && d3==PIM && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;} if ( d2==KS && d1==PIM && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;} if ( d2==KS && d3==PIM && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;} if ( d3==KS && d1==PIM && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;} if ( d3==KS && d2==PIM && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;} if ( d1==KP && d2==PIM && d3==PIM ) { _flag=1; _d1=0; _d2=1; _d3=2;} if ( d2==KP && d1==PIM && d3==PIM ) { _flag=1; _d1=1; _d2=0; _d3=2;} if ( d3==KP && d1==PIM && d2==PIM ) { _flag=1; _d1=2; _d2=0; _d3=1;} } if ( parnum == DSP ) { if ( d1==KM && d2==KP && d3==PIP ) { _flag=6; _d1=0; _d2=1; _d3=2; } if ( d1==KM && d3==KP && d2==PIP ) { _flag=6; _d1=0; _d2=2; _d3=1; } if ( d2==KM && d1==KP && d3==PIP ) { _flag=6; _d1=1; _d2=0; _d3=2; } if ( d2==KM && d3==KP && d1==PIP ) { _flag=6; _d1=1; _d2=2; _d3=0; } if ( d3==KM && d1==KP && d2==PIP ) { _flag=6; _d1=2; _d2=0; _d3=1; } if ( d3==KM && d2==KP && d1==PIP ) { _flag=6; _d1=2; _d2=1; _d3=0; } } if ( parnum == DSM ) { if ( d1==KP && d2==KM && d3==PIM ) { _flag=6; _d1=0; _d2=1; _d3=2; } if ( d1==KP && d3==KM && d2==PIM ) { _flag=6; _d1=0; _d2=2; _d3=1; } if ( d2==KP && d1==KM && d3==PIM ) { _flag=6; _d1=1; _d2=0; _d3=2; } if ( d2==KP && d3==KM && d1==PIM ) { _flag=6; _d1=1; _d2=2; _d3=0; } if ( d3==KP && d1==KM && d2==PIM ) { _flag=6; _d1=2; _d2=0; _d3=1; } if ( d3==KP && d2==KM && d1==PIM ) { _flag=6; _d1=2; _d2=1; _d3=0; } } if ( _flag==6) { _kkpi_params.push_back(EvtFlatteParam(MPI, MPI, 0.406)); _kkpi_params.push_back(EvtFlatteParam(MKP, MKP, 0.800)); } if ( _flag==0) { report(ERROR,"EvtGen") << "EvtDDaltiz: Invalid mode."<<endl; assert(0); } }

void EvtDDalitz::initProbMax (   )  [virtual]

Reimplemented from EvtDecayBase. { //probmax different for different modes! if ( _flag==1 ) {setProbMax(9.6);} if ( _flag==2 ) {setProbMax(147.9);} if ( _flag==3 ) {setProbMax(5000.0);} if ( _flag==4 ) {setProbMax(3000.0);} if ( _flag==5 ) {setProbMax(10000000.0);} if ( _flag==6 ) {setProbMax(50000.0);} }

void EvtDecayAmp::makeDecay (  EvtParticle *  p  )  [virtual, inherited]

Implements EvtDecayBase. { int ntimes=10000; int i,more; EvtSpinDensity rho; double prob,prob_max; _amp2.init(p->getId(),getNDaug(),getDaugs()); // report(INFO,"EvtGen") << "Decaying " << EvtPDL::name(p->getId()) << endl; do{ _daugsDecayedByParentModel=false; _weight = 1.0; decay(p); rho=_amp2.getSpinDensity(); prob=p->getSpinDensityForward().NormalizedProb(rho); if (prob<0.0) { report(ERROR,"EvtGen")<<"Negative prob:"<<p->getId().getId() <<" "<<p->getChannel()<<endl; report(ERROR,"EvtGen") << "rho_forward:"<<endl; report(ERROR,"EvtGen") << p->getSpinDensityForward(); report(ERROR,"EvtGen") << "rho decay:"<<endl; report(ERROR,"EvtGen") << rho <<endl; } if (prob!=prob) { report(DEBUG,"EvtGen") << "Forward density matrix:"<<endl; report(ERROR,"EvtGen") << p->getSpinDensityForward(); report(DEBUG,"EvtGen") << "Decay density matrix:"<<endl; report(ERROR,"EvtGen") << rho; report(DEBUG,"EvtGen") << "prob:"<<prob<<endl; report(DEBUG,"EvtGen") << "Particle:" <<EvtPDL::name(p->getId()).c_str()<<endl; report(DEBUG,"EvtGen") << "channel :"<<p->getChannel()<<endl; report(DEBUG,"EvtGen") << "Momentum:" << p->getP4() << " " << p->mass() << endl; if( p->getParent()!=0){ report(DEBUG,"EvtGen") << "parent:" <<EvtPDL::name( p->getParent()->getId()).c_str()<<endl; report(DEBUG,"EvtGen") << "parent channel :" <<p->getParent()->getChannel()<<endl; int i; report(DEBUG,"EvtGen") << "parent daughters :"; for (i=0;i<p->getParent()->getNDaug();i++){ report(DEBUG,"") << EvtPDL::name( p->getParent()->getDaug(i)->getId()).c_str() << " "; } report(DEBUG,"") << endl; report(DEBUG,"EvtGen") << "daughters :"; for (i=0;i<p->getNDaug();i++){ report(DEBUG,"") << EvtPDL::name( p->getDaug(i)->getId()).c_str() << " "; } report(DEBUG,"") << endl; report(DEBUG,"EvtGen") << "daughter momenta :" << endl;; for (i=0;i<p->getNDaug();i++){ report(DEBUG,"") << p->getDaug(i)->getP4() << " " << p->getDaug(i)->mass(); report(DEBUG,"") << endl; } } } prob/=_weight; prob_max = getProbMax(prob); p->setDecayProb(prob/prob_max); //report(INFO,"EvtGen") << "Prob,prob_max,weight:"<<prob<<" "<<prob_max<<" "<<_weight<<endl; more=prob<EvtRandom::Flat(prob_max); ntimes--; }while(ntimes&&more); //report(INFO,"EvtGen") << "Done\n"; if (ntimes==0){ report(DEBUG,"EvtGen") << "Tried accept/reject:10000" <<" times, and rejected all the times!"<<endl; report(DEBUG,"")<<p->getSpinDensityForward()<<endl; report(DEBUG,"EvtGen") << "Is therefore accepting the last event!"<<endl; report(DEBUG,"EvtGen") << "Decay of particle:"<< EvtPDL::name(p->getId()).c_str()<<"(channel:"<< p->getChannel()<<") with mass "<<p->mass()<<endl; int ii; for(ii=0;ii<p->getNDaug();ii++){ report(DEBUG,"EvtGen") <<"Daughter "<<ii<<":"<< EvtPDL::name(p->getDaug(ii)->getId()).c_str()<<" with mass "<< p->getDaug(ii)->mass()<<endl; } } EvtSpinDensity rho_list[10]; rho_list[0]=p->getSpinDensityForward(); EvtAmp ampcont; if (_amp2._pstates!=1){ ampcont=_amp2.contract(0,p->getSpinDensityForward()); // J2BB2 bugging here } else{ ampcont=_amp2; } // it may be that the parent decay model has already // done the decay - this should be rare and the // model better know what it is doing.. if ( !daugsDecayedByParentModel() ){ // report(INFO,"EvtGen") << "Found " << p->getNDaug() << " daughters\n"; for(i=0;i<p->getNDaug();i++){ rho.SetDim(_amp2.dstates[i]); if (_amp2.dstates[i]==1) { rho.Set(0,0,EvtComplex(1.0,0.0)); } else{ rho=ampcont.contract(_amp2._dnontrivial[i],_amp2); } if (!rho.Check()) { report(ERROR,"EvtGen") << "-------start error-------"<<endl; report(ERROR,"EvtGen")<<"forward rho failed Check:"<< EvtPDL::name(p->getId()).c_str()<<" "<<p->getChannel()<<" "<<i<<endl; report(ERROR,"EvtGen")<<"Parent:"<<EvtPDL::name(p->getParent()->getId()).c_str()<<endl; report(ERROR,"EvtGen")<<"GrandParent:"<<EvtPDL::name(p->getParent()->getParent()->getId()).c_str()<<endl; report(ERROR,"EvtGen")<<"GrandGrandParent:"<<EvtPDL::name(p->getParent()->getParent()->getParent()->getId()).c_str()<<endl; report(ERROR,"EvtGen") << rho; int ii; _amp2.dump(); for(ii=0;ii<i+1;ii++){ report(ERROR,"EvtGen") << rho_list[ii]; } report(ERROR,"EvtGen") << "-------Done with error-------"<<endl; } p->getDaug(i)->setSpinDensityForward(rho); p->getDaug(i)->decay(); rho_list[i+1]= p->getDaug(i)->getSpinDensityBackward(); if (_amp2.dstates[i]!=1){ ampcont=ampcont.contract(_amp2._dnontrivial[i],rho_list[i+1]); } } p->setSpinDensityBackward(_amp2.getBackwardSpinDensity(rho_list)); if (!p->getSpinDensityBackward().Check()) { report(ERROR,"EvtGen")<<"rho_backward failed Check"<< p->getId().getId()<<" "<<p->getChannel()<<endl; report(ERROR,"EvtGen") << p->getSpinDensityBackward(); } } if (getPHOTOS() || EvtRadCorr::alwaysRadCorr()) { EvtRadCorr::doRadCorr(p); } }

bool EvtDecayBase::matchingDecay (  const EvtDecayBase &  other  )  const [virtual, inherited]

{ if ( _ndaug != other._ndaug) return false; if ( _parent != other._parent) return false; std::vector<int> useDs; for ( unsigned int i=0; i<_ndaug; i++) useDs.push_back(0); for ( unsigned int i=0; i<_ndaug; i++) { bool foundIt=false; for ( unsigned int j=0; j<_ndaug; j++) { if ( useDs[j] == 1 ) continue; if ( _daug[i] == other._daug[j] && _daug[i].getAlias() == other._daug[j].getAlias()) { foundIt=true; useDs[j]=1; break; } } if ( foundIt==false) return false; } for ( unsigned int i=0; i<_ndaug; i++) if ( useDs[i]==0) return false; return true; }

void EvtDecayBase::noProbMax (   )  [inherited]

{ defaultprobmax=0; }

virtual int EvtDecayBase::nRealDaughters (   )  [inline, virtual, inherited]

Reimplemented in EvtBtoKD3P, and EvtVSSBMixCPT. { return _ndaug;}

void EvtDecayBase::printSummary (   )  [inherited]

{ int i; if (ntimes_prob>0) { report(INFO,"EvtGen") << "Calls="<<ntimes_prob<<" eff:"<< sum_prob/(probmax*ntimes_prob)<<" frac. max:"<<max_prob/probmax; report(INFO,"") <<" probmax:"<<probmax<<" max:"<<max_prob<<" : "; } report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> "; for(i=0;i<_ndaug;i++){ report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " "; } report(INFO,"") << " ("<<_modelname.c_str()<<"):"<< endl; }

double EvtDecayBase::resetProbMax (  double  prob  )  [inherited]

{ report(INFO,"EvtGen") << "Reseting prob max\n"; report(INFO,"EvtGen") << "prob > probmax:("<<prob<<">"<<probmax<<")"; report(INFO,"") << "("<<_modelname.c_str()<<")"; report(INFO,"") << EvtPDL::getStdHep(_parent)<<"->"; for( int i=0;i<_ndaug;i++){ report(INFO,"") << EvtPDL::getStdHep(_daug[i]) << " "; } report(INFO,"") << endl; probmax = 0.0; defaultprobmax = 0; ntimes_prob = 0; return prob; }

void EvtDecayBase::saveDecayInfo (  EvtId  ipar, int  ndaug, EvtId *  daug, int  narg, std::vector< std::string > &  args, std::string  name, double  brfr )  [inherited]

{ int i; _brfr=brfr; _ndaug=ndaug; _narg=narg; _parent=ipar; _dsum=0; if (_ndaug>0) { _daug=new EvtId [_ndaug]; for(i=0;i<_ndaug;i++){ _daug[i]=daug[i]; _dsum+=daug[i].getAlias(); } } else{ _daug=0; } if (_narg>0) { _args=new std::string[_narg+1]; for(i=0;i<_narg;i++){ _args[i]=args[i]; } } else{ _args = 0; } _modelname=name; this->init(); this->initProbMax(); if (_chkCharge){ this->checkQ(); } if (defaultprobmax){ report(INFO,"EvtGen") << "No default probmax for "; report(INFO,"") << "("<<_modelname.c_str()<<") "; report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> "; for(i=0;i<_ndaug;i++){ report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " "; } report(INFO,"") << endl; report(INFO,"") << "This is fine for development, but must be provided for production."<<endl; report(INFO,"EvtGen") << "Never fear though - the decay will use the \n"; report(INFO,"EvtGen") << "500 iterations to build up a good probmax \n"; report(INFO,"EvtGen") << "before accepting a decay. "<<endl; } }

void EvtDecayBase::setPHOTOS (   )  [inline, inherited]

{_photos=1;}

void EvtDecayBase::setProbMax (  double  prbmx  )  [inherited]

{ defaultprobmax=0; probmax=prbmx; }

void EvtDecayBase::setSummary (   )  [inline, inherited]

{_summary=1;}

void EvtDecayBase::setVerbose (   )  [inline, inherited]

{_verbose=1;}

void EvtDecayAmp::setWeight (  double  weight  )  [inline, inherited]

{_weight=weight;}

int EvtDecayBase::summary (   )  [inline, inherited]

{return _summary; }

int EvtDecayBase::verbose (   )  [inline, inherited]

{return _verbose; }

void EvtDecayAmp::vertex (  int *  i1, const EvtComplex &  amp )  [inline, inherited]

sets the amplitudes calculated in the decay objects { _amp2.vertex(i1,amp);}

void EvtDecayAmp::vertex (  int  i1, int  i2, int  i3, const EvtComplex &  amp )  [inline, inherited]

sets the amplitudes calculated in the decay objects {_amp2.vertex(i1,i2,i3,amp);}

void EvtDecayAmp::vertex (  int  i1, int  i2, const EvtComplex &  amp )  [inline, inherited]

sets the amplitudes calculated in the decay objects {_amp2.vertex(i1,i2,amp);}

void EvtDecayAmp::vertex (  int  i1, const EvtComplex &  amp )  [inline, inherited]

sets the amplitudes calculated in the decay objects {_amp2.vertex(i1,amp);}

void EvtDecayAmp::vertex (  const EvtComplex &  amp  )  [inline, inherited]

sets the amplitudes calculated in the decay objects {_amp2.vertex(amp);}

---

Member Data Documentation

EvtAmp EvtDecayAmp::_amp2 [protected, inherited]

int EvtDDalitz::_d1 [private]

int EvtDDalitz::_d2 [private]

int EvtDDalitz::_d3 [private]

bool EvtDecayBase::_daugsDecayedByParentModel [protected, inherited]

int EvtDDalitz::_flag [private]

vector<EvtFlatteParam> EvtDDalitz::_kkpi_params [private]

---

The documentation for this class was generated from the following files:

• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/Generator/BesEvtGen/BesEvtGen-00-01-11/src/EvtGen/EvtGenModels/EvtDDalitz.hh
• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/Generator/BesEvtGen/BesEvtGen-00-01-11/src/EvtGen/EvtGenModels/EvtDDalitz.cc

---