+++ /dev/null
-#include "EvtGenBase/EvtPatches.hh"
-#include "EvtGenBase/EvtReport.hh"
-#include "EvtGenBase/EvtSpinAmp.hh"
-#include <stdlib.h>
-
-using std::endl;
-
-std::ostream&
-operator<<( std::ostream& os, const EvtSpinAmp& amp )
-{
- vector<int> index = amp.iterinit();
-
- os << ":";
- do {
- os<<"<";
- for(size_t i=0; i<index.size()-1; ++i) {
- os<<index[i];
- }
- os<<index[index.size()-1]<<">"<<amp(index)<<":";
- } while( amp.iterate( index ) );
-
- return os;
-}
-
-EvtSpinAmp operator*( const EvtComplex& real, const EvtSpinAmp& cont )
-{
- EvtSpinAmp ret( cont );
-
- for( size_t i=0; i<ret._elem.size(); ++i ) {
- ret._elem[i] *= real;
- }
-
- return ret;
-}
-
-EvtSpinAmp operator*( const EvtSpinAmp& cont, const EvtComplex& real )
-{
- return real*cont;
-}
-
-EvtSpinAmp operator/( const EvtSpinAmp& cont, const EvtComplex& real )
-{
- EvtSpinAmp ret( cont );
-
- for( size_t i=0; i<ret._elem.size(); ++i ) {
- ret._elem[i] /= real;
- }
-
- return ret;
-}
-
-vector<unsigned int> EvtSpinAmp::calctwospin( const vector<EvtSpinType::spintype>& type ) const
-{
- vector<unsigned int> twospin;
-
- for( size_t i=0; i<type.size(); ++i ) {
- twospin.push_back( EvtSpinType::getSpin2( type[i] ) );
- }
-
- return twospin;
-}
-
-EvtSpinAmp::EvtSpinAmp( const vector<EvtSpinType::spintype>& type )
-{
- int num = 1;
- _type = type;
- _twospin=calctwospin( type );
-
- for( size_t i=0; i<_twospin.size(); ++i )
- num*=_twospin[i]+1;
-
- _elem=vector<EvtComplex>( num );
-}
-
-EvtSpinAmp::EvtSpinAmp( const vector<EvtSpinType::spintype>& type, const EvtComplex & val )
-{
- int num = 1;
- _type = type;
- _twospin=calctwospin( type );
-
- for( size_t i=0; i<_twospin.size(); ++i )
- num*=_twospin[i]+1;
-
- _elem=vector<EvtComplex>( num, val );
-}
-
-EvtSpinAmp::EvtSpinAmp( const vector<EvtSpinType::spintype>& type,
- const vector<EvtComplex>& elem )
-{
- size_t num = 1;
-
- _type = type;
- _twospin=calctwospin( type );
- _elem=elem;
-
- for( size_t i=0; i<_twospin.size(); ++i ){
- num*=(_twospin[i]+1);
- }
-
- if(_elem.size() != num ) {
- report(ERROR,"EvtGen")<<"Wrong number of elements input:"
- <<_elem.size()<<" vs. "<<num<<endl;
- ::abort();
- }
-
-}
-
-EvtSpinAmp::EvtSpinAmp( const EvtSpinAmp& copy )
-{
- _twospin = copy._twospin;
- _elem = copy._elem;
- _type = copy._type;
-}
-
-void EvtSpinAmp::checktwospin( const vector<unsigned int>& twospin ) const
-{
- if( _twospin == twospin )
- return;
-
- report( ERROR, "EvtGen" )
- <<"Dimension or order of tensors being operated on does not match"
- <<endl;
- ::abort();
-}
-
-void EvtSpinAmp::checkindexargs( const vector<int>& index ) const
-{
- if( index.size()==0 ) {
- report(ERROR,"EvtGen") << "EvtSpinAmp can't handle no indices" << endl;
- ::abort();
- }
-
- if( index.size() != _twospin.size() ) {
- report( ERROR, "EvtGen" ) << "Rank of EvtSpinAmp index does not match: "
- <<_twospin.size()<<" expected "<<index.size()<<" input."<<endl;
- ::abort();
- }
-
- for( size_t i=0; i<_twospin.size(); ++i ) {
- if( static_cast<int>(_twospin[i])>=abs(index[i]) && static_cast<int>(_twospin[i])%2==index[i]%2 )
- continue;
- report(ERROR,"EvtGen")<<"EvtSpinAmp index out of range" << endl;
- report(ERROR,"EvtGen")<<" Index: ";
- for(size_t j=0; j<_twospin.size(); ++j )
- report(ERROR," ")<<_twospin[j];
-
- report(ERROR, " ")<<endl<<" Index: ";
- for(size_t j=0; j<index.size(); ++j )
- report(ERROR," ")<<index[j];
- report(ERROR, " ")<<endl;
- ::abort();
- }
-}
-
-int EvtSpinAmp::findtrueindex( const vector<int>& index ) const
-{
- int trueindex = 0;
-
- for( size_t i = index.size()-1; i>0; --i ) {
- trueindex += (index[i]+_twospin[i])/2;
- trueindex *= _twospin[i-1]+1;
- }
-
- trueindex += (index[0]+_twospin[0])/2;
-
- return trueindex;
-}
-
-EvtComplex & EvtSpinAmp::operator()( const vector<int>& index )
-{
- checkindexargs( index );
-
- size_t trueindex = findtrueindex(index);
- if(trueindex >= _elem.size()) {
- report(ERROR,"EvtGen")<<"indexing error "<<trueindex<<" "<<_elem.size()<<endl;
- for(size_t i=0; i<_twospin.size(); ++i) {
- report(ERROR,"")<<_twospin[i]<<" ";
- }
- report(ERROR,"")<<endl;
-
- for(size_t i=0; i<index.size(); ++i) {
- report(ERROR,"")<<index[i]<<" ";
- }
- report(ERROR,"")<<endl;
-
- ::abort();
- }
-
- return _elem[trueindex];
-}
-
-const EvtComplex & EvtSpinAmp::operator()( const vector<int>& index ) const
-{
- checkindexargs( index );
-
- size_t trueindex = findtrueindex(index);
- if(trueindex >= _elem.size()) {
- report(ERROR,"EvtGen")<<"indexing error "<<trueindex<<" "<<_elem.size()<<endl;
- for(size_t i=0; i<_twospin.size(); ++i) {
- report(ERROR,"")<<_twospin[i]<<" ";
- }
- report(ERROR,"")<<endl;
-
- for(size_t i=0; i<index.size(); ++i) {
- report(ERROR,"")<<index[i]<<" ";
- }
- report(ERROR,"")<<endl;
-
- ::abort();
- }
-
- return _elem[trueindex];
-}
-
-EvtComplex & EvtSpinAmp::operator()( int i, ... )
-{
- va_list ap;
- vector<int> index( _twospin.size() );
-
- va_start(ap, i);
-
- index[0]=i;
- for(size_t n=1; n<_twospin.size(); ++n )
- index[n]=va_arg( ap, int );
-
- va_end(ap);
-
- return (*this)( index );
-}
-
-const EvtComplex & EvtSpinAmp::operator()( int i, ... ) const
-{
- vector<int> index( _twospin.size() );
- va_list ap;
-
- va_start(ap, i);
-
- index[0]=i;
- for(size_t n=1; n<_twospin.size(); ++n )
- index[n]=va_arg( ap, int );
-
- va_end(ap);
-
- return (*this)( index );
-}
-
-EvtSpinAmp& EvtSpinAmp::operator=( const EvtSpinAmp& cont )
-{
- _twospin=cont._twospin;
- _elem=cont._elem;
- _type=cont._type;
-
- return *this;
-}
-
-EvtSpinAmp EvtSpinAmp::operator+( const EvtSpinAmp & cont ) const
-{
- checktwospin( cont._twospin );
-
- EvtSpinAmp ret( cont );
- for( size_t i=0; i<ret._elem.size(); ++i ) {
- ret._elem[i]+=_elem[i];
- }
-
- return ret;
-}
-
-EvtSpinAmp& EvtSpinAmp::operator+=( const EvtSpinAmp&
- cont )
-{
- checktwospin( cont._twospin );
-
- for( size_t i=0; i<_elem.size(); ++i )
- _elem[i]+=cont._elem[i];
-
- return *this;
-}
-
-EvtSpinAmp EvtSpinAmp::operator-( const EvtSpinAmp & cont ) const
-{
- checktwospin( cont._twospin );
-
- EvtSpinAmp ret( *this );
- for( size_t i=0; i<ret._elem.size(); ++i )
- ret._elem[i]-=cont._elem[i];
-
- return ret;
-}
-
-EvtSpinAmp& EvtSpinAmp::operator-=( const EvtSpinAmp& cont )
-{
- checktwospin( cont._twospin );
-
- for( size_t i=0; i<_elem.size(); ++i )
- _elem[i]-=cont._elem[i];
-
- return *this;
-}
-
-// amp = amp1 * amp2
-EvtSpinAmp EvtSpinAmp::operator*( const EvtSpinAmp & amp2 ) const
-{
- vector<int> index(rank()+amp2.rank());
- vector<int> index1(rank()), index2(amp2.rank());
- EvtSpinAmp amp;
-
- amp._twospin=_twospin;
- amp._type=_type;
-
- for( size_t i=0; i<amp2._twospin.size(); ++i ) {
- amp._twospin.push_back( amp2._twospin[i] );
- amp._type.push_back( amp2._type[i] );
- }
-
- amp._elem = vector<EvtComplex>( _elem.size() * amp2._elem.size() );
-
- for( size_t i=0; i<index1.size(); ++i )
- index[i]=index1[i]=-_twospin[i];
-
- for( size_t i=0; i<index2.size(); ++i )
- index[i+rank()]=index2[i]=-amp2._twospin[i];
-
- while(true) {
- amp( index ) = (*this)( index1 )*amp2( index2 );
- if(!amp.iterate( index )) break;
-
- for( size_t i=0; i<index1.size(); ++i )
- index1[i]=index[i];
-
- for( size_t i=0; i<index2.size(); ++i )
- index2[i]=index[i+rank()];
- }
-
- return amp;
-}
-
-EvtSpinAmp& EvtSpinAmp::operator*=( const EvtSpinAmp& cont )
-{
- EvtSpinAmp ret = (*this)*cont;
- *this=ret;
- return *this;
-}
-
-EvtSpinAmp& EvtSpinAmp::operator*=( const EvtComplex& real )
-{
- for( size_t i=0; i<_elem.size(); ++i )
- _elem[i] *= real;
-
- return *this;
-}
-
-EvtSpinAmp& EvtSpinAmp::operator/=( const EvtComplex& real )
-{
- for( size_t i=0; i<_elem.size(); ++i )
- _elem[i] /= real;
-
- return *this;
-}
-
-vector<int> EvtSpinAmp::iterinit() const
-{
- vector<int> init( _twospin.size() );
-
- for( size_t i=0; i<_twospin.size(); ++i )
- init[i]=-_twospin[i];
-
- return init;
-}
-
-bool EvtSpinAmp::iterate( vector<int>& index ) const
-{
- int last = _twospin.size() - 1;
-
- index[0]+=2;
- for( size_t j=0; static_cast<int>(j)<last; ++j ) {
- if( index[j] > static_cast<int>(_twospin[j]) ) {
- index[j] = -_twospin[j];
- index[j+1]+=2;
- }
- }
-
- return (abs(index[last]))<=((int)_twospin[last]);
-}
-
-// Test whether a particular index is an allowed one (specifically to deal with
-// photons and possibly neutrinos)
-bool EvtSpinAmp::allowed( const vector<int>& index ) const
-{
- if( index.size() != _type.size() ) {
- report(ERROR,"EvtGen")
- <<"Wrong dimensino index input to allowed."<<endl;
- ::abort();
- }
-
- for(size_t i=0; i<index.size(); ++i) {
- switch(_type[i]) {
- case EvtSpinType::PHOTON:
- if(abs(index[i])!=2) return false;
- break;
- case EvtSpinType::NEUTRINO:
- report(ERROR,"EvtGen")
- <<"EvMultibody currently cannot handle neutrinos."<<endl;
- ::abort();
- default:
- break;
- }
- }
-
- return true;
-}
-
-bool EvtSpinAmp::iterateallowed( vector<int>& index ) const
-{
- while(true) {
- if(!iterate( index ))
- return false;
- if(allowed( index ))
- return true;
- }
-}
-
-vector<int> EvtSpinAmp::iterallowedinit() const
-{
- vector<int> init = iterinit();
- while(!allowed(init)) {
- iterate(init);
- }
-
- return init;
-}
-
-void EvtSpinAmp::intcont( size_t a, size_t b )
-{
-
- if(rank()<=2) {
- report(ERROR,"EvtGen")<<"EvtSpinAmp can't handle no indices" << endl;
- ::abort();
- }
-
- size_t newrank=rank()-2;
-
- if(_twospin[a]!=_twospin[b]) {
- report(ERROR,"EvtGen")
- <<"Contaction called on indices of different dimension"
- <<endl;
- report(ERROR,"EvtGen")
- <<"Called on "<<_twospin[a]<<" and "<<_twospin[b]
- <<endl;
- ::abort();
- }
-
- vector<int> newtwospin( newrank );
- vector<EvtSpinType::spintype> newtype( newrank );
-
- for( size_t i=0, j=0; i<_twospin.size(); ++i ){
- if(i==a || i==b) continue;
-
- newtwospin[j] = _twospin[i];
- newtype[j] = _type[i];
- ++j;
- }
-
- EvtSpinAmp newamp( newtype );
- vector<int> index( rank() ), newindex = newamp.iterinit();
-
- for( size_t i=0; i<newrank; ++i )
- newindex[i] = -newtwospin[i];
-
- while(true) {
-
- for( size_t i=0, j=0; i<rank(); ++i ) {
- if(i==a || i==b) continue;
- index[i]=newindex[j];
- ++j;
- }
-
- index[b]=index[a]=-_twospin[a];
- newamp(newindex) = (*this)(index);
- for( size_t i=-_twospin[a]+2; i<=_twospin[a]; i+=2 ) {
- index[b]=index[a]=i;
- newamp(newindex) += (*this)(index);
- }
-
- if(!newamp.iterate(newindex)) break;
- }
-
- *this=newamp;
-}
-
-// In A.extcont(B), a is the index in A and b is the index in B - note that this
-// routine can be extremely improved!
-void EvtSpinAmp::extcont( const EvtSpinAmp& cont, int a, int b )
-{
- EvtSpinAmp ret = (*this)*cont;
- ret.intcont( a, rank()+b );
-
- *this=ret;
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff