00001 #include "EvtGenBase/EvtPatches.hh"
00002 #include "EvtGenBase/EvtReport.hh"
00003 #include "EvtGenBase/EvtSpinAmp.hh"
00004 #include <cstdlib>
00005
00006 using std::endl;
00007
00008 std::ostream&
00009 operator<<( std::ostream& os, const EvtSpinAmp& amp )
00010 {
00011 vector<int> index = amp.iterinit();
00012
00013 os << ":";
00014 do {
00015 os<<"<";
00016 for(int i=0; i<index.size()-1; ++i) {
00017 os<<index[i];
00018 }
00019 os<<index[index.size()-1]<<">"<<amp(index)<<":";
00020 } while( amp.iterate( index ) );
00021
00022 return os;
00023 }
00024
00025 EvtSpinAmp operator*( const EvtComplex& real, const EvtSpinAmp& cont )
00026 {
00027 EvtSpinAmp ret( cont );
00028
00029 for( int i=0; i<ret._elem.size(); ++i ) {
00030 ret._elem[i] *= real;
00031 }
00032
00033 return ret;
00034 }
00035
00036 EvtSpinAmp operator*( const EvtSpinAmp& cont, const EvtComplex& real )
00037 {
00038 return real*cont;
00039 }
00040
00041 EvtSpinAmp operator/( const EvtSpinAmp& cont, const EvtComplex& real )
00042 {
00043 EvtSpinAmp ret( cont );
00044
00045 for( int i=0; i<ret._elem.size(); ++i ) {
00046 ret._elem[i] /= real;
00047 }
00048
00049 return ret;
00050 }
00051
00052 vector<int> EvtSpinAmp::calctwospin( const vector<EvtSpinType::spintype>& type ) const
00053 {
00054 vector<int> twospin;
00055
00056 for( int i=0; i<type.size(); ++i ) {
00057 twospin.push_back( EvtSpinType::getSpin2( type[i] ) );
00058 }
00059
00060 return twospin;
00061 }
00062
00063 EvtSpinAmp::EvtSpinAmp( const vector<EvtSpinType::spintype>& type )
00064 {
00065 int num = 1;
00066 _type = type;
00067 _twospin=calctwospin( type );
00068
00069 for( int i=0; i<_twospin.size(); ++i )
00070 num*=_twospin[i]+1;
00071
00072 _elem=vector<EvtComplex>( num );
00073 }
00074
00075 EvtSpinAmp::EvtSpinAmp( const vector<EvtSpinType::spintype>& type, const EvtComplex & val )
00076 {
00077 int num = 1;
00078 _type = type;
00079 _twospin=calctwospin( type );
00080
00081 for( int i=0; i<_twospin.size(); ++i )
00082 num*=_twospin[i]+1;
00083
00084 _elem=vector<EvtComplex>( num, val );
00085 }
00086
00087 EvtSpinAmp::EvtSpinAmp( const vector<EvtSpinType::spintype>& type,
00088 const vector<EvtComplex>& elem )
00089 {
00090 int num = 1;
00091
00092 _type = type;
00093 _twospin=calctwospin( type );
00094 _elem=elem;
00095
00096 for( int i=0; i<_twospin.size(); ++i )
00097 num*=_twospin[i]+1;
00098
00099 if(_elem.size() != num ) {
00100 report(ERROR,"EvtGen")<<"Wrong number of elements input:"
00101 <<_elem.size()<<" vs. "<<num<<endl;
00102 ::abort();
00103 }
00104
00105 }
00106
00107 EvtSpinAmp::EvtSpinAmp( const EvtSpinAmp& copy )
00108 {
00109 _twospin = copy._twospin;
00110 _elem = copy._elem;
00111 _type = copy._type;
00112 }
00113
00114 void EvtSpinAmp::checktwospin( const vector<int>& twospin ) const
00115 {
00116 if( _twospin == twospin )
00117 return;
00118
00119 report( ERROR, "EvtGen" )
00120 <<"Dimension or order of tensors being operated on does not match"
00121 <<endl;
00122 ::abort();
00123 }
00124
00125 void EvtSpinAmp::checkindexargs( const vector<int>& index ) const
00126 {
00127 if( index.size()==0 ) {
00128 report(ERROR,"EvtGen") << "EvtSpinAmp can't handle no indices" << endl;
00129 ::abort();
00130 }
00131
00132 if( index.size() != _twospin.size() ) {
00133 report( ERROR, "EvtGen" ) << "Rank of EvtSpinAmp index does not match: "
00134 <<_twospin.size()<<" expected "<<index.size()<<" input."<<endl;
00135 ::abort();
00136 }
00137
00138 for( int i=0; i<_twospin.size(); ++i ) {
00139 if( _twospin[i]>=abs(index[i]) && _twospin[i]%2==index[i]%2 )
00140 continue;
00141 report(ERROR,"EvtGen")<<"EvtSpinAmp index out of range" << endl;
00142 report(ERROR,"EvtGen")<<" Index: ";
00143 for(int j=0; j<_twospin.size(); ++j )
00144 report(ERROR," ")<<_twospin[j];
00145
00146 report(ERROR, " ")<<endl<<" Index: ";
00147 for(int j=0; j<index.size(); ++j )
00148 report(ERROR," ")<<index[j];
00149 report(ERROR, " ")<<endl;
00150 ::abort();
00151 }
00152 }
00153
00154 int EvtSpinAmp::findtrueindex( const vector<int>& index ) const
00155 {
00156 int trueindex = 0;
00157
00158 for( int i = index.size()-1; i>0; --i ) {
00159 trueindex += (index[i]+_twospin[i])/2;
00160 trueindex *= _twospin[i-1]+1;
00161 }
00162
00163 trueindex += (index[0]+_twospin[0])/2;
00164
00165 return trueindex;
00166 }
00167
00168 EvtComplex & EvtSpinAmp::operator()( const vector<int>& index )
00169 {
00170 checkindexargs( index );
00171
00172 int trueindex = findtrueindex(index);
00173 if(trueindex >= _elem.size()) {
00174 report(ERROR,"EvtGen")<<"indexing error "<<trueindex<<" "<<_elem.size()<<endl;
00175 for(int i=0; i<_twospin.size(); ++i) {
00176 report(ERROR,"")<<_twospin[i]<<" ";
00177 }
00178 report(ERROR,"")<<endl;
00179
00180 for(int i=0; i<index.size(); ++i) {
00181 report(ERROR,"")<<index[i]<<" ";
00182 }
00183 report(ERROR,"")<<endl;
00184
00185 ::abort();
00186 }
00187
00188 return _elem[trueindex];
00189 }
00190
00191 const EvtComplex & EvtSpinAmp::operator()( const vector<int>& index ) const
00192 {
00193 checkindexargs( index );
00194
00195 int trueindex = findtrueindex(index);
00196 if(trueindex >= _elem.size()) {
00197 report(ERROR,"EvtGen")<<"indexing error "<<trueindex<<" "<<_elem.size()<<endl;
00198 for(int i=0; i<_twospin.size(); ++i) {
00199 report(ERROR,"")<<_twospin[i]<<" ";
00200 }
00201 report(ERROR,"")<<endl;
00202
00203 for(int i=0; i<index.size(); ++i) {
00204 report(ERROR,"")<<index[i]<<" ";
00205 }
00206 report(ERROR,"")<<endl;
00207
00208 ::abort();
00209 }
00210
00211 return _elem[trueindex];
00212 }
00213
00214 EvtComplex & EvtSpinAmp::operator()( int i, ... )
00215 {
00216 va_list ap;
00217 vector<int> index( _twospin.size() );
00218
00219 va_start(ap, i);
00220
00221 index[0]=i;
00222 for(int n=1; n<_twospin.size(); ++n )
00223 index[n]=va_arg( ap, int );
00224
00225 va_end(ap);
00226
00227 return (*this)( index );
00228 }
00229
00230 const EvtComplex & EvtSpinAmp::operator()( int i, ... ) const
00231 {
00232 vector<int> index( _twospin.size() );
00233 va_list ap;
00234
00235 va_start(ap, i);
00236
00237 index[0]=i;
00238 for(int n=1; n<_twospin.size(); ++n )
00239 index[n]=va_arg( ap, int );
00240
00241 va_end(ap);
00242
00243 return (*this)( index );
00244 }
00245
00246 EvtSpinAmp& EvtSpinAmp::operator=( const EvtSpinAmp& cont )
00247 {
00248 _twospin=cont._twospin;
00249 _elem=cont._elem;
00250 _type=cont._type;
00251
00252 return *this;
00253 }
00254
00255 EvtSpinAmp EvtSpinAmp::operator+( const EvtSpinAmp & cont ) const
00256 {
00257 checktwospin( cont._twospin );
00258
00259 EvtSpinAmp ret( cont );
00260 for( int i=0; i<ret._elem.size(); ++i ) {
00261 ret._elem[i]+=_elem[i];
00262 }
00263
00264 return ret;
00265 }
00266
00267 EvtSpinAmp& EvtSpinAmp::operator+=( const EvtSpinAmp&
00268 cont )
00269 {
00270 checktwospin( cont._twospin );
00271
00272 for( int i=0; i<_elem.size(); ++i )
00273 _elem[i]+=cont._elem[i];
00274
00275 return *this;
00276 }
00277
00278 EvtSpinAmp EvtSpinAmp::operator-( const EvtSpinAmp & cont ) const
00279 {
00280 checktwospin( cont._twospin );
00281
00282 EvtSpinAmp ret( *this );
00283 for( int i=0; i<ret._elem.size(); ++i )
00284 ret._elem[i]-=cont._elem[i];
00285
00286 return ret;
00287 }
00288
00289 EvtSpinAmp& EvtSpinAmp::operator-=( const EvtSpinAmp& cont )
00290 {
00291 checktwospin( cont._twospin );
00292
00293 for( int i=0; i<_elem.size(); ++i )
00294 _elem[i]-=cont._elem[i];
00295
00296 return *this;
00297 }
00298
00299
00300 EvtSpinAmp EvtSpinAmp::operator*( const EvtSpinAmp & amp2 ) const
00301 {
00302 vector<int> index(rank()+amp2.rank());
00303 vector<int> index1(rank()), index2(amp2.rank());
00304 EvtSpinAmp amp;
00305
00306 amp._twospin=_twospin;
00307 amp._type=_type;
00308
00309 for( int i=0; i<amp2._twospin.size(); ++i ) {
00310 amp._twospin.push_back( amp2._twospin[i] );
00311 amp._type.push_back( amp2._type[i] );
00312 }
00313
00314 amp._elem = vector<EvtComplex>( _elem.size() * amp2._elem.size() );
00315
00316 for( int i=0; i<index1.size(); ++i )
00317 index[i]=index1[i]=-_twospin[i];
00318
00319 for( int i=0; i<index2.size(); ++i )
00320 index[i+rank()]=index2[i]=-amp2._twospin[i];
00321
00322 while(true) {
00323 amp( index ) = (*this)( index1 )*amp2( index2 );
00324 if(!amp.iterate( index )) break;
00325
00326 for( int i=0; i<index1.size(); ++i )
00327 index1[i]=index[i];
00328
00329 for( int i=0; i<index2.size(); ++i )
00330 index2[i]=index[i+rank()];
00331 }
00332
00333 return amp;
00334 }
00335
00336 EvtSpinAmp& EvtSpinAmp::operator*=( const EvtSpinAmp& cont )
00337 {
00338 EvtSpinAmp ret = (*this)*cont;
00339 *this=ret;
00340 return *this;
00341 }
00342
00343 EvtSpinAmp& EvtSpinAmp::operator*=( const EvtComplex& real )
00344 {
00345 for( int i=0; i<_elem.size(); ++i )
00346 _elem[i] *= real;
00347
00348 return *this;
00349 }
00350
00351 EvtSpinAmp& EvtSpinAmp::operator/=( const EvtComplex& real )
00352 {
00353 for( int i=0; i<_elem.size(); ++i )
00354 _elem[i] /= real;
00355
00356 return *this;
00357 }
00358
00359 vector<int> EvtSpinAmp::iterinit() const
00360 {
00361 vector<int> init( _twospin.size() );
00362
00363 for( int i=0; i<_twospin.size(); ++i )
00364 init[i]=-_twospin[i];
00365
00366 return init;
00367 }
00368
00369 bool EvtSpinAmp::iterate( vector<int>& index ) const
00370 {
00371 int last = _twospin.size() - 1;
00372
00373 index[0]+=2;
00374 for( int j=0; j<last; ++j ) {
00375 if( index[j] > _twospin[j] ) {
00376 index[j] = -_twospin[j];
00377 index[j+1]+=2;
00378 }
00379 }
00380
00381 return abs(index[last])<=_twospin[last];
00382 }
00383
00384
00385
00386 bool EvtSpinAmp::allowed( const vector<int>& index ) const
00387 {
00388 if( index.size() != _type.size() ) {
00389 report(ERROR,"EvtGen")
00390 <<"Wrong dimensino index input to allowed."<<endl;
00391 ::abort();
00392 }
00393
00394 for(int i=0; i<index.size(); ++i) {
00395 switch(_type[i]) {
00396 case EvtSpinType::PHOTON:
00397 if(abs(index[i])!=2) return false;
00398 break;
00399 case EvtSpinType::NEUTRINO:
00400 report(ERROR,"EvtGen")
00401 <<"EvMultibody currently cannot handle neutrinos."<<endl;
00402 ::abort();
00403 default:
00404 break;
00405 }
00406 }
00407
00408 return true;
00409 }
00410
00411 bool EvtSpinAmp::iterateallowed( vector<int>& index ) const
00412 {
00413 while(true) {
00414 if(!iterate( index ))
00415 return false;
00416 if(allowed( index ))
00417 return true;
00418 }
00419 }
00420
00421 vector<int> EvtSpinAmp::iterallowedinit() const
00422 {
00423 vector<int> init = iterinit();
00424 while(!allowed(init)) {
00425 iterate(init);
00426 }
00427
00428 return init;
00429 }
00430
00431 void EvtSpinAmp::intcont( int a, int b )
00432 {
00433 int newrank=rank()-2;
00434 if(newrank<=0) {
00435 report(ERROR,"EvtGen")<<"EvtSpinAmp can't handle no indices" << endl;
00436 ::abort();
00437 }
00438
00439 if(_twospin[a]!=_twospin[b]) {
00440 report(ERROR,"EvtGen")
00441 <<"Contaction called on indices of different dimension"
00442 <<endl;
00443 report(ERROR,"EvtGen")
00444 <<"Called on "<<_twospin[a]<<" and "<<_twospin[b]
00445 <<endl;
00446 ::abort();
00447 }
00448
00449 vector<int> newtwospin( newrank );
00450 vector<EvtSpinType::spintype> newtype( newrank );
00451
00452 for( int i=0, j=0; i<_twospin.size(); ++i ){
00453 if(i==a || i==b) continue;
00454
00455 newtwospin[j] = _twospin[i];
00456 newtype[j] = _type[i];
00457 ++j;
00458 }
00459
00460 EvtSpinAmp newamp( newtype );
00461 vector<int> index( rank() ), newindex = newamp.iterinit();
00462
00463 for( int i=0; i<newrank; ++i )
00464 newindex[i] = -newtwospin[i];
00465
00466 while(true) {
00467
00468 for( int i=0, j=0; i<rank(); ++i ) {
00469 if(i==a || i==b) continue;
00470 index[i]=newindex[j];
00471 ++j;
00472 }
00473
00474 index[b]=index[a]=-_twospin[a];
00475 newamp(newindex) = (*this)(index);
00476 for( int i=-_twospin[a]+2; i<=_twospin[a]; i+=2 ) {
00477 index[b]=index[a]=i;
00478 newamp(newindex) += (*this)(index);
00479 }
00480
00481 if(!newamp.iterate(newindex)) break;
00482 }
00483
00484 *this=newamp;
00485 }
00486
00487
00488
00489 void EvtSpinAmp::extcont( const EvtSpinAmp& cont, int a, int b )
00490 {
00491 EvtSpinAmp ret = (*this)*cont;
00492 ret.intcont( a, rank()+b );
00493
00494 *this=ret;
00495 }