[u/mrichter/AliRoot.git] / TEvtGen / EvtGenBase / EvtMTree.cxx

#include <stdio.h>
#include <stdlib.h>
#include <algorithm>

#include "EvtGenBase/EvtParticle.hh"

#include "EvtGenBase/EvtMTree.hh"
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtKine.hh"
#include "EvtGenBase/EvtReport.hh"

// Make sure to include Lineshapes here
#include "EvtGenBase/EvtMTrivialLS.hh"
#include "EvtGenBase/EvtMBreitWigner.hh"

// Make sure to include Parametrizations
#include "EvtGenBase/EvtMHelAmp.hh"

using std::endl;

EvtMTree::EvtMTree( const EvtId * idtbl, unsigned int ndaug )
{
    for( size_t i=0; i<ndaug; ++i ) {
        _lbltbl.push_back( EvtPDL::name( idtbl[i] ) );
    }
}

EvtMTree::~EvtMTree()
{
    for(size_t i=0; i<_root.size(); ++i) delete _root[i];
}

bool EvtMTree::parsecheck( char arg, const string& chars )
{
    bool ret = false;

    for(size_t i=0; i<chars.size(); ++i) {
        ret = ret || (chars[i]==arg);
    }

    return ret;
}

vector<EvtMNode *> EvtMTree::makeparticles( const string& strid ) 
{
    vector<EvtMNode *> particles;
    vector<int> labels;
   
    for( size_t i = 0; i<_lbltbl.size(); ++i ) {
        if( _lbltbl[i] == strid ) labels.push_back( i );
    }
    
    if( labels.size() == 0 ) {
        report(ERROR,"EvtGen")<<"Error unknown particle label "<<strid<<endl;
        ::abort();
    }

    for( size_t i = 0; i<labels.size(); ++i )
        particles.push_back( new EvtMParticle( labels[i], EvtPDL::getId( strid ) ) );

    return particles;
}

EvtMRes * EvtMTree::makeresonance( const EvtId& id, const string& ls,
        const vector<string>& lsarg, const string& type,
        const vector<EvtComplex>& amps, const vector<EvtMNode *>& children )
{
    EvtMRes * resonance = NULL;
    EvtMLineShape * lineshape = NULL;

    if( ls=="BREITWIGNER" ) {
        lineshape = new EvtMBreitWigner( id, lsarg );
    } else if( ls=="TRIVIAL" ) {
        lineshape = new EvtMTrivialLS( id, lsarg );
    } else {
        report(ERROR,"EvtGen")<<"Lineshape "<<lineshape
                              <<" not recognized."<<endl;
        ::abort();
    }

    if( type=="HELAMP" ) {
        resonance = new EvtMHelAmp( id, lineshape, children, amps );
    } else {
        report(ERROR,"EvtGen")<<"Model "<<type<<" not recognized."<<endl;
        ::abort();
    }

    lineshape->setres( resonance );

    return resonance;
}

void EvtMTree::parseerror( bool flag, ptype& c_iter, ptype& c_begin, 
        ptype& c_end )
{ 
    if(!flag) return;

    string error;
    
    while( c_begin != c_end ) {
        if(c_begin == c_iter) {
            error+='_';
            error+=*c_begin;
            error+='_';
        } else 
            error+=*c_begin;

        ++c_begin;
    }

    report(ERROR,"EvtGen")<<"Parse error at: "<<error<<endl;
    ::abort();
}

string EvtMTree::parseId( ptype& c_iter, ptype& c_begin, ptype& c_end ) 
{
    string strid;

    while(c_iter != c_end) {
        parseerror(parsecheck(*c_iter, ")[],"), c_iter, c_begin, c_end);
        if( *c_iter == '(' ) {
            ++c_iter;
            return strid;
        }

        strid += *c_iter;
        ++c_iter;
    }

    return strid;
}

