#include "AMEGIC++/DipoleSubtraction/Single_LOProcess_MHV.H" #include "ATOOLS/Org/Run_Parameter.H" #include "AMEGIC++/DipoleSubtraction/Single_LOProcess_MHV.H" #include "AMEGIC++/Amplitude/FullAmplitude_MHV_Base.H" #include "MODEL/Main/Running_AlphaS.H" #include "PHASIC++/Main/Process_Integrator.H" #include "PHASIC++/Main/Phase_Space_Handler.H" #include "AMEGIC++/Phasespace/Phase_Space_Generator.H" #include "BEAM/Main/Beam_Spectra_Handler.H" #include "PDF/Main/ISR_Handler.H" #include "ATOOLS/Org/Run_Parameter.H" #include "ATOOLS/Org/Shell_Tools.H" #include "ATOOLS/Org/MyStrStream.H" #include using namespace AMEGIC; using namespace MODEL; using namespace PHASIC; using namespace PDF; using namespace BEAM; using namespace ATOOLS; using namespace std; /*------------------------------------------------------------------------------- Constructors ------------------------------------------------------------------------------- */ Single_LOProcess_MHV::Single_LOProcess_MHV(const Process_Info &pi, BEAM::Beam_Spectra_Handler *const beam, PDF::ISR_Handler *const isr, const ATOOLS::sbt::subtype& st) : Single_LOProcess(pi, beam, isr, st) { m_ownamps = false; m_emitgluon = false; } Single_LOProcess_MHV::~Single_LOProcess_MHV() { if (p_hel) {delete p_hel; p_hel=0;} if (p_BS) {delete p_BS; p_BS=0;} if (p_shand) {delete p_shand;p_shand=0;} if (p_ampl) {delete p_ampl; p_ampl=0;} #ifndef Basic_Sfuncs_In_MHV if (p_momlist) {delete p_momlist; p_momlist=0;} #endif if (p_MHVamp && m_ownamps) {delete p_MHVamp; p_MHVamp=0;} } /*------------------------------------------------------------------------------ Initializing libraries, amplitudes, etc. ------------------------------------------------------------------------------*/ int Single_LOProcess_MHV::InitAmplitude(Amegic_Model * model,Topology* top, vector & links, vector & errs,int checkloopmap) { m_type = 21; if (!model->p_model->CheckFlavours(m_nin,m_nout,&m_flavs.front())) return 0; model->p_model->GetCouplings(m_cpls); m_partonlistqcd.clear(); m_partonlistqed.clear(); if (m_stype&sbt::qcd) { for (size_t i=0;iSetk0(s_gauge); //////////////////////////////////////////////// #ifdef Basic_Sfuncs_In_MHV p_momlist = p_BS; #else p_momlist = new MomentumList(m_nin,m_nout); #endif Flavour* fl = &m_flavs.front(); int *plist = new int[m_nin+m_nout]; for (size_t i=0;ip_model,&m_cpls,m_nin+m_nout,plist,p_momlist,m_ownamps,127,127); delete [] plist; ////////////////////////////////////////////// p_shand = new String_Handler(m_gen_str,p_BS,model->p_model->GetCouplings()); p_ampl = new Amplitude_Handler(m_nin+m_nout,fl,p_b,p_pinfo,model,top,m_maxcpl,m_mincpl, m_pinfo.m_ntchan,m_pinfo.m_mtchan, &m_cpls,p_BS,p_shand,m_print_graphs,0,true,m_ptypename+"/"+m_libname); if (p_ampl->GetGraphNumber()==0) { msg_Tracking()<<"Single_LOProcess_MHV::InitAmplitude : No diagrams for "<PossibleConfigsExist(m_maxcpl,m_mincpl)) { msg_Tracking()<<"Single_LOProcess::InitAmplitude : No possible combinations exist for "< cplmap; for (size_t j=0;jType()) { cplmap.clear(); if (m_allowmap && FlavCompare(links[j]) && p_ampl->CompareAmplitudes(links[j]->GetAmplitudeHandler(),m_sfactor,cplmap)) { if (p_hel->Compare(links[j]->GetHelicity(),m_nin+m_nout)) { m_sfactor = sqr(m_sfactor); msg_Tracking()<<"AMEGIC::Single_Process_MHV::InitAmplitude : Found compatible process for "<Name()<Result()*m_sfactor; InitFlavmap(p_partner); Minimize(); return 1; } } } p_ampl->FillPointlist(); p_BS->Initialize(); switch (Tests()) { case 1 : if (p_partner==this) links.push_back(this); return 1; case -3: return 0; default : msg_Error()<<"ERROR in Single_Fin_Process_MHV::InitAmplitude : "< & links, vector & errs, std::vector* epol,std::vector * pfactors) { m_type = 11; if (!model->p_model->CheckFlavours(m_nin,m_nout,&m_flavs.front())) return 0; model->p_model->GetCouplings(m_cpls); int cnt=0; vector fi_tags; m_pinfo.m_fi.GetTags(fi_tags); if (fi_tags.size()!=m_nout) THROW(fatal_error, "Internal error."); for (size_t i(0);iSetk0(s_gauge); p_epol = epol; //////////////////////////////////////////////// #ifdef Basic_Sfuncs_In_MHV p_momlist = p_BS; #else p_momlist = new MomentumList(m_nin,m_nout); #endif int *plist = new int[m_nin+m_nout]; for (size_t i=0;ip_model,&m_cpls,m_nin+m_nout,plist,p_momlist,m_ownamps,m_emit,m_spect); delete [] plist; ////////////////////////////////////////////// p_shand = new String_Handler(m_gen_str,p_BS,model->p_model->GetCouplings()); p_ampl = new Amplitude_Handler(m_nin+m_nout,&m_flavs.front(),p_b,p_pinfo,model,top,m_maxcpl,m_mincpl, m_pinfo.m_ntchan,m_pinfo.m_mtchan, &m_cpls,p_BS,p_shand,m_print_graphs,0,true,m_ptypename+"/"+m_libname); if (p_ampl->GetGraphNumber()==0) { msg_Tracking()<<"Single_LOProcess_MHV::InitAmplitude : No diagrams for "<PossibleConfigsExist(m_maxcpl,m_mincpl)) { msg_Tracking()<<"Single_LOProcess::InitAmplitude : No possible combinations exist for "< cplmap; for (size_t j=0;jType()) { cplmap.clear(); if (m_allowmap && FlavCompare(links[j]) && p_ampl->CompareAmplitudes(links[j]->GetAmplitudeHandler(),m_sfactor,cplmap)) { if (p_hel->Compare(links[j]->GetHelicity(),m_nin+m_nout)) { Single_LOProcess_MHV *pp=dynamic_cast(links[j]); if (m_emit!=pp->m_emit || m_spect!=pp->m_spect || p_sub->m_ijt!=pp->p_sub->m_ijt || p_sub->m_kt!=pp->p_sub->m_kt || p_sub->m_i!=pp->p_sub->m_i || p_sub->m_j!=pp->p_sub->m_j || p_sub->m_k!=pp->p_sub->m_k) continue; m_sfactor = sqr(m_sfactor); msg_Tracking()<<"AMEGIC::Single_Process_MHV::InitAmplitude : Found compatible process for "<Name()<Result()*m_sfactor; InitFlavmap(p_partner); Minimize(); return 1; } } } p_ampl->FillPointlist(); p_BS->Initialize(); int tr=Tests(pfactors); // PRINT_INFO("Tests Result: "< * pfactors) { int number = 1; int gauge_test = 1; int string_test = 1; /* --------------------------------------------------- The reference result for momenta moms --------------------------------------------------- */ string testname = string(""); if (FoundMappingFile(testname,m_pslibname)) { if (testname != string("")) { gauge_test = string_test = 0; } } double M2 = 0.; double helvalue; if (gauge_test) { #ifdef Basic_Sfuncs_In_MHV p_BS->Setk0(0); p_BS->CalcEtaMu(p_testmoms); #else p_momlist->PutMomenta(p_testmoms); #endif msg_Tracking()<<"Single_LOProcess_MHV::Tests for "<SetSqMatrix((*pfactors)[1],CalculateHelicityPhase(p_testmoms)); for (size_t i=0;iMaxHel();i++) { if (p_hel->On(i) && p_hel->GetEPol(i)==90) { helvalue = p_MHVamp->MSquare((*p_hel)[i],p_BS)*p_hel->PolarizationFactor(i); M2 += helvalue; } } M2 *= p_MHVamp->ParticlesNorm(); m_iresult = M2; } /* --------------------------------------------------- First test : gauge test --------------------------------------------------- */ #ifdef Basic_Sfuncs_In_MHV p_BS->Setk0(s_gauge); p_BS->CalcEtaMu(p_testmoms); #else p_momlist->PutMomenta(p_testmoms); #endif number++; double M2g = 0.; double * M_doub = new double[p_hel->MaxHel()]; for (size_t i=0; iMaxHel(); ++i) M_doub[i]=0.; if (m_emitgluon) p_MHVamp->SetSqMatrix((*pfactors)[1],CalculateHelicityPhase(p_testmoms)); for (size_t i=0; iMaxHel(); ++i) { if (p_hel->On(i) && p_hel->GetEPol(i)==90) { M_doub[i] = p_MHVamp->MSquare((*p_hel)[i],p_BS)*p_hel->PolarizationFactor(i); M2g += M_doub[i]; } } M2g *= p_MHVamp->ParticlesNorm(); m_iresult = M2g; p_BS->StartPrecalc(); if (gauge_test) { if (!ATOOLS::IsEqual(M2,M2g)) { msg_Out()<<"WARNING: Gauge test not satisfied: " <MaxHel();i++) { if (p_hel->On(i)) { for (size_t j=i+1;jMaxHel();j++) { if (p_hel->On(j)) { #ifdef FuckUp_Helicity_Mapping if (ATOOLS::IsEqual(M_doub[i],M_doub[j])) { p_hel->SwitchOff(j); p_hel->SetPartner(i,j); p_hel->IncMultiplicity(i); } #endif } } } } delete[] M_doub; return 1; } Complex Single_LOProcess_MHV::CalculateHelicityPhase(const ATOOLS::Vec4D * mom) { #ifndef Basic_Sfuncs_In_MHV msg_Error()<<"Must use Basic_Sfuncs to generate dipole subtraction terms!"<Getk0())); Complex eta13=2.*csqrt(p*p_BS->Getk0())*csqrt(q*p_BS->Getk0()); Complex eip1 = Complex(0.,1.)*(-p_BS->CalcS(q,e1)*conj(p_BS->CalcS(e1,p))+mu12*eta13)/p_BS->CalcS(q,p); return eip1; } double Single_LOProcess_MHV::operator()(const ATOOLS::Vec4D_Vector &labmom,const ATOOLS::Vec4D *mom, std::vector * pfactors,std::vector* epol,const int mode) { if (p_partner!=this) { return m_lastxs = p_partner->operator()(labmom,mom,pfactors,epol,mode)*m_sfactor; } p_int->SetMomenta(labmom); p_scale->CalculateScale(labmom,m_cmode); #ifdef Basic_Sfuncs_In_MHV p_BS->CalcEtaMu((ATOOLS::Vec4D*)mom); #else p_momlist->PutMomenta(mom); #endif double M2(0.); if (m_emitgluon) p_MHVamp->SetSqMatrix((*pfactors)[1],CalculateHelicityPhase(mom)); for (size_t i=0;iMaxHel();i++) { if (p_hel->On(i) && p_hel->GetEPol(i)==90) { double mh=p_MHVamp->MSquare((*p_hel)[i],p_BS) * p_hel->Multiplicity(i) * p_hel->PolarizationFactor(i); M2 += mh; } } m_lastxs = M2*p_MHVamp->ParticlesNorm(); return m_lastxs; } void Single_LOProcess_MHV::Calc_AllXS(const ATOOLS::Vec4D_Vector &labmom, const ATOOLS::Vec4D *mom, std::vector > &dsijqcd, std::vector > &dsijew, const int mode) { p_int->SetMomenta(labmom); p_scale->CalculateScale(labmom,m_cmode); #ifdef Basic_Sfuncs_In_MHV p_BS->CalcEtaMu((ATOOLS::Vec4D*)mom); #else p_momlist->PutMomenta(mom); #endif dsijqcd[0][0]=0.; for (size_t i=0;iMaxHel();i++) { if (p_hel->On(i)) { double fac = p_hel->Multiplicity(i) * p_hel->PolarizationFactor(i) * p_MHVamp->ParticlesNorm(); p_MHVamp->CalculateAmps((*p_hel)[i],p_BS); dsijqcd[0][0] += p_MHVamp->MSquare(0,0)*fac; for (size_t i=0;iMSquare(m_partonlistqcd[i],m_partonlistqcd[k])*fac; } } } } // PRINT_INFO(Name()); // cout<<"Result 00: "<