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