TEvtGen/EvtGenBase/EvtMHelAmp.cpp

   1 #include "EvtGenBase/EvtPatches.hh"
   2 #include "EvtGenBase/EvtMHelAmp.hh"
   3 #include "EvtGenBase/EvtKine.hh"
   4 #include "EvtGenBase/EvtReport.hh"
   5 #include <stdlib.h>
   6
   7 using std::endl;
   8
   9 EvtMHelAmp::EvtMHelAmp( const EvtId& id, EvtMLineShape * lineshape,
  10         const vector<EvtMNode* >& children, const vector<EvtComplex>& elem )
  11 {
  12     _id = id;
  13     _twospin = EvtSpinType::getSpin2( EvtPDL::getSpinType( id ) );
  14     _parent = NULL;
  15     _lineshape = lineshape;
  16
  17     _elem = elem;
  18
  19     vector<EvtSpinType::spintype> type;
  20     for(size_t i=0; i<children.size(); ++i ) {
  21         _children.push_back( children[i] );
  22         type.push_back( children[i]->getspintype() );
  23         const vector<int> &res = children[i]->getresonance();
  24         for(size_t j=0; j<res.size(); ++j )
  25             _resonance.push_back( res[j] );
  26         children[i]->setparent( this );
  27     }
  28
  29     // XXX New code - bugs could appear here XXX
  30     _amp = EvtSpinAmp( type );
  31     vector<int> index = _amp.iterinit();
  32     size_t i = 0;
  33     do {
  34         if( !_amp.allowed(index) )
  35            _amp( index ) = 0.0;
  36         else if( abs(index[0] - index[1]) > _twospin )
  37             _amp( index ) = 0.0;
  38         else {
  39             _amp( index ) = elem[i];
  40             ++i;
  41         }
  42     } while( _amp.iterate( index ) );
  43     if(elem.size() != i) {
  44         report(ERROR,"EvtGen")
  45             <<"Wrong number of elements input in helicity amplitude."<<endl;
  46         ::abort();
  47     }
  48
  49     if( children.size() > 2 ) {
  50         report(ERROR,"EvtGen")
  51             <<"Helicity amplitude formalism can only handle two body resonances"
  52             <<endl;
  53         ::abort();
  54     }
  55 }
  56
  57 EvtSpinAmp EvtMHelAmp::amplitude( const vector<EvtVector4R> &
  58         product ) const
  59 {
  60     EvtVector4R d = _children[0]->get4vector(product);
  61     double phi, theta;
  62
  63     if( _parent == NULL ) {
  64
  65         // This means that we're calculating the first level and we need to just
  66         // calculate the polar and azymuthal angles daughters in rest frame of
  67         // this (root) particle (this is automatic).
  68         phi = atan2( d.get(1), d.get(2) );
  69         theta = acos( d.get(3)/d.d3mag() );
  70
  71     } else {
  72
  73         // We have parents therefore calculate things in correct coordinate
  74         // system
  75         EvtVector4R p = _parent->get4vector(product);
  76         EvtVector4R q = get4vector(product);
  77
  78         // See if we have a grandparent - if no then the z-axis is defined by
  79         // the z-axis of the root particle
  80         EvtVector4R g = _parent->getparent()==NULL ?
  81             EvtVector4R(0.0, 0.0, 0.0, 1.0) :
  82             _parent->getparent()->get4vector(product);
  83
  84         theta = acos(EvtDecayAngle(p, q, d));
  85         phi = EvtDecayAnglePhi( g, p, q, d );
  86
  87     }
  88
  89     vector<EvtSpinType::spintype> types( 3 );
  90     types[0] = getspintype();
  91     types[1] = _children[0]->getspintype();
  92     types[2] = _children[1]->getspintype();
  93     EvtSpinAmp amp( types, EvtComplex(0.0, 0.0) );
  94     vector<int> index = amp.iterallowedinit();
  95
  96     do {
  97         if( abs(index[1]-index[2]) > _twospin ) continue;
  98         amp(index) +=
  99             conj(wignerD(_twospin,index[0],index[1]-index[2],phi,theta,0.0)) *
 100             _amp(index[1],index[2]);
 101     } while(amp.iterateallowed(index));
 102
 103     EvtSpinAmp amp0 = _children[0]->amplitude(product);
 104     EvtSpinAmp amp1 = _children[1]->amplitude(product);
 105
 106     amp.extcont( amp0, 1, 0 );
 107     amp.extcont( amp1, 1, 0 );
 108
 109     amp *= sqrt( ( _twospin + 1 ) / ( 2 * EvtConst::twoPi ) ) *
 110         _children[0]->line(product) * _children[1]->line(product);
 111
 112     return amp;
 113 }
 114
 115 EvtMNode * EvtMHelAmp::duplicate() const
 116 {
 117     vector<EvtMNode *> children;
 118
 119     for(size_t i=0; i<_children.size(); ++i ) {
 120         children.push_back( _children[i]->duplicate() );
 121     }
 122
 123     EvtMLineShape * lineshape = _lineshape->duplicate();
 124     EvtMHelAmp * ret = new EvtMHelAmp( _id, lineshape, children, _elem );
 125     lineshape->setres( ret );
 126
 127     return ret;
 128 }