HepLib
QGRAF.cpp
Go to the documentation of this file.
1 
6 #include "QGRAF.h"
7 
8 namespace HepLib::QGRAF {
9 
10  namespace {
11  string di_("di");
12  string li_("li");
13  string ti_("ti");
14  string fi_("fi");
15  string ci_("ci");
16  string ai_("ai");
17  string rdi_("rdi");
18  string rli_("rli");
19  string rti_("rti");
20  string rfi_("rfi");
21  string rci_("rci");
22  string rai_("rai");
23  inline string n2s(ex fn) {
24  int n = ex_to<numeric>(fn).to_int();
25  return (n<0 ? "m" : "") + to_string(abs(n));
26  }
27  }
28 
29  REGISTER_FUNCTION(Propagator, do_not_evalf_params())
30  REGISTER_FUNCTION(InField, do_not_evalf_params())
31  REGISTER_FUNCTION(OutField, do_not_evalf_params())
32 
33  #ifndef DOXYGEN_SKIP
34 
35  unsigned Field2_SERIAL::serial = GiNaC::function::register_new(function_options("Field",2).do_not_evalf_params().overloaded(2));
36  unsigned Field3_SERIAL::serial = GiNaC::function::register_new(function_options("Field",3).do_not_evalf_params().overloaded(2));
37  unsigned Vertex2_SERIAL::serial = GiNaC::function::register_new(function_options("Vertex",2).do_not_evalf_params().overloaded(5));
38  unsigned Vertex3_SERIAL::serial = GiNaC::function::register_new(function_options("Vertex",3).do_not_evalf_params().overloaded(5));
39  unsigned Vertex4_SERIAL::serial = GiNaC::function::register_new(function_options("Vertex",4).do_not_evalf_params().overloaded(5));
40  unsigned Vertex5_SERIAL::serial = GiNaC::function::register_new(function_options("Vertex",5).do_not_evalf_params().overloaded(5));
41  unsigned Vertex6_SERIAL::serial = GiNaC::function::register_new(function_options("Vertex",6).do_not_evalf_params().overloaded(5));
42 
43  #endif
44 
45  Index DI(ex fn) { return Index(di_+n2s(fn),Index::Type::VD); }
46  Index LI(ex fn) { return Index(li_+n2s(fn),Index::Type::VD); }
47  Index TI(ex fn) { return Index(ti_+n2s(fn),Index::Type::CF); }
48  Index FI(ex fn) { return Index(fi_+n2s(fn),Index::Type::CF); }
49  Index CI(ex fn) { return Index(ci_+n2s(fn),Index::Type::CA); }
50  Index AI(ex fn) { return Index(ai_+n2s(fn),Index::Type::CA); }
51  Index RDI(ex fn) { return Index(rdi_+n2s(fn),Index::Type::VD); }
52  Index RLI(ex fn) { return Index(rli_+n2s(fn),Index::Type::VD); }
53  Index RTI(ex fn) { return Index(rti_+n2s(fn),Index::Type::CF); }
54  Index RFI(ex fn) { return Index(rfi_+n2s(fn),Index::Type::CF); }
55  Index RCI(ex fn) { return Index(rci_+n2s(fn),Index::Type::CA); }
56  Index RAI(ex fn) { return Index(rai_+n2s(fn),Index::Type::CA); }
57 
63  lst Process::Amplitudes(symtab st) {
64  auto rc = system("rm -f qgraf.dat qgraf.out qgraf.sty qgraf.mod");
65  std::ofstream style;
66  style.open("qgraf.sty", ios::out);
67  style << Style << endl;
68  style.close();
69 
70  std::ofstream model;
71  style.open("qgraf.mod", ios::out);
72  style << Model << endl;
73  style.close();
74 
75  std::ofstream ofs;
76  ofs.open("qgraf.dat", ios::out);
77  ofs << "output='qgraf.out';" << endl;
78  ofs << "style='qgraf.sty';" << endl;
79  ofs << "model='qgraf.mod';" << endl;
80  ofs << "in=" << In << ";" << endl;
81  ofs << "out=" << Out << ";" << endl;
82  ofs << "loops=" << Loops << ";" << endl;
83  ofs << "loop_momentum=" << LoopPrefix << ";" << endl;
84  ofs << "options=" << Options << ";" << endl;
85  for(auto vs : Others) ofs << vs << ";" << endl;
86  ofs.close();
87 
88  if(Debug) rc = system((InstallPrefix+"/bin/qgraf").c_str());
89  else rc = system((InstallPrefix+"/bin/qgraf > /dev/null").c_str());
90 
91  ifstream ifs("qgraf.out");
92  string ostr((istreambuf_iterator<char>(ifs)), (istreambuf_iterator<char>()));
93  ifs.close();
94 
95  if(!Debug) {
96  if(access("qgraf.dat",F_OK)!=-1) remove("qgraf.dat");
97  if(access("qgraf.mod",F_OK)!=-1) remove("qgraf.mod");
98  if(access("qgraf.sty",F_OK)!=-1) remove("qgraf.sty");
99  if(access("qgraf.out",F_OK)!=-1) remove("qgraf.out");
100  }
101 
102  const char* rm_chars = " \t\v\r\n,";
103  if(!ostr.empty()) {
104  ostr.erase(0, ostr.find_first_not_of(rm_chars));
105  ostr.erase(ostr.find_last_not_of(rm_chars)+1);
106  }
107  ostr = "{" + ostr + "}";
108 
109  Parser amp(st);
110  ex ret = amp.Read(ostr);
111 
112  return ex_to<lst>(ret);
113  }
114 
120  lst TopoLines(const ex & amp) {
121  map<ex,int,ex_is_less> v2id, fid2vid;
122  map<int,ex> vid2fs; // fileds in the vertex
123  int cid = 0;
124  lst lines;
125  for(const_preorder_iterator i = amp.preorder_begin(); i != amp.preorder_end(); ++i) {
126  ex e = (*i);
127  if(isFunction(e,"OutField")) {
128  lines.append(lst{ e.op(1), iWF(e.op(1)), lst{e.op(0)}, e.op(2) });
129  } else if(isFunction(e,"InField")) {
130  lines.append(lst{ iWF(e.op(1)), e.op(1), lst{e.op(0)}, e.op(2) });
131  } else if(isFunction(e, "Propagator")) {
132  lines.append(lst{ iWF(e.op(0).op(1)), iWF(e.op(1).op(1)), lst{e.op(0).op(0), e.op(1).op(0)}, e.op(2) });
133  } else if(isFunction(e, "Vertex")) {
134  auto itr = v2id.find(e);
135  if(itr==v2id.end()) {
136  cid++;
137  v2id.emplace(make_pair(e,cid));
138  lst fs;
139  for(auto f : e) {
140  ex fid = f.op(1);
141  fid2vid[fid] = cid;
142  fs.append(f.op(0));
143  }
144  vid2fs[cid] = fs;
145  }
146  }
147  }
148  lines = ex_to<lst>(MapFunction([&vid2fs,&fid2vid](const ex &e, MapFunction &self)->ex{
149  if(!e.has(iWF(w))) return e;
150  else if(e.match(iWF(w))) {
151  auto vid = fid2vid[e.op(0)];
152  return lst{vid, vid2fs[vid]};
153  } else return e.map(self);
154  })(lines));
155 
156  return lines.sort();
157  }
158 
165  void DrawPDF(const lst & amps, string fn) {
166  int id=0, rc;
167  exvector amp_vec;
168  for(auto item : amps) amp_vec.push_back(item);
169  string tex_path = to_string(getpid()) + "_TeX/";
170  if(!dir_exists(tex_path)) rc = system(("mkdir -p "+tex_path).c_str());
171  int limit = 300;
172 
173  GiNaC_Parallel(amp_vec.size(), [&amp_vec,tex_path](int idx)->ex {
174  auto amp = amp_vec[idx];
175  ofstream out(tex_path+to_string(idx)+".tex");
176  out << "\\documentclass[tikz]{standalone}" << endl;
177  out << "\\usepackage{tikz-feynman}" << endl;
178  out << "\\tikzfeynmanset{compat=1.1.0}" << endl;
179  out << "\\begin{document}" << endl;
180  out << "\\feynmandiagram{" << endl;
181  auto lines = TopoLines(amp);
182 
183  exmap bend_map;
184  std::map<ex,int,ex_is_less> vtex_map; // vertex option, only once
185  for(auto l : lines) {
186  lst ll = lst{l.op(0), l.op(1)};
187  bool isExt = (is_a<numeric>(ll.op(0)) && ll.op(0)<0) || (is_a<numeric>(ll.op(1)) && ll.op(1)<0);
188  ll.sort();
189  bend_map[ll] = bend_map[ll] + 1;
190 
191  auto fidL = (is_a<numeric>(l.op(1)) ? l.op(1) : l.op(1).op(0));
192  out << "\"" << fidL << "\"";
193  if(fidL<0) {
194  out << "[particle=";
195  if(InOutTeX[fidL].length()>0) out << InOutTeX[fidL];
196  else out << fidL;
197  out << "]";
198  } else if(vtex_map[l.op(1).op(1)]==0) {
199  out << VerTeX[l.op(1).op(1)];
200  vtex_map[l.op(1).op(1)]=1;
201  }
202 
203  out << " --[";
204  auto f = l.op(2).op(0);
205  if(LineTeX[f].length()>0) {
206  if(!isExt) out << LineTeX[f];
207  else {
208  auto cpos = LineTeX[f].find(", edge");
209  if(cpos>0) out << LineTeX[f].substr(0,cpos);
210  else out << LineTeX[f];
211  }
212  }
213  if(bend_map[ll]>2) out << ",half right";
214  else if(bend_map[ll]>1) out << ",half left";
215  if(is_zero(l.op(0)-l.op(1))) out << ",loop,distance=2cm";
216  out << "]";
217 
218  auto fidR = (is_a<numeric>(l.op(0)) ? l.op(0) : l.op(0).op(0));
219  out << " \"" << fidR << "\"";
220  if(fidR<0) {
221  out << "[particle=";
222  if(InOutTeX[fidR].length()>0) out << InOutTeX[fidR];
223  else out<< fidR;
224  out << "]";
225  } else if(vtex_map[l.op(0).op(1)]==0) {
226  out << VerTeX[l.op(0).op(1)];
227  vtex_map[l.op(0).op(1)]=1;
228  }
229  out << ";" << endl;
230  }
231  out << "};" << endl;
232  out << "\\end{document}" << endl;
233  out.close();
234  auto rc = system(("cd "+tex_path+" && echo X | lualatex " + to_string(idx) + " 1>/dev/null").c_str());
235  return 0;
236  }, "TeX");
237 
238  ofstream out(tex_path+"diagram.tex");
239  out << "\\let\\mypdfximage\\pdfximage" << endl;
240  out << "\\def\\pdfximage{\\immediate\\mypdfximage}" << endl;
241  out << "\\documentclass{standalone}" << endl;
242  out << "\\usepackage{graphicx}" << endl;
243  out << "\\usepackage{adjustbox}" << endl;
244  out << "\\begin{document}" << endl;
245  out << "\\begin{adjustbox}{valign=T,width=\\textwidth}" << endl;
246  out << "\\begin{tabular}{|cc|cc|cc|cc|}" << endl;
247  out << "\\hline" << endl;
248  int total = amps.nops();
249  int namps = total;
250  if((total%4)!=0) total = (total/4+1)*4;
251  for(int i=0 ; i<total; i++) {
252 
253  if((i!=0) && (i+1!=total) && (i%limit)==0) {
254  out << "\\end{tabular}" << endl << endl;
255  out << "\\begin{tabular}{|cc|cc|cc|cc|}" << endl;
256  out << "\\hline" << endl;
257  }
258 
259  out << "{\\tiny " << i+1 << "}&" << endl;
260  if(i<namps) {
261  out << "\\includegraphics[keepaspectratio,";
262  out << "height=0.22\\textwidth,";
263  out << "width=0.22\\textwidth]";
264  out << "{"<<i<<".pdf}" << endl;
265  }
266  if((i+1)%4==0) out << "\\\\ \\hline";
267  else out << "&";
268  }
269  out << "\\end{tabular}" << endl;
270  out << "\\end{adjustbox}" << endl;
271  out << "\\end{document}" << endl;
272  out.close();
273  if(Debug) rc = system(("cd "+tex_path+" && pdflatex diagram && mv diagram.pdf ../"+fn).c_str());
274  else rc = system(("cd "+tex_path+" && echo X | pdflatex diagram 1>/dev/null && mv diagram.pdf ../"+fn).c_str());
275  if(!Debug) rc = system(("rm -r "+tex_path).c_str());
276  }
277 
285  vector<lst> ShrinkCut(ex amp, lst prop, int n) {
286  vector<lst> ret;
287  if(prop.nops()<1) throw Error("ShrinkCut: no cut provided!");
288  auto tls = TopoLines(amp);
289  vector<int> cls_vec;
290  for(int i=0; i<tls.nops(); i++) {
291  auto pi = tls.op(i).op(2);
292  if(pi.nops()<2) continue;
293  if(!is_a<lst>(prop.op(0))) { // prop : lst{ g, g }
294  if(is_zero(pi.op(0)-prop.op(0)) && is_zero(pi.op(1)-prop.op(1))) cls_vec.push_back(i);
295  } else {
296  for(auto iprop : prop) {
297  if(is_zero(pi.op(0)-iprop.op(0)) && is_zero(pi.op(1)-iprop.op(1))) cls_vec.push_back(i);
298  }
299  }
300  }
301  if(cls_vec.size()<n) return ret;
302 
303  Combinations(cls_vec.size(), n, [n,&ret,cls_vec,tls](const int * is)->void {
304  int cls[n];
305  for(int i=0; i<n; i++) cls[i] = cls_vec[is[i]];
306 
307  // cut each line into 2 half-lines, labeled with 0
308  auto tls2 = tls;
309  for(auto ci : cls) {
310  auto ol = tls2.op(ci);
311  tls2.let_op(ci) = lst{ol.op(0), 0, lst{ol.op(2).op(0)}, ol.op(3)};
312  tls2.append(lst{0, ol.op(1), lst{ol.op(2).op(1)}, ol.op(3)} );
313  }
314 
315  // shrink internal lines
316  int last = 0;
317  int ntls2 = tls2.nops();
318  while(true) {
319  ex lp = 0;
320  for(int i=last; i<ntls2; i++) {
321  auto li = tls2.op(i);
322  if(is_zero(li) || li.op(2).nops()<2) continue;
323  last = i;
324  lp = li;
325  tls2.let_op(last) = 0;
326  break;
327  }
328  if(is_zero(lp)) break;
329 
330  for(int i=0; i<ntls2; i++) {
331  if(is_zero(tls2.op(i))) continue;
332  if(is_zero(tls2.op(i).op(0)-lp.op(0))) tls2.let_op(i).let_op(0) = lp.op(1);
333  if(is_zero(tls2.op(i).op(1)-lp.op(0))) tls2.let_op(i).let_op(1) = lp.op(1);
334  }
335  }
336 
337  // final connected parts
338  map<int, lst> con_map;
339  for(auto li : tls2) {
340  if(is_zero(li)) continue;
341  ex key, val;
342  ex fiL = li.op(0);
343  if(is_a<lst>(fiL)) fiL = fiL.op(0);
344  ex fiR = li.op(1);
345  if(is_a<lst>(fiR)) fiR = fiR.op(0);
346  if(fiL>0 && fiR<0) {
347  con_map[ex_to<numeric>(fiL).to_int()].append(fiR);
348  } else if(fiR>0 && fiL<0) {
349  con_map[ex_to<numeric>(fiR).to_int()].append(fiL);
350  } else {
351  if(fiL>0 && is_zero(fiR)) key = fiL;
352  else if(fiR>0 && is_zero(fiL)) key = fiR;
353  else throw Error("ShrinkCut: unexpcected point reached.");
354  val = li.op(2).op(0);
355  con_map[ex_to<numeric>(key).to_int()].append(val);
356  }
357  }
358 
359  lst item;
360  for(auto kv : con_map) item.append(kv.second.sort());
361  ret.push_back(item);
362  });
363 
364  return ret;
365  }
366 
373  bool HasLoop(ex amp, lst prop) {
374  auto tls = TopoLines(amp);
375  int ntls = tls.nops();
376  // shrink internal lines of prop
377  int last = 0;
378  while(true) {
379  ex lp = 0;
380  for(int i=last; i<ntls; i++) {
381  auto li = tls.op(i);
382  if(is_zero(li) || li.op(2).nops()<2) continue; // external line
383  auto cpi = li.op(2);
384  bool m1 = !(is_zero(cpi.op(0)-prop.op(0)) && is_zero(cpi.op(1)-prop.op(1)));
385  bool m2 = !(is_zero(cpi.op(0)-prop.op(1)) && is_zero(cpi.op(1)-prop.op(0)));
386  if(m1 && m2) continue; // other propagators
387  if(is_zero(li.op(0)-li.op(1))) return true; // a loop found
388  last = i;
389  lp = li;
390  tls.let_op(last) = 0;
391  break;
392  }
393  if(is_zero(lp)) break;
394 
395  for(int i=0; i<ntls; i++) {
396  if(is_zero(tls.op(i))) continue;
397  if(is_zero(tls.op(i).op(0)-lp.op(0))) tls.let_op(i).let_op(0) = lp.op(1);
398  if(is_zero(tls.op(i).op(1)-lp.op(0))) tls.let_op(i).let_op(1) = lp.op(1);
399  }
400  }
401  return false;
402  }
403 
410  ex QuarkPropagator(const ex & e, const ex & m) {
411  auto fi1 = e.op(0).op(1);
412  auto fi2 = e.op(1).op(1);
413  auto mom = e.op(2);
414  return I * SP(TI(fi1),TI(fi2)) * Matrix(GAS(mom)+GAS(1)*m, DI(fi1),DI(fi2)) / (SP(mom)-m*m);
415  }
416 
423  ex LeptonPropagator(const ex & e, const ex & m) {
424  auto fi1 = e.op(0).op(1);
425  auto fi2 = e.op(1).op(1);
426  auto mom = e.op(2);
427  return I * Matrix(GAS(mom)+GAS(1)*m, DI(fi1),DI(fi2)) / (SP(mom)-m*m);
428  }
429 
436  ex GluonPropagator(const ex & e, const ex & xi) {
437  auto fi1 = e.op(0).op(1);
438  auto fi2 = e.op(1).op(1);
439  auto mom = e.op(2);
440  return (-I) * SP(CI(fi1),CI(fi2)) * ( SP(LI(fi1),LI(fi2)) / SP(mom) - (1-xi)*SP(mom,LI(fi1))*SP(mom,LI(fi2))/pow(SP(mom),2) );
441  }
442 
449  ex GluonGhostPropagator(const ex & e, const ex & xi, const ex & eta_G) {
450  auto fi1 = e.op(0).op(1);
451  auto fi2 = e.op(1).op(1);
452  auto mom = e.op(2);
453  return I * eta_G * SP(CI(fi1),CI(fi2)) / SP(mom);
454  }
455 
463  ex AZWPropagator(const ex & e, const ex & m, const ex & xi) {
464  auto fi1 = e.op(0).op(1);
465  auto fi2 = e.op(1).op(1);
466  auto mom = e.op(2);
467  return -I/(SP(mom)-m*m) * (SP(LI(fi1),LI(fi2))-(1-xi)*SP(LI(fi1))*SP(LI(fi2))/(SP(mom)-xi*m*m));
468  }
469 
477  ex AZWGhostPropagator(const ex & e, const ex & m, const ex & xi, const ex & eta_G) {
478  auto fi1 = e.op(0).op(1);
479  auto fi2 = e.op(1).op(1);
480  auto mom = e.op(2);
481  return eta_G * I / (SP(mom)-xi*m*m);
482  }
483 
491  ex ScalarPropagator(const ex & e, const ex & m, const ex & xi) {
492  auto fi1 = e.op(0).op(1);
493  auto fi2 = e.op(1).op(1);
494  auto mom = e.op(2);
495  return I / (SP(mom)-xi*m*m);
496  }
497 
504  ex QuarkGluonVertex(const ex & e, const ex & eta_s) {
505  auto fi1 = e.op(0).op(1);
506  auto fi2 = e.op(1).op(1);
507  auto fi3 = e.op(2).op(1);
508  return -I*eta_s*gs*Matrix(GAS(LI(fi3)),DI(fi1),DI(fi2))*SUNT(CI(fi3),TI(fi1),TI(fi2));
509  }
510 
517  ex Gluon3Vertex(const ex & e, const ex & eta_s) {
518  auto fi1 = e.op(0).op(1);
519  auto fi2 = e.op(1).op(1);
520  auto fi3 = e.op(2).op(1);
521  auto mom1 = e.op(0).op(2);
522  auto mom2 = e.op(1).op(2);
523  auto mom3 = e.op(2).op(2);
524  return -eta_s*gs*SUNF(CI(fi1),CI(fi2),CI(fi3))*(
525  SP(mom1-mom2,LI(fi3))*SP(LI(fi1),LI(fi2)) +
526  SP(mom2-mom3,LI(fi1))*SP(LI(fi2),LI(fi3)) +
527  SP(mom3-mom1,LI(fi2))*SP(LI(fi3),LI(fi1))
528  );
529  }
530 
537  ex Gluon4Vertex(const ex & e, const ex & eta_s) {
538  auto fi1 = e.op(0).op(1);
539  auto fi2 = e.op(1).op(1);
540  auto fi3 = e.op(2).op(1);
541  auto fi4 = e.op(3).op(1);
542  return -I*gs*gs*(
543  SUNF4(CI(fi1),CI(fi2),CI(fi3),CI(fi4))*(SP(LI(fi1),LI(fi3))*SP(LI(fi2),LI(fi4))-SP(LI(fi1),LI(fi4))*SP(LI(fi2),LI(fi3))) +
544  SUNF4(CI(fi1),CI(fi3),CI(fi2),CI(fi4))*(SP(LI(fi1),LI(fi2))*SP(LI(fi3),LI(fi4))-SP(LI(fi1),LI(fi4))*SP(LI(fi2),LI(fi3))) +
545  SUNF4(CI(fi1),CI(fi4),CI(fi2),CI(fi3))*(SP(LI(fi1),LI(fi2))*SP(LI(fi4),LI(fi3))-SP(LI(fi1),LI(fi3))*SP(LI(fi4),LI(fi2)))
546  );
547  }
548 
556  ex GhostGluonVertex(const ex & e, const ex & eta_s, const ex & eta_G) {
557  auto fi1 = e.op(0).op(1);
558  auto fi2 = e.op(1).op(1);
559  auto fi3 = e.op(2).op(1);
560  auto mom1 = e.op(0).op(2);
561  return eta_s*eta_G*gs*SUNF(CI(fi1),CI(fi2),CI(fi3))*SP(mom1,LI(fi3));
562  }
563 
571  ex QuarkAVertex(const ex & e, const ex & eq, const ex & eta_e) {
572  auto fi1 = e.op(0).op(1);
573  auto fi2 = e.op(1).op(1);
574  auto fi3 = e.op(2).op(1);
575  return -I*eta_e*eq*Matrix(GAS(LI(fi3)),DI(fi1),DI(fi2))*SP(TI(fi1),TI(fi2));
576  }
577 
585  ex LeptonAVertex(const ex & e, const ex & eq, const ex & eta_e) {
586  auto fi1 = e.op(0).op(1);
587  auto fi2 = e.op(1).op(1);
588  auto fi3 = e.op(2).op(1);
589  return -I*eta_e*eq*Matrix(GAS(LI(fi3)),DI(fi1),DI(fi2));
590  }
591 
598  ex IndexL2R(ex e, bool all) {
599  static MapFunction map([all](const ex &e, MapFunction &self)->ex {
600  if(!Index::has(e)) return e;
601  else if(is_a<Index>(e)) {
602  auto idx = ex_to<Index>(e);
603  auto nstr = idx.name.get_name();
604  if(!all && nstr.rfind("lim",0)==0) return e;
605  else if(!all && nstr.rfind("dim",0)==0) return e;
606  else if(!all && nstr.rfind("cim",0)==0) return e;
607  else if(!all && nstr.rfind("tim",0)==0) return e;
608  else if(nstr.rfind("li",0)==0 || nstr.rfind("di",0)==0 || nstr.rfind("ci",0)==0 || nstr.rfind("ti",0)==0) return Index("r"+nstr, idx.type);
609  else return e;
610  }
611  else return e.map(self);
612  });
613  return map(e);
614  }
615 
616  ex IndexCC(ex e, bool all) {
617  return IndexL2R(e,all);
618  }
619 
626  ex GluonSumL(int qi, bool color) {
627  if(color) return -SP(LI(qi), RLI(qi)) * SP(CI(qi), RCI(qi));
628  else return -SP(LI(qi), RLI(qi));
629  }
630 
637  ex GluonSumR(int qi, bool color) {
638  if(color) return -SP(RLI(qi),LI(qi)) * SP(RCI(qi),CI(qi));
639  else return -SP(RLI(qi),LI(qi));
640  }
641 
650  ex QuarkSumL(int qi, ex p, ex m, bool color) {
651  if(color) return Matrix(GAS(p)+m*GAS(1), RDI(qi), DI(qi)) * SP(RTI(qi), TI(qi));
652  else return Matrix(GAS(p)+m*GAS(1), RDI(qi), DI(qi));
653  }
654 
663  ex QuarkSumR(int qi, ex p, ex m, bool color) {
664  if(color) return Matrix(GAS(p)+m*GAS(1), DI(qi), RDI(qi)) * SP(TI(qi),RTI(qi));
665  else return Matrix(GAS(p)+m*GAS(1), DI(qi), RDI(qi));
666  }
667 
676  ex AntiQuarkSumL(int qi, ex p, ex m, bool color) {
677  if(color) return Matrix(GAS(p)-m*GAS(1), DI(qi), RDI(qi)) * SP(TI(qi), RTI(qi));
678  else return Matrix(GAS(p)-m*GAS(1), DI(qi), RDI(qi));
679  }
680 
689  ex AntiQuarkSumR(int qi, ex p, ex m, bool color) {
690  if(color) return Matrix(GAS(p)-m*GAS(1), RDI(qi), DI(qi)) * SP(RTI(qi), TI(qi));
691  else return Matrix(GAS(p)-m*GAS(1), RDI(qi), DI(qi));
692  }
693 
699  ex GhostSumL(int qi) {
700  return SP(CI(qi), RCI(qi));
701  }
702 
708  ex GhostSumR(int qi) {
709  return SP(RCI(qi),CI(qi));
710  }
711 
717  ex AntiGhostSumL(int qi) {
718  return -SP(CI(qi), RCI(qi));
719  }
720 
726  ex AntiGhostSumR(int qi) {
727  return -SP(RCI(qi),CI(qi));
728  }
729 
736  ex J1SumL(int qi, ex p) {
737  if(is_zero(p)) return -SP(LI(qi),RLI(qi));
738  else return -SP(LI(qi),RLI(qi)) + SP(p,LI(qi)) * SP(p,RLI(qi)) / SP(p);
739  }
740 
747  ex J1SumR(int qi, ex p) {
748  if(is_zero(p)) return -SP(RLI(qi),LI(qi));
749  else return -SP(RLI(qi),LI(qi)) + SP(p,RLI(qi)) * SP(p,LI(qi)) / SP(p);
750  }
751 
752  // https://arxiv.org/abs/1209.6213v2
753  lst Models::FeynRulesSM(const lst & amps, const ex & xi) {
754  lst ret;
755  for(auto item : amps) ret.append(FeynRulesSM(item,xi));
756  return ret;
757  }
758  ex Models::FeynRulesSM(const ex & amp, const ex & xi) {
759  if(is_a<lst>(amp)) return FeynRulesSM(ex_to<lst>(amp), xi);
760 
761  static Symbol CW("CW"); // cos(theta)
762  static Symbol SW("SW"); // sin(theta)
763  //static Symbol C2W("C2W"); // cos(2theta)
764  static ex CW2 = CW*CW;
765  static ex SW2 = SW*SW;
766  static ex C2W = CW2-SW2;
767 
768  static Symbol U("U"); // U-quark
769  static Symbol Ubar("Ubar"); // anti U-quark
770  static Symbol D("D"); // D-quark
771  static Symbol Dbar("Dbar"); // anti D-quark
772  static Symbol C("C"); // C-quark
773  static Symbol Cbar("Cbar"); // anti C-quark
774  static Symbol S("S"); // S-quark
775  static Symbol Sbar("Sbar"); // anti S-quark
776  static Symbol T("T"); // T-quark
777  static Symbol Tbar("Tbar"); // anti T-quark
778  static Symbol B("B"); // B-quark
779  static Symbol Bbar("Bbar"); // anti B-quark
780 
781  static Symbol g("g"); // gluon
782  static Symbol gh("gh"); // gluon ghost
783  static Symbol ghbar("ghbar"); // anti gluon ghost
784 
785  static Symbol A("A"); // photon
786  static Symbol Wm("Wm"); // W-
787  static Symbol Wp("Wp"); // W+
788  static Symbol Z("Z"); // Z
789 
790  static Symbol ghA("ghA"); // photon ghost
791  static Symbol ghAbar("ghAbar"); // anti photon ghost
792  static Symbol ghWm("ghWm"); // W- ghost
793  static Symbol ghWmbar("ghWmbar"); // anti W- ghost
794  static Symbol ghWp("ghWp"); // W+ ghost
795  static Symbol ghWpbar("ghWpbar"); // anti W+ ghost
796  static Symbol ghZ("ghZ"); // Z ghost
797  static Symbol ghZbar("ghZbar"); // anti Z ghost
798 
799  static Symbol em("em"); // e-
800  static Symbol ep("ep"); // e+
801  static Symbol ne("ne"); // e-neutrino
802  static Symbol nebar("nebar"); // anti e-neutrino
803  static Symbol mum("mum"); // mu-
804  static Symbol mup("mup"); // mu+
805  static Symbol nmu("nmu"); // mu-neutrino
806  static Symbol nmubar("nmubar"); // anti mu-neutrino
807  static Symbol taum("taum"); // tau-
808  static Symbol taup("taup"); // tau+
809  static Symbol ntau("ntau"); // tau-neutrino
810  static Symbol ntaubar("ntaubar"); // anti tau-neutrino
811 
812  static Symbol chi("chi"); // Z goldstone
813  static Symbol phim("phim"); // W- goldstone
814  static Symbol phip("phip"); // W+ goldstone
815  static Symbol H("H"); // higgs
816 
817  static auto M = [](const string & si)->ex {
818  return Symbol("M"+si);
819  };
820  static ex MW = M("W");
821  static ex MZ = M("Z");
822  static ex MH = M("H");
823 
824  static Symbol EL("e"); // charge e
825 
826  static ex eta_s = -1;
827  // Peskin's book convention
828  static ex eta = -1;
829  static ex eta_prime = -1;
830  static ex eta_Z = 1;
831  static ex eta_theta = 1;
832  static ex eta_Y = 1;
833  static ex eta_e = -1;
834  static ex eta_G = 1;
835 
836  static ex GEW = eta_e*eta*eta_theta * EL/SW;
837  static ex EL2 = EL*EL;
838  static ex MW2 = MW*MW;
839  static ex MZ2 = MZ*MZ;
840  static ex MH2 = MH*MH;
841  static ex GEW2 = GEW*GEW;
842  static ex sqrt2 = sqrt(ex(2));
843 
844  map<string,ex> T3;
845  T3["U"] = T3["C"] = T3["T"] = ex(1)/2;
846  T3["D"] = T3["S"] = T3["B"] = -ex(1)/2;
847  T3["ne"] = T3["nmu"] = T3["ntau"] = ex(1)/2;
848  T3["em"] = T3["mum"] = T3["taum"] = -ex(1)/2;
849  T3["ep"] = T3["mup"] = T3["taup"] = -ex(1)/2; // for sure
850 
851  map<string,ex> Q;
852  Q["U"] = Q["C"] = Q["T"] = ex(2)/3;
853  Q["D"] = Q["S"] = Q["B"] = -ex(1)/3;
854  Q["ne"] = Q["nmu"] = Q["ntau"] = 0;
855  Q["em"] = Q["mum"] = Q["taum"] = -1;
856  Q["ep"] = Q["mup"] = Q["taup"] = -1; // for sure
857 
858  auto is_qbar_q = [&](const ex pi1, const ex pi2)->bool {
859  return (pi1==Ubar && pi2==U) || (pi1==Dbar && pi2==D) || (pi1==Cbar && pi2==C) || (pi1==Sbar && pi2==S) || (pi1==Tbar && pi2==T) || (pi1==Bbar && pi2==B);
860  };
861  auto is_lbar_l = [&](const ex pi1, const ex pi2)->bool {
862  return (pi1==ep && pi2==em) || (pi1==mup && pi2==mum) || (pi1==taup && pi2==taum);
863  };
864  auto is_lbar_n = [&](const ex pi1, const ex pi2)->bool {
865  return (pi1==ep && pi2==ne) || (pi1==mup && pi2==nmu) || (pi1==taup && pi2==ntau);
866  };
867  auto is_nbar_n = [&](const ex pi1, const ex pi2)->bool {
868  return (pi1==nebar && pi2==ne) || (pi1==nmubar && pi2==nmu) || (pi1==ntaubar && pi2==ntau);
869  };
870  auto is_nbar_l = [&](const ex pi1, const ex pi2)->bool {
871  return (pi1==nebar && pi2==em) || (pi1==nmubar && pi2==mum) || (pi1==ntaubar && pi2==taum);
872  };
873  auto is_uqbar = [&](const ex pi)->bool {
874  return (pi==Ubar) || (pi==Cbar) || (pi==Tbar);
875  };
876  auto is_uq = [&](const ex pi)->bool {
877  return (pi==U) || (pi==C) || (pi==T);
878  };
879  auto is_dqbar = [&](const ex pi)->bool {
880  return (pi==Dbar) || (pi==Sbar) || (pi==Bbar);
881  };
882  auto is_dq = [&](const ex pi)->bool {
883  return (pi==D) || (pi==S) || (pi==B);
884  };
885  auto is_uqbar_dq = [&](const ex pi1, const ex pi2)->bool {
886  return is_uqbar(pi1) && is_dq(pi2);
887  };
888  auto is_dqbar_uq = [&](const ex pi1, const ex pi2)->bool {
889  return is_dqbar(pi1) && is_uq(pi2);
890  };
891  auto CKM = [&](const string & si1, const string & si2)->ex {
892  return Symbol("V"+si1+si2);
893  };
894 
895  auto gfV = [&](const string & si)->ex {
896  return T3[si]/2-Q[si]*SW2;
897  };
898  auto gfA = [&](const string & si)->ex {
899  return T3[si]/2;
900  };
901 
902  auto fr = MapFunction([&](const ex &e, MapFunction &self)->ex {
903  if(isFunction(e,"OutField") || isFunction(e,"InField")) return 1;
904  else if(isFunction(e, "Propagator")) {
905  auto pi = e.op(0).op(0);
906  auto si = ex2str(pi);
907  if(pi==U || pi==D || pi==C || pi==S || pi==T || pi==B) {
908  return QuarkPropagator(e, M(si));
909  } else if(pi==g) {
910  return GluonPropagator(e, xi);
911  } else if(pi==gh) {
912  return GluonGhostPropagator(e, xi, eta_G);
913  } else if(pi==A) {
914  return AZWPropagator(e, 0, xi);
915  } else if(pi==Z) {
916  return AZWPropagator(e, MZ, xi);
917  } else if(pi==Wm) {
918  return AZWPropagator(e, MW, xi);
919  } else if(pi==ghA) {
920  return AZWGhostPropagator(e, 0, xi, eta_G);
921  } else if(pi==ghZ) {
922  return AZWGhostPropagator(e, MZ, xi, eta_G);
923  } else if(pi==ghWm || pi==ghWp) {
924  return AZWGhostPropagator(e, MW, xi, eta_G);
925  } else if(pi==em || pi==mum || pi==taum) {
926  return LeptonPropagator(e, M(si));
927  } else if(pi==ne || pi==nmu || pi==ntau) {
928  return LeptonPropagator(e, 0);
929  } else if(pi==H) {
930  return ScalarPropagator(e, MH, 1);
931  } else if(pi==chi) {
932  return ScalarPropagator(e, MZ, xi);
933  } else if(pi==phim) {
934  return ScalarPropagator(e, MW, xi);
935  } else {
936  cout << endl << e << endl;
937  throw Error("Propagator Not defined!");
938  }
939  } else if(isFunction(e, "Vertex")) {
940  auto pi1 = e.op(0).op(0);
941  auto pi2 = e.op(1).op(0);
942  auto pi3 = e.op(2).op(0);
943  auto fi1 = e.op(0).op(1);
944  auto fi2 = e.op(1).op(1);
945  auto fi3 = e.op(2).op(1);
946  auto mom1 = e.op(0).op(2);
947  auto mom2 = e.op(1).op(2);
948  auto mom3 = e.op(2).op(2);
949  int nn = e.nops();
950 
951  auto si1 = ex2str(pi1);
952  auto si2 = ex2str(pi2);
953  string_replace_all(si1, "bar", "");
954  string_replace_all(si2, "bar", "");
955  auto si = si2;
956 
957  if(nn==3) {
958  if(pi1==ghbar && pi2==gh && pi3==g) {
959  // ghbar-gh-g
960  return GhostGluonVertex(e, eta_s, eta_G);
961  } else if(pi1==g && pi2==g && pi3==g) {
962  // g^3
963  return Gluon3Vertex(e, eta_s);
964  } else if((pi3==g) && is_qbar_q(pi1, pi2)) {
965  // qbar-q-g
966  return QuarkGluonVertex(e, eta_s);
967  } else if((pi3==A) && is_qbar_q(pi1, pi2)) {
968  // qbar-q-A
969  return QuarkAVertex(e, EL*Q[si], eta_e);
970  } else if((pi3==A) && is_lbar_l(pi1, pi2)) {
971  // lbar-l-A
972  return LeptonAVertex(e, EL*Q[si], eta_e);
973  } else if((pi3==Z) && is_qbar_q(pi1, pi2)) {
974  // qbar-q-Z
975  return -I*eta*eta_Z*GEW/CW*(gfV(si)*Matrix(GAS(LI(fi3)),DI(fi1),DI(fi2))-gfA(si)*Matrix(GAS(LI(fi3))*GAS(5),DI(fi1),DI(fi2)))*SP(TI(fi1),TI(fi2));
976  } else if((pi3==Z) && (is_lbar_l(pi1, pi2) || is_nbar_n(pi1, pi2))) {
977  // lbar-l-Z & nbar-n-Z
978  return -I*eta*eta_Z*GEW/CW*(gfV(si)*Matrix(GAS(LI(fi3)),DI(fi1),DI(fi2))-gfA(si)*Matrix(GAS(LI(fi3))*GAS(5),DI(fi1),DI(fi2)));
979  } else if(pi1==Wp && pi2==Wm && pi3==A) {
980  // Wp-Wm-A
981  return -I*eta_e*EL*(SP(LI(fi1),LI(fi2))*SP(mom2-mom1,LI(fi3))+SP(LI(fi2),LI(fi3))*SP(mom3-mom2,LI(fi1))+SP(LI(fi3),LI(fi1))*SP(mom1-mom3,LI(fi2)));
982  } else if(pi1==Wp && pi2==Wm && pi3==Z) {
983  // Wp-Wm-Z
984  return -I*eta*eta_Z*GEW*CW*(SP(LI(fi1),LI(fi2))*SP(mom2-mom1,LI(fi3))+SP(LI(fi2),LI(fi3))*SP(mom3-mom2,LI(fi1))+SP(LI(fi3),LI(fi1))*SP(mom1-mom3,LI(fi2)));
985  } else if(is_uqbar(pi1) && is_dq(pi2) && pi3==Wp) {
986  // uqbar-dq-Wp
987  auto ckm = CKM(si1, si2);
988  return -I*eta*GEW/sqrt2*Matrix(GAS(LI(fi3))-GAS(LI(fi3))*GAS(5),DI(fi1),DI(fi2))/2*ckm*SP(TI(fi1),TI(fi2));
989  } else if(is_dqbar(pi1) && is_uq(pi2) && pi3==Wm) {
990  // dqbar-uq-Wm
991  auto ckm = CKM(si1, si2);
992  return -I*eta*GEW/sqrt2*Matrix(GAS(LI(fi3))-GAS(LI(fi3))*GAS(5),DI(fi1),DI(fi2))/2*ckm*SP(TI(fi1),TI(fi2));
993  } else if( (is_nbar_l(pi1, pi2) && pi3==Wp) || (is_lbar_n(pi1, pi2) && pi3==Wm) ) {
994  // nbar-l-Wp & lbar-n-Wm
995  return -I*eta*GEW/sqrt2*Matrix(GAS(LI(fi3))-GAS(LI(fi3))*GAS(5),DI(fi1),DI(fi2))/2;
996  } else if(is_qbar_q(pi1, pi2) && (pi3==H)) {
997  // qbar-q-H
998  return -I*GEW/2*M(si)/MW*Matrix(GAS(1),DI(fi1),DI(fi2))*SP(TI(fi1),TI(fi2));
999  } else if( (is_lbar_l(pi1, pi2) || is_nbar_n(pi1, pi2)) && (pi3==H)) {
1000  // lbar-l-H & nbar-n-H
1001  return -I*GEW/2*M(si)/MW*Matrix(GAS(1),DI(fi1),DI(fi2));
1002  } else if(is_qbar_q(pi1, pi2) && (pi3==chi)) {
1003  // qbar-q-chi
1004  return -GEW*T3[si]*M(si)/MW*Matrix(GAS(5),DI(fi1),DI(fi2))*SP(TI(fi1),TI(fi2));
1005  } else if( (is_lbar_l(pi1, pi2) || is_nbar_n(pi1, pi2)) && (pi3==chi)) {
1006  // lbar-l-chi & nbar-n-chi
1007  return -GEW*T3[si]*M(si)/MW*Matrix(GAS(5),DI(fi1),DI(fi2));
1008  } else if(is_uqbar_dq(pi1, pi2) && pi3==phip) {
1009  // uqbar-dq-phip
1010  auto ckm = CKM(si1,si2);
1011  return I*GEW/sqrt2*(M(si1)/MW*Matrix(GAS(1)-GAS(5),DI(fi1),DI(fi2))/2-M(si2)/MW*Matrix(GAS(1)+GAS(5),DI(fi1),DI(fi2))/2)*ckm*SP(TI(fi1),TI(fi2));
1012  } else if(is_dqbar_uq(pi1, pi2) && pi3==phim) {
1013  // dqbar-uq-phim
1014  auto ckm = CKM(si1,si2);
1015  return I*GEW/sqrt2*(M(si2)/MW*Matrix(GAS(1)+GAS(5),DI(fi1),DI(fi2))/2-M(si1)/MW*Matrix(GAS(1)-GAS(5),DI(fi1),DI(fi2))/2)*ckm*SP(TI(fi1),TI(fi2));
1016  } else if(is_nbar_l(pi1, pi2) && pi3==phip) {
1017  // nbar-l-phip
1018  return -I*GEW/sqrt2*M(si2)/MW*Matrix(GAS(1)+GAS(5),DI(fi1),DI(fi2))/2; // Mnl = 0
1019  } else if(is_lbar_n(pi1, pi2) && pi3==phim) {
1020  // lbar-n-phim
1021  return -I*GEW/sqrt2*M(si1)/MW*Matrix(GAS(1)-GAS(5),DI(fi1),DI(fi2))/2; // Mnl = 0
1022  } else if(pi1==A && pi2==phip && pi3==phim) {
1023  // A-phip-phim
1024  return -I*eta_e*EL*SP(mom2-mom3,LI(fi1));
1025  } else if(pi1==Z && pi2==phip && pi3==phim) {
1026  // Z-phip-phim
1027  return -I*eta*eta_Z*GEW*C2W/(2*CW)*SP(mom2-mom3,LI(fi1));
1028  } else if(pi1==Wp && pi2==phim && pi3==H) {
1029  // Wp-phim-H
1030  return I/2*eta*GEW*SP(mom2-mom3,LI(fi1));
1031  } else if(pi1==Wm && pi2==phip && pi3==H) {
1032  // Wm-phip-H
1033  return -I/2*eta*GEW*SP(mom2-mom3,LI(fi1));
1034  } else if( ((pi1==Wp && pi2==phim) || (pi1==Wm && pi2==phip)) && pi3==chi) {
1035  // Wp-phim-chi & Wm-phip-chi
1036  return -eta*GEW/2*SP(mom2-mom3,LI(fi1));
1037  } else if(pi1==Z && pi2==chi && pi3==H) {
1038  // Z-chi-H
1039  return -eta*eta_Z*GEW/(2*CW)*SP(mom2-mom3,LI(fi1));
1040  } else if( ((pi1==phim && pi2==Wp) || (pi1==phip && pi2==Wm)) && pi3==A) {
1041  // phim-Wp-A & phip-Wm-A
1042  return I*eta_e*eta*EL*MW*SP(LI(fi2),LI(fi3));
1043  } else if( ((pi1==phim && pi2==Wp) || (pi1==phip && pi2==Wm)) && pi3==Z) {
1044  // phim-Wp-Z & phip-Wm-Z
1045  return -I*eta_Z*GEW*MZ*SW2*SP(LI(fi2),LI(fi3));
1046  } else if(pi1==H && pi2==Wp && pi3==Wm) {
1047  // H-Wp-Wm
1048  return I*GEW*MW*SP(LI(fi2),LI(fi3));
1049  } else if(pi1==H && pi2==Z && pi3==Z) {
1050  // H-Z-Z
1051  return I*GEW/CW*MZ*SP(LI(fi2),LI(fi3));
1052  } else if(pi1==H && pi2==phim && pi3==phip) {
1053  // H-phim-phip
1054  return -I/2*GEW*MH2/MW;
1055  } else if(pi1==H && pi2==H && pi3==H) {
1056  // H-H-H
1057  return -3/ex(2)*I*GEW*MH2/MW;
1058  } else if(pi1==H && pi2==chi && pi3==chi) {
1059  // H-chi-chi
1060  return -I/2*GEW*MH2/MW;
1061  } else if(pi1==ghWpbar && pi2==ghWp && pi3==A) {
1062  // ghWpbar-ghWp-A
1063  return I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1064  } else if(pi1==ghWmbar && pi2==ghWm && pi3==A) {
1065  // ghWmbar-ghWm-A
1066  return -I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1067  } else if(pi1==ghWpbar && pi2==ghWp && pi3==Z) {
1068  // ghWpbar-ghWp-Z
1069  return I*eta_G*eta*eta_Z*GEW*CW*SP(mom1,LI(fi3));
1070  } else if(pi1==ghWmbar && pi2==ghWm && pi3==Z) {
1071  // ghWmbar-ghWm-Z
1072  return -I*eta_G*eta*eta_Z*GEW*CW*SP(mom1,LI(fi3));
1073  } else if(pi1==ghWpbar && pi2==ghZ && pi3==Wp) {
1074  // ghWpbar-ghZ-Wp
1075  return -I*eta_G*eta*eta_Z*GEW*CW*SP(mom1,LI(fi3));
1076  } else if(pi1==ghWmbar && pi2==ghZ && pi3==Wm) {
1077  // ghWmbar-ghZ-Wm
1078  return I*eta_G*eta*eta_Z*GEW*CW*SP(mom1,LI(fi3));
1079  } else if(pi1==ghWpbar && pi2==ghA && pi3==Wp) {
1080  // ghWpbar-ghA-Wp
1081  return -I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1082  } else if(pi1==ghWmbar && pi2==ghA && pi3==Wm) {
1083  // ghWmbar-ghA-Wm
1084  return I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1085  } else if(pi1==ghZbar && pi2==ghWp && pi3==Wm) {
1086  // ghZbar-ghWp-Wm
1087  return -I*eta_G*eta*GEW*CW*SP(mom1,LI(fi3));
1088  } else if(pi1==ghZbar && pi2==ghWm && pi3==Wp) {
1089  // ghZbar-ghWm-Wp
1090  return I*eta_G*eta*GEW*CW*SP(mom1,LI(fi3));
1091  } else if(pi1==ghAbar && pi2==ghWp && pi3==Wm) {
1092  // ghAbar-ghWp-Wm
1093  return -I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1094  } else if(pi1==ghAbar && pi2==ghWm && pi3==Wp) {
1095  // ghAbar-ghWm-Wp
1096  return I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1097  } else if(pi1==ghWpbar && pi2==ghWp && pi3==chi) {
1098  // ghWpbar-ghWp-chi
1099  return eta_G*GEW/2*xi*MW;
1100  } else if(pi1==ghWmbar && pi2==ghWm && pi3==chi) {
1101  // ghWmbar-ghWm-chi
1102  return -eta_G*GEW/2*xi*MW;
1103  } else if( ((pi1==ghWpbar && pi2==ghWp) ||(pi1==ghWmbar && pi2==ghWm)) && pi3==H) {
1104  // ghWpbar-ghWp-H & ghWmbar-ghWm-H
1105  return -I/2*eta_G*GEW*xi*MW;
1106  } else if(pi1==ghZbar && pi2==ghZ && pi3==H) {
1107  // ghZbar-ghZ-H
1108  return -eta_G*I*GEW/(2*CW)*xi*MZ;
1109  } else if( pi1==ghZbar && ((pi2==ghWp && pi3==phim) || (pi2==ghWm && pi3==phip)) ) {
1110  // ghZbar-ghWp-phim & ghZbar-ghWm-phip
1111  return I/2*eta_G*eta_Z*GEW*xi*MZ;
1112  } else if( ((pi1==ghWpbar && pi3==phip) || (pi1==ghWmbar && pi3==phim)) && pi2==ghZ) {
1113  // ghWpbar-ghZ-phip & ghWmbar-ghZ-phim
1114  return -I*eta_G*eta_Z*GEW*C2W/(2*CW)*xi*MW;
1115  } else if( ((pi1==ghWpbar && pi3==phip) || (pi1==ghWmbar && pi3==phim)) && pi2==ghA) {
1116  // ghWpbar-ghA-phip & ghWmbar-ghA-phim
1117  return -I*eta_G*eta_e*eta*EL*xi*MW;
1118  } else {
1119  cout << endl << e << endl;
1120  throw Error("Vertex Not defined!");
1121  }
1122  } else if(nn==4) {
1123  auto pi4 = e.op(3).op(0);
1124  auto fi4 = e.op(3).op(1);
1125  auto mom4 = e.op(3).op(2);
1126 
1127  if(pi1==g && pi2==g && pi3==g && pi4==g) {
1128  // g^4
1129  return Gluon4Vertex(e, eta_s);
1130  } else if(pi1==Wp && pi2==Wm && pi3==A && pi4==A) {
1131  // Wp-Wm-A-A
1132  return -I*EL2*(2*SP(LI(fi1),LI(fi2))*SP(LI(fi3),LI(fi4))-SP(LI(fi1),LI(fi3))*SP(LI(fi2),LI(fi4))-SP(LI(fi1),LI(fi4))*SP(LI(fi2),LI(fi3)));
1133  } else if(pi1==Wp && pi2==Wm && pi3==Z && pi4==Z) {
1134  // Wp-Wm-Z-Z
1135  return -I*GEW2*CW2*(2*SP(LI(fi1),LI(fi2))*SP(LI(fi3),LI(fi4))-SP(LI(fi1),LI(fi3))*SP(LI(fi2),LI(fi4))-SP(LI(fi1),LI(fi4))*SP(LI(fi2),LI(fi3)));
1136  } else if(pi1==Wp && pi2==Wm && pi3==A && pi4==Z) {
1137  // Wp-Wm-A-Z
1138  return -I*eta_e*eta*eta_Z*EL*GEW*CW*(2*SP(LI(fi1),LI(fi2))*SP(LI(fi3),LI(fi4))-SP(LI(fi1),LI(fi3))*SP(LI(fi2),LI(fi4))-SP(LI(fi1),LI(fi4))*SP(LI(fi2),LI(fi3)));
1139  } else if(pi1==Wp && pi2==Wp && pi3==Wm && pi4==Wm) {
1140  // Wp-Wp-Wm-Wm
1141  return I*GEW2*(2*SP(LI(fi1),LI(fi2))*SP(LI(fi3),LI(fi4))-SP(LI(fi1),LI(fi3))*SP(LI(fi2),LI(fi4))-SP(LI(fi1),LI(fi4))*SP(LI(fi2),LI(fi3)));
1142  } else if(pi1==Wp && pi2==Wm && pi3==H && pi4==H) {
1143  // Wp-Wm-H-H
1144  return I/2*GEW2*SP(LI(fi1),LI(fi2));
1145  } else if(pi1==Wp && pi2==Wm && pi3==chi && pi4==chi) {
1146  // Wp-Wm-chi-chi
1147  return I/2*GEW2*SP(LI(fi1),LI(fi2));
1148  } else if(pi1==Z && pi2==Z && pi3==H && pi4==H) {
1149  // Z-Z-H-H
1150  return I/2*GEW2/CW2*SP(LI(fi1),LI(fi2));
1151  } else if(pi1==Z && pi2==Z && pi3==chi && pi4==chi) {
1152  // Z-Z-chi-chi
1153  return I/2*GEW2/CW2*SP(LI(fi1),LI(fi2));
1154  } else if(pi1==A && pi2==A && pi3==phip && pi4==phim) {
1155  // A-A-phip-phim
1156  return 2*I*EL2*SP(LI(fi1),LI(fi2));
1157  } else if(pi1==Z && pi2==Z && pi3==phip && pi4==phim) {
1158  // Z-Z-phip-phim
1159  return I/2*pow(GEW*C2W/CW,2)*SP(LI(fi1),LI(fi2));
1160  } else if(pi1==Wp && pi2==Wm && pi3==phip && pi4==phim) {
1161  // Wp-Wm-phip-phim
1162  return I/2*GEW2*SP(LI(fi1),LI(fi2));
1163  } else if( ((pi1==Wp && pi3==phim) || (pi1==Wm && pi3==phip)) && pi2==Z && pi4==H) {
1164  // Wp-Z-phim-H & Wm-Z-phip-H
1165  return -I*eta_Z*GEW2*SW2/(2*CW)*SP(LI(fi1),LI(fi2));
1166  } else if(pi1==Wm && pi2==Z && pi3==phip && pi4==chi) {
1167  // Wm-Z-phip-chi
1168  return -eta_Z*GEW2*SW2/(2*CW)*SP(LI(fi1),LI(fi2));
1169  } else if(pi1==Wp && pi2==Z && pi3==phim && pi4==chi) {
1170  // Wp-Z-phim-chi
1171  return eta_Z*GEW2*SW2/(2*CW)*SP(LI(fi1),LI(fi2));
1172  } else if( ((pi1==Wm && pi3==phip) || (pi1==Wp && pi3==phim)) && pi2==A && pi4==H) {
1173  // Wm-A-phip-H & Wp-A-phim-H
1174  return I/2*eta_e*eta*EL*GEW*SP(LI(fi1),LI(fi2));
1175  } else if(pi1==Wp && pi2==A && pi3==phim && pi4==chi) {
1176  // Wp-A-phim-chi
1177  return -1/ex(2)*eta_e*eta*EL*GEW*SP(LI(fi1),LI(fi2));
1178  } else if(pi1==Wm && pi2==A && pi3==phip && pi4==chi) {
1179  // Wm-A-phip-chi
1180  return 1/ex(2)*eta_e*eta*EL*GEW*SP(LI(fi1),LI(fi2));
1181  } else if(pi1==Z && pi2==A && pi3==phip && pi4==phim) {
1182  // Z-A-phip-phim
1183  return I*eta_e*eta*eta_Z*EL*GEW*C2W/CW*SP(LI(fi1),LI(fi2));
1184  } else if(pi1==phim && pi2==phip && pi3==phim && pi4==phip) {
1185  // phim-phip-phim-phip
1186  return -I/2*GEW2*MH2/MW2;
1187  } else if(pi1==H && pi2==H && pi3==phim && pi4==phip) {
1188  // H-H-phim-phip
1189  return -I/4*GEW2*MH2/MW2;
1190  } else if(pi1==chi && pi2==chi && pi3==phim && pi4==phip) {
1191  // chi-chi-phim-phip
1192  return -I/4*GEW2*MH2/MW2;
1193  } else if(pi1==H && pi2==H && pi3==H && pi4==H) {
1194  // H-H-H-H
1195  return -3/ex(4)*I*GEW2*MH2/MW2;
1196  } else if(pi1==H && pi2==H && pi3==chi && pi4==chi) {
1197  // H-H-chi-chi
1198  return -I/4*GEW2*MH2/MW2;
1199  } else if(pi1==chi && pi2==chi && pi3==chi && pi4==chi) {
1200  // chi-chi-chi-chi
1201  return -3/ex(4)*I*GEW2*MH2/MW2;
1202  } else {
1203  cout << endl << e << endl;
1204  throw Error("Vertex Not defined!");
1205  }
1206  } else {
1207  cout << endl << e << endl;
1208  throw Error("Vertex Not defined!");
1209  }
1210  } else return e.map(self);
1211  return e;
1212  });
1213  return fr(amp);
1214  }
1215 
1216 
1217  lst Models::FeynRulesQCD(const lst & amps, const ex & xi) {
1218  lst ret;
1219  for(auto item : amps) ret.append(FeynRulesQCD(item,xi));
1220  return ret;
1221  }
1222  ex Models::FeynRulesQCD(const ex & amp, const ex & xi) {
1223  if(is_a<lst>(amp)) return FeynRulesQCD(ex_to<lst>(amp), xi);
1224 
1225  static Symbol A("A"), Q("Q"), Qbar("Qbar"), q("q"), qbar("qbar"), g("g"), gh("gh"), ghbar("ghbar");
1226  static Symbol m("m"), eq("eq"), eQ("eQ");
1227 
1228  auto fr = MapFunction([&](const ex &e, MapFunction &self)->ex {
1229  if(isFunction(e,"OutField") || isFunction(e,"InField")) return 1;
1230  else if(isFunction(e, "Propagator")) {
1231  if(e.op(0).op(0)==q) {
1232  return QuarkPropagator(e, 0);
1233  } else if(e.op(0).op(0)==Q) {
1234  return QuarkPropagator(e, m);
1235  } else if(e.op(0).op(0)==g) {
1236  return GluonPropagator(e, xi);
1237  } else if(e.op(0).op(0)==gh) {
1238  return GluonGhostPropagator(e, xi);
1239  } else {
1240  cout << "expr: " << e << endl;
1241  throw Error("Propagator Not Found!");
1242  }
1243  } else if(isFunction(e, "Vertex")) {
1244  if(e.nops()==3 && ((e.op(0).op(0)==qbar && e.op(1).op(0)==q) || (e.op(0).op(0)==Qbar && e.op(1).op(0)==Q)) && (e.op(2).op(0)==g || e.op(2).op(0)==A) ) {
1245  // qbar-q-g or Qbar-Q-g or g -> A
1246  if(e.op(2).op(0)==g) return QuarkGluonVertex(e);
1247  else {
1248  auto fi1 = e.op(0).op(1);
1249  auto fi2 = e.op(1).op(1);
1250  auto fi3 = e.op(2).op(1);
1251  if(e.op(1).op(0)==q) QuarkAVertex(e, eq);
1252  else if(e.op(1).op(0)==Q) QuarkAVertex(e, eQ);
1253  else throw Error("Vertex Error.");
1254  }
1255  } else if(e.nops()==3 && e.op(0).op(0)==ghbar && e.op(1).op(0)==gh) {
1256  // ghbar-gh-g
1257  return GhostGluonVertex(e);
1258  } else if(e.nops()==3 && e.op(0).op(0)==g && e.op(1).op(0)==g && e.op(2).op(0)==g) {
1259  // g^3
1260  return Gluon3Vertex(e);
1261  } else if(e.nops()==4 && e.op(0).op(0)==g && e.op(1).op(0)==g && e.op(2).op(0)==g && e.op(3).op(0)==g) {
1262  // g^4
1263  return Gluon4Vertex(e);
1264  } else {
1265  cout << "expr: " << e << endl;
1266  throw Error("Vertex Not Found!");
1267  }
1268  } else return e.map(self);
1269  return e;
1270  });
1271  return fr(amp);
1272  }
1273 
1274 }
1275 
QGRAF.h
QGRAF header file.
HepLib::Error
class used to wrap error message
Definition: BASIC.h:242
HepLib::Index
class for index object
Definition: HEP.h:104
HepLib::Index::has
static bool has(const ex &e)
HepLib::MapFunction
class to wrap map_function of GiNaC
Definition: BASIC.h:632
HepLib::Parser
class to parse for string or file, helpful with Symbol class
Definition: BASIC.h:645
HepLib::Parser::Read
ex Read(const string &instr, bool s2S=true)
Definition: Functions.cpp:111
HepLib::QGRAF::Process::Style
static string Style
Definition: QGRAF.h:77
HepLib::QGRAF::Process::Others
vector< string > Others
Definition: QGRAF.h:84
HepLib::QGRAF::Process::Amplitudes
lst Amplitudes(symtab st)
generte the Amplitudes
Definition: QGRAF.cpp:63
HepLib::SUNF4
class for SUNF4 object
Definition: HEP.h:278
HepLib::SUNF
class for SUNF object
Definition: HEP.h:233
HepLib::SUNT
class for SUNT object
Definition: HEP.h:189
HepLib::Symbol
class extended to GiNaC symbol class, represent a positive symbol
Definition: BASIC.h:113
HepLib::QGRAF
namespace for generating Feynman diagrams or amplitudes.
Definition: QGRAF.cpp:8
HepLib::QGRAF::AntiQuarkSumL
ex AntiQuarkSumL(int qi, ex p, ex m, bool color)
polarization sum for anti-quark
Definition: QGRAF.cpp:676
HepLib::QGRAF::HasLoop
bool HasLoop(ex amp, lst prop)
check a amplitude has a loop w.r.t. propapagtor
Definition: QGRAF.cpp:373
HepLib::QGRAF::AntiGhostSumL
ex AntiGhostSumL(int qi)
polarization sum for anti-ghost, note that we will add extra - here
Definition: QGRAF.cpp:717
HepLib::QGRAF::IndexL2R
ex IndexL2R(ex e, bool all)
Change Index from left to right, only affect li/di/ci/ti, external index start with dim/lim/cim/tim w...
Definition: QGRAF.cpp:598
HepLib::QGRAF::GluonPropagator
ex GluonPropagator(const ex &e, const ex &xi)
Propagator for gluon.
Definition: QGRAF.cpp:436
HepLib::QGRAF::LeptonAVertex
ex LeptonAVertex(const ex &e, const ex &eq, const ex &eta_e)
l-lbar-A vertex
Definition: QGRAF.cpp:585
HepLib::QGRAF::Gluon3Vertex
ex Gluon3Vertex(const ex &e, const ex &eta_s)
g-g-g vertex
Definition: QGRAF.cpp:517
HepLib::QGRAF::GhostSumL
ex GhostSumL(int qi)
polarization sum for ghost
Definition: QGRAF.cpp:699
HepLib::QGRAF::RDI
Index RDI(ex fn)
Definition: QGRAF.cpp:51
HepLib::QGRAF::AntiGhostSumR
ex AntiGhostSumR(int qi)
polarization sum for anti-ghost, note that we will add extra - here
Definition: QGRAF.cpp:726
HepLib::QGRAF::DrawPDF
void DrawPDF(const lst &amps, string fn)
generate Feynman diagrams for the amplitudes, in PDF format
Definition: QGRAF.cpp:165
HepLib::QGRAF::QuarkGluonVertex
ex QuarkGluonVertex(const ex &e, const ex &eta_s)
q-qbar-g vertex
Definition: QGRAF.cpp:504
HepLib::QGRAF::ScalarPropagator
ex ScalarPropagator(const ex &e, const ex &m, const ex &xi)
Propagator for scalar.
Definition: QGRAF.cpp:491
HepLib::QGRAF::J1SumL
ex J1SumL(int qi, ex p)
polarization sum for total angular momentum
Definition: QGRAF.cpp:736
HepLib::QGRAF::RTI
Index RTI(ex fn)
Definition: QGRAF.cpp:53
HepLib::QGRAF::LeptonPropagator
ex LeptonPropagator(const ex &e, const ex &m)
Propagator for lepton.
Definition: QGRAF.cpp:423
HepLib::QGRAF::ShrinkCut
vector< lst > ShrinkCut(ex amp, lst prop, int n)
cut the amplitude, and return the connected parts, may have many different cuts
Definition: QGRAF.cpp:285
HepLib::QGRAF::AntiQuarkSumR
ex AntiQuarkSumR(int qi, ex p, ex m, bool color)
polarization sum for anti-quark
Definition: QGRAF.cpp:689
HepLib::QGRAF::RCI
Index RCI(ex fn)
Definition: QGRAF.cpp:55
HepLib::QGRAF::J1SumR
ex J1SumR(int qi, ex p)
polarization sum for total angular momentum
Definition: QGRAF.cpp:747
HepLib::QGRAF::CI
Index CI(ex fn)
Definition: QGRAF.cpp:49
HepLib::QGRAF::DI
Index DI(ex fn)
Definition: QGRAF.cpp:45
HepLib::QGRAF::RLI
Index RLI(ex fn)
Definition: QGRAF.cpp:52
HepLib::QGRAF::AZWGhostPropagator
ex AZWGhostPropagator(const ex &e, const ex &m, const ex &xi, const ex &eta_G)
Propagator for A/Z/W ghost.
Definition: QGRAF.cpp:477
HepLib::QGRAF::Gluon4Vertex
ex Gluon4Vertex(const ex &e, const ex &eta_s)
g-g-g-g vertex
Definition: QGRAF.cpp:537
HepLib::QGRAF::IndexCC
ex IndexCC(ex e, bool all)
Definition: QGRAF.cpp:616
HepLib::QGRAF::GhostGluonVertex
ex GhostGluonVertex(const ex &e, const ex &eta_s, const ex &eta_G)
ghbar-gh-g vertex
Definition: QGRAF.cpp:556
HepLib::QGRAF::RAI
Index RAI(ex fn)
Definition: QGRAF.cpp:56
HepLib::QGRAF::QuarkPropagator
ex QuarkPropagator(const ex &e, const ex &m)
Propagator for quark.
Definition: QGRAF.cpp:410
HepLib::QGRAF::AI
Index AI(ex fn)
Definition: QGRAF.cpp:50
HepLib::QGRAF::QuarkSumR
ex QuarkSumR(int qi, ex p, ex m, bool color)
polarization sum for quark
Definition: QGRAF.cpp:663
HepLib::QGRAF::QuarkSumL
ex QuarkSumL(int qi, ex p, ex m, bool color)
polarization sum for quark
Definition: QGRAF.cpp:650
HepLib::QGRAF::GhostSumR
ex GhostSumR(int qi)
polarization sum for ghost
Definition: QGRAF.cpp:708
HepLib::QGRAF::LI
Index LI(ex fn)
Definition: QGRAF.cpp:46
HepLib::QGRAF::AZWPropagator
ex AZWPropagator(const ex &e, const ex &m, const ex &xi)
Propagator for A/Z/W boson.
Definition: QGRAF.cpp:463
HepLib::QGRAF::QuarkAVertex
ex QuarkAVertex(const ex &e, const ex &eq, const ex &eta_e)
qbar-q-A vertex
Definition: QGRAF.cpp:571
HepLib::QGRAF::TI
Index TI(ex fn)
Definition: QGRAF.cpp:47
HepLib::QGRAF::GluonGhostPropagator
ex GluonGhostPropagator(const ex &e, const ex &xi, const ex &eta_G)
Propagator for gluon ghost.
Definition: QGRAF.cpp:449
HepLib::QGRAF::FI
Index FI(ex fn)
Definition: QGRAF.cpp:48
HepLib::QGRAF::RFI
Index RFI(ex fn)
Definition: QGRAF.cpp:54
HepLib::QGRAF::TopoLines
lst TopoLines(const ex &amp)
generate the topological lines for the amplitude
Definition: QGRAF.cpp:120
HepLib::QGRAF::GluonSumR
ex GluonSumR(int qi, bool color)
polarization sum for gluon
Definition: QGRAF.cpp:637
HepLib::QGRAF::GluonSumL
ex GluonSumL(int qi, bool color)
polarization sum for gluon
Definition: QGRAF.cpp:626
HepLib::SD::do_not_evalf_params
do_not_evalf_params().expl_derivative_func(zd1D).derivative_func(zp1D)) REGISTER_FUNCTION(FTX
HepLib::gs
const Symbol gs
HepLib::Combinations
void Combinations(int n, int m, std::function< void(const int *)> f)
Definition: Functions.cpp:184
HepLib::ep
const Symbol ep
HepLib::GAS
ex GAS(const ex &expr, unsigned rl)
function similar to GAD/GSD in FeynClac
Definition: DGamma.cpp:246
HepLib::REGISTER_FUNCTION
REGISTER_FUNCTION(Matrix, do_not_evalf_params().print_func< FCFormat >(&Matrix_fc_print).conjugate_func(mat_conj).set_return_type(return_types::commutative)) bool IsZero(const ex &e)
Definition: Basic.cpp:715
HepLib::Debug
bool Debug
Definition: Init.cpp:141
HepLib::w
ex w
Definition: Init.cpp:90
HepLib::vs
const Symbol vs
HepLib::isFunction
bool isFunction(const ex &e, string func_name)
Definition: BASIC.h:658
HepLib::CA
const Symbol CA
HepLib::CF
const Symbol CF
HepLib::InstallPrefix
string InstallPrefix
Definition: Init.cpp:159
HepLib::string_replace_all
void string_replace_all(string &str, const string &from, const string &to)
Definition: Functions.cpp:148
HepLib::ex2str
string ex2str(const ex &expr)
convert ex to output string, the defalut printer format will be used
Definition: BASIC.cpp:715
HepLib::dir_exists
bool dir_exists(string dir)
Definition: BASIC.h:297
HepLib::GiNaC_Parallel
exvector GiNaC_Parallel(int ntotal, std::function< ex(int)> f, const string &key)
GiNaC Parallel Evaluation using fork.
Definition: BASIC.cpp:260
HepLib::SP
ex SP(const ex &a, bool use_map=false)
Definition: Pair.cpp:166
Generated by doxygen 1.9.1