Commit 366b2d06 authored by Frank Siegert's avatar Frank Siegert

Completion of colour treatment in hard decay handler, colour is now integrated...

Completion of colour treatment in hard decay handler, colour is now integrated out explicitly with the Color_Integrator. Fixes MC4BSM setup.
parent bcb8878f
......@@ -4,11 +4,13 @@
#include "METOOLS/Explicit/Vertex.H"
#include "MODEL/Main/Model_Base.H"
#include "MODEL/Interaction_Models/Single_Vertex.H"
#include "PHASIC++/Main/Color_Integrator.H"
#include <assert.h>
using namespace EXTRAXS;
using namespace ATOOLS;
using namespace METOOLS;
using namespace PHASIC;
using namespace std;
Comix1to3::Comix1to3(const vector<Flavour>& flavs, const Flavour& prop,
......@@ -103,6 +105,22 @@ Comix1to3::Comix1to3(const vector<Flavour>& flavs, const Flavour& prop,
m_antifcur->InitPols(pols2);
m_antifcur->HM().resize(m_n);
for (size_t i(0);i<m_n;++i) m_antifcur->HM()[i]=i;
p_ci=new Color_Integrator();
Idx_Vector cids(4,0);
Int_Vector acts(4,0), types(4,0);
for (size_t i(0);i<flavs.size();++i) {
cids[i]=i;
acts[i]=flavs[i].Strong();
if (acts[i]) {
if (flavs[i].StrongCharge()==8) types[i]=0;
else if (flavs[i].IsAnti()) types[i]=(i==0?1:-1);
else types[i]=(i==0?-1:1);
}
}
if (!p_ci->ConstructRepresentations(cids,types,acts)) {
THROW(fatal_error, "Internal error.");
}
}
Comix1to3::~Comix1to3()
......@@ -155,9 +173,10 @@ Vertex* Comix1to3::GetVertex(Current* cur1, Current* cur2, Current* prop) {
void Comix1to3::Calculate(const ATOOLS::Vec4D_Vector& momenta, bool anti) {
DEBUG_FUNC(momenta.size());
p_ci->GeneratePoint();
if (anti) {
for (size_t i(0);i<m_anticur.size();++i) {
m_anticur[i]->ConstructJ(i==0?-momenta[i]:momenta[i],0,0,0,0);
m_anticur[i]->ConstructJ(i==0?-momenta[i]:momenta[i],0,p_ci->I()[i],p_ci->J()[i],0);
m_anticur[i]->Print();
}
m_antiscur->Evaluate();
......@@ -165,7 +184,7 @@ void Comix1to3::Calculate(const ATOOLS::Vec4D_Vector& momenta, bool anti) {
}
else {
for (size_t i(0);i<m_cur.size();++i) {
m_cur[i]->ConstructJ(i==0?-momenta[i]:momenta[i],0,0,0,0);
m_cur[i]->ConstructJ(i==0?-momenta[i]:momenta[i],0,p_ci->I()[i],p_ci->J()[i],0);
m_cur[i]->Print();
}
m_scur->Evaluate();
......@@ -181,6 +200,10 @@ void Comix1to3::Calculate(const ATOOLS::Vec4D_Vector& momenta, bool anti) {
else {
m_fcur->Contract<double>(*m_cur.front(),fill,*this,0);
}
for (size_t i=0; i<size(); ++i) {
(*this)[i] *= sqrt(p_ci->GlobalWeight());
}
}
size_t Comix1to3::NHel(const Flavour& fl)
......
......@@ -5,6 +5,10 @@ namespace METOOLS {
class Vertex;
}
namespace PHASIC {
class Color_Integrator;
}
namespace EXTRAXS {
class Comix1to3 : public METOOLS::Spin_Amplitudes {
std::vector<METOOLS::Current*> m_cur;
......@@ -19,6 +23,8 @@ namespace EXTRAXS {
ATOOLS::Flavour m_prop;
PHASIC::Color_Integrator* p_ci;
size_t NHel(const ATOOLS::Flavour& fl);
METOOLS::Vertex* GetVertex(METOOLS::Current* cur1,
METOOLS::Current* cur2,
......
......@@ -4,10 +4,12 @@
#include "METOOLS/Explicit/Vertex.H"
#include "MODEL/Main/Model_Base.H"
#include "MODEL/Interaction_Models/Single_Vertex.H"
#include "PHASIC++/Main/Color_Integrator.H"
using namespace EXTRAXS;
using namespace ATOOLS;
using namespace METOOLS;
using namespace PHASIC;
using namespace std;
Comix1to2::Comix1to2(const vector<Flavour>& flavs) :
......@@ -68,6 +70,22 @@ Comix1to2::Comix1to2(const vector<Flavour>& flavs) :
m_antifcur->Print();
m_antifcur->HM().resize(m_n);
for (size_t i(0);i<m_n;++i) m_antifcur->HM()[i]=i;
p_ci=new Color_Integrator();
Idx_Vector cids(3,0);
Int_Vector acts(3,0), types(3,0);
for (size_t i(0);i<flavs.size();++i) {
cids[i]=i;
acts[i]=flavs[i].Strong();
if (acts[i]) {
if (flavs[i].StrongCharge()==8) types[i]=0;
else if (flavs[i].IsAnti()) types[i]=(i==0?1:-1);
else types[i]=(i==0?-1:1);
}
}
if (!p_ci->ConstructRepresentations(cids,types,acts)) {
THROW(fatal_error, "Internal error.");
}
}
......@@ -108,16 +126,17 @@ Vertex* Comix1to2::GetVertex(Current* cur1, Current* cur2, Current* prop) {
void Comix1to2::Calculate(const ATOOLS::Vec4D_Vector& momenta, bool anti) {
p_ci->GeneratePoint();
if (anti) {
for (size_t i(0);i<m_anticur.size();++i) {
m_anticur[i]->ConstructJ(i==0?-momenta[i]:momenta[i],0,0,0,0);
m_anticur[i]->ConstructJ(i==0?-momenta[i]:momenta[i],0,p_ci->I()[i],p_ci->J()[i],0);
m_anticur[i]->Print();
}
m_antifcur->Evaluate();
}
else {
for (size_t i(0);i<m_cur.size();++i) {
m_cur[i]->ConstructJ(i==0?-momenta[i]:momenta[i],0,0,0,0);
m_cur[i]->ConstructJ(i==0?-momenta[i]:momenta[i],0,p_ci->I()[i],p_ci->J()[i],0);
m_cur[i]->Print();
}
m_fcur->Evaluate();
......@@ -131,6 +150,10 @@ void Comix1to2::Calculate(const ATOOLS::Vec4D_Vector& momenta, bool anti) {
else {
m_fcur->Contract<double>(*m_cur.front(),fill,*this,0);
}
for (size_t i=0; i<size(); ++i) {
(*this)[i] *= sqrt(p_ci->GlobalWeight());
}
}
size_t Comix1to2::NHel(const Flavour& fl)
......
......@@ -5,6 +5,10 @@ namespace METOOLS {
class Vertex;
}
namespace PHASIC {
class Color_Integrator;
}
namespace EXTRAXS {
class Comix1to2 : public METOOLS::Spin_Amplitudes {
......@@ -18,6 +22,8 @@ namespace EXTRAXS {
std::vector<size_t> m_nhel;
PHASIC::Color_Integrator* p_ci;
size_t NHel(const ATOOLS::Flavour& fl);
METOOLS::Vertex* GetVertex(METOOLS::Current* cur1,
METOOLS::Current* cur2,
......
......@@ -5,7 +5,6 @@
#include "HADRONS++/Current_Library/Current_Base.H"
#include "HADRONS++/PS_Library/HD_PS_Base.H"
#include "PHASIC++/Decays/Decay_Table.H"
#include "PHASIC++/Decays/Color_Function_Decay.H"
#include "PHASIC++/Channels/Multi_Channel.H"
#include "ATOOLS/Org/MyStrStream.H"
#include "ATOOLS/Math/Vector.H"
......@@ -135,8 +134,7 @@ bool Hadron_Decay_Channel::Initialise(GeneralModel startmd)
vector<int> decayindices(n);
for(int i=0;i<n;i++) decayindices[i]=i;
HD_ME_Base* me=new Generic(m_physicalflavours,decayindices,"Generic");
PHASIC::Color_Function_Decay* col=new PHASIC::Color_Function_Decay();
AddDiagram(me, col);
AddDiagram(me);
AddPSChannel( string("Isotropic"), 1., m_startmd);
msg_Tracking()<<"Calculating width for "<<Name()<<":\n";
CalculateWidth();
......@@ -212,8 +210,7 @@ void Hadron_Decay_Channel::ProcessME( vector<vector<string> > me_svv,
Complex factor = Complex(ToType<double>(ip.Interprete(me_svv[i][0])),
ToType<double>(ip.Interprete(me_svv[i][1])));
me->SetFactor(factor);
PHASIC::Color_Function_Decay* col=new PHASIC::Color_Function_Decay();
AddDiagram(me, col);
AddDiagram(me);
nr_of_mes++;
}
if(me_svv[i].size()==4) {
......@@ -258,8 +255,7 @@ void Hadron_Decay_Channel::ProcessME( vector<vector<string> > me_svv,
me->SetCurrent1(current1);
me->SetCurrent2(current2);
me->SetFactor(factor);
PHASIC::Color_Function_Decay* col=new PHASIC::Color_Function_Decay();
AddDiagram(me, col);
AddDiagram(me);
nr_of_mes++;
}
}
......@@ -270,8 +266,7 @@ void Hadron_Decay_Channel::ProcessME( vector<vector<string> > me_svv,
vector<int> decayindices(n);
for(int i=0;i<n;i++) decayindices[i]=i;
HD_ME_Base* me=new Generic(m_physicalflavours,decayindices,"Generic");
PHASIC::Color_Function_Decay* col=new PHASIC::Color_Function_Decay();
AddDiagram(me, col);
AddDiagram(me);
}
}
......@@ -432,7 +427,7 @@ bool Hadron_Decay_Channel::SetColorFlow(ATOOLS::Blob* blob)
Particle_Vector outparts=blob->GetOutParticles();
if(m_diagrams.size()>0) {
// try if the matrix element knows how to set the color flow
HD_ME_Base* firstme=(HD_ME_Base*) m_diagrams[0].first;
HD_ME_Base* firstme=(HD_ME_Base*) m_diagrams[0];
bool anti=blob->InParticle(0)->Flav().IsAnti();
if(firstme->SetColorFlow(outparts,n_q,n_g,anti)) return true;
}
......@@ -559,7 +554,7 @@ void Hadron_Decay_Channel::LatexOutput(std::ostream& f, double totalwidth)
}
f<<"\\\\"<<endl;
if((m_diagrams.size()>0 &&
((HD_ME_Base*) m_diagrams[0].first)->Name()!="Generic")) {
((HD_ME_Base*) m_diagrams[0])->Name()!="Generic")) {
sprintf( helpstr, "%.4f", IWidth()/totalwidth*100. );
f<<" & "<<helpstr;
if( IDeltaWidth() > 0. ) {
......@@ -569,7 +564,7 @@ void Hadron_Decay_Channel::LatexOutput(std::ostream& f, double totalwidth)
f<<" \\% ";
}
for(size_t i=0;i<m_diagrams.size();i++) {
HD_ME_Base* me=(HD_ME_Base*) m_diagrams[i].first;
HD_ME_Base* me=(HD_ME_Base*) m_diagrams[i];
if(me->Name()=="Current_ME") {
Current_ME* cme=(Current_ME*) me;
f<<"\\verb;"<<cme->GetCurrent1()->Name()
......
......@@ -31,10 +31,10 @@ Amplitude2_Tensor::Amplitude2_Tensor(const std::vector<ATOOLS::Particle*>& parts
}
}
Amplitude2_Tensor::Amplitude2_Tensor(const vector<ATOOLS::Particle*>& parts,
size_t level,
const PHASIC::DiagColVec& diagrams,
const vector<vector<Complex> >& cmatrix,
const std::vector<Spin_Amplitudes*>& diagrams,
vector<int>& spin_i,
vector<int>& spin_j) :
p_next(NULL), m_value(-1.0,0.0), p_part(NULL), m_nhel(0)
......@@ -45,8 +45,8 @@ Amplitude2_Tensor::Amplitude2_Tensor(const vector<ATOOLS::Particle*>& parts,
m_value=Complex(0.0, 0.0);
for (size_t i(0); i<diagrams.size(); ++i) {
for (size_t j(0); j<diagrams.size(); ++j) {
m_value+=cmatrix[i][j]*diagrams[i].first->Get(spin_i)*
conj(diagrams[j].first->Get(spin_j));
m_value+=diagrams[i]->Get(spin_i)*
conj(diagrams[j]->Get(spin_j));
}
}
}
......@@ -61,23 +61,23 @@ Amplitude2_Tensor::Amplitude2_Tensor(const vector<ATOOLS::Particle*>& parts,
spin_i[level]=(i%m_nhel);
spin_j[level]=(i/m_nhel);
(*p_next)[i]=new Amplitude2_Tensor(parts, level+1,
diagrams, cmatrix, spin_i, spin_j);
diagrams, spin_i, spin_j);
}
}
}
Amplitude2_Tensor::Amplitude2_Tensor(const vector<ATOOLS::Particle*>& parts,
const vector<int>& permutation,
size_t level,
const vector<Spin_Amplitudes>& diagrams,
const vector<vector<Complex> >& cmatrix,
vector<int>& spin_i, vector<int>& spin_j) :
p_next(NULL), m_value(-1.0,0.0), p_part(NULL), m_nhel(0)
{
if (level>parts.size()) THROW(fatal_error, "Internal error 1");
if (level==parts.size() || parts[level]->RefFlav().IsStable()) {
m_value=ContractRemaining(parts,permutation,level,diagrams,cmatrix,
m_value=ContractRemaining(parts,permutation,level,diagrams,
spin_i,spin_j, 1.0);
}
else {
......@@ -91,17 +91,17 @@ Amplitude2_Tensor::Amplitude2_Tensor(const vector<ATOOLS::Particle*>& parts,
spin_i[level]=(i%m_nhel);
spin_j[level]=(i/m_nhel);
(*p_next)[i]=new Amplitude2_Tensor(parts, permutation, level+1,
diagrams, cmatrix, spin_i, spin_j);
diagrams, spin_i, spin_j);
}
}
}
Complex Amplitude2_Tensor::ContractRemaining
(const std::vector<ATOOLS::Particle*>& parts,
const vector<int>& permutation,
size_t level,
const vector<Spin_Amplitudes>& diagrams,
const std::vector<std::vector<Complex> >& cmatrix,
vector<int>& spin_i, vector<int>& spin_j, double factor) const
{
if (level>parts.size()) THROW(fatal_error, "Internal error 1");
......@@ -116,8 +116,8 @@ Complex Amplitude2_Tensor::ContractRemaining
}
for (size_t i(0); i<diagrams.size(); ++i) {
for (size_t j(0); j<diagrams.size(); ++j) {
ret+=cmatrix[i][j]*diagrams[i].Get(spin_i_perm)*
conj(diagrams[j].Get(spin_j_perm))*factor;
ret+=diagrams[i].Get(spin_i_perm)*
conj(diagrams[j].Get(spin_j_perm))*factor;
}
}
}
......@@ -129,7 +129,7 @@ Complex Amplitude2_Tensor::ContractRemaining
spin_i[level]=i;
spin_j[level]=i;
ret+=ContractRemaining(parts, permutation, level+1,
diagrams, cmatrix, spin_i, spin_j, newfactor);
diagrams, spin_i, spin_j, newfactor);
}
}
return ret;
......
......@@ -7,13 +7,6 @@
#include "METOOLS/Main/Spin_Structure.H"
#include "ATOOLS/Phys/Particle.H"
namespace PHASIC {
class Color_Function_Decay;
typedef
std::vector<std::pair<METOOLS::Spin_Amplitudes*,Color_Function_Decay*> >
DiagColVec;
}
namespace METOOLS {
class Decay_Matrix;
......@@ -28,7 +21,6 @@ namespace METOOLS {
const std::vector<int> &permutation,
size_t level,
const std::vector<Spin_Amplitudes> &diagrams,
const std::vector<std::vector<Complex> > &cmatrix,
std::vector<int> &spin_i, std::vector<int> &spin_j,
double factor) const;
......@@ -36,15 +28,13 @@ namespace METOOLS {
Amplitude2_Tensor(const std::vector<ATOOLS::Particle*>& parts, size_t level);
Amplitude2_Tensor(const std::vector<ATOOLS::Particle*>& parts,
size_t level,
const PHASIC::DiagColVec& diagrams,
const std::vector<std::vector<Complex> >& colormatrix,
const std::vector<Spin_Amplitudes*>& diagrams,
std::vector<int>& spin_i,
std::vector<int>& spin_j);
Amplitude2_Tensor(const std::vector<ATOOLS::Particle*>& parts,
const std::vector<int>& permutation,
size_t level,
const std::vector<Spin_Amplitudes>& diagrams,
const std::vector<std::vector<Complex> >& cmatrix,
std::vector<int>& spin_i, std::vector<int>& spin_j);
~Amplitude2_Tensor();
......
#include "PHASIC++/Decays/Color_Function_Decay.H"
#include "MODEL/Interaction_Models/Color_Function.H"
#include "ATOOLS/Phys/Color.H"
#include "ATOOLS/Org/Exception.H"
#include "ATOOLS/Org/MyStrStream.H"
#include <algorithm>
using namespace PHASIC;
using namespace ATOOLS;
using namespace std;
size_t Color_Function_Decay::m_nid(100);
Color_Function_Decay::Color_Function_Decay() : m_max(-1)
{
push_back(make_pair("1",vector<int>()));
}
Color_Function_Decay::Color_Function_Decay(const MODEL::Color_Function& c1) :
m_max(-1)
{
push_back(make_pair(string(""),vector<int>()));
std::string& name = back().first;
std::vector<int>& inds = back().second;
switch(c1.Type()) {
case MODEL::cf::T:
name="T";
for (size_t i(0); i<3; ++i)
inds.push_back(c1.ParticleArg(i));
break;
case MODEL::cf::F:
name="F";
for (size_t i(0); i<3; ++i)
inds.push_back(c1.ParticleArg(i));
break;
case MODEL::cf::D:
name="D";
for (size_t i(0); i<2; ++i)
inds.push_back(c1.ParticleArg(i));
break;
case MODEL::cf::G:
name="G";
for (size_t i(0); i<2; ++i)
inds.push_back(c1.ParticleArg(i));
break;
case MODEL::cf::None:
name="1";
break;
default:
THROW(not_implemented,"Encountered unknown type of color function: "+
ToString(c1.Type()));
}
for (size_t i(0); i<back().second.size(); ++i) {
if (back().second[i]>m_max) m_max=back().second[i];
}
}
void Color_Function_Decay::Conjugate()
{
for (size_t i=0; i<size(); ++i) {
if (at(i).first=="T") {
int help=at(i).second[1];
at(i).second[1]=at(i).second[2];
at(i).second[2]=help;
}
else if (at(i).first=="D") {
int help=at(i).second[0];
at(i).second[0]=at(i).second[1];
at(i).second[1]=help;
}
else if (at(i).first=="G") {
int help=at(i).second[0];
at(i).second[0]=at(i).second[1];
at(i).second[1]=help;
}
}
}
std::string Color_Function_Decay::String() const
{
std::string ret;
for (size_t i(0); i<size(); ++i) {
if (i>0) ret+="*";
if (at(i).first=="G") {
ret+="2*T["+ToString(at(i).second[0])+","
+ToString(m_nid)+","+ToString(m_nid+1)+"]"
+"*T["+ToString(at(i).second[1])+","
+ToString(m_nid+1)+","+ToString(m_nid)+"]";
m_nid+=2;
}
else {
ret+=at(i).first;
for (size_t j(0); j<at(i).second.size(); ++j) {
if (j==0) ret+="["+ToString(at(i).second[j]);
else ret+=","+ToString(at(i).second[j]);
if (j==at(i).second.size()-1) ret+="]";
}
}
}
return ret;
}
void Color_Function_Decay::BumpIndices(size_t i, int bump)
{
for (size_t j(0); j<at(i).second.size(); ++j) at(i).second[j]+=bump;
for (size_t j(0); j<at(i).second.size(); ++j) {
if (at(i).second[j]>m_max) m_max=at(i).second[j];
}
}
std::vector<int> Color_Function_Decay::Multiply(const Color_Function_Decay& c)
{
DEBUG_FUNC(String()<<" * "<<c.String());
// returns a vector of what has been added to the indices of each color
// function in c to make the indices unique.
std::vector<int> ret(c.size());
for (size_t i(0); i<c.size(); ++i) {
push_back(make_pair(c[i].first, c[i].second));
// make indices unique
ret[i]=m_max+1;
BumpIndices(size()-1, ret[i]);
}
return ret;
}
void Color_Function_Decay::Contract(int i1, int i2)
{
DEBUG_FUNC(i1<<" "<<i2);
int has_i1(HasIndex(i1)), has_i2(HasIndex(i2));
if (has_i1==0 || has_i2==0) return;
else if (has_i1==1 && has_i2==1) {
for (size_t i(0); i<size(); ++i) {
for (size_t j(0); j<at(i).second.size(); ++j) {
if (at(i).second[j]==i2) at(i).second[j]=i1;
}
}
m_internal.push_back(i1);
}
else
THROW(fatal_error, "Internal error.");
}
void Color_Function_Decay::ReplaceIndex(int i1, int i2)
{
for (size_t i(0); i<size(); ++i) {
for (size_t j(0); j<at(i).second.size(); ++j) {
if (at(i).second[j]==i1) at(i).second[j]=i2;
}
}
}
int Color_Function_Decay::HasIndex(int i1) const
{
int ret(0);
for (size_t i(0); i<size(); ++i) {
for (size_t j(0); j<at(i).second.size(); ++j) {
if (at(i).second[j]==i1) ++ret;
}
}
return ret;
}
Complex Color_Function_Decay::Contract(const Color_Function_Decay& c)
{
DEBUG_FUNC(String()+" * "+c.String());
DEBUG_VAR(this->m_max);
DEBUG_VAR(c.m_max);
double max=std::max(m_max, c.m_max);
Color_Function_Decay c1(c);
for (size_t i(0); i<c1.Internal().size(); ++i) {
DEBUG_VAR(i);
DEBUG_VAR(c1.m_internal[i]);
if (c1.m_internal[i]<=max) continue;
c1.ReplaceIndex(c1.m_internal[i], max+1+i);
}
c1.Conjugate();
DEBUG_VAR(String());
DEBUG_VAR(c1.String());
Expression expr(String()+"*"+c1.String());
if (msg_LevelIsDebugging()) expr.Print();
if (expr.Evaluate()) {
DEBUG_VAR(expr.Result());
return expr.Result();
}
else {
THROW(fatal_error, "Color factor evaluation failed: "+String()+
" * "+c1.String());
}
return Complex(1.0,0.0);
}
namespace PHASIC {
std::ostream &operator<<(std::ostream &os,const Color_Function_Decay &c)
{
os<<c.String();
return os;
}
}
#ifndef Color_Function_Decay_h
#define Color_Function_Decay_h
#include "ATOOLS/Math/MyComplex.H"
#include <string>
#include <vector>
namespace MODEL {
class Color_Function;
}
namespace PHASIC {
class Color_Function_Decay :
public std::vector<std::pair<std::string,std::vector<int> > >
{
std::vector<int> m_internal;
int m_max;
static size_t m_nid;
public:
Color_Function_Decay();
Color_Function_Decay(const MODEL::Color_Function& c1);
void Conjugate();
std::string String() const;
std::vector<int> Multiply(const Color_Function_Decay& c);
void Contract(int c1, int c2);
void ReplaceIndex(int c1, int c2);
void BumpIndices(size_t i, int bump);
int HasIndex(int i) const;
Complex Contract(const Color_Function_Decay& c);
inline const std::vector<int>& Internal() const { return m_internal; }
friend std::ostream &operator<<(std::ostream &os,
const Color_Function_Decay &);
};
}
#endif
......@@ -9,13 +9,11 @@
#include "METOOLS/Main/Spin_Structure.H"
#include "METOOLS/SpinCorrelations/Amplitude2_Tensor.H"
#include "METOOLS/SpinCorrelations/Spin_Density.H"
#include "PHASIC++/Decays/Color_Function_Decay.H"
#include <algorithm>
using namespace PHASIC;
using namespace ATOOLS;
using namespace METOOLS;
using namespace MODEL;
using namespace std;
Decay_Channel::Decay_Channel(const Flavour & _flin,
......@@ -30,8 +28,7 @@ Decay_Channel::Decay_Channel(const Flavour & _flin,
Decay_Channel::~Decay_Channel()
{
for (size_t i(0); i<m_diagrams.size(); ++i) {
delete m_diagrams[i].first;
delete m_diagrams[i].second;
delete m_diagrams[i];
}
if (p_channels) delete p_channels;
if (p_amps) delete p_amps;
......@@ -76,28 +73,8 @@ void Decay_Channel::AddDecayProduct(const ATOOLS::Flavour& flout,
m_minmass += p_ms->Mass(flout);
}
void Decay_Channel::AddDiagram(METOOLS::Spin_Amplitudes* amp,
Color_Function_Decay* col) {
DEBUG_FUNC(*col);
m_diagrams.push_back(make_pair(amp, col));
size_t index=m_diagrams.size()-1;
DEBUG_INFO("Add one column to all existing rows");
if (m_colormatrix.size()!=index) THROW(fatal_error,"Wrong size of cols.");
for (size_t i=0; i<m_colormatrix.size(); ++i) {
if (m_colormatrix[i].size()!=index) THROW(fatal_error,"Wrong size of cols");
DEBUG_INFO("Contracting "<<*m_diagrams[i].second<<" with "<<*col);
m_colormatrix[i].push_back(m_diagrams[i].second->Contract(*col));
DEBUG_VAR(m_colormatrix[i].back());
}
DEBUG_INFO("Add an additional row");
m_colormatrix.resize(m_colormatrix.size()+1);
for (size_t i=0; i<m_diagrams.size(); ++i) {
DEBUG_INFO("Contracting "<<*col<<" with "<<*m_diagrams[i].second);
m_colormatrix.back().push_back(col->Contract(*m_diagrams[i].second));
DEBUG_VAR(m_colormatrix.back().back());
}
void Decay_Channel::AddDiagram(METOOLS::Spin_Amplitudes* amp) {
m_diagrams.push_back(amp);
}
void Decay_Channel::AddChannel(PHASIC::Single_Channel* chan)
......@@ -127,11 +104,6 @@ namespace PHASIC {
if (dc.Active()<1) {
os<<" [disabled]";
}
if (msg_LevelIsTracking()) {
for (size_t i(0); i<dc.GetDiagrams().size(); ++i) {
os<<" "<<setw(10)<<*dc.GetDiagrams()[i].second;
}
}
return os;
}
}
......@@ -241,8 +213,7 @@ void Decay_Channel::CalculateWidth()
long int n=0;
int opt=0;
double value, oldvalue=0., sum=0., sum2=0., result=1., disc;
bool simple=false;
double value, sum=0., sum2=0., result=1., disc;
m_ideltawidth=1.0;
std::vector<Vec4D> momenta(1+NOut());
......@@ -256,14 +227,11 @@ void Decay_Channel::CalculateWidth()
if (value>m_max) {
m_max = value;
}
if (value!=0. && value==oldvalue) { simple = true; break; }