Fix for definitions for CINT
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenBase / EvtMTree.cpp
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 }