Converting TEvtGen to native cmake
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenBase / EvtMHelAmp.cpp
CommitLineData
da0e9ce3 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
7using std::endl;
8
9EvtMHelAmp::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
57EvtSpinAmp 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
115EvtMNode * 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}