00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #ifdef WIN32
00016 #pragma warning( disable : 4786 )
00017
00018 #endif
00019
00020 #include "EvtGenModels/EvtHypNonLepton.hh"
00021 #include "EvtGenBase/EvtParticle.hh"
00022 #include "EvtGenBase/EvtPDL.hh"
00023 #include "EvtGenBase/EvtDiracSpinor.hh"
00024 #include "EvtGenBase/EvtTensor4C.hh"
00025 #include "EvtGenBase/EvtVector4C.hh"
00026 #include "EvtGenBase/EvtVector4R.hh"
00027 #include "EvtGenBase/EvtComplex.hh"
00028 #include "EvtGenBase/EvtGammaMatrix.hh"
00029
00030
00031 EvtHypNonLepton::~EvtHypNonLepton() {}
00032
00033 EvtDecayBase* EvtHypNonLepton::clone(){
00034 return new EvtHypNonLepton;
00035 }
00036
00037 void EvtHypNonLepton::getName(std::string& model_name){
00038 model_name= "HypNonLepton";
00039 }
00040
00041 void EvtHypNonLepton::init(){
00042
00043 if(getNArg()<2 || getNArg()>3){
00044 report(ERROR,"EvtGen") << " ERROR: EvtHypNonLepton generator expected 2 or 3 arguments but found: " << getNArg() << std::endl;
00045 report(INFO ,"EvtGen") << " 1. Decay asymmetry parameter - alpha" << std::endl;
00046 report(INFO ,"EvtGen") << " 2. Parameter phi - in degrees (not radians)" << std::endl;
00047 report(INFO ,"EvtGen") << " 3. Note on every x-th decay" << std::endl;
00048 ::abort();
00049 }
00050
00051 if(getNDaug()!=2){
00052 report(ERROR,"EvtGen") << " ERROR: EvtHypNonLepton generator expected 2 daughters but found: " << getNDaug() << std::endl;
00053 ::abort();
00054 }
00055
00056
00057 if(EvtSpinType::getSpin2(EvtPDL::getSpinType(getParentId()))!=1){
00058 report(ERROR,"EvtGen") << " ERROR: EvtHypNonLepton generator expected dirac parent particle, but found " << EvtSpinType::getSpin2(EvtPDL::getSpinType(getParentId())) << " spin degrees of freedom" << std::endl;
00059 ::abort();
00060 }
00061 if(EvtSpinType::getSpin2(EvtPDL::getSpinType(getDaug(0)))!=1){
00062 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;
00063 ::abort();
00064 }
00065 if(EvtSpinType::getSpin2(EvtPDL::getSpinType(getDaug(1)))!=0){
00066 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;
00067 ::abort();
00068 }
00069
00070
00071 m_alpha = getArg(0);
00072 m_phi = getArg(1)*EvtConst::pi/180;
00073 if(getNArg()==3) m_noTries = static_cast<long>(getArg(2));
00074 else m_noTries = 0;
00075
00076
00077 double p,M,m1,m2;
00078 double p_to_s,beta,delta,gamma;
00079
00080 M = EvtPDL::getMass(getParentId());
00081 m1 = EvtPDL::getMass(getDaug(0));
00082 m2 = EvtPDL::getMass(getDaug(1));
00083
00084 if(m1+m2>=M){
00085 report(ERROR,"EvtGen") << " ERROR: EvtHypNonLepton found impossible decay: " << M << " --> " << m1 << " + " << m2 << " GeV\n" << std::endl;
00086 ::abort();
00087 }
00088
00089 p = sqrt(M*M-(m1+m2)*(m1+m2))*sqrt(M*M-(m1-m2)*(m1-m2))/2./M;
00090
00091 beta = sqrt(1.-m_alpha*m_alpha)*sin(m_phi);
00092 delta = -atan2(beta,m_alpha);
00093 gamma = sqrt(1.-m_alpha*m_alpha-beta*beta);
00094 p_to_s = sqrt((1.-gamma)/(1.+gamma));
00095
00096 m_B_to_A = p_to_s*(m1+sqrt(p*p+m1*m1))/p*EvtComplex(cos(delta),sin(delta));
00097
00098 }
00099
00100 void EvtHypNonLepton::initProbMax(){
00101
00102 double maxProb,m1,m2,M,p;
00103
00104 M=EvtPDL::getMass(getParentId());
00105 m1=EvtPDL::getMass(getDaug(0));
00106 m2=EvtPDL::getMass(getDaug(1));
00107
00108 if(m1+m2>=M){
00109 report(ERROR,"EvtGen") << " ERROR: EvtHypNonLepton found impossible decay: " << M << " --> " << m1 << " + " << m2 << " GeV\n" << std::endl;
00110 ::abort();
00111 }
00112
00113 p=sqrt(M*M-(m1+m2)*(m1+m2))*sqrt(M*M-(m1-m2)*(m1-m2))/2/M;
00114 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));
00115
00116
00117 setProbMax(maxProb);
00118 report(INFO,"EvtGen") << " EvtHypNonLepton set up maximum probability to " << maxProb << std::endl;
00119
00120 }
00121
00122 void EvtHypNonLepton::decay(EvtParticle* parent){
00123
00124 parent->initializePhaseSpace(getNDaug(),getDaugs());
00125 calcAmp(&_amp2,parent);
00126
00127 }
00128
00129 void EvtHypNonLepton::calcAmp(EvtAmp *amp,EvtParticle *parent){
00130
00131 static long noTries=0;
00132 int i;
00133 EvtComplex Matrix[2][2],B_to_A;
00134
00135
00136
00137
00138 for(i=0;i<4;i++){
00139
00140 Matrix[i/2][i%2] = EvtLeptonSCurrent(parent->sp(i/2),parent->getDaug(0)->spParent(i%2));
00141
00142 Matrix[i/2][i%2] -= m_B_to_A*EvtLeptonPCurrent(parent->sp(i/2),parent->getDaug(0)->spParent(i%2));
00143
00144
00145
00146
00147 amp->vertex(i/2,i%2,Matrix[i/2][i%2]);
00148 }
00149
00150 if(m_noTries>0) if(!((++noTries)%m_noTries)) report(DEBUG,"EvtGen") << " EvtHypNonLepton already finished " << noTries << " matrix element calculations" << std::endl;
00151 }