#include "AMEGIC++/Phasespace/Channel_Generator3V.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_Generator3V::Channel_Generator3V(int _nin,int _nout, Point * _plist,int _aid) : Channel_Generator_Base(_nin,_nout,_plist) { IdentifyProps(plist); m_idstr = string(""); GenerateTopos(plist); m_aid = _aid; } Channel_Generator3V::~Channel_Generator3V() { for (size_t i=0;ileft; 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."<middle = NULL; ph->left = CopyPoints(p->left); ph->right = CopyPoints(p->right); ph->m = 1; } return ph; } Point* Channel_Generator3V::TransformTS(Point* p) { Point **props = new Point*[tcount+1]; Point **propt = new Point*[tcount+1]; int counts = 0; SetProps(p,props,propt,counts); Point* ph = new Point(*p); Point* ps = ph; m_pclist.push_back(ph); ph->m = 1; ph->right = CopyPoints(propt[tcount]); 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)) )) { ph->left = new Point(*propt[tcount-1]); m_pclist.push_back(ph->left); ph = ph->left; for (int i=0;im = 0; if (i%2==0) { ph->left = new Point(*propt[i/2]); m_pclist.push_back(ph->left); ph->right = CopyPoints(props[tcount-i/2]); ph = ph->left; } else { ph->right = new Point(*propt[tcount-(i+1)/2-1]); m_pclist.push_back(ph->right); ph->left = CopyPoints(props[(i-1)/2]); ph = ph->right; } } ph->m = 0; ph->left = CopyPoints(props[(tcount-1)/2]); ph->right = CopyPoints(props[(tcount+1)/2]); } else { ph->left = new Point(*propt[0]); m_pclist.push_back(ph->left); ph = ph->left; for (int i=0;im = 0; if (i%2==0) { ph->right = new Point(*propt[tcount-i/2-1]); m_pclist.push_back(ph->right); ph->left = CopyPoints(props[i/2]); ph = ph->right; } else { ph->left = new Point(*propt[(i+1)/2]); m_pclist.push_back(ph->left); ph->right = CopyPoints(props[tcount-(i-1)/2]); ph = ph->left; } } ph->m = 0; ph->left = CopyPoints(props[tcount/2]); ph->right = CopyPoints(props[tcount/2+1]); } delete[] props; delete[] propt; return ps; } void Channel_Generator3V::MRPScan() {} int Channel_Generator3V::MakeChannel(int& echflag,int n,string& path,string& pID) { if (m_idstr==string("")) CreateChannelID(echflag); //add Channel char name[22]; snprintf(name,22,"C%i_%i",nout,n); string filename = rpa->gen.Variable("SHERPA_CPP_PATH")+string("/Process/Amegic/")+path+string("/")+ string(name)+string(".C"); // cout<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 = s[\"TCHANNEL_ALPHA\"].Get();"<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."<m==0) { help+=string("ZS_"); } else { if (!ATOOLS::IsZero(ph->fl.Mass())) { help+=string("ZR")+ToString(ph->fl.Kfcode())+string("_"); } else help+=string("ZS_"); } { int pt; double mass = PMassSum(ph,&pt); if (pt>1) mass*=pow(1.5,(double)pt-1.); help+=ToString((int)mass); m_idc.push_back(help); } case 0: case 1: sf<<" Vec4D p"<m!=0 && !ATOOLS::IsZero(ph->fl.Mass())) { sf<<" type = 1;"<fl.Kfcode()<<")).Mass();"<fl.Kfcode()<<")).Width();"<1) mass*=pow(1.5,(double)pt-1.); sf<<" type = 2;"< pin0; vector pin1; GenerateDecayChain(flag,ph,rannum,sf,pin0,pin1); } void Channel_Generator3V::GenerateDecayChain(int flag,Point* p,int& rannum,ofstream& sf, vector pin0, vector pin1) { if (p->left==0) return; // int sa = AntennaS(p); // if (sa>2) return QCDAntenna(flag,p,rannum,sf,sa); if (p->m==0) { int hi = p->fl.Kfcode(); string tmstr; if (flag>=0 && !ATOOLS::IsZero(p->fl.Mass())) { tmstr = string("tmass")+ToString(p->number); sf<<" double "<0) { for (size_t i=0;i0) help1 = string("p0_") + Order(help1); else help1 = string("p[0]"); if (help2.length()>1) AddToVariables(flag,string("0_")+pin0sum,help1+string("-p")+help2,1,sf); else AddToVariables(flag,string("0_")+pin0sum,help1+string("-p[")+GetMassIndex(help2)+string("]"),1,sf); } if (pin1.size()>0) { for (size_t i=0;i0) help1 = string("p1_") + Order(help1); else help1 = string("p[1]"); if (help2.length()>1) AddToVariables(flag,string("1_")+pin1sum,help1+string("-p")+help2,1,sf); else AddToVariables(flag,string("1_")+pin1sum,help1+string("-p[")+GetMassIndex(help2)+string("]"),1,sf); } pin0sum = string("0_") + pin0sum; pin1sum = string("1_") + pin1sum; string pout0sum = Order(LinkedMasses(p->left)); string pout1sum = Order(LinkedMasses(p->right)); string sctmax("m_ctmax"); string sctmin("m_ctmin"); string his(""); switch (flag) { case -11: if(!ATOOLS::IsZero(p->fl.Mass())) his=ToString(hi); m_idc.push_back(string("TC")+his+ string("_")+pin0sum+string("_")+pin1sum+ string("_")+pout0sum+string("_")+pout1sum); break; case 0: sf<<" CE.TChannelMomenta("; if (pin0.size()==0) sf<<"p[0]"; else sf<<"p"<fl.Mass())) his=ToString(hi); string idh = string("TC")+his+ string("_")+pin0sum+string("_")+pin1sum+string("_")+pout0sum+string("_")+pout1sum; sf<<" if (m_k"<left->m==0) pin1.push_back(Order(LinkedMasses(p->right))); if (p->right->m==0) pin0.push_back(Order(LinkedMasses(p->left))); } else { 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; bool first = p->prev->number==0; // 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,rannum,sf,pin0,pin1); GenerateDecayChain(flag,p->right,rannum,sf,pin0,pin1); } bool Channel_Generator3V::QCDAntenna(int flag,Point* p,int& rannum,ofstream& sf,int n) { string lm = LinkedMasses(p->left); string rm = LinkedMasses(p->right); string mummy = Order(lm+rm); if (flag<0) { m_idc.push_back(string("AP_")+mummy); return 1; } if (flag>9 || flag==-1) { flag -= 10; } switch(flag) { case 0: sf <<" double s0"<scut["<scut["<left==0) { string m = LinkedMasses(p); AddToVariables(flag,m,string("p_ms[")+GetMassIndex(m)+string("]"),0,sf); return; } string lm,rm; string mummy,prt; string clm = Order(LinkedMasses(clmp)); lm = Order(LinkedMasses(p->left)); rm = Order(LinkedMasses(p->right)); mummy = Order(lm+rm); for (size_t i=0;izwf==0) { if (prt.length()>1) { Point* po=clmp->left; if (po==p) po=clmp->right; if (po->zwf) { if (!CheckVariables(flag,prt+string("_min"),0)) { if (!CheckVariables(flag,prt,0)) { sf<<" Flavour fl"<fl).Kfcode()<<"));"<=0 && CheckVariables(flag,prt,0)) sf<<" s"<< prt <<"_min = s"<< prt <<";"<fl).Kfcode()<<"));"<0) { for (size_t i=0;ifl).Kfcode())+string("_")+mummy); AddToVariables(flag,mummy+string("_max"),"s"+mummy,0,sf); AddToVariables(flag,mummy+string("_min"),"s"+mummy,0,sf); sclmp = p; } if (rm.length()>1 && lm.length()>1) sclmp = p; else if (clmp->left==p && LinkedMasses(clmp->right).length()>1) sclmp = p; else if (clmp->right==p && LinkedMasses(clmp->left).length()>1) sclmp = p; //cout<<":GenerateMassChain: right: "<right)<right,sclmp,rannum,sf); //cout<<":GenerateMassChain: left: "<left)<left,sclmp,rannum,sf); if (p->zwf!=0) return; if (clmp==p) return; // if (sclm!=mummy) { // if (prt.length()>1) { // if (!CheckVariables(flag,prt+string("_min"),0)) CalcSmin(flag,"min",prt,sf,0); // else if (flag>=0 && CheckVariables(flag,prt,0)) sf<<" s"<< prt <<"_min = s"<< prt <<";"<2 && flag>=0) { sf <<" s"<< mummy <<"_min = Max(s"<< mummy <<"_min,sqr(sqrt(s"<< lm <<")+sqrt(s"<< rm <<")));"<fl.Width()*p->fl.Mass()); if (p->m>0 && !ATOOLS::IsZero(res) && Massive(p->fl)) maxpole = 1./res; int hi = 4; if (mummy.length()>2) hi = 2; if (maxpole>0.) { hi = (p->fl).Kfcode(); if (flag>=0) sf<<" Flavour fl"<0.) { m_idc.push_back(string("MP")+ToString(hi)+string("_")+mummy); } else m_idc.push_back(string("MTH_")+Order(mummy)); break; case 0: sf<<" Vec4D p"<0.) { if (p->zwf==0) { sf<<" double s"<< mummy <<" = CE.MassivePropMomenta(fl"<0) { for (size_t i=0;i0.) { if (p->zwf==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; propt[count] = p->right; } else { props[count] = p->right; propt[count] = p->left; } return; } } count++; SetProps(propt[count-1],props,propt,count); } void Channel_Generator3V::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(")"); if (!CheckVariables(flag,Order(psum)+string("_min"),0)) { AddToVariables(flag,Order(psum) + string("_min"),s,0,sf); } else sf<<" s"<1) { AddToVariables(flag,Order(lm) + string("_") + string(min), string("cuts->GetscutAmegic(std::string(\"") + Order(lm) + string("\"))"),0,sf); } else { AddToVariables(flag,Order(lm) + string("_") + string(min), string("p_ms[") + GetMassIndex(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) return string(""); if (p->fl.Mass()>PMassSum(p->left,0)+PMassSum(p->right,0)) { return string("Flavour((kf_code)(")+ToString((p->fl).Kfcode())+string(")).Mass()"); } string h1 = GetFlMass(p->left); string h2 = GetFlMass(p->right); if (h1.length()==0) return h2; if (h2.length()==0) return h1; return h1+string("+")+h2; } string Channel_Generator3V::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); } void Channel_Generator3V::IdentifyProps(Point* _plist) { InitT(_plist); Point* endp; tcount = 0; _plist->prev = 0; BackLinks(_plist,endp); Point* p = endp; if (p->prev != _plist && p->prev!=NULL) { for (;;) { p = p->prev; p->t = 1; tcount++; if (p->prev == _plist) break; } } } void Channel_Generator3V::BackLinks(Point* p,Point* &endp) { if ((p->left==0) && (p->right==0)) { if (p->b == -1) endp = p; return; } p->t = 0; p->left->prev = p; p->right->prev = p; BackLinks(p->left,endp); BackLinks(p->right,endp); } void Channel_Generator3V::InitT(Point* p) { if (p->fl.Mass()>rpa->gen.Ecms()) m_valid=0; p->t = 0; if (p->left==0) return; InitT(p->left); InitT(p->right); } string Channel_Generator3V::IString(int i) { MyStrStream sstr; sstr<>istr; return istr; } string Channel_Generator3V::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_Generator3V::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 if (rhs!=string("")) { declarations[name]=rhs; if (type == 0) sf<<" double s"; else sf<<" Vec4D p"; sf<fl.Strong() || p->fl.IsMassive() || p->m!=1) return 0; if (p->left==0) return 1; int ls = AntennaS(p->left); int rs = AntennaS(p->right); if (ls==0 || rs==0) return 0; return ls+rs; } namespace AMEGIC { class Compare_String { public: int operator()(const string & a, const string & b) { return (aleft) return 0.; double m = 0.; if (p->m>0 && p->fl.IsMassive()) { m = p->fl.Mass(); } double mc = PMassSum(p->left,&tptl) + PMassSum(p->right,&tptr); if (pt) { if (m>mc) { *pt=1; } else if (tptl+tptr>0) { *pt=Max(tptl,tptr)+1; } } return Max(m,mc); }