Line data    Source code

       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           0 : EvtMHelAmp::EvtMHelAmp( const EvtId& id, EvtMLineShape * lineshape, 
      10             :         const vector<EvtMNode* >& children, const vector<EvtComplex>& elem ) 
      11           0 : {
      12           0 :     _id = id;
      13           0 :     _twospin = EvtSpinType::getSpin2( EvtPDL::getSpinType( id ) );
      14           0 :     _parent = NULL;
      15           0 :     _lineshape = lineshape;
      16             :     
      17           0 :     _elem = elem;
      18             : 
      19           0 :     vector<EvtSpinType::spintype> type;
      20           0 :     for(size_t i=0; i<children.size(); ++i ) {
      21           0 :         _children.push_back( children[i] );
      22           0 :         type.push_back( children[i]->getspintype() );
      23           0 :         const vector<int> &res = children[i]->getresonance();
      24           0 :         for(size_t j=0; j<res.size(); ++j )
      25           0 :             _resonance.push_back( res[j] );
      26           0 :         children[i]->setparent( this );
      27             :     }
      28             :    
      29             :     // XXX New code - bugs could appear here XXX
      30           0 :     _amp = EvtSpinAmp( type );
      31           0 :     vector<int> index = _amp.iterinit();
      32             :     size_t i = 0;
      33           0 :     do {
      34           0 :         if( !_amp.allowed(index) )
      35           0 :            _amp( index ) = 0.0;
      36           0 :         else if( abs(index[0] - index[1]) > _twospin )
      37           0 :             _amp( index ) = 0.0;
      38             :         else {
      39           0 :             _amp( index ) = elem[i];
      40           0 :             ++i;
      41             :         }
      42           0 :     } while( _amp.iterate( index ) );
      43           0 :     if(elem.size() != i) {
      44           0 :         report(Severity::Error,"EvtGen")
      45           0 :             <<"Wrong number of elements input in helicity amplitude."<<endl;
      46           0 :         ::abort();
      47             :     }
      48             : 
      49           0 :     if( children.size() > 2 ) {
      50           0 :         report(Severity::Error,"EvtGen")
      51           0 :             <<"Helicity amplitude formalism can only handle two body resonances"
      52           0 :             <<endl;
      53           0 :         ::abort();
      54             :     }
      55           0 : }
      56             : 
      57             : EvtSpinAmp EvtMHelAmp::amplitude( const vector<EvtVector4R> &
      58             :         product ) const
      59             : {
      60           0 :     EvtVector4R d = _children[0]->get4vector(product);
      61             :     double phi, theta;
      62             :     
      63           0 :     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           0 :         phi = atan2( d.get(1), d.get(2) ); 
      69           0 :         theta = acos( d.get(3)/d.d3mag() );
      70             : 
      71           0 :     } else {
      72             : 
      73             :         // We have parents therefore calculate things in correct coordinate
      74             :         // system
      75           0 :         EvtVector4R p = _parent->get4vector(product);
      76           0 :         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           0 :         EvtVector4R g = _parent->getparent()==NULL ?
      81           0 :             EvtVector4R(0.0, 0.0, 0.0, 1.0) :
      82           0 :             _parent->getparent()->get4vector(product);
      83             : 
      84           0 :         theta = acos(EvtDecayAngle(p, q, d));
      85           0 :         phi = EvtDecayAnglePhi( g, p, q, d );
      86             : 
      87           0 :     }
      88             : 
      89           0 :     vector<EvtSpinType::spintype> types( 3 );
      90           0 :     types[0] = getspintype();
      91           0 :     types[1] = _children[0]->getspintype();
      92           0 :     types[2] = _children[1]->getspintype();
      93           0 :     EvtSpinAmp amp( types, EvtComplex(0.0, 0.0) );
      94           0 :     vector<int> index = amp.iterallowedinit();
      95             : 
      96             :     do {
      97           0 :         if( abs(index[1]-index[2]) > _twospin ) continue;
      98           0 :         amp(index) +=
      99           0 :             conj(wignerD(_twospin,index[0],index[1]-index[2],phi,theta,0.0)) *
     100           0 :             _amp(index[1],index[2]); 
     101           0 :     } while(amp.iterateallowed(index));
     102             : 
     103           0 :     EvtSpinAmp amp0 = _children[0]->amplitude(product);
     104           0 :     EvtSpinAmp amp1 = _children[1]->amplitude(product);
     105             : 
     106           0 :     amp.extcont( amp0, 1, 0 );
     107           0 :     amp.extcont( amp1, 1, 0 );
     108             : 
     109           0 :     amp *= sqrt( ( _twospin + 1 ) / ( 2 * EvtConst::twoPi ) ) *
     110           0 :         _children[0]->line(product) * _children[1]->line(product);
     111             : 
     112             :     return amp;
     113           0 : }
     114             : 
     115             : EvtMNode * EvtMHelAmp::duplicate() const
     116             : {
     117           0 :     vector<EvtMNode *> children;
     118             : 
     119           0 :     for(size_t i=0; i<_children.size(); ++i ) {
     120           0 :         children.push_back( _children[i]->duplicate() );
     121             :     }
     122             :     
     123           0 :     EvtMLineShape * lineshape = _lineshape->duplicate();
     124           0 :     EvtMHelAmp * ret = new EvtMHelAmp( _id, lineshape, children, _elem );
     125           0 :     lineshape->setres( ret );
     126             : 
     127           0 :     return ret;
     128           0 : }