string EvtMTree::parseKey( ptype& c_iter, ptype& c_begin, ptype& c_end ) 
{
    string key;

    while( *c_iter != ',' ) {
        parseerror(c_iter==c_end || parsecheck(*c_iter, "()[]"),
            c_iter, c_begin, c_end);
        key += *c_iter;
        ++c_iter;
    }

    ++c_iter;

    parseerror(c_iter == c_end, c_iter, c_begin, c_end);
    
    return key;
}

vector<string> EvtMTree::parseArg( ptype &c_iter, ptype &c_begin, ptype &c_end )
{
    vector<string> arg;

    if( *c_iter != '[' ) return arg;
    ++c_iter;

    string temp;
    while(true) {
        parseerror( c_iter == c_end || parsecheck(*c_iter, "[()"),
                c_iter, c_begin, c_end );

        if( *c_iter == ']' ) {
            ++c_iter;
            if(temp.size() > 0) arg.push_back( temp );
            break;
        }

        if( *c_iter == ',') {
            arg.push_back( temp );
            temp.clear();
            ++c_iter;
            continue;
        }

        temp += *c_iter;
        ++c_iter;
    }
    parseerror(c_iter == c_end || *c_iter != ',', c_iter, c_begin, c_end);
    ++c_iter;

    return arg;
}

vector<EvtComplex> EvtMTree::parseAmps( ptype &c_iter, 
        ptype &c_begin, ptype &c_end )
{
    vector<string> parg = parseArg( c_iter, c_begin, c_end );
    parseerror( parg.size() == 0, c_iter, c_begin, c_end );

    // Get parametrization amplitudes
    vector<string>::iterator amp_iter = parg.begin();
    vector<string>::iterator amp_end = parg.end();
    vector<EvtComplex> amps;

    while( amp_iter != amp_end ) {
        const char * nptr;
        char * endptr = NULL;
        double amp=0.0, phase=0.0;

        nptr = (*amp_iter).c_str();
        amp = strtod(nptr, &endptr);
        parseerror( nptr==endptr, c_iter, c_begin, c_end );

        ++amp_iter;
        parseerror( amp_iter == amp_end, c_iter, c_begin, c_end );

        nptr = (*amp_iter).c_str();
        phase = strtod(nptr, &endptr);
        parseerror( nptr==endptr, c_iter, c_begin, c_end );

        amps.push_back( amp*exp(EvtComplex(0.0, phase)) );

        ++amp_iter;
    }

    return amps;
}

vector<EvtMNode *> EvtMTree::duplicate( const vector<EvtMNode *>& list ) const
{
    vector<EvtMNode *> newlist;

    for(size_t i=0; i<list.size(); ++i)
        newlist.push_back( list[i]->duplicate() );

    return newlist;
}

// XXX Warning it is unsafe to use cl1 after a call to this function XXX
vector< vector<EvtMNode * > > EvtMTree::unionChildren( const string& nodestr,
        vector< vector<EvtMNode * > >& cl1 ) 
{
    vector<EvtMNode *> cl2 = parsenode( nodestr, false );
    vector< vector<EvtMNode * > > cl;
    
    if( cl1.size() == 0 ) {
        for( size_t i=0; i<cl2.size(); ++i ) {
            vector<EvtMNode *> temp(1, cl2[i]);
            cl.push_back( temp );
        }

        return cl;
    }

    for( size_t i=0; i<cl1.size(); ++i ) {
        for( size_t j=0; j<cl2.size(); ++j ) {
            vector<EvtMNode *> temp;
            temp = duplicate( cl1[i] );
            temp.push_back( cl2[j]->duplicate() );

            cl.push_back( temp );
        }
    }
 
    for(size_t i=0; i<cl1.size(); ++i)
        for(size_t j=0; j<cl1[i].size(); ++j)
            delete cl1[i][j];
    for(size_t i=0; i<cl2.size(); ++i)
        delete (cl2[i]);

    return cl;
}

