TEvtGen/EvtGenBase/EvtSpinAmp.cpp

   1 #include "EvtGenBase/EvtPatches.hh"
   2 #include "EvtGenBase/EvtReport.hh"
   3 #include "EvtGenBase/EvtSpinAmp.hh"
   4 #include <stdlib.h>
   5
   6 using std::endl;
   7
   8 std::ostream&
   9 operator<<( std::ostream& os, const EvtSpinAmp& amp )
  10 {
  11     vector<int> index = amp.iterinit();
  12
  13     os << ":";
  14     do {
  15         os<<"<";
  16         for(size_t i=0; i<index.size()-1; ++i) {
  17             os<<index[i];
  18         }
  19         os<<index[index.size()-1]<<">"<<amp(index)<<":";
  20     } while( amp.iterate( index ) );
  21
  22     return os;
  23 }
  24
  25 EvtSpinAmp operator*( const EvtComplex& real, const EvtSpinAmp& cont )
  26 {
  27     EvtSpinAmp ret( cont );
  28
  29     for( size_t i=0; i<ret._elem.size(); ++i ) {
  30         ret._elem[i] *= real;
  31     }
  32
  33     return ret;
  34 }
  35
  36 EvtSpinAmp operator*( const EvtSpinAmp& cont, const EvtComplex& real )
  37 {
  38     return real*cont;
  39 }
  40
  41 EvtSpinAmp operator/( const EvtSpinAmp& cont, const EvtComplex& real )
  42 {
  43     EvtSpinAmp ret( cont );
  44
  45     for( size_t i=0; i<ret._elem.size(); ++i ) {
  46         ret._elem[i] /= real;
  47     }
  48
  49     return ret;
  50 }
  51
  52 vector<unsigned int> EvtSpinAmp::calctwospin( const vector<EvtSpinType::spintype>& type  ) const
  53 {
  54     vector<unsigned int> twospin;
  55
  56     for( size_t i=0; i<type.size(); ++i ) {
  57         twospin.push_back( EvtSpinType::getSpin2( type[i] ) );
  58     }
  59
  60     return twospin;
  61 }
  62
  63 EvtSpinAmp::EvtSpinAmp( const vector<EvtSpinType::spintype>& type )
  64 {
  65     int num = 1;
  66     _type = type;
  67     _twospin=calctwospin( type );
  68
  69     for( size_t i=0; i<_twospin.size(); ++i )
  70         num*=_twospin[i]+1;
  71
  72     _elem=vector<EvtComplex>( num );
  73 }
  74
  75 EvtSpinAmp::EvtSpinAmp( const vector<EvtSpinType::spintype>& type, const EvtComplex & val )
  76 {
  77     int num = 1;
  78     _type = type;
  79     _twospin=calctwospin( type );
  80
  81     for( size_t i=0; i<_twospin.size(); ++i )
  82         num*=_twospin[i]+1;
  83
  84     _elem=vector<EvtComplex>( num, val );
  85 }
  86
  87 EvtSpinAmp::EvtSpinAmp( const vector<EvtSpinType::spintype>& type,
  88         const vector<EvtComplex>& elem )
  89 {
  90     size_t num = 1;
  91
  92     _type = type;
  93     _twospin=calctwospin( type );
  94     _elem=elem;
  95
  96     for( size_t i=0; i<_twospin.size(); ++i ){
  97         num*=(_twospin[i]+1);
  98     }
  99
 100     if(_elem.size() != num ) {
 101         report(ERROR,"EvtGen")<<"Wrong number of elements input:"
 102             <<_elem.size()<<" vs. "<<num<<endl;
 103        ::abort();
 104     }
 105
 106 }
 107
 108 EvtSpinAmp::EvtSpinAmp( const EvtSpinAmp& copy )
 109 {
 110     _twospin = copy._twospin;
 111     _elem = copy._elem;
 112     _type = copy._type;
 113 }
 114
 115 void EvtSpinAmp::checktwospin( const vector<unsigned int>& twospin ) const
 116 {
 117     if( _twospin == twospin )
 118         return;
 119
 120     report( ERROR, "EvtGen" )
 121         <<"Dimension or order of tensors being operated on does not match"
 122         <<endl;
 123     ::abort();
 124 }
 125
 126 void EvtSpinAmp::checkindexargs( const vector<int>& index ) const
 127 {
 128     if( index.size()==0 ) {
 129         report(ERROR,"EvtGen") << "EvtSpinAmp can't handle no indices" << endl;
 130         ::abort();
 131     }
 132
 133     if( index.size() != _twospin.size() ) {
 134         report( ERROR, "EvtGen" ) << "Rank of EvtSpinAmp index does not match: "
 135             <<_twospin.size()<<" expected "<<index.size()<<" input."<<endl;
 136         ::abort();
 137     }
 138
 139     for( size_t i=0; i<_twospin.size(); ++i ) {
 140         if( static_cast<int>(_twospin[i])>=abs(index[i]) && static_cast<int>(_twospin[i])%2==index[i]%2 )
 141             continue;
 142         report(ERROR,"EvtGen")<<"EvtSpinAmp index out of range" << endl;
 143         report(ERROR,"EvtGen")<<" Index: ";
 144         for(size_t j=0; j<_twospin.size(); ++j )
 145             report(ERROR," ")<<_twospin[j];
 146
 147         report(ERROR, " ")<<endl<<" Index: ";
 148         for(size_t j=0; j<index.size(); ++j )
 149             report(ERROR," ")<<index[j];
 150         report(ERROR, " ")<<endl;
 151         ::abort();
 152     }
 153 }
 154
 155 int EvtSpinAmp::findtrueindex( const vector<int>& index ) const
 156 {
 157     int trueindex = 0;
 158
 159     for( size_t i = index.size()-1; i>0; --i ) {
 160         trueindex += (index[i]+_twospin[i])/2;
 161         trueindex *= _twospin[i-1]+1;
 162     }
 163
 164     trueindex += (index[0]+_twospin[0])/2;
 165
 166     return trueindex;
 167 }
 168
 169 EvtComplex & EvtSpinAmp::operator()( const vector<int>& index )
 170 {
 171     checkindexargs( index );
 172
 173     size_t trueindex = findtrueindex(index);
 174     if(trueindex >= _elem.size()) {
 175         report(ERROR,"EvtGen")<<"indexing error "<<trueindex<<" "<<_elem.size()<<endl;
 176         for(size_t i=0; i<_twospin.size(); ++i) {
 177             report(ERROR,"")<<_twospin[i]<<" ";
 178         }
 179         report(ERROR,"")<<endl;
 180
 181         for(size_t i=0; i<index.size(); ++i) {
 182             report(ERROR,"")<<index[i]<<" ";
 183         }
 184         report(ERROR,"")<<endl;
 185
 186         ::abort();
 187     }
 188
 189     return _elem[trueindex];
 190 }
 191
 192 const EvtComplex & EvtSpinAmp::operator()( const vector<int>& index ) const
 193 {
 194     checkindexargs( index );
 195
 196     size_t trueindex = findtrueindex(index);
 197     if(trueindex >= _elem.size()) {
 198         report(ERROR,"EvtGen")<<"indexing error "<<trueindex<<" "<<_elem.size()<<endl;
 199         for(size_t i=0; i<_twospin.size(); ++i) {
 200             report(ERROR,"")<<_twospin[i]<<" ";
 201         }
 202         report(ERROR,"")<<endl;
 203
 204         for(size_t i=0; i<index.size(); ++i) {
 205             report(ERROR,"")<<index[i]<<" ";
 206         }
 207         report(ERROR,"")<<endl;
 208
 209         ::abort();
 210     }
 211
 212     return _elem[trueindex];
 213 }
 214
 215 EvtComplex & EvtSpinAmp::operator()( int i, ... )
 216 {
 217     va_list ap;
 218     vector<int> index( _twospin.size() );
 219
 220     va_start(ap, i);
 221
 222     index[0]=i;
 223     for(size_t n=1; n<_twospin.size(); ++n )
 224         index[n]=va_arg( ap, int );
 225
 226     va_end(ap);
 227
 228     return (*this)( index );
 229 }
 230
 231 const EvtComplex & EvtSpinAmp::operator()( int i, ... ) const
 232 {
 233     vector<int> index( _twospin.size() );
 234     va_list ap;
 235
 236     va_start(ap, i);
 237
 238     index[0]=i;
 239     for(size_t n=1; n<_twospin.size(); ++n )
 240         index[n]=va_arg( ap, int );
 241
 242     va_end(ap);
 243
 244     return (*this)( index );
 245 }
 246
 247 EvtSpinAmp& EvtSpinAmp::operator=( const EvtSpinAmp& cont )
 248 {
 249     _twospin=cont._twospin;
 250     _elem=cont._elem;
 251     _type=cont._type;
 252
 253     return *this;
 254 }
 255
 256 EvtSpinAmp EvtSpinAmp::operator+( const EvtSpinAmp & cont ) const
 257 {
 258     checktwospin( cont._twospin );
 259
 260     EvtSpinAmp ret( cont );
 261     for( size_t i=0; i<ret._elem.size(); ++i ) {
 262         ret._elem[i]+=_elem[i];
 263     }
 264
 265     return ret;
 266 }
 267
 268 EvtSpinAmp& EvtSpinAmp::operator+=( const EvtSpinAmp&
 269         cont )
 270 {
 271     checktwospin( cont._twospin );
 272
 273     for( size_t i=0; i<_elem.size(); ++i )
 274         _elem[i]+=cont._elem[i];
 275
 276     return *this;
 277 }
 278
 279 EvtSpinAmp EvtSpinAmp::operator-( const EvtSpinAmp & cont ) const
 280 {
 281     checktwospin( cont._twospin );
 282
 283     EvtSpinAmp ret( *this );
 284     for( size_t i=0; i<ret._elem.size(); ++i )
 285         ret._elem[i]-=cont._elem[i];
 286
 287     return ret;
 288 }
 289
 290 EvtSpinAmp& EvtSpinAmp::operator-=( const EvtSpinAmp& cont )
 291 {
 292     checktwospin( cont._twospin );
 293
 294     for( size_t i=0; i<_elem.size(); ++i )
 295         _elem[i]-=cont._elem[i];
 296
 297     return *this;
 298 }
 299
 300 // amp = amp1 * amp2
 301 EvtSpinAmp EvtSpinAmp::operator*( const EvtSpinAmp & amp2 ) const
 302 {
 303     vector<int> index(rank()+amp2.rank());
 304     vector<int> index1(rank()), index2(amp2.rank());
 305     EvtSpinAmp amp;
 306
 307     amp._twospin=_twospin;
 308     amp._type=_type;
 309
 310     for( size_t i=0; i<amp2._twospin.size(); ++i ) {
 311         amp._twospin.push_back( amp2._twospin[i] );
 312         amp._type.push_back( amp2._type[i] );
 313     }
 314
 315     amp._elem = vector<EvtComplex>( _elem.size() * amp2._elem.size() );
 316
 317     for( size_t i=0; i<index1.size(); ++i )
 318         index[i]=index1[i]=-_twospin[i];
 319
 320     for( size_t i=0; i<index2.size(); ++i )
 321         index[i+rank()]=index2[i]=-amp2._twospin[i];
 322
 323     while(true) {
 324         amp( index ) = (*this)( index1 )*amp2( index2 );
 325         if(!amp.iterate( index )) break;
 326
 327         for( size_t i=0; i<index1.size(); ++i )
 328             index1[i]=index[i];
 329
 330         for( size_t i=0; i<index2.size(); ++i )
 331             index2[i]=index[i+rank()];
 332     }
 333
 334     return amp;
 335 }
 336
 337 EvtSpinAmp& EvtSpinAmp::operator*=( const EvtSpinAmp& cont )
 338 {
 339     EvtSpinAmp ret = (*this)*cont;
 340     *this=ret;
 341     return *this;
 342 }
 343
 344 EvtSpinAmp& EvtSpinAmp::operator*=( const EvtComplex& real )
 345 {
 346     for( size_t i=0; i<_elem.size(); ++i )
 347         _elem[i] *= real;
 348
 349     return *this;
 350 }
 351
 352 EvtSpinAmp& EvtSpinAmp::operator/=( const EvtComplex& real )
 353 {
 354     for( size_t i=0; i<_elem.size(); ++i )
 355         _elem[i] /= real;
 356
 357     return *this;
 358 }
 359
 360 vector<int> EvtSpinAmp::iterinit() const
 361 {
 362     vector<int> init( _twospin.size() );
 363
 364     for( size_t i=0; i<_twospin.size(); ++i )
 365         init[i]=-_twospin[i];
 366
 367     return init;
 368 }
 369
 370 bool EvtSpinAmp::iterate( vector<int>& index ) const
 371 {
 372     int last = _twospin.size() - 1;
 373
 374     index[0]+=2;
 375     for( size_t j=0; static_cast<int>(j)<last; ++j ) {
 376         if( index[j] > static_cast<int>(_twospin[j]) ) {
 377             index[j] = -_twospin[j];
 378             index[j+1]+=2;
 379         }
 380     }
 381
 382     return (abs(index[last]))<=((int)_twospin[last]);
 383 }
 384
 385 // Test whether a particular index is an allowed one (specifically to deal with
 386 // photons and possibly neutrinos)
 387 bool EvtSpinAmp::allowed( const vector<int>& index ) const
 388 {
 389     if( index.size() != _type.size() ) {
 390         report(ERROR,"EvtGen")
 391             <<"Wrong dimensino index input to allowed."<<endl;
 392         ::abort();
 393     }
 394
 395     for(size_t i=0; i<index.size(); ++i) {
 396         switch(_type[i]) {
 397             case EvtSpinType::PHOTON:
 398                 if(abs(index[i])!=2) return false;
 399                 break;
 400             case EvtSpinType::NEUTRINO:
 401                 report(ERROR,"EvtGen")
 402                     <<"EvMultibody currently cannot handle neutrinos."<<endl;
 403                 ::abort();
 404             default:
 405               break;
 406         }
 407     }
 408
 409     return true;
 410 }
 411
 412 bool EvtSpinAmp::iterateallowed( vector<int>& index ) const
 413 {
 414     while(true) {
 415         if(!iterate( index ))
 416             return false;
 417         if(allowed( index ))
 418             return true;
 419     }
 420 }
 421
 422 vector<int> EvtSpinAmp::iterallowedinit() const
 423 {
 424     vector<int> init = iterinit();
 425     while(!allowed(init)) {
 426         iterate(init);
 427     }
 428
 429     return init;
 430 }
 431
 432 void EvtSpinAmp::intcont( size_t a, size_t b )
 433 {
 434
 435     if(rank()<=2) {
 436       report(ERROR,"EvtGen")<<"EvtSpinAmp can't handle no indices" << endl;
 437       ::abort();
 438     }
 439
 440     size_t newrank=rank()-2;
 441
 442     if(_twospin[a]!=_twospin[b]) {
 443         report(ERROR,"EvtGen")
 444             <<"Contaction called on indices of different dimension"
 445             <<endl;
 446         report(ERROR,"EvtGen")
 447             <<"Called on "<<_twospin[a]<<" and "<<_twospin[b]
 448             <<endl;
 449         ::abort();
 450     }
 451
 452     vector<int> newtwospin( newrank );
 453     vector<EvtSpinType::spintype> newtype( newrank );
 454
 455     for( size_t i=0, j=0; i<_twospin.size(); ++i ){
 456         if(i==a || i==b) continue;
 457
 458         newtwospin[j] = _twospin[i];
 459         newtype[j] = _type[i];
 460         ++j;
 461     }
 462
 463     EvtSpinAmp newamp( newtype );
 464     vector<int> index( rank() ), newindex = newamp.iterinit();
 465
 466     for( size_t i=0; i<newrank; ++i )
 467         newindex[i] = -newtwospin[i];
 468
 469     while(true) {
 470
 471         for( size_t i=0, j=0; i<rank(); ++i ) {
 472             if(i==a || i==b) continue;
 473             index[i]=newindex[j];
 474             ++j;
 475         }
 476
 477         index[b]=index[a]=-_twospin[a];
 478         newamp(newindex) = (*this)(index);
 479         for( size_t i=-_twospin[a]+2; i<=_twospin[a]; i+=2 ) {
 480             index[b]=index[a]=i;
 481             newamp(newindex) += (*this)(index);
 482         }
 483
 484         if(!newamp.iterate(newindex)) break;
 485     }
 486
 487     *this=newamp;
 488 }
 489
 490 // In A.extcont(B), a is the index in A and b is the index in B - note that this
 491 // routine can be extremely improved!
 492 void EvtSpinAmp::extcont( const EvtSpinAmp& cont, int a, int b )
 493 {
 494     EvtSpinAmp ret = (*this)*cont;
 495     ret.intcont( a, rank()+b );
 496
 497     *this=ret;
 498 }