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 }