vector< vector<EvtMNode * > > EvtMTree::parseChildren( ptype &c_iter, 
        ptype &c_begin, ptype &c_end ) 
{
    bool test = true;
    int pcount=0;
    string nodestr;
    vector< vector<EvtMNode * > > children;

    parseerror(c_iter == c_end || *c_iter != '[', c_iter, c_begin, c_end );
    ++c_iter;

    while( test ) {
        parseerror( c_iter==c_end || pcount < 0, c_iter, c_begin, c_end );

        switch( *c_iter ) {
            case ')':
                --pcount;
                nodestr += *c_iter;
                break;
            case '(':
                ++pcount;
                nodestr += *c_iter;
                break;
            case ']':
                if( pcount==0 ) {
                    children = unionChildren( nodestr, children );
                    test=false;
                } else {
                    nodestr += *c_iter;
                }
                break;
            case ',':
                if( pcount==0 ) {
                    children = unionChildren( nodestr, children );
                    nodestr.clear();
                } else {
                    nodestr += *c_iter;
                }
                break;
            default:
                nodestr += *c_iter;
                break;
        }

        ++c_iter;
    }

    return children;
}
    
vector<EvtMNode *> EvtMTree::parsenode( const string& args, bool rootnode )
{
    ptype c_iter, c_begin, c_end;

    c_iter=c_begin=args.begin();
    c_end = args.end();

    string strid = parseId( c_iter, c_begin, c_end );

    // Case 1: Particle
    if( c_iter == c_end ) return makeparticles( strid );

    // Case 2: Resonance - parse further
    EvtId id = EvtPDL::getId(strid);
    parseerror(EvtId( -1, -1 )==id, c_iter, c_begin, c_end);
    
    string ls;
    vector<string> lsarg;

    if( rootnode ) {
        ls = "TRIVIAL";
    } else {
        // Get lineshape (e.g. BREITWIGNER)
        ls = parseKey( c_iter, c_begin, c_end );
        lsarg = parseArg( c_iter, c_begin, c_end );
    }

    // Get resonance parametrization type (e.g. HELAMP)
    string type = parseKey( c_iter, c_begin, c_end );
    vector<EvtComplex> amps = parseAmps( c_iter, c_begin, c_end );

    // Children
    vector<vector<EvtMNode * > > children = parseChildren( c_iter, c_begin,
            c_end );

    report(ERROR,"EvtGen")<<children.size()<<endl;
    vector<EvtMNode *> resonances;
    for(size_t i=0; i<children.size(); ++i ) {
        resonances.push_back(makeresonance(id,ls,lsarg,type,amps,children[i]));
    }

    parseerror(c_iter == c_end || *c_iter!=')', c_iter, c_begin, c_end);

    return resonances;
}

bool EvtMTree::validTree( const EvtMNode * root ) const
{
    vector<int> res = root->getresonance();
    vector<bool> check(res.size(), false);

    for( size_t i=0; i<res.size(); ++i) {
        check[res[i]] = true;
    }

    bool ret = true;

    for( size_t i=0; i<check.size(); ++i ) {
        ret = ret&&check[i];
    }

    return ret;
}

void EvtMTree::addtree( const string& str )
{
    vector<EvtMNode *> roots = parsenode( str, true );
    _norm = 0;

    for( size_t i=0; i<roots.size(); ++i ) {
        if( validTree( roots[i] ) ) {
            _root.push_back( roots[i] );
            _norm = _norm + 1;
        } else
            delete roots[i];
    }
    
    _norm = 1.0/sqrt(_norm);
}
EvtSpinAmp EvtMTree::getrotation( EvtParticle * p ) const
{
    // Set up the rotation matrix for the root particle (for now)
    EvtSpinDensity sd = p->rotateToHelicityBasis();
    EvtSpinType::spintype type = EvtPDL::getSpinType(_root[0]->getid());
    int twospin = EvtSpinType::getSpin2(type);
    
    vector<EvtSpinType::spintype> types(2, type);
    EvtSpinAmp rot( types, EvtComplex(0.0, 0.0) );
    vector<int> index = rot.iterallowedinit();
    do {
        rot(index) = sd.get((index[0]+twospin)/2,(index[1]+twospin)/2);
    } while( rot.iterateallowed( index ) );

    return rot;
}

