HepLib
Loading...
Searching...
No Matches
QGRAF.cpp
Go to the documentation of this file.
1
6#include "QGRAF.h"
7
8namespace 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)) * GMat(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 * GMat(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*GMat(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*GMat(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*GMat(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 GMat(GAS(p)+m*GAS(1), RDI(qi), DI(qi)) * SP(RTI(qi), TI(qi));
652 else return GMat(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 GMat(GAS(p)+m*GAS(1), DI(qi), RDI(qi)) * SP(TI(qi),RTI(qi));
665 else return GMat(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 GMat(GAS(p)-m*GAS(1), DI(qi), RDI(qi)) * SP(TI(qi), RTI(qi));
678 else return GMat(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 GMat(GAS(p)-m*GAS(1), RDI(qi), DI(qi)) * SP(RTI(qi), TI(qi));
691 else return GMat(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)*GMat(GAS(LI(fi3)),DI(fi1),DI(fi2))-gfA(si)*GMat(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)*GMat(GAS(LI(fi3)),DI(fi1),DI(fi2))-gfA(si)*GMat(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*GMat(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*GMat(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*GMat(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*GMat(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*GMat(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*GMat(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*GMat(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*GMat(GAS(1)-GAS(5),DI(fi1),DI(fi2))/2-M(si2)/MW*GMat(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*GMat(GAS(1)+GAS(5),DI(fi1),DI(fi2))/2-M(si1)/MW*GMat(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*GMat(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*GMat(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::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:257
HepLib::Debug
bool Debug
Definition Init.cpp:144
HepLib::w
ex w
Definition Init.cpp:93
HepLib::vs
const Symbol vs
HepLib::isFunction
bool isFunction(const ex &e, string func_name)
Definition BASIC.h:658
HepLib::InstallPrefix
string InstallPrefix
Definition Init.cpp:162
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::REGISTER_FUNCTION
REGISTER_FUNCTION(GMat, do_not_evalf_params().print_func< FCFormat >(&GMat_fc_print).conjugate_func(mat_conj).set_return_type(return_types::commutative)) bool IsZero(const ex &e)
Definition Basic.cpp:1002
HepLib::SP
ex SP(const ex &a, bool use_map=false)
Definition Pair.cpp:166
Generated by doxygen 1.9.8