00001 #include "EvtGenBase/EvtPatches.hh"
00002 #include "EvtGenBase/EvtParticle.hh"
00003 #include "EvtGenBase/EvtGenKine.hh"
00004 #include "EvtGenBase/EvtPDL.hh"
00005 #include "EvtGenBase/EvtReport.hh"
00006 #include "EvtGenBase/EvtResonance.hh"
00007 #include "EvtGenBase/EvtResonance2.hh"
00008 #include "EvtGenModels/EvtMultibody.hh"
00009 #include "EvtGenBase/EvtConst.hh"
00010 #include "EvtGenBase/EvtdFunction.hh"
00011 #include "EvtGenBase/EvtKine.hh"
00012
00013 EvtMultibody::~EvtMultibody()
00014 {
00015 if( _decayTree != NULL ) delete _decayTree;
00016 _decayTree=NULL;
00017 }
00018
00019 void EvtMultibody::getName(std::string& model_name)
00020 {
00021 model_name = "D_MULTIBODY";
00022 }
00023
00024 EvtDecayBase* EvtMultibody::clone()
00025 {
00026 return new EvtMultibody;
00027 }
00028
00029 void EvtMultibody::init()
00030 {
00031 int N = getNArg();
00032
00033 _decayTree = new EvtMTree( getDaugs(), getNDaug() );
00034
00035 for(int i=0; i<N-1; ++i) {
00036 if(getArgStr( i )=="RESONANCE") {
00037 _decayTree->addtree( getArgStr( ++i ) );
00038 } else {
00039 report(ERROR,"EvtGen")
00040 << "Syntax error at " << getArgStr( i ) << std::endl;
00041 ::abort();
00042 }
00043 }
00044 }
00045
00046
00047
00048
00049 void EvtMultibody::initProbMax()
00050 {
00051
00052 }
00053
00054 void EvtMultibody::decay( EvtParticle *p )
00055 {
00056
00057 p->initializePhaseSpace(getNDaug(),getDaugs());
00058 std::vector<EvtVector4R> product;
00059
00060 for(int i=0; i<getNDaug(); ++i)
00061 product.push_back(p->getDaug(i)->getP4Lab());
00062
00063 EvtSpinAmp amp = _decayTree->amplitude( product );
00064
00065 int * ilist = new int[amp.rank()];
00066
00067
00068 EvtSpinDensity R = p->rotateToHelicityBasis();
00069
00070 EvtSpinType::spintype type = EvtPDL::getSpinType(getParentId());
00071 int twospin = EvtSpinType::getSpin2( EvtPDL::getSpinType(getParentId()) );
00072
00073 std::vector<EvtSpinType::spintype> types(2, type);
00074 EvtSpinAmp newamp( types, EvtComplex(0.0, 0.0) );
00075 std::vector<int> index = newamp.iterallowedinit();
00076 do {
00077 newamp(index) = R.Get((index[0]+twospin)/2,(index[1]+twospin)/2);
00078 } while( newamp.iterateallowed( index ) );
00079
00080 newamp.extcont(amp, 1, 0);
00081 amp=newamp;
00082
00083 index = amp.iterallowedinit();
00084 std::vector<int> spins = amp.dims();
00085
00086 do {
00087 for( int i=0; i<index.size(); ++i ) {
00088 ilist[i]=index[i]+spins[i];
00089 }
00090
00091 vertex( ilist, amp( index ) );
00092 } while( amp.iterateallowed( index ) );
00093
00094 delete [] ilist;
00095 }