EvtSpinAmp EvtMTree::amplitude( EvtParticle * p ) const
{
    vector<EvtVector4R> product;
    for(size_t i=0; i<p->getNDaug(); ++i)
        product.push_back(p->getDaug(i)->getP4Lab());

    if( _root.size() == 0 ) {
        report(ERROR, "EvtGen")<<"No decay tree present."<<endl;
        ::abort();
    }
    
    EvtSpinAmp amp = _root[0]->amplitude( product );
    for( size_t i=1; i<_root.size(); ++i ) {
        // Assume that helicity amplitude is returned 
        amp += _root[i]->amplitude( product );
    }
    amp = _norm*amp;

    //ryd
    return amp;

    // Do Rotation to Proper Frame
    EvtSpinAmp newamp = getrotation( p );
    newamp.extcont(amp, 1, 0);

    return newamp;
}
Commit	Line	Data
da0e9ce3	1	#include <stdio.h>
	2	#include <stdlib.h>
	3	#include <algorithm>
	4
	5	#include "EvtGenBase/EvtParticle.hh"
	6
	7	#include "EvtGenBase/EvtMTree.hh"
	8	#include "EvtGenBase/EvtConst.hh"
	9	#include "EvtGenBase/EvtKine.hh"
	10	#include "EvtGenBase/EvtReport.hh"
	11
	12	// Make sure to include Lineshapes here
	13	#include "EvtGenBase/EvtMTrivialLS.hh"
	14	#include "EvtGenBase/EvtMBreitWigner.hh"
	15
	16	// Make sure to include Parametrizations
	17	#include "EvtGenBase/EvtMHelAmp.hh"
	18
	19	using std::endl;
	20
	21	EvtMTree::EvtMTree( const EvtId * idtbl, unsigned int ndaug )
	22	{
	23	for( size_t i=0; i<ndaug; ++i ) {
	24	_lbltbl.push_back( EvtPDL::name( idtbl[i] ) );
	25	}
	26	}
	27
	28	EvtMTree::~EvtMTree()
	29	{
	30	for(size_t i=0; i<_root.size(); ++i) delete _root[i];
	31	}
	32
	33	bool EvtMTree::parsecheck( char arg, const string& chars )
	34	{
	35	bool ret = false;
	36
	37	for(size_t i=0; i<chars.size(); ++i) {
	38	ret = ret \|\| (chars[i]==arg);
	39	}
	40
	41	return ret;
	42	}
	43
	44	vector<EvtMNode *> EvtMTree::makeparticles( const string& strid )
	45	{
	46	vector<EvtMNode *> particles;
	47	vector<int> labels;
	48
	49	for( size_t i = 0; i<_lbltbl.size(); ++i ) {
	50	if( _lbltbl[i] == strid ) labels.push_back( i );
	51	}
	52
	53	if( labels.size() == 0 ) {
	54	report(ERROR,"EvtGen")<<"Error unknown particle label "<<strid<<endl;
	55	::abort();
	56	}
	57
	58	for( size_t i = 0; i<labels.size(); ++i )
	59	particles.push_back( new EvtMParticle( labels[i], EvtPDL::getId( strid ) ) );
	60
	61	return particles;
	62	}
	63
	64	EvtMRes * EvtMTree::makeresonance( const EvtId& id, const string& ls,
65	const vector<string>& lsarg, const string& type,
66	const vector<EvtComplex>& amps, const vector<EvtMNode *>& children )
67	{
68	EvtMRes * resonance = NULL;
69	EvtMLineShape * lineshape = NULL;
70
71	if( ls=="BREITWIGNER" ) {
72	lineshape = new EvtMBreitWigner( id, lsarg );
73	} else if( ls=="TRIVIAL" ) {
74	lineshape = new EvtMTrivialLS( id, lsarg );
75	} else {
76	report(ERROR,"EvtGen")<<"Lineshape "<<lineshape
77	<<" not recognized."<<endl;
78	::abort();
79	}
80
81	if( type=="HELAMP" ) {
82	resonance = new EvtMHelAmp( id, lineshape, children, amps );
83	} else {
84	report(ERROR,"EvtGen")<<"Model "<<type<<" not recognized."<<endl;
85	::abort();
86	}
87
88	lineshape->setres( resonance );
89
90	return resonance;
91	}
92
93	void EvtMTree::parseerror( bool flag, ptype& c_iter, ptype& c_begin,
94	ptype& c_end )
95	{
96	if(!flag) return;
97
98	string error;
99
100	while( c_begin != c_end ) {
101	if(c_begin == c_iter) {
102	error+='_';
103	error+=*c_begin;
104	error+='_';
105	} else
106	error+=*c_begin;
107
108	++c_begin;
109	}
110
111	report(ERROR,"EvtGen")<<"Parse error at: "<<error<<endl;
112	::abort();
113	}
114
115	string EvtMTree::parseId( ptype& c_iter, ptype& c_begin, ptype& c_end )
116	{
117	string strid;
118
119	while(c_iter != c_end) {
120	parseerror(parsecheck(*c_iter, ")[],"), c_iter, c_begin, c_end);
121	if( *c_iter == '(' ) {
122	++c_iter;
123	return strid;
124	}
125
126	strid += *c_iter;
127	++c_iter;
128	}
129
130	return strid;
131	}
132
133	string EvtMTree::parseKey( ptype& c_iter, ptype& c_begin, ptype& c_end )
134	{
135	string key;
136
137	while( *c_iter != ',' ) {
138	parseerror(c_iter==c_end \|\| parsecheck(*c_iter, "()[]"),
139	c_iter, c_begin, c_end);
140	key += *c_iter;
141	++c_iter;
142	}
143
144	++c_iter;
145
146	parseerror(c_iter == c_end, c_iter, c_begin, c_end);
147
148	return key;
149	}
150
151	vector<string> EvtMTree::parseArg( ptype &c_iter, ptype &c_begin, ptype &c_end )
152	{
153	vector<string> arg;
154
155	if( *c_iter != '[' ) return arg;
156	++c_iter;
157
158	string temp;
159	while(true) {
160	parseerror( c_iter == c_end \|\| parsecheck(*c_iter, "[()"),
161	c_iter, c_begin, c_end );
162
163	if( *c_iter == ']' ) {
164	++c_iter;
165	if(temp.size() > 0) arg.push_back( temp );
166	break;
167	}
168
169	if( *c_iter == ',') {
170	arg.push_back( temp );
171	temp.clear();
172	++c_iter;
173	continue;
174	}
175
176	temp += *c_iter;
177	++c_iter;
178	}
179	parseerror(c_iter == c_end \|\| *c_iter != ',', c_iter, c_begin, c_end);
180	++c_iter;
181
182	return arg;
183	}
184
185	vector<EvtComplex> EvtMTree::parseAmps( ptype &c_iter,
186	ptype &c_begin, ptype &c_end )
187	{
188	vector<string> parg = parseArg( c_iter, c_begin, c_end );
189	parseerror( parg.size() == 0, c_iter, c_begin, c_end );
190
191	// Get parametrization amplitudes
192	vector<string>::iterator amp_iter = parg.begin();
193	vector<string>::iterator amp_end = parg.end();
194	vector<EvtComplex> amps;
195
196	while( amp_iter != amp_end ) {
197	const char * nptr;
198	char * endptr = NULL;
199	double amp=0.0, phase=0.0;
200
201	nptr = (*amp_iter).c_str();
202	amp = strtod(nptr, &endptr);
203	parseerror( nptr==endptr, c_iter, c_begin, c_end );
204
205	++amp_iter;
206	parseerror( amp_iter == amp_end, c_iter, c_begin, c_end );
207
208	nptr = (*amp_iter).c_str();
209	phase = strtod(nptr, &endptr);
210	parseerror( nptr==endptr, c_iter, c_begin, c_end );
211
212	amps.push_back( amp*exp(EvtComplex(0.0, phase)) );
213
214	++amp_iter;
215	}
216
217	return amps;
218	}
219
220	vector<EvtMNode > EvtMTree::duplicate( const vector<EvtMNode >& list ) const
221	{
222	vector<EvtMNode *> newlist;
223
224	for(size_t i=0; i<list.size(); ++i)
225	newlist.push_back( list[i]->duplicate() );
226
227	return newlist;
228	}
229
230	// XXX Warning it is unsafe to use cl1 after a call to this function XXX
231	vector< vector<EvtMNode * > > EvtMTree::unionChildren( const string& nodestr,
232	vector< vector<EvtMNode * > >& cl1 )
233	{
234	vector<EvtMNode *> cl2 = parsenode( nodestr, false );
235	vector< vector<EvtMNode * > > cl;
236
237	if( cl1.size() == 0 ) {
238	for( size_t i=0; i<cl2.size(); ++i ) {
239	vector<EvtMNode *> temp(1, cl2[i]);
240	cl.push_back( temp );
241	}
242
243	return cl;
244	}
245
246	for( size_t i=0; i<cl1.size(); ++i ) {
247	for( size_t j=0; j<cl2.size(); ++j ) {
248	vector<EvtMNode *> temp;
249	temp = duplicate( cl1[i] );
250	temp.push_back( cl2[j]->duplicate() );
251
252	cl.push_back( temp );
253	}
254	}
255
256	for(size_t i=0; i<cl1.size(); ++i)
257	for(size_t j=0; j<cl1[i].size(); ++j)
258	delete cl1[i][j];
259	for(size_t i=0; i<cl2.size(); ++i)
260	delete (cl2[i]);
261
262	return cl;
263	}
264
265	vector< vector<EvtMNode * > > EvtMTree::parseChildren( ptype &c_iter,
266	ptype &c_begin, ptype &c_end )
267	{
268	bool test = true;
269	int pcount=0;
270	string nodestr;
271	vector< vector<EvtMNode * > > children;
272
273	parseerror(c_iter == c_end \|\| *c_iter != '[', c_iter, c_begin, c_end );
274	++c_iter;
275
276	while( test ) {
277	parseerror( c_iter==c_end \|\| pcount < 0, c_iter, c_begin, c_end );
278
279	switch( *c_iter ) {
280	case ')':
281	--pcount;
282	nodestr += *c_iter;
283	break;
284	case '(':
285	++pcount;
286	nodestr += *c_iter;
287	break;
288	case ']':
289	if( pcount==0 ) {
290	children = unionChildren( nodestr, children );
291	test=false;
292	} else {
293	nodestr += *c_iter;
294	}
295	break;
296	case ',':
297	if( pcount==0 ) {
298	children = unionChildren( nodestr, children );
299	nodestr.clear();
300	} else {
301	nodestr += *c_iter;
302	}
303	break;
304	default:
305	nodestr += *c_iter;
306	break;
307	}
308
309	++c_iter;
310	}
311
312	return children;
313	}
314
315	vector<EvtMNode *> EvtMTree::parsenode( const string& args, bool rootnode )
316	{
317	ptype c_iter, c_begin, c_end;
318
319	c_iter=c_begin=args.begin();
320	c_end = args.end();
321
322	string strid = parseId( c_iter, c_begin, c_end );
323
324	// Case 1: Particle
325	if( c_iter == c_end ) return makeparticles( strid );
326
327	// Case 2: Resonance - parse further
328	EvtId id = EvtPDL::getId(strid);
329	parseerror(EvtId( -1, -1 )==id, c_iter, c_begin, c_end);
330
331	string ls;
332	vector<string> lsarg;
333
334	if( rootnode ) {
335	ls = "TRIVIAL";
336	} else {
337	// Get lineshape (e.g. BREITWIGNER)
338	ls = parseKey( c_iter, c_begin, c_end );
339	lsarg = parseArg( c_iter, c_begin, c_end );
340	}
341
342	// Get resonance parametrization type (e.g. HELAMP)
343	string type = parseKey( c_iter, c_begin, c_end );
344	vector<EvtComplex> amps = parseAmps( c_iter, c_begin, c_end );
345
346	// Children
347	vector<vector<EvtMNode * > > children = parseChildren( c_iter, c_begin,
348	c_end );
349
350	report(ERROR,"EvtGen")<<children.size()<<endl;
351	vector<EvtMNode *> resonances;
352	for(size_t i=0; i<children.size(); ++i ) {
353	resonances.push_back(makeresonance(id,ls,lsarg,type,amps,children[i]));
354	}
355
356	parseerror(c_iter == c_end \|\| *c_iter!=')', c_iter, c_begin, c_end);
357
358	return resonances;
359	}
360
361	bool EvtMTree::validTree( const EvtMNode * root ) const
362	{
363	vector<int> res = root->getresonance();
364	vector<bool> check(res.size(), false);
365
366	for( size_t i=0; i<res.size(); ++i) {
367	check[res[i]] = true;
368	}
369
370	bool ret = true;
371
372	for( size_t i=0; i<check.size(); ++i ) {
373	ret = ret&&check[i];
374	}
375
376	return ret;
377	}
378
379	void EvtMTree::addtree( const string& str )
380	{
381	vector<EvtMNode *> roots = parsenode( str, true );
382	_norm = 0;
383
384	for( size_t i=0; i<roots.size(); ++i ) {
385	if( validTree( roots[i] ) ) {
386	_root.push_back( roots[i] );
387	_norm = _norm + 1;
388	} else
389	delete roots[i];
390	}
391
392	_norm = 1.0/sqrt(_norm);
393	}
394	EvtSpinAmp EvtMTree::getrotation( EvtParticle * p ) const
395	{
396	// Set up the rotation matrix for the root particle (for now)
397	EvtSpinDensity sd = p->rotateToHelicityBasis();
398	EvtSpinType::spintype type = EvtPDL::getSpinType(_root[0]->getid());
399	int twospin = EvtSpinType::getSpin2(type);
400
401	vector<EvtSpinType::spintype> types(2, type);
402	EvtSpinAmp rot( types, EvtComplex(0.0, 0.0) );
403	vector<int> index = rot.iterallowedinit();
404	do {
405	rot(index) = sd.get((index[0]+twospin)/2,(index[1]+twospin)/2);
406	} while( rot.iterateallowed( index ) );
407
408	return rot;
409	}
410
411	EvtSpinAmp EvtMTree::amplitude( EvtParticle * p ) const
412	{
413	vector<EvtVector4R> product;
414	for(size_t i=0; i<p->getNDaug(); ++i)
415	product.push_back(p->getDaug(i)->getP4Lab());
416
417	if( _root.size() == 0 ) {
418	report(ERROR, "EvtGen")<<"No decay tree present."<<endl;
419	::abort();
420	}
421
422	EvtSpinAmp amp = _root[0]->amplitude( product );
423	for( size_t i=1; i<_root.size(); ++i ) {
424	// Assume that helicity amplitude is returned
425	amp += _root[i]->amplitude( product );
426	}
427	amp = _norm*amp;
428
429	//ryd
430	return amp;
431
432	// Do Rotation to Proper Frame
433	EvtSpinAmp newamp = getrotation( p );
434	newamp.extcont(amp, 1, 0);
435
436	return newamp;
437	}