00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <algorithm>
00004 #include <string>
00005
00006 #include "EvtGenBase/EvtMTree.hh"
00007 #include "EvtGenBase/EvtConst.hh"
00008 #include "EvtGenBase/EvtKine.hh"
00009 #include "EvtGenBase/EvtReport.hh"
00010
00011
00012 #include "EvtGenBase/EvtMTrivialLS.hh"
00013 #include "EvtGenBase/EvtMBreitWigner.hh"
00014
00015
00016 #include "EvtGenBase/EvtMHelAmp.hh"
00017
00018 using std::endl;
00019
00020 EvtMTree::EvtMTree( const EvtId * idtbl, int ndaug )
00021 {
00022 for( int i=0; i<ndaug; ++i ) {
00023 _lbltbl.push_back( EvtPDL::name( idtbl[i] ) );
00024 }
00025 }
00026
00027 EvtMTree::~EvtMTree()
00028 {
00029 for(int i=0; i<_root.size(); ++i) delete _root[i];
00030 }
00031
00032 bool EvtMTree::parsecheck( char arg, const string& chars )
00033 {
00034 bool ret = false;
00035
00036 for(int i=0; i<chars.size(); ++i) {
00037 ret = ret || (chars[i]==arg);
00038 }
00039
00040 return ret;
00041 }
00042
00043 vector<EvtMNode *> EvtMTree::makeparticles( const string& strid )
00044 {
00045 vector<EvtMNode *> particles;
00046 vector<int> labels;
00047
00048 for( int i = 0; i<_lbltbl.size(); ++i ) {
00049 if( _lbltbl[i] == strid ) labels.push_back( i );
00050 }
00051
00052 if( labels.size() == 0 ) {
00053 report(ERROR,"EvtGen")<<"Error unknown particle label "<<strid<<endl;
00054 ::abort();
00055 }
00056
00057 for( int i = 0; i<labels.size(); ++i )
00058 particles.push_back( new EvtMParticle( labels[i], EvtPDL::getId( strid ) ) );
00059
00060 return particles;
00061 }
00062
00063 EvtMRes * EvtMTree::makeresonance( const EvtId& id, const string& ls,
00064 const vector<string>& lsarg, const string& type,
00065 const vector<EvtComplex>& amps, const vector<EvtMNode *>& children )
00066 {
00067 EvtMRes * resonance = NULL;
00068 EvtMLineShape * lineshape = NULL;
00069
00070 if( ls=="BREITWIGNER" ) {
00071 lineshape = new EvtMBreitWigner( id, lsarg );
00072 } else if( ls=="TRIVIAL" ) {
00073 lineshape = new EvtMTrivialLS( id, lsarg );
00074 } else {
00075 report(ERROR,"EvtGen")<<"Lineshape "<<lineshape
00076 <<" not recognized."<<endl;
00077 ::abort();
00078 }
00079
00080 if( type=="HELAMP" ) {
00081 resonance = new EvtMHelAmp( id, lineshape, children, amps );
00082 } else {
00083 report(ERROR,"EvtGen")<<"Model "<<type<<" not recognized."<<endl;
00084 ::abort();
00085 }
00086
00087 lineshape->setres( resonance );
00088
00089 return resonance;
00090 }
00091
00092 void EvtMTree::parseerror( bool flag, ptype& c_iter, ptype& c_begin,
00093 ptype& c_end )
00094 {
00095 if(!flag) return;
00096
00097 string error;
00098
00099 while( c_begin != c_end ) {
00100 if(c_begin == c_iter) {
00101 error+='_';
00102 error+=*c_begin;
00103 error+='_';
00104 } else
00105 error+=*c_begin;
00106
00107 ++c_begin;
00108 }
00109
00110 report(ERROR,"EvtGen")<<"Parse error at: "<<error<<endl;
00111 ::abort();
00112 }
00113
00114 string EvtMTree::parseId( ptype& c_iter, ptype& c_begin, ptype& c_end )
00115 {
00116 string strid;
00117
00118 while(c_iter != c_end) {
00119 parseerror(parsecheck(*c_iter, ")[],"), c_iter, c_begin, c_end);
00120 if( *c_iter == '(' ) {
00121 ++c_iter;
00122 return strid;
00123 }
00124
00125 strid += *c_iter;
00126 ++c_iter;
00127 }
00128
00129 return strid;
00130 }
00131
00132 string EvtMTree::parseKey( ptype& c_iter, ptype& c_begin, ptype& c_end )
00133 {
00134 string key;
00135
00136 while( *c_iter != ',' ) {
00137 parseerror(c_iter==c_end || parsecheck(*c_iter, "()[]"),
00138 c_iter, c_begin, c_end);
00139 key += *c_iter;
00140 ++c_iter;
00141 }
00142
00143 ++c_iter;
00144
00145 parseerror(c_iter == c_end, c_iter, c_begin, c_end);
00146
00147 return key;
00148 }
00149
00150 vector<string> EvtMTree::parseArg( ptype &c_iter, ptype &c_begin, ptype &c_end )
00151 {
00152 vector<string> arg;
00153
00154 if( *c_iter != '[' ) return arg;
00155 ++c_iter;
00156
00157 string temp;
00158 while(true) {
00159 parseerror( c_iter == c_end || parsecheck(*c_iter, "[()"),
00160 c_iter, c_begin, c_end );
00161
00162 if( *c_iter == ']' ) {
00163 ++c_iter;
00164 if(temp.size() > 0) arg.push_back( temp );
00165 break;
00166 }
00167
00168 if( *c_iter == ',') {
00169 arg.push_back( temp );
00170 temp.erase();
00171 ++c_iter;
00172 continue;
00173 }
00174
00175 temp += *c_iter;
00176 ++c_iter;
00177 }
00178 parseerror(c_iter == c_end || *c_iter != ',', c_iter, c_begin, c_end);
00179 ++c_iter;
00180
00181 return arg;
00182 }
00183
00184 vector<EvtComplex> EvtMTree::parseAmps( ptype &c_iter,
00185 ptype &c_begin, ptype &c_end )
00186 {
00187 vector<string> parg = parseArg( c_iter, c_begin, c_end );
00188 parseerror( parg.size() == 0, c_iter, c_begin, c_end );
00189
00190
00191 vector<string>::iterator amp_iter = parg.begin();
00192 vector<string>::iterator amp_end = parg.end();
00193 vector<EvtComplex> amps;
00194
00195 while( amp_iter != amp_end ) {
00196 const char * nptr;
00197 char * endptr = NULL;
00198 double amp=0.0, phase=0.0;
00199
00200 nptr = (*amp_iter).c_str();
00201 amp = strtod(nptr, &endptr);
00202 parseerror( nptr==endptr, c_iter, c_begin, c_end );
00203
00204 ++amp_iter;
00205 parseerror( amp_iter == amp_end, c_iter, c_begin, c_end );
00206
00207 nptr = (*amp_iter).c_str();
00208 phase = strtod(nptr, &endptr);
00209 parseerror( nptr==endptr, c_iter, c_begin, c_end );
00210
00211 amps.push_back( amp*exp(EvtComplex(0.0, phase)) );
00212
00213 ++amp_iter;
00214 }
00215
00216 return amps;
00217 }
00218
00219 vector<EvtMNode *> EvtMTree::duplicate( const vector<EvtMNode *>& list ) const
00220 {
00221 vector<EvtMNode *> newlist;
00222
00223 for(int i=0; i<list.size(); ++i)
00224 newlist.push_back( list[i]->duplicate() );
00225
00226 return newlist;
00227 }
00228
00229
00230
00231 vector< vector<EvtMNode * > > EvtMTree::unionChildren( const string& nodestr,
00232 vector< vector<EvtMNode * > >& cl1 )
00233 {
00234 vector<EvtMNode *> cl2 = parsenode( nodestr, false );
00235 vector< vector<EvtMNode * > > cl;
00236
00237 if( cl1.size() == 0 ) {
00238 for( int i=0; i<cl2.size(); ++i ) {
00239 vector<EvtMNode *> temp(1, cl2[i]);
00240 cl.push_back( temp );
00241 }
00242
00243 return cl;
00244 }
00245
00246 for( int i=0; i<cl1.size(); ++i ) {
00247 for( int j=0; j<cl2.size(); ++j ) {
00248 vector<EvtMNode *> temp;
00249 temp = duplicate( cl1[i] );
00250 temp.push_back( cl2[j]->duplicate() );
00251
00252 cl.push_back( temp );
00253 }
00254 }
00255
00256 for(int i=0; i<cl1.size(); ++i)
00257 for(int j=0; j<cl1[i].size(); ++j)
00258 delete cl1[i][j];
00259 for(int i=0; i<cl2.size(); ++i)
00260 delete (cl2[i]);
00261
00262 return cl;
00263 }
00264
00265 vector< vector<EvtMNode * > > EvtMTree::parseChildren( ptype &c_iter,
00266 ptype &c_begin, ptype &c_end )
00267 {
00268 bool test = true;
00269 int pcount=0;
00270 string nodestr;
00271 vector< vector<EvtMNode * > > children;
00272
00273 parseerror(c_iter == c_end || *c_iter != '[', c_iter, c_begin, c_end );
00274 ++c_iter;
00275
00276 while( test ) {
00277 parseerror( c_iter==c_end || pcount < 0, c_iter, c_begin, c_end );
00278
00279 switch( *c_iter ) {
00280 case ')':
00281 --pcount;
00282 nodestr += *c_iter;
00283 break;
00284 case '(':
00285 ++pcount;
00286 nodestr += *c_iter;
00287 break;
00288 case ']':
00289 if( pcount==0 ) {
00290 children = unionChildren( nodestr, children );
00291 test=false;
00292 } else {
00293 nodestr += *c_iter;
00294 }
00295 break;
00296 case ',':
00297 if( pcount==0 ) {
00298 children = unionChildren( nodestr, children );
00299 nodestr.erase();
00300 } else {
00301 nodestr += *c_iter;
00302 }
00303 break;
00304 default:
00305 nodestr += *c_iter;
00306 break;
00307 }
00308
00309 ++c_iter;
00310 }
00311
00312 return children;
00313 }
00314
00315 vector<EvtMNode *> EvtMTree::parsenode( const string& args, bool rootnode )
00316 {
00317 ptype c_iter, c_begin, c_end;
00318
00319 c_iter=c_begin=args.begin();
00320 c_end = args.end();
00321
00322 string strid = parseId( c_iter, c_begin, c_end );
00323
00324
00325 if( c_iter == c_end ) return makeparticles( strid );
00326
00327
00328 EvtId id = EvtPDL::getId(strid);
00329 parseerror(EvtId( -1, -1 )==id, c_iter, c_begin, c_end);
00330
00331 string ls;
00332 vector<string> lsarg;
00333
00334 if( rootnode ) {
00335 ls = "TRIVIAL";
00336 } else {
00337
00338 ls = parseKey( c_iter, c_begin, c_end );
00339 lsarg = parseArg( c_iter, c_begin, c_end );
00340 }
00341
00342
00343 string type = parseKey( c_iter, c_begin, c_end );
00344 vector<EvtComplex> amps = parseAmps( c_iter, c_begin, c_end );
00345
00346
00347 vector<vector<EvtMNode * > > children = parseChildren( c_iter, c_begin,
00348 c_end );
00349
00350 report(ERROR,"EvtGen")<<children.size()<<endl;
00351 vector<EvtMNode *> resonances;
00352 for(int i=0; i<children.size(); ++i ) {
00353 resonances.push_back(makeresonance(id,ls,lsarg,type,amps,children[i]));
00354 }
00355
00356 parseerror(c_iter == c_end || *c_iter!=')', c_iter, c_begin, c_end);
00357
00358 return resonances;
00359 }
00360
00361 bool EvtMTree::validTree( const EvtMNode * root ) const
00362 {
00363 vector<int> res = root->getresonance();
00364 vector<bool> check(res.size(), false);
00365
00366 for( int i=0; i<res.size(); ++i) {
00367 check[res[i]] = true;
00368 }
00369
00370 bool ret = true;
00371
00372 for( int i=0; i<check.size(); ++i ) {
00373 ret = ret&&check[i];
00374 }
00375
00376 return ret;
00377 }
00378
00379 void EvtMTree::addtree( const string& str )
00380 {
00381 vector<EvtMNode *> roots = parsenode( str, true );
00382 _norm = 0;
00383
00384 for( int i=0; i<roots.size(); ++i ) {
00385 if( validTree( roots[i] ) ) {
00386 _root.push_back( roots[i] );
00387 _norm = _norm + 1;
00388 } else
00389 delete roots[i];
00390 }
00391
00392 _norm = 1.0/sqrt(_norm);
00393 }
00394
00395 EvtSpinAmp EvtMTree::amplitude( const vector<EvtVector4R>&
00396 product) const
00397 {
00398 if( _root.size() == 0 ) {
00399 report(ERROR, "EvtGen")<<"No decay tree present."<<endl;
00400 ::abort();
00401 }
00402
00403 EvtSpinAmp amp = _root[0]->amplitude( product );
00404 for( int i=1; i<_root.size(); ++i ) {
00405
00406
00407
00408 amp += _root[i]->amplitude( product );
00409 }
00410
00411 return _norm*amp;
00412 }