1 //-----------------------------------------------------------------------------------------------
3 // Module: EvtHypNonLepton.cpp
5 // Desription: Routine to implement Hyperon(s=1/2) -> Baryon(s=1/2) + Scalar decays accroding to
6 // Review Of Particle Physics 2004, Phys.Lett.B, Vol.592, p.864
8 // Modification history:
10 // 09/02/2009 PR Corrected Delta sign
11 // 20/02/2005 PR Module created according to PHSP and Lb2Lll model
13 //-----------------------------------------------------------------------------------------------
16 #pragma warning( disable : 4786 )
17 // Disable anoying warning about symbol size
20 #include "EvtGenModels/EvtHypNonLepton.hh"
21 #include "EvtGenBase/EvtParticle.hh"
22 #include "EvtGenBase/EvtPDL.hh"
23 #include "EvtGenBase/EvtDiracSpinor.hh"
24 #include "EvtGenBase/EvtTensor4C.hh"
25 #include "EvtGenBase/EvtVector4C.hh"
26 #include "EvtGenBase/EvtVector4R.hh"
27 #include "EvtGenBase/EvtComplex.hh"
28 #include "EvtGenBase/EvtGammaMatrix.hh"
31 EvtHypNonLepton::~EvtHypNonLepton() {}
33 EvtDecayBase* EvtHypNonLepton::clone(){
34 return new EvtHypNonLepton;
37 std::string EvtHypNonLepton::getName(){
38 return "HypNonLepton";
41 void EvtHypNonLepton::init(){
43 if(getNArg()<2 || getNArg()>3){ // alpha phi gamma delta
44 report(ERROR,"EvtGen") << " ERROR: EvtHypNonLepton generator expected 2 or 3 arguments but found: " << getNArg() << std::endl;
45 report(INFO ,"EvtGen") << " 1. Decay asymmetry parameter - alpha" << std::endl;
46 report(INFO ,"EvtGen") << " 2. Parameter phi - in degrees (not radians)" << std::endl;
47 report(INFO ,"EvtGen") << " 3. Note on every x-th decay" << std::endl;
51 if(getNDaug()!=2){ // Check that there are 2 daughters only
52 report(ERROR,"EvtGen") << " ERROR: EvtHypNonLepton generator expected 2 daughters but found: " << getNDaug() << std::endl;
56 // Check particles spins
57 if(EvtSpinType::getSpin2(EvtPDL::getSpinType(getParentId()))!=1){
58 report(ERROR,"EvtGen") << " ERROR: EvtHypNonLepton generator expected dirac parent particle, but found " << EvtSpinType::getSpin2(EvtPDL::getSpinType(getParentId())) << " spin degrees of freedom" << std::endl;
61 if(EvtSpinType::getSpin2(EvtPDL::getSpinType(getDaug(0)))!=1){
62 report(ERROR,"EvtGen") << " ERROR: EvtHypNonLepton generator expected the first child to be dirac particle, but found " << EvtSpinType::getSpin2(EvtPDL::getSpinType(getDaug(0))) << " spin degrees of freedom" << std::endl;
65 if(EvtSpinType::getSpin2(EvtPDL::getSpinType(getDaug(1)))!=0){
66 report(ERROR,"EvtGen") << " ERROR: EvtHypNonLepton generator expected the second child to be scalar particle, but found " << EvtSpinType::getSpin2(EvtPDL::getSpinType(getDaug(1))) << " spin degrees of freedom" << std::endl;
70 // Read all parameters
72 m_phi = getArg(1)*EvtConst::pi/180;
73 if(getNArg()==3) m_noTries = static_cast<long>(getArg(2));
76 // calculate additional parameters
78 double p_to_s,beta,delta,gamma;
80 M = EvtPDL::getMass(getParentId());
81 m1 = EvtPDL::getMass(getDaug(0));
82 m2 = EvtPDL::getMass(getDaug(1));
85 report(ERROR,"EvtGen") << " ERROR: EvtHypNonLepton found impossible decay: " << M << " --> " << m1 << " + " << m2 << " GeV\n" << std::endl;
89 p = sqrt(M*M-(m1+m2)*(m1+m2))*sqrt(M*M-(m1-m2)*(m1-m2))/2./M;
91 beta = sqrt(1.-m_alpha*m_alpha)*sin(m_phi);
92 delta = -atan2(beta,m_alpha);
93 gamma = sqrt(1.-m_alpha*m_alpha-beta*beta);
94 p_to_s = sqrt((1.-gamma)/(1.+gamma));
96 m_B_to_A = p_to_s*(m1+sqrt(p*p+m1*m1))/p*EvtComplex(cos(delta),sin(delta));
100 void EvtHypNonLepton::initProbMax(){
102 double maxProb,m1,m2,M,p;
104 M=EvtPDL::getMass(getParentId());
105 m1=EvtPDL::getMass(getDaug(0));
106 m2=EvtPDL::getMass(getDaug(1));
109 report(ERROR,"EvtGen") << " ERROR: EvtHypNonLepton found impossible decay: " << M << " --> " << m1 << " + " << m2 << " GeV\n" << std::endl;
113 p=sqrt(M*M-(m1+m2)*(m1+m2))*sqrt(M*M-(m1-m2)*(m1-m2))/2/M;
114 maxProb=16*M*(sqrt(p*p+m1*m1)+m1+abs(m_B_to_A)*abs(m_B_to_A)*(sqrt(p*p+m1*m1)-m1));
115 //maxProb *= G_F*M_pi*M_pi;
118 report(INFO,"EvtGen") << " EvtHypNonLepton set up maximum probability to " << maxProb << std::endl;
122 void EvtHypNonLepton::decay(EvtParticle* parent){
124 parent->initializePhaseSpace(getNDaug(),getDaugs());
125 calcAmp(&_amp2,parent);
129 void EvtHypNonLepton::calcAmp(EvtAmp *amp,EvtParticle *parent){
131 static long noTries=0;
133 EvtComplex Matrix[2][2],B_to_A;
139 //std::cout << "--------------------------------------------------" << std::endl;
140 Matrix[i/2][i%2] = EvtLeptonSCurrent(parent->sp(i/2),parent->getDaug(0)->spParent(i%2));
141 //std::cout << "Matrix = " << Matrix[i/2][i%2] << std::endl;
142 Matrix[i/2][i%2] -= m_B_to_A*EvtLeptonPCurrent(parent->sp(i/2),parent->getDaug(0)->spParent(i%2));
143 //std::cout << "Matrix = " << Matrix[i/2][i%2] << std::endl;
144 //Matrix[i/2][i%2] *= G_F*M_pi*M_pi;
145 //std::cout << "Matrix = " << Matrix[i/2][i%2] << std::endl;
146 //std::cout << "--------------------------------------------------" << std::endl;
147 amp->vertex(i/2,i%2,Matrix[i/2][i%2]);
150 if(m_noTries>0) if(!((++noTries)%m_noTries)) report(DEBUG,"EvtGen") << " EvtHypNonLepton already finished " << noTries << " matrix element calculations" << std::endl;