#include "AMEGIC++/Amplitude/FullAmplitude_External.H" #include "AMEGIC++/Amplitude/Zfunctions/Basic_Sfuncs.H" #include "AMEGIC++/Main/ColorSC.H" #include "AMEGIC++/Main/Helicity.H" #include "PHASIC++/Process/Tree_ME2_Base.H" #include "ATOOLS/Phys/Color.H" #include "ATOOLS/Org/Exception.H" #include "ATOOLS/Org/Run_Parameter.H" #include "ATOOLS/Org/MyStrStream.H" #include "ATOOLS/Math/MathTools.H" using namespace AMEGIC; using namespace PHASIC; using namespace MODEL; using namespace ATOOLS; FullAmplitude_External::FullAmplitude_External (const Process_Info &info,Model_Base *model,Coupling_Map *const cpls, Helicity *hel,const size_t &emit,const size_t &spect): p_model(model), m_emit(emit), m_spect(spect) { DEBUG_FUNC(""); p_calc = Tree_ME2_Base::GetME2(info); if (p_calc==NULL) return; msg_Debugging()< "<SetCouplings(*cpls); m_oqcd=p_calc->OrderQCD(); m_oew=p_calc->OrderEW(); m_nin=info.m_ii.NExternal(); m_fls=info.ExtractFlavours(); for (size_t i(0);iGetFlavourHelicityMap(); BuildHelicityMap(hel); BuildColorMatrix(); } FullAmplitude_External::~FullAmplitude_External() { delete p_calc; } double FullAmplitude_External::Calc(const ATOOLS::Vec4D *_p) { DEBUG_FUNC(""); Vec4D_Vector p(m_fls.size()); for (size_t i(0);iCalc(p)); #ifdef DEBUG__External for (size_t i(0);iNAmps();++id) { const std::vector &s(p_calc->GetAmplitudes(id)); for (size_t k=0;kGetPhase(id)*std::conj(amps[j]); amp*=m_colfs[m_emit][m_spect][i][j]; msg_Debugging()<0?'+':'-') <<", id = "<GetPhase(id)<<"\n"; if (m_pmap[h][k]<0) { if (m_pmap[h][l]<0) ampsqmm+=amp; else ampsqmp+=amp; } else { if (m_pmap[h][l]>0) ampsqpp+=amp; else ampsqpm+=amp; } } } } Complex sum(ampsqpp); if (m_emit==m_spect) { msg_Debugging()<<"AA* = "< "<GetHelicityPhase(pijt,eps1); } double FullAmplitude_External::MSquare (const size_t &h,const size_t &ci,const size_t &cj) { DEBUG_FUNC("h="<"<NAmps();++id) { const std::vector &s(p_calc->GetAmplitudes(id)); for (size_t k=0;kGetPhase(id)*std::conj(amps[j]); amp*=m_colfs[ci][cj][i][j]; msg_Debugging()<0?'+':'-') <<", id = "<GetPhase(id)<<"\n"; sum+=amp; } } } msg_Debugging()<<"AA* = "<MaxHel()); m_pmap.resize(hel->MaxHel()); for (size_t k=0;k chel(m_fls.size()); for (size_t j=0;jMaxHel();++i) { int epol(hel->GetEPol(i)==90?1:-1); std::vector rhel((*hel)[i],&(*hel)[i][chel.size()]); for (size_t l=0;l=90) epol=rhel[l]=chel[l]; if (chel==rhel) { msg_Debugging()<<" map "< "< " < "< &ids,std::vector &cid, Flavour_Vector &cfl,int &nsub,int &psub,int &swap) const { for (size_t i(0);i fid(cid); for (int oswap(-1);oswap(fid[i],fid[i-1]); if (m_fls[fid[i]].IsFermion() && m_fls[fid[i]]==m_fls[fid[i-1]]) ++swap; } } } void FullAmplitude_External::BuildColorMatrix() { DEBUG_FUNC(m_amap.size()); m_colfs.resize(m_fls.size()); for (size_t i(0);i (m_amap.size(),Complex(0.0,0.0))); for (size_t a=0;a &cp(m_amap[a].m_perm); int nsuba, psuba, swapa; Flavour_Vector fla; std::vector ida, perma(cp.begin(),cp.end()); GetPermutation(perma,ida,fla,nsuba,psuba,swapa); msg_Debugging()<<"A: "< &cp(m_amap[b].m_perm); int nsubb, psubb, swapb; Flavour_Vector flb; std::vector idb, permb(cp.begin(),cp.end()); GetPermutation(permb,idb,flb,nsubb,psubb,swapb); msg_Debugging()<<"B: "<TR()); size_t ad(expression.AIndex()); size_t lf(expression.FIndex()), fl(lf), qc(0), tf(0); for (size_t i=0;i0) lf=expression.FIndex(); expression.push_back(Delta::New(ida[i]+1,lf)); if (ci!=cj) { if (ci==ida[i]) { msg_Debugging()<<"|> 3-insertion on ci = "<0 && tf==1) { if (fla[i-1].StrongCharge()!=-3) THROW(fatal_error,"Invalid permutation"); msg_Debugging()<<"|> connect quark lines i = "< connect quark lines i = "< ~3-insertion on ci = "< 8-insertion on ci = "<0) lf=expression.FIndex(); expression.push_back(Delta::New(lf,idb[i]+1)); if (ci!=cj) { if (cj==idb[i]) { msg_Debugging()<<"<| 3-insertion on cj = "<0 && tf==1) { if (flb[i-1].StrongCharge()!=-3) THROW(fatal_error,"Invalid permutation"); msg_Debugging()<<"<| connect quark lines i = "<