TEvtGen/EvtGen/EvtGenModels/EvtHypNonLepton.cpp

   1 //-----------------------------------------------------------------------------------------------
   2 //
   3 // Module: EvtHypNonLepton.cpp
   4 //
   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
   7 //
   8 // Modification history:
   9 //
  10 //  09/02/2009  PR   Corrected Delta sign
  11 //  20/02/2005  PR   Module created according to PHSP and Lb2Lll model
  12 //
  13 //-----------------------------------------------------------------------------------------------
  14
  15 #ifdef WIN32
  16 #pragma warning( disable : 4786 )
  17 // Disable anoying warning about symbol size
  18 #endif
  19
  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"
  29
  30
  31 EvtHypNonLepton::~EvtHypNonLepton() {}
  32
  33 EvtDecayBase* EvtHypNonLepton::clone(){
  34   return new EvtHypNonLepton;
  35 }
  36
  37 std::string EvtHypNonLepton::getName(){
  38   return "HypNonLepton";
  39 }
  40
  41 void EvtHypNonLepton::init(){
  42
  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;
  48     ::abort();
  49   }
  50
  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;
  53     ::abort();
  54   }
  55
  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;
  59     ::abort();
  60   }
  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;
  63     ::abort();
  64   }
  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;
  67     ::abort();
  68   }
  69
  70   // Read all parameters
  71   m_alpha = getArg(0);
  72   m_phi   = getArg(1)*EvtConst::pi/180;
  73   if(getNArg()==3) m_noTries = static_cast<long>(getArg(2));
  74   else             m_noTries = 0;
  75
  76   // calculate additional parameters
  77   double p,M,m1,m2;
  78   double p_to_s,beta,delta,gamma;
  79
  80   M  = EvtPDL::getMass(getParentId());
  81   m1 = EvtPDL::getMass(getDaug(0));
  82   m2 = EvtPDL::getMass(getDaug(1));
  83
  84   if(m1+m2>=M){
  85     report(ERROR,"EvtGen") << " ERROR: EvtHypNonLepton found impossible decay: " << M << " --> " << m1 << " + " << m2 << " GeV\n" << std::endl;
  86     ::abort();
  87   }
  88
  89   p = sqrt(M*M-(m1+m2)*(m1+m2))*sqrt(M*M-(m1-m2)*(m1-m2))/2./M;
  90
  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));
  95
  96   m_B_to_A = p_to_s*(m1+sqrt(p*p+m1*m1))/p*EvtComplex(cos(delta),sin(delta));
  97
  98 }
  99
 100 void EvtHypNonLepton::initProbMax(){
 101
 102   double maxProb,m1,m2,M,p;
 103
 104   M=EvtPDL::getMass(getParentId());
 105   m1=EvtPDL::getMass(getDaug(0));
 106   m2=EvtPDL::getMass(getDaug(1));
 107
 108   if(m1+m2>=M){
 109     report(ERROR,"EvtGen") << " ERROR: EvtHypNonLepton found impossible decay: " << M << " --> " << m1 << " + " << m2 << " GeV\n" << std::endl;
 110     ::abort();
 111   }
 112
 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;
 116
 117   setProbMax(maxProb);
 118   report(INFO,"EvtGen") << " EvtHypNonLepton set up maximum probability to " << maxProb << std::endl;
 119
 120 }
 121
 122 void EvtHypNonLepton::decay(EvtParticle* parent){
 123
 124   parent->initializePhaseSpace(getNDaug(),getDaugs());
 125   calcAmp(&_amp2,parent);
 126
 127 }
 128
 129 void EvtHypNonLepton::calcAmp(EvtAmp *amp,EvtParticle *parent){
 130
 131   static long noTries=0;
 132   int i;
 133   EvtComplex Matrix[2][2];
 134
 135   //G_F  = 1.16637e-5;
 136   //M_pi = 0.13957;
 137
 138   for(i=0;i<4;i++){
 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]);
 148   }
 149
 150   if(m_noTries>0) if(!((++noTries)%m_noTries)) report(DEBUG,"EvtGen") << " EvtHypNonLepton already finished " << noTries << " matrix element calculations" << std::endl;
 151 }