#include "AMEGIC++/Phasespace/Channel_Generator_Decays.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 using namespace AMEGIC; using namespace PHASIC; using namespace ATOOLS; using namespace std; Channel_Generator_Decays::Channel_Generator_Decays(int _nin,int _nout, Point * _plist,int _aid) : Channel_Generator_Base(_nin,_nout,_plist) { m_idstr = string(""); m_aid = _aid; } Channel_Generator_Decays::~Channel_Generator_Decays() { } int Channel_Generator_Decays::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,"CD%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_Decays"<0) { chf <<" Info_Key "; bool first=true; for (size_t i=0; iMPISync(); }"<Optimize(); } "<EndOptimize(); } "<WriteOut(pId); } "<ReadIn(pId); } "<GeneratePoint(_ran);"<GenerateWeight(p_rans);" << endl; chf << " if (wt!=0.) wt = vw/wt/pow(2.*M_PI," << nout << "*3.-4.);" << endl; chf<AddPoint(Value,rans);"<1 process !!!"<fl<<" -> { "<left->fl<<" "<right->fl<<" }. Abort the run."<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) { if (flag==0 || flag==10) AddToVariables(flag,lm[i],string("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"<GetscutAmegic(std::string(\"") + Order(lm) + string("\"))"),0,sf); } string Channel_Generator_Decays::LinkedMasses(Point* p) { if (p->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); } string Channel_Generator_Decays::IString(int i) { MyStrStream sstr; sstr<>istr; return istr; } string Channel_Generator_Decays::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_Decays::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); }