]>
Commit | Line | Data |
---|---|---|
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 | ||
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 | } |