00001 #include "EvtGenBase/EvtPatches.hh"
00002 #include "EvtGenBase/EvtMHelAmp.hh"
00003 #include "EvtGenBase/EvtKine.hh"
00004 #include <cstdlib>
00005
00006 using std::endl;
00007
00008 EvtMHelAmp::EvtMHelAmp( const EvtId& id, EvtMLineShape * lineshape,
00009 const vector<EvtMNode* >& children, const vector<EvtComplex>& elem )
00010 {
00011 _id = id;
00012 _twospin = EvtSpinType::getSpin2( EvtPDL::getSpinType( id ) );
00013 _parent = NULL;
00014 _lineshape = lineshape;
00015
00016 _elem = elem;
00017
00018 vector<EvtSpinType::spintype> type;
00019 for( int i=0; i<children.size(); ++i ) {
00020 _children.push_back( children[i] );
00021 type.push_back( children[i]->getspintype() );
00022 const vector<int> &res = children[i]->getresonance();
00023 for( int j=0; j<res.size(); ++j )
00024 _resonance.push_back( res[j] );
00025 children[i]->setparent( this );
00026 }
00027
00028
00029 _amp = EvtSpinAmp( type );
00030 vector<int> index = _amp.iterinit();
00031 int i = 0;
00032 do {
00033 if( !_amp.allowed(index) )
00034 _amp( index ) = 0.0;
00035 else if( abs(index[0] - index[1]) > _twospin )
00036 _amp( index ) = 0.0;
00037 else {
00038 _amp( index ) = elem[i];
00039 ++i;
00040 }
00041 } while( _amp.iterate( index ) );
00042 if(elem.size() != i) {
00043 report(ERROR,"EvtGen")
00044 <<"Wrong number of elements input in helicity amplitude."<<endl;
00045 ::abort();
00046 }
00047
00048 if( children.size() > 2 ) {
00049 report(ERROR,"EvtGen")
00050 <<"Helicity amplitude formalism can only handle two body resonances"
00051 <<endl;
00052 ::abort();
00053 }
00054 }
00055
00056 EvtSpinAmp EvtMHelAmp::amplitude( const vector<EvtVector4R> &
00057 product ) const
00058 {
00059 EvtVector4R d = _children[0]->get4vector(product);
00060 double phi, theta;
00061
00062 if( _parent == NULL ) {
00063
00064
00065
00066
00067 phi = atan2( d.get(1), d.get(2) );
00068 theta = acos( d.get(3)/d.d3mag() );
00069
00070 } else {
00071
00072
00073
00074 EvtVector4R p = _parent->get4vector(product);
00075 EvtVector4R q = get4vector(product);
00076
00077
00078
00079 EvtVector4R g = _parent->getparent()==NULL ?
00080 EvtVector4R(0.0, 0.0, 0.0, 1.0) :
00081 _parent->getparent()->get4vector(product);
00082
00083 theta = acos(EvtDecayAngle(p, q, d));
00084 phi = EvtDecayAnglePhi( g, p, q, d );
00085
00086 }
00087
00088 vector<EvtSpinType::spintype> types( 3 );
00089 types[0] = getspintype();
00090 types[1] = _children[0]->getspintype();
00091 types[2] = _children[1]->getspintype();
00092 EvtSpinAmp amp( types, EvtComplex(0.0, 0.0) );
00093 vector<int> index = amp.iterallowedinit();
00094
00095 do {
00096 if( abs(index[1]-index[2]) > _twospin ) continue;
00097 amp(index) +=
00098 conj(wignerD(_twospin,index[0],index[1]-index[2],phi,theta,0.0)) *
00099 _amp(index[1],index[2]);
00100 } while(amp.iterateallowed(index));
00101
00102 EvtSpinAmp amp0 = _children[0]->amplitude(product);
00103 EvtSpinAmp amp1 = _children[1]->amplitude(product);
00104
00105 amp.extcont( amp0, 1, 0 );
00106 amp.extcont( amp1, 1, 0 );
00107
00108 amp *= sqrt( ( _twospin + 1 ) / ( 2 * EvtConst::twoPi ) ) *
00109 _children[0]->line(product) * _children[1]->line(product);
00110
00111 return amp;
00112 }
00113
00114 EvtMNode * EvtMHelAmp::duplicate() const
00115 {
00116 vector<EvtMNode *> children;
00117
00118 for( int i=0; i<_children.size(); ++i ) {
00119 children.push_back( _children[i]->duplicate() );
00120 }
00121
00122 EvtMLineShape * lineshape = _lineshape->duplicate();
00123 EvtMHelAmp * ret = new EvtMHelAmp( _id, lineshape, children, _elem );
00124 lineshape->setres( ret );
00125
00126 return ret;
00127 }