#include "AMEGIC++/Phasespace/Channel_Generator.H" #include "AMEGIC++/Main/Topology.H" #include "ATOOLS/Phys/Flavour.H" #include "AMEGIC++/Main/Point.H" #include "ATOOLS/Org/Message.H" #include "ATOOLS/Math/MathTools.H" #include "ATOOLS/Org/MyStrStream.H" #include "ATOOLS/Org/Run_Parameter.H" #include "ATOOLS/Org/My_MPI.H" #include #include #include using namespace AMEGIC; using namespace PHASIC; using namespace ATOOLS; using namespace std; Channel_Generator::Channel_Generator(int _nin,int _nout, Point * _plist,int _aid) : Channel_Generator_Base(_nin,_nout,_plist) { IdentifyProps(plist); m_idstr = string(""); m_aid = _aid; } Channel_Generator::~Channel_Generator() { } int Channel_Generator::MakeChannel(int& echflag,int n,string& path,string& pID) { if (m_idstr==string("")) CreateChannelID(echflag); int oldflag = echflag; extrachannelflag = echflag; //add Channel char name[22]; snprintf(name,22,"C%i_%i",nout,n); if (echflag!=0) snprintf(name,22,"%s%c",name,'a'+extrachannelflag-1); string filename = rpa->gen.Variable("SHERPA_CPP_PATH")+string("/Process/Amegic/")+path+string("/")+ string(name)+string(".C"); int rannum = 0; ofstream chf; chf.open((filename).c_str()); chf<<"// Channel_Generator"<0) chf <<" double m_alpha,m_ctmax,m_ctmin;"<0) { chf <<" Info_Key "; bool first=true; for (size_t i=0; iMPISync(); }"<Optimize(); } "<EndOptimize(); } "<WriteOut(pId); } "<ReadIn(pId); } "<GeneratePoint(_ran);"<GenerateWeight(p_rans);"<0) { chf <<" m_alpha = .5;"<AddPoint(Value,p_rans);"<left; if (ph->left==0) { ph = p->right; if (ph->left==0 && p->middle) ph = p->middle; } if (ph == 0) { msg_Error()<1 process. " <fl<<" -> { "<left->fl<<" "<right->fl<<" }." <<" Aborting."<fl.Mass())) { help+=string("ZR")+ToString(ph->fl.Kfcode())+string("_"); } else help+=string("ZS_"); help+=ToString((int)PMassSum(ph,0)); m_idc.push_back(help); } if (flag<2) { m = LinkedMasses(p->left); if (m.length()<2) m = LinkedMasses(p->right); AddToVariables(flag,m,string("p[0] + p[1]"),1,sf); AddToVariables(flag,m,string("dabs(p")+Order(m)+string(".Abs2())"),0,sf); flag += 10; if (!StepS(flag,p->left,rannum,sf,flav,maxnumb)) { if (!StepS(flag,p->right,rannum,sf,flav,maxnumb)) { if (p->middle == 0) { msg_Error()<1 process. " <fl<<" -> { "<left->fl<<" "<right->fl<<" }." <<" Aborting."<left; if (ph->left==0) { ph = p->right; if (ph->left==0 && p->middle) ph = p->middle; if (ph->left==0) { msg_Error()<1 process. " <fl<<" -> { "<left->fl<<" "<right->fl<<" }." <<" Aborting."<fl.Kfcode()<<")).Mass();"<fl.Kfcode()<<")).Width();"<left==0) return 0; Point* l = p->left; Point* r = p->right; string lm = LinkedMasses(l); string rm = LinkedMasses(r); string mummy = Order(lm+rm); string moml,momr,momlalt,momralt; //Minima if (l->left==0) moml = string("p[") + GetMassIndex(lm) + string("]"); else moml = string("p") + Order(lm); if (r->left==0) momr = string("p[") + GetMassIndex(rm) + string("]"); else momr = string("p") + Order(rm); momlalt = string("p")+Order(mummy)+string("-")+momr; momralt = string("p")+Order(mummy)+string("-")+moml; Point** _plist = new Point*[2]; _plist[0] = p->left; _plist[1] = p->right; // Count resonating props, massless s-channels cannot be resonant.... for (short int i=0;i<2;i++) { if ((_plist[i]->left !=0) && (_plist[i]->fl.IsMassive())) { flav[maxnumb] = _plist[i]->fl; maxnumb++; } } GenerateMasses(flag,_plist,2,rannum,sf); delete[] _plist; bool first = 0; if (flag>9 || flag==-1) { first = 1; flag -= 10; } // Check for decay type. if ((!first) && (l->left==0) && (l->fl.IsVector()) && (!(l->fl.IsMassive())) && (r->fl.IsFermion()) && m_aid) { switch (flag) { case -11: m_idc.push_back(string("AI_")+Order(lm)+string("_")+Order(rm)); break; case 0: sf<<" CE.Anisotropic2Momenta(p"<left==0) && (r->fl.IsVector()) && (!(r->fl.IsMassive())) && (l->fl.IsFermion()) && m_aid) { //anisotropic decay for left outgoing massless vectors switch (flag) { case -11: m_idc.push_back(string("AI_")+Order(rm)+string("_")+Order(lm)); break; case 0: sf<<" CE.Anisotropic2Momenta(p"<number) < (r->number)) { switch (flag) { case -11: m_idc.push_back(string("I_")+Order(lm)+string("_")+Order(rm)); break; case 0: sf<<" CE.Isotropic2Momenta(p"<left==0) return; Point **props = new Point*[tcount+1]; Point **propt = new Point*[tcount+1]; int counts = 0; SetProps(p,props,propt,counts); if (flag==2) { sf<<" type = 2;"<fl.Kfcode()<<")).Mass() + "; } sf<<"Flavour((kf_code)("<fl.Kfcode()<<")).Mass();"<fl;maxnumb++; } m = Order(m); AddToVariables(flag,m,string("p[0] + p[1]"),1,sf); AddToVariables(flag,m,string("dabs(p")+Order(m)+string(".Abs2())"),0,sf); if (flag==-1||flag==-11) { double sms=0.; for (int i=0;i pin0;pin0.push_back(string("0")); vector pin1;pin1.push_back(string("1")); vector pout0; vector pout1; int count = 0; if (props[tcount]->number<99 && (props[tcount]->fl.IsVector() || props[0]->number>99 || !(props[0]->fl.Strong() || (props[0]->fl.IsLepton() && props[0]->fl.IntCharge()!=0)) )) { for (short int i=0;i pin0, vector pin1, vector pout0,vector pout1) { int hi = propt[count]->fl.Kfcode(); string tmstr; if (flag>=0 && !ATOOLS::IsZero(propt[count]->fl.Mass())) { tmstr = string("tmass")+ToString(count); sf<<" double "<=0) { //mins if (pout0.size()>1) CalcTSmin(flag,pout0,sf); if (pout1.size()>1) CalcTSmin(flag,pout1,sf); //maxs if (pout0.size()>1) { string s = string("sqr(sqrt(dabs((p")+Order(pin0sum)+string("+p")+ Order(pin1sum)+string(").Abs2()))-"); if (pout1.size()==1) { if (pout1[0].length()==1) s += string("sqrt(p_ms[")+GetMassIndex(pout1[0])+string("])"); else s += string("sqrt(s")+Order(pout1[0])+string(")"); } else s += string("sqrt(s")+Order(pout1sum)+string("_min)"); s += ")"; AddToVariables(flag,Order(pout0sum) + string("_max"),s,0,sf); //sf<<" double s"<1) { switch (flag) { case -11: if (pout1.size()==1 && pin1.size()==1 && pout1[0].length()==1 && extrachannelflag==1) { m_idc.push_back(string("LLP_")+Order(pout0sum)); extrachannelflag = 0; } else { m_idc.push_back(string("MlP_")+Order(pout0sum)); } break; case 0: sf<<" Vec4D p"<>i; } if (pout1.size()==1 && pin1.size()==1 && pout1[0].length()==1 && extrachannelflag==1) { sf<<" = CE.LLPropMomenta(0.99,1.001*sqr(rpa->gen.Ecms()),s"<gen.Ecms()),s"<1) { if (flag>=0) { string s = string("sqr(sqrt(dabs((p")+Order(pin0sum)+string("+p")+ Order(pin1sum)+string(").Abs2()))-"); if (pout0.size()==1) { if (pout0[0].size()==1) s += string("sqrt(p_ms[")+GetMassIndex(pout0[0])+string("])"); else s += string("sqrt(s")+Order(pout0[0])+string(")"); } else s += string("sqrt(s")+Order(pout1sum)+string(")"); s += string(")"); AddToVariables(flag,pout1sum + string("_max"),s,0,sf); //sf<<" double s"<fl.Mass())) his = ToString(hi); m_idc.push_back(string("TC")+his+ string("_")+Order(pin0sum)+string("_")+Order(pin1sum)+ string("_")+Order(pout0sum)+string("_")+Order(pout1sum)); break; case 0: sf<<" CE.TChannelMomenta(p"<fl.Mass())) his = ToString(hi); string idh = string("TC")+his+ string("_")+Order(pin0sum)+string("_")+Order(pin1sum)+ string("_")+Order(pout0sum)+string("_")+Order(pout1sum); //sf<<" std::cout<<\""<1 || pout1.size()>1) { if (pout0.size()>1) { //pout1 -> pin1 if (pout1.size()==1) { pin1.push_back(pout1[0]); string pin1sum = string("1_"); for (size_t i=1;i1) { //pout0 -> pin0 if (pout0.size()==1) { pin0.push_back(pout0[0]); string pin0sum = string("0_"); for (size_t i=1;ileft==0) { if (flag==0 || flag==10) AddToVariables(flag,lm[i],string("p_ms[")+GetMassIndex(lm[i])+string("]"),0,sf); //sf<<" double s"<0) { //CalcSmin(flag,"restmin",help,sf,0); //sum_s_i = string("-sqrt(s")+Order(help)+string("_restmin)"); CalcSmin(flag,"min",help,sf,0); sum_s_i = string("-sqrt(s")+Order(help)+string("_min)"); } int hit; double maxpole; double res; Flavour flav; string smax; for (;;) { hit = -1; maxpole = -1.; for (short int j=0;jfl; res = ATOOLS::sqr(flav.Width()*flav.Mass()); if (!ATOOLS::IsZero(res) && Massive(flav)) { if (1./res>maxpole) { maxpole = 1./res; hit = j; } } else { if (hit==-1) hit = j; } } } if (hit==-1) break; smax = string("sqr(sqrt(s")+Order(mummy)+string(")")+sum_s_i; for (short int j=0;j0.) { m_idc.push_back(string("MP")+ToString(hi)+string("_")+Order(lm[hit])); } else m_idc.push_back(string("MlP_")+Order(lm[hit])); break; case 0: case 10: sf<<" Vec4D "<0.) { sf<<" s"<0.) { sf<<" wt *= CE.MassivePropWeight(fl"<left==0) return; if (p->right->t) { props[count] = p->left; propt[count] = p->right; } else { if (p->left->t) { props[count] = p->right; propt[count] = p->left; } else { if (p->right->b == -1 && p->right->number<99) props[count] = p->left; else props[count] = p->right; return; } } count++; SetProps(propt[count-1],props,propt,count); } void Channel_Generator::CalcTSmin(int flag,vector& p,ofstream& sf) { string help; for (size_t i=0;i0) { if (help.length()==psum.length()) { CalcSmin(flag,"min",help,sf,0); return; } else CalcSmin(flag,"min",help,sf,0); s = string("sqr(sqrt(s")+Order(help)+string("_min)"); } else s = string("sqr("); for (size_t i=0;i1) s += string("+sqrt(s") + Order(p[i]) + string(")"); } s += string(")"); sf<<" double s"<GetscutAmegic(std::string(\"") + Order(lm) + string("\"))"),0,sf); /* string s(""); for (short int i=0;iscut[")+lm[i]+string("][")+lm[j]+string("]+"); } if (lm.length()==1) s += string("p_ms[")+lm[0]+string("]"); else { if (lm.length()==2) s +=string("0"); else { s += string("(2-")+IString(lm.length())+string(")*(p_ms[")+lm[0]+string("]"); for (short int i=1;i2 && p!=0) { AddToVariables(flag,Order(lm) + string("_") + string(min) + string("1"),s,0,sf); //sf<<" double s"<left==0) { char help[4]; snprintf(help,4,"%i",0); if (p->number<10) help[0]=p->number+48; else help[0]=p->number+55; return string(help); } return LinkedMasses(p->left)+LinkedMasses(p->right); } void Channel_Generator::IdentifyProps(Point* _plist) { InitT(_plist); Point* endp; tcount = 0; _plist->prev = 0; BackLinks(_plist,endp); Point* p = endp; if (p->prev != _plist) { for (;;) { p = p->prev; p->t = 1; tcount++; if (p->prev == _plist) break; } } } void Channel_Generator::BackLinks(Point* p,Point* &endp) { if ((p->left==0) && (p->right==0)) { if (p->b == -1) endp = p; return; } assert(p->left); assert(p->right); p->t = 0; p->left->prev = p; p->right->prev = p; BackLinks(p->left,endp); BackLinks(p->right,endp); } void Channel_Generator::InitT(Point* p) { p->t = 0; if (p->left==0) return; InitT(p->left); InitT(p->right); } string Channel_Generator::IString(int i) { MyStrStream sstr; sstr<>istr; return istr; } string Channel_Generator::Order(string s) { int beg = s.find("_"); if (beg!=-1) { return Order(s.substr(0,beg)) + string("_") + Order(s.substr(beg+1)); } if (s[0]>85 || s[0]<='0') return s; for (size_t i=0;is[j]) { char help = s[i]; s[i] = s[j]; s[j] = help; } } return s; } void Channel_Generator::AddToVariables(int flag,const string& lhs,const string& rhs,const int& type, ofstream& sf) { if (flag<0) return; string lhso = Order(lhs); std::string name; if (type ==0) name=string("s")+lhso; else name=string("p")+lhso; Decls::const_iterator cit=declarations.find(name); if (cit==declarations.end()) { // daoes not exist declarations[name]=rhs; if (type == 0) sf<<" double s"; else sf<<" Vec4D p"; sf<* kfs) { if (!p->left) return 0.; double m = 0.; if (p->fl.IsMassive()) { m = p->fl.Mass(); if (kfs) kfs->push_back(p->fl.Kfcode()); } return m + PMassSum(p->left,kfs) + PMassSum(p->right,kfs); }