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),d); }
46 Index LI(ex fn) { return Index(li_+n2s(fn),d); }
47 Index TI(ex fn) { return Index(ti_+n2s(fn),NF); }
48 Index FI(ex fn) { return Index(fi_+n2s(fn),NF); }
49 Index CI(ex fn) { return Index(ci_+n2s(fn),NA); }
50 Index AI(ex fn) { return Index(ai_+n2s(fn),NA); }
51 Index RDI(ex fn) { return Index(rdi_+n2s(fn),d); }
52 Index RLI(ex fn) { return Index(rli_+n2s(fn),d); }
53 Index RTI(ex fn) { return Index(rti_+n2s(fn),NF); }
54 Index RFI(ex fn) { return Index(rfi_+n2s(fn),NF); }
55 Index RCI(ex fn) { return Index(rci_+n2s(fn),NA); }
56 Index RAI(ex fn) { return Index(rai_+n2s(fn),NA); }
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) {
86 if(vs.find(";") != std::string::npos) ofs << vs << endl;
87 else ofs << vs << ";" << endl;
88 }
89 ofs.close();
90
91 if(Debug) rc = system((InstallPrefix+"/bin/qgraf").c_str());
92 else rc = system((InstallPrefix+"/bin/qgraf > /dev/null").c_str());
93
94 ifstream ifs("qgraf.out");
95 string ostr((istreambuf_iterator<char>(ifs)), (istreambuf_iterator<char>()));
96 ifs.close();
97
98 if(!Debug) {
99 if(access("qgraf.dat",F_OK)!=-1) remove("qgraf.dat");
100 if(access("qgraf.mod",F_OK)!=-1) remove("qgraf.mod");
101 if(access("qgraf.sty",F_OK)!=-1) remove("qgraf.sty");
102 if(access("qgraf.out",F_OK)!=-1) remove("qgraf.out");
103 }
104
105 const char* rm_chars = " \t\v\r\n,";
106 if(!ostr.empty()) {
107 ostr.erase(0, ostr.find_first_not_of(rm_chars));
108 ostr.erase(ostr.find_last_not_of(rm_chars)+1);
109 }
110 ostr = "{" + ostr + "}";
111
112 Parser amp(st);
113 ex ret = amp.Read(ostr);
114
115 return ex_to<lst>(ret);
116 }
117
123 lst TopoLines(const ex & amp) {
124 map<ex,int,ex_is_less> v2id, fid2vid;
125 map<int,ex> vid2fs; // fileds in the vertex
126 int cid = 0;
127 lst lines;
128 for(const_preorder_iterator i = amp.preorder_begin(); i != amp.preorder_end(); ++i) {
129 ex e = (*i);
130 if(isFunction(e,"OutField")) {
131 lines.append(lst{ e.op(1), iWF(e.op(1)), lst{e.op(0)}, e.op(2) });
132 } else if(isFunction(e,"InField")) {
133 lines.append(lst{ iWF(e.op(1)), e.op(1), lst{e.op(0)}, e.op(2) });
134 } else if(isFunction(e, "Propagator")) {
135 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) });
136 } else if(isFunction(e, "Vertex")) {
137 auto itr = v2id.find(e);
138 if(itr==v2id.end()) {
139 cid++;
140 v2id.emplace(make_pair(e,cid));
141 lst fs;
142 for(auto f : e) {
143 ex fid = f.op(1);
144 fid2vid[fid] = cid;
145 fs.append(f.op(0));
146 }
147 vid2fs[cid] = fs;
148 }
149 }
150 }
151 lines = ex_to<lst>(MapFunction([&vid2fs,&fid2vid](const ex &e, MapFunction &self)->ex{
152 if(!e.has(iWF(w))) return e;
153 else if(e.match(iWF(w))) {
154 auto vid = fid2vid[e.op(0)];
155 return lst{vid, vid2fs[vid]};
156 } else return e.map(self);
157 })(lines));
158
159 return lines.sort();
160 }
161
169 void DrawPDF(const lst & amps, string fn, int nr, bool single_page) {
170 int id=0, rc;
171 exvector amp_vec;
172 for(auto item : amps) amp_vec.push_back(item);
173 string tex_path = to_string(getpid()) + "_TeX/";
174 if(!dir_exists(tex_path)) rc = system(("mkdir -p "+tex_path).c_str());
175
176 GiNaC_Parallel(amp_vec.size(), [&amp_vec,single_page,tex_path](int idx)->ex {
177 auto amp = amp_vec[idx];
178 ofstream out(tex_path+to_string(idx)+".tex");
179 out << "\\documentclass[tikz]{standalone}" << endl;
180 out << "\\usepackage{tikz-feynman}" << endl;
181 out << "\\tikzfeynmanset{compat=1.1.0}" << endl;
182 out << "\\begin{document}" << endl;
183 out << "\\feynmandiagram{" << endl;
184 auto lines = TopoLines(amp);
185
186 map<ex,int,ex_is_less> tot_bend_map;
187 for(auto l : lines) {
188 lst ll = lst{l.op(0), l.op(1)};
189 ll.sort();
190 tot_bend_map[ll]++;
191 }
192
193 map<ex,int,ex_is_less> bend_map;
194 std::map<ex,int,ex_is_less> vtex_map; // vertex option, only once
195 for(auto l : lines) {
196 lst ll = lst{l.op(0), l.op(1)};
197 bool isExt = (is_a<numeric>(ll.op(0)) && ll.op(0)<0) || (is_a<numeric>(ll.op(1)) && ll.op(1)<0);
198 lst oll = ll;
199 ll.sort();
200 bool ll_same = oll.is_equal(ll);
201 bend_map[ll]++;
202
203 auto fidL = (is_a<numeric>(l.op(1)) ? l.op(1) : l.op(1).op(0));
204 out << "\"" << fidL << "\"";
205 if(fidL<0) {
206 out << "[particle=";
207 if(InOutTeX[fidL].length()>0) out << InOutTeX[fidL];
208 else out << fidL;
209 out << "]";
210 } else if(vtex_map[l.op(1).op(1)]==0) {
211 out << VerTeX[l.op(1).op(1)];
212 vtex_map[l.op(1).op(1)]=1;
213 }
214
215 out << " --[";
216 auto f = l.op(2).op(0);
217 if(LineTeX[f].length()>0) {
218 if(!isExt) out << LineTeX[f];
219 else {
220 auto cpos = LineTeX[f].find(", edge");
221 if(cpos>0) out << LineTeX[f].substr(0,cpos);
222 else out << LineTeX[f];
223 }
224 }
225 if(tot_bend_map[ll]>1) {
226 if((tot_bend_map[ll]%2)==0) { // even case
227 if((bend_map[ll]%2)==0) {
228 if(ll_same) out << ",half right";
229 else out << ",half left";
230 } else {
231 if(ll_same) out << ",half left";
232 else out << ",half right";
233 }
234 } else if(bend_map[ll]>2) { // odd case
235 if((bend_map[ll]%2)==0) {
236 if(ll_same) out << ",half right";
237 else out << ",half left";
238 } else {
239 if(ll_same) out << ",half left";
240 else out << ",half right";
241 }
242 }
243 }
244 if(is_zero(l.op(0)-l.op(1))) out << ",loop,distance=2cm";
245 out << "]";
246
247 auto fidR = (is_a<numeric>(l.op(0)) ? l.op(0) : l.op(0).op(0));
248 out << " \"" << fidR << "\"";
249 if(fidR<0) {
250 out << "[particle=";
251 if(InOutTeX[fidR].length()>0) out << InOutTeX[fidR];
252 else out<< fidR;
253 out << "]";
254 } else if(vtex_map[l.op(0).op(1)]==0) {
255 out << VerTeX[l.op(0).op(1)];
256 vtex_map[l.op(0).op(1)]=1;
257 }
258 out << ";" << endl;
259 }
260 out << "};" << endl;
261 out << "\\end{document}" << endl;
262 out.close();
263 auto rc = system(("cd "+tex_path+" && echo X | lualatex " + to_string(idx) + " 1>/dev/null").c_str());
264 return 0;
265 }, "TeX");
266
267 ofstream out(tex_path+"diagram.tex");
268 if(single_page) {
269 out << "\\let\\mypdfximage\\pdfximage" << endl;
270 out << "\\def\\pdfximage{\\immediate\\mypdfximage}" << endl;
271 out << "\\documentclass{standalone}" << endl;
272 } else {
273 out << "\\documentclass{article}" << endl;
274 }
275 out << "\\usepackage{graphicx}" << endl;
276 out << "\\usepackage{adjustbox}" << endl;
277 out << "\\usepackage{longtable}" << endl;
278 out << "\\usepackage{scalefnt}" << endl;
279 out << "\\begin{document}" << endl;
280 //out << "\\scalefont{" << .5/nr << "}" << endl;
281 if(single_page) {
282 out << "\\begin{adjustbox}{valign=T,width=\\textwidth}" << endl;
283 out << "\\begin{tabular}{|"; for(int i=0; i<nr; i++) out << "cc|"; out << "}" << endl;
284 } else {
285 out << "\\setlength\\LTleft{-2.5cm} \\setlength\\LTright{0pt}" << endl;
286 out << "\\begin{longtable}{|"; for(int i=0; i<nr; i++) out << "cc|"; out << "}" << endl;
287 }
288 out << "\\hline" << endl;
289 int total = amps.nops();
290 int namps = total;
291 if((total%nr)!=0) total = (total/nr+1)*nr;
292 for(int i=0; i<total; i++) {
293
294 //int limit = 300;
295 //if((i!=0) && (i+1!=total)) {
296 // out << "\\end{tabular}" << endl << endl;
297 // out << "\\begin{tabular}{|"; for(int i=0; i<nr; i++) out << "cc|"; out << "}" << endl;
298 // out << "\\hline" << endl;
299 //}
300
301 out << "{\\tiny " << i+1 << "}&" << endl;
302 if(i<namps) {
303 out << "\\includegraphics[keepaspectratio,";
304 out << "height=" << 1.0/nr << "\\textwidth,";
305 out << "width=" << 1.0/nr << "\\textwidth]";
306 out << "{"<<i<<".pdf}" << endl;
307 }
308 if((i+1)%nr==0) out << "\\\\ \\hline";
309 else out << "&";
310 }
311 if(single_page) {
312 out << "\\end{tabular}" << endl;
313 out << "\\end{adjustbox}" << endl;
314 } else {
315 out << "\\end{longtable}" << endl;
316 }
317 out << "\\end{document}" << endl;
318 out.close();
319 if(Debug) rc = system(("ulimit -n unlimited > /dev/null; cd "+tex_path+" && pdflatex diagram && mv diagram.pdf ../"+fn).c_str());
320 else rc = system(("ulimit -n unlimited > /dev/null; cd "+tex_path+" && echo X | pdflatex diagram 1>/dev/null && mv diagram.pdf ../"+fn).c_str());
321
322 if(!file_exists(fn)) {
323 cout << "DrawPDF failed, files not removed, please check the TeX files!" << endl;
324 cout << "latex command: pdflatex diagram.tex" << endl;
325 cout << "if there are many pdf files, one needs ulimint -n <large number>." << endl;
326 exit(1);
327 } else if(!Debug) {
328 rc = system(("rm -r "+tex_path).c_str());
329 }
330 }
331
339 vector<lst> ShrinkCut(ex amp, lst prop, int n) {
340 vector<lst> ret;
341 if(prop.nops()<1) throw Error("ShrinkCut: no cut provided!");
342 auto tls = TopoLines(amp);
343 vector<int> cls_vec;
344 for(int i=0; i<tls.nops(); i++) {
345 auto pi = tls.op(i).op(2);
346 if(pi.nops()<2) continue;
347 if(!is_a<lst>(prop.op(0))) { // prop : lst{ g, g }
348 if(is_zero(pi.op(0)-prop.op(0)) && is_zero(pi.op(1)-prop.op(1))) cls_vec.push_back(i);
349 } else {
350 for(auto iprop : prop) {
351 if(is_zero(pi.op(0)-iprop.op(0)) && is_zero(pi.op(1)-iprop.op(1))) cls_vec.push_back(i);
352 }
353 }
354 }
355 if(cls_vec.size()<n) return ret;
356
357 Combinations(cls_vec.size(), n, [n,&ret,cls_vec,tls](const int * is)->void {
358 int cls[n];
359 for(int i=0; i<n; i++) cls[i] = cls_vec[is[i]];
360
361 // cut each line into 2 half-lines, labeled with 0
362 auto tls2 = tls;
363 for(auto ci : cls) {
364 auto ol = tls2.op(ci);
365 tls2[ci] = lst{ol.op(0), 0, lst{ol.op(2).op(0)}, ol.op(3)};
366 tls2.append(lst{0, ol.op(1), lst{ol.op(2).op(1)}, ol.op(3)} );
367 }
368
369 // shrink internal lines
370 int last = 0;
371 int ntls2 = tls2.nops();
372 while(true) {
373 ex lp = 0;
374 for(int i=last; i<ntls2; i++) {
375 auto li = tls2.op(i);
376 if(is_zero(li) || li.op(2).nops()<2) continue;
377 last = i;
378 lp = li;
379 tls2[last] = 0;
380 break;
381 }
382 if(is_zero(lp)) break;
383
384 for(int i=0; i<ntls2; i++) {
385 if(is_zero(tls2.op(i))) continue;
386 if(is_zero(tls2.op(i).op(0)-lp.op(0))) tls2[i][0] = lp.op(1);
387 if(is_zero(tls2.op(i).op(1)-lp.op(0))) tls2[i][1] = lp.op(1);
388 }
389 }
390
391 // final connected parts
392 map<int, lst> con_map;
393 for(auto li : tls2) {
394 if(is_zero(li)) continue;
395 ex key, val;
396 ex fiL = li.op(0);
397 if(is_a<lst>(fiL)) fiL = fiL.op(0);
398 ex fiR = li.op(1);
399 if(is_a<lst>(fiR)) fiR = fiR.op(0);
400 if(fiL>0 && fiR<0) {
401 con_map[ex_to<numeric>(fiL).to_int()].append(fiR);
402 } else if(fiR>0 && fiL<0) {
403 con_map[ex_to<numeric>(fiR).to_int()].append(fiL);
404 } else {
405 if(fiL>0 && is_zero(fiR)) key = fiL;
406 else if(fiR>0 && is_zero(fiL)) key = fiR;
407 else throw Error("ShrinkCut: unexpcected point reached.");
408 val = li.op(2).op(0);
409 con_map[ex_to<numeric>(key).to_int()].append(val);
410 }
411 }
412
413 lst item;
414 for(auto kv : con_map) item.append(kv.second.sort());
415 ret.push_back(item);
416 });
417
418 return ret;
419 }
420
427 bool HasLoop(ex amp, lst prop) {
428 auto tls = TopoLines(amp);
429 int ntls = tls.nops();
430 // shrink internal lines of prop
431 int last = 0;
432 while(true) {
433 ex lp = 0;
434 for(int i=last; i<ntls; i++) {
435 auto li = tls.op(i);
436 if(is_zero(li) || li.op(2).nops()<2) continue; // external line
437 auto cpi = li.op(2);
438 bool m1 = !(is_zero(cpi.op(0)-prop.op(0)) && is_zero(cpi.op(1)-prop.op(1)));
439 bool m2 = !(is_zero(cpi.op(0)-prop.op(1)) && is_zero(cpi.op(1)-prop.op(0)));
440 if(m1 && m2) continue; // other propagators
441 if(is_zero(li.op(0)-li.op(1))) return true; // a loop found
442 last = i;
443 lp = li;
444 tls[last] = 0;
445 break;
446 }
447 if(is_zero(lp)) break;
448
449 for(int i=0; i<ntls; i++) {
450 if(is_zero(tls.op(i))) continue;
451 if(is_zero(tls.op(i).op(0)-lp.op(0))) tls[i][0] = lp.op(1);
452 if(is_zero(tls.op(i).op(1)-lp.op(0))) tls[i][1] = lp.op(1);
453 }
454 }
455 return false;
456 }
457
464 ex QuarkPropagator(const ex & e, const ex & m) {
465 auto fi1 = e.op(0).op(1);
466 auto fi2 = e.op(1).op(1);
467 auto mom = e.op(2);
468 return I * SP(TI(fi1),TI(fi2)) * GMat(GAS(mom)+GAS(1)*m, DI(fi1),DI(fi2)) / (SP(mom)-m*m);
469 }
470
477 ex LeptonPropagator(const ex & e, const ex & m) {
478 auto fi1 = e.op(0).op(1);
479 auto fi2 = e.op(1).op(1);
480 auto mom = e.op(2);
481 return I * GMat(GAS(mom)+GAS(1)*m, DI(fi1),DI(fi2)) / (SP(mom)-m*m);
482 }
483
490 ex GluonPropagator(const ex & e, const ex & xi) {
491 auto fi1 = e.op(0).op(1);
492 auto fi2 = e.op(1).op(1);
493 auto mom = e.op(2);
494 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) );
495 }
496
503 ex GluonGhostPropagator(const ex & e, const ex & xi, const ex & eta_G) {
504 auto fi1 = e.op(0).op(1);
505 auto fi2 = e.op(1).op(1);
506 auto mom = e.op(2);
507 return I * eta_G * SP(CI(fi1),CI(fi2)) / SP(mom);
508 }
509
517 ex AZWPropagator(const ex & e, const ex & m, const ex & xi) {
518 auto fi1 = e.op(0).op(1);
519 auto fi2 = e.op(1).op(1);
520 auto mom = e.op(2);
521 return -I/(SP(mom)-m*m) * (SP(LI(fi1),LI(fi2))-(1-xi)*SP(LI(fi1))*SP(LI(fi2))/(SP(mom)-xi*m*m));
522 }
523
531 ex AZWGhostPropagator(const ex & e, const ex & m, const ex & xi, const ex & eta_G) {
532 auto fi1 = e.op(0).op(1);
533 auto fi2 = e.op(1).op(1);
534 auto mom = e.op(2);
535 return eta_G * I / (SP(mom)-xi*m*m);
536 }
537
545 ex ScalarPropagator(const ex & e, const ex & m, const ex & xi) {
546 auto fi1 = e.op(0).op(1);
547 auto fi2 = e.op(1).op(1);
548 auto mom = e.op(2);
549 return I / (SP(mom)-xi*m*m);
550 }
551
558 ex QuarkGluonVertex(const ex & e, const ex & eta_s) {
559 auto fi1 = e.op(0).op(1);
560 auto fi2 = e.op(1).op(1);
561 auto fi3 = e.op(2).op(1);
562 return -I*eta_s*gs*GMat(GAS(LI(fi3)),DI(fi1),DI(fi2))*SUNT(CI(fi3),TI(fi1),TI(fi2));
563 }
564
571 ex Gluon3Vertex(const ex & e, const ex & eta_s) {
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 auto mom1 = e.op(0).op(2);
576 auto mom2 = e.op(1).op(2);
577 auto mom3 = e.op(2).op(2);
578 return -eta_s*gs*SUNF(CI(fi1),CI(fi2),CI(fi3))*(
579 SP(mom1-mom2,LI(fi3))*SP(LI(fi1),LI(fi2)) +
580 SP(mom2-mom3,LI(fi1))*SP(LI(fi2),LI(fi3)) +
581 SP(mom3-mom1,LI(fi2))*SP(LI(fi3),LI(fi1))
582 );
583 }
584
591 ex Gluon4Vertex(const ex & e, const ex & eta_s) {
592 auto fi1 = e.op(0).op(1);
593 auto fi2 = e.op(1).op(1);
594 auto fi3 = e.op(2).op(1);
595 auto fi4 = e.op(3).op(1);
596 return -I*gs*gs*(
597 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))) +
598 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))) +
599 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)))
600 );
601 }
602
610 ex GhostGluonVertex(const ex & e, const ex & eta_s, const ex & eta_G) {
611 auto fi1 = e.op(0).op(1);
612 auto fi2 = e.op(1).op(1);
613 auto fi3 = e.op(2).op(1);
614 auto mom1 = e.op(0).op(2);
615 return eta_s*eta_G*gs*SUNF(CI(fi1),CI(fi2),CI(fi3))*SP(mom1,LI(fi3));
616 }
617
625 ex QuarkAVertex(const ex & e, const ex & eq, const ex & eta_e) {
626 auto fi1 = e.op(0).op(1);
627 auto fi2 = e.op(1).op(1);
628 auto fi3 = e.op(2).op(1);
629 return -I*eta_e*eq*GMat(GAS(LI(fi3)),DI(fi1),DI(fi2))*SP(TI(fi1),TI(fi2));
630 }
631
639 ex LeptonAVertex(const ex & e, const ex & eq, const ex & eta_e) {
640 auto fi1 = e.op(0).op(1);
641 auto fi2 = e.op(1).op(1);
642 auto fi3 = e.op(2).op(1);
643 return -I*eta_e*eq*GMat(GAS(LI(fi3)),DI(fi1),DI(fi2));
644 }
645
652 ex IndexL2R(ex e, bool all) {
653 static MapFunction map([all](const ex &e, MapFunction &self)->ex {
654 if(!Index::has(e)) return e;
655 else if(is_a<Index>(e)) {
656 auto idx = ex_to<Index>(e);
657 auto nstr = idx.name.get_name();
658 if(!all && nstr.rfind("lim",0)==0) return e;
659 else if(!all && nstr.rfind("dim",0)==0) return e;
660 else if(!all && nstr.rfind("cim",0)==0) return e;
661 else if(!all && nstr.rfind("tim",0)==0) return e;
662 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.dim);
663 else return e;
664 }
665 else return e.map(self);
666 });
667 return map(e);
668 }
669
670 ex IndexCC(ex e, bool all) {
671 return IndexL2R(e,all);
672 }
673
680 ex GluonSumL(int qi, bool color) {
681 if(color) return -SP(LI(qi), RLI(qi)) * SP(CI(qi), RCI(qi));
682 else return -SP(LI(qi), RLI(qi));
683 }
684
691 ex GluonSumR(int qi, bool color) {
692 if(color) return -SP(RLI(qi),LI(qi)) * SP(RCI(qi),CI(qi));
693 else return -SP(RLI(qi),LI(qi));
694 }
695
704 ex QuarkSumL(int qi, ex p, ex m, bool color) {
705 if(color) return GMat(GAS(p)+m*GAS(1), RDI(qi), DI(qi)) * SP(RTI(qi), TI(qi));
706 else return GMat(GAS(p)+m*GAS(1), RDI(qi), DI(qi));
707 }
708
717 ex QuarkSumR(int qi, ex p, ex m, bool color) {
718 if(color) return GMat(GAS(p)+m*GAS(1), DI(qi), RDI(qi)) * SP(TI(qi),RTI(qi));
719 else return GMat(GAS(p)+m*GAS(1), DI(qi), RDI(qi));
720 }
721
730 ex AntiQuarkSumL(int qi, ex p, ex m, bool color) {
731 if(color) return GMat(GAS(p)-m*GAS(1), DI(qi), RDI(qi)) * SP(TI(qi), RTI(qi));
732 else return GMat(GAS(p)-m*GAS(1), DI(qi), RDI(qi));
733 }
734
743 ex AntiQuarkSumR(int qi, ex p, ex m, bool color) {
744 if(color) return GMat(GAS(p)-m*GAS(1), RDI(qi), DI(qi)) * SP(RTI(qi), TI(qi));
745 else return GMat(GAS(p)-m*GAS(1), RDI(qi), DI(qi));
746 }
747
753 ex GhostSumL(int qi) {
754 return SP(CI(qi), RCI(qi));
755 }
756
762 ex GhostSumR(int qi) {
763 return SP(RCI(qi),CI(qi));
764 }
765
771 ex AntiGhostSumL(int qi) {
772 return -SP(CI(qi), RCI(qi));
773 }
774
780 ex AntiGhostSumR(int qi) {
781 return -SP(RCI(qi),CI(qi));
782 }
783
790 ex J1SumL(int qi, ex p) {
791 if(is_zero(p)) return -SP(LI(qi),RLI(qi));
792 else return -SP(LI(qi),RLI(qi)) + SP(p,LI(qi)) * SP(p,RLI(qi)) / SP(p);
793 }
794
795 ex J1SumTr(int si, const ex & p, const ex & n) {
796 ex mu = LI(si), nu = RLI(si);
797 return -SP(mu,nu) + (SP(p,mu)*SP(n,nu)+SP(p,nu)*SP(n,mu))/SP(p,n) - SP(p)*SP(n,mu)*SP(n,nu)/pow(SP(n,p),2);
798 }
799
806 ex J1SumR(int qi, ex p) {
807 if(is_zero(p)) return -SP(RLI(qi),LI(qi));
808 else return -SP(RLI(qi),LI(qi)) + SP(p,RLI(qi)) * SP(p,LI(qi)) / SP(p);
809 }
810
811 // https://arxiv.org/abs/1209.6213v2
812 lst Models::FeynRulesSM(const lst & amps, const ex & xi) {
813 lst ret;
814 for(auto item : amps) ret.append(FeynRulesSM(item,xi));
815 return ret;
816 }
817 ex Models::FeynRulesSM(const ex & amp, const ex & xi) {
818 if(is_a<lst>(amp)) return FeynRulesSM(ex_to<lst>(amp), xi);
819
820 static Symbol CW("CW"); // cos(theta)
821 static Symbol SW("SW"); // sin(theta)
822 //static Symbol C2W("C2W"); // cos(2theta)
823 static ex CW2 = CW*CW;
824 static ex SW2 = SW*SW;
825 static ex C2W = CW2-SW2;
826
827 static Symbol U("U"); // U-quark
828 static Symbol Ubar("Ubar"); // anti U-quark
829 static Symbol D("D"); // D-quark
830 static Symbol Dbar("Dbar"); // anti D-quark
831 static Symbol C("C"); // C-quark
832 static Symbol Cbar("Cbar"); // anti C-quark
833 static Symbol S("S"); // S-quark
834 static Symbol Sbar("Sbar"); // anti S-quark
835 static Symbol T("T"); // T-quark
836 static Symbol Tbar("Tbar"); // anti T-quark
837 static Symbol B("B"); // B-quark
838 static Symbol Bbar("Bbar"); // anti B-quark
839
840 static Symbol g("g"); // gluon
841 static Symbol gh("gh"); // gluon ghost
842 static Symbol ghbar("ghbar"); // anti gluon ghost
843
844 static Symbol A("A"); // photon
845 static Symbol Wm("Wm"); // W-
846 static Symbol Wp("Wp"); // W+
847 static Symbol Z("Z"); // Z
848
849 static Symbol ghA("ghA"); // photon ghost
850 static Symbol ghAbar("ghAbar"); // anti photon ghost
851 static Symbol ghWm("ghWm"); // W- ghost
852 static Symbol ghWmbar("ghWmbar"); // anti W- ghost
853 static Symbol ghWp("ghWp"); // W+ ghost
854 static Symbol ghWpbar("ghWpbar"); // anti W+ ghost
855 static Symbol ghZ("ghZ"); // Z ghost
856 static Symbol ghZbar("ghZbar"); // anti Z ghost
857
858 static Symbol em("em"); // e-
859 static Symbol ep("ep"); // e+
860 static Symbol ne("ne"); // e-neutrino
861 static Symbol nebar("nebar"); // anti e-neutrino
862 static Symbol mum("mum"); // mu-
863 static Symbol mup("mup"); // mu+
864 static Symbol nmu("nmu"); // mu-neutrino
865 static Symbol nmubar("nmubar"); // anti mu-neutrino
866 static Symbol taum("taum"); // tau-
867 static Symbol taup("taup"); // tau+
868 static Symbol ntau("ntau"); // tau-neutrino
869 static Symbol ntaubar("ntaubar"); // anti tau-neutrino
870
871 static Symbol chi("chi"); // Z goldstone
872 static Symbol phim("phim"); // W- goldstone
873 static Symbol phip("phip"); // W+ goldstone
874 static Symbol H("H"); // higgs
875
876 static auto M = [](const string & si)->ex {
877 return Symbol("M"+si);
878 };
879 static ex MW = M("W");
880 static ex MZ = M("Z");
881 static ex MH = M("H");
882
883 static Symbol EL("e"); // charge e
884
885 static ex eta_s = -1;
886 // Peskin's book convention
887 static ex eta = -1;
888 static ex eta_prime = -1;
889 static ex eta_Z = 1;
890 static ex eta_theta = 1;
891 static ex eta_Y = 1;
892 static ex eta_e = -1;
893 static ex eta_G = 1;
894
895 static ex GEW = eta_e*eta*eta_theta * EL/SW;
896 static ex EL2 = EL*EL;
897 static ex MW2 = MW*MW;
898 static ex MZ2 = MZ*MZ;
899 static ex MH2 = MH*MH;
900 static ex GEW2 = GEW*GEW;
901 static ex sqrt2 = sqrt(ex(2));
902
903 map<string,ex> T3;
904 T3["U"] = T3["C"] = T3["T"] = ex(1)/2;
905 T3["D"] = T3["S"] = T3["B"] = -ex(1)/2;
906 T3["ne"] = T3["nmu"] = T3["ntau"] = ex(1)/2;
907 T3["em"] = T3["mum"] = T3["taum"] = -ex(1)/2;
908 T3["ep"] = T3["mup"] = T3["taup"] = -ex(1)/2; // for sure
909
910 map<string,ex> Q;
911 Q["U"] = Q["C"] = Q["T"] = ex(2)/3;
912 Q["D"] = Q["S"] = Q["B"] = -ex(1)/3;
913 Q["ne"] = Q["nmu"] = Q["ntau"] = 0;
914 Q["em"] = Q["mum"] = Q["taum"] = -1;
915 Q["ep"] = Q["mup"] = Q["taup"] = -1; // for sure
916
917 auto is_qbar_q = [&](const ex pi1, const ex pi2)->bool {
918 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);
919 };
920 auto is_lbar_l = [&](const ex pi1, const ex pi2)->bool {
921 return (pi1==ep && pi2==em) || (pi1==mup && pi2==mum) || (pi1==taup && pi2==taum);
922 };
923 auto is_lbar_n = [&](const ex pi1, const ex pi2)->bool {
924 return (pi1==ep && pi2==ne) || (pi1==mup && pi2==nmu) || (pi1==taup && pi2==ntau);
925 };
926 auto is_nbar_n = [&](const ex pi1, const ex pi2)->bool {
927 return (pi1==nebar && pi2==ne) || (pi1==nmubar && pi2==nmu) || (pi1==ntaubar && pi2==ntau);
928 };
929 auto is_nbar_l = [&](const ex pi1, const ex pi2)->bool {
930 return (pi1==nebar && pi2==em) || (pi1==nmubar && pi2==mum) || (pi1==ntaubar && pi2==taum);
931 };
932 auto is_uqbar = [&](const ex pi)->bool {
933 return (pi==Ubar) || (pi==Cbar) || (pi==Tbar);
934 };
935 auto is_uq = [&](const ex pi)->bool {
936 return (pi==U) || (pi==C) || (pi==T);
937 };
938 auto is_dqbar = [&](const ex pi)->bool {
939 return (pi==Dbar) || (pi==Sbar) || (pi==Bbar);
940 };
941 auto is_dq = [&](const ex pi)->bool {
942 return (pi==D) || (pi==S) || (pi==B);
943 };
944 auto is_uqbar_dq = [&](const ex pi1, const ex pi2)->bool {
945 return is_uqbar(pi1) && is_dq(pi2);
946 };
947 auto is_dqbar_uq = [&](const ex pi1, const ex pi2)->bool {
948 return is_dqbar(pi1) && is_uq(pi2);
949 };
950 auto CKM = [&](const string & si1, const string & si2)->ex {
951 return Symbol("V"+si1+si2);
952 };
953
954 auto gfV = [&](const string & si)->ex {
955 return T3[si]/2-Q[si]*SW2;
956 };
957 auto gfA = [&](const string & si)->ex {
958 return T3[si]/2;
959 };
960
961 auto fr = MapFunction([&](const ex &e, MapFunction &self)->ex {
962 if(isFunction(e,"OutField") || isFunction(e,"InField")) return 1;
963 else if(isFunction(e, "Propagator")) {
964 auto pi = e.op(0).op(0);
965 auto si = ex2str(pi);
966 if(pi==U || pi==D || pi==C || pi==S || pi==T || pi==B) {
967 return QuarkPropagator(e, M(si));
968 } else if(pi==g) {
969 return GluonPropagator(e, xi);
970 } else if(pi==gh) {
971 return GluonGhostPropagator(e, xi, eta_G);
972 } else if(pi==A) {
973 return AZWPropagator(e, 0, xi);
974 } else if(pi==Z) {
975 return AZWPropagator(e, MZ, xi);
976 } else if(pi==Wm) {
977 return AZWPropagator(e, MW, xi);
978 } else if(pi==ghA) {
979 return AZWGhostPropagator(e, 0, xi, eta_G);
980 } else if(pi==ghZ) {
981 return AZWGhostPropagator(e, MZ, xi, eta_G);
982 } else if(pi==ghWm || pi==ghWp) {
983 return AZWGhostPropagator(e, MW, xi, eta_G);
984 } else if(pi==em || pi==mum || pi==taum) {
985 return LeptonPropagator(e, M(si));
986 } else if(pi==ne || pi==nmu || pi==ntau) {
987 return LeptonPropagator(e, 0);
988 } else if(pi==H) {
989 return ScalarPropagator(e, MH, 1);
990 } else if(pi==chi) {
991 return ScalarPropagator(e, MZ, xi);
992 } else if(pi==phim) {
993 return ScalarPropagator(e, MW, xi);
994 } else {
995 cout << endl << e << endl;
996 throw Error("Propagator Not defined!");
997 }
998 } else if(isFunction(e, "Vertex")) {
999 auto pi1 = e.op(0).op(0);
1000 auto pi2 = e.op(1).op(0);
1001 auto pi3 = e.op(2).op(0);
1002 auto fi1 = e.op(0).op(1);
1003 auto fi2 = e.op(1).op(1);
1004 auto fi3 = e.op(2).op(1);
1005 auto mom1 = e.op(0).op(2);
1006 auto mom2 = e.op(1).op(2);
1007 auto mom3 = e.op(2).op(2);
1008 int nn = e.nops();
1009
1010 auto si1 = ex2str(pi1);
1011 auto si2 = ex2str(pi2);
1012 string_replace_all(si1, "bar", "");
1013 string_replace_all(si2, "bar", "");
1014 auto si = si2;
1015
1016 if(nn==3) {
1017 if(pi1==ghbar && pi2==gh && pi3==g) {
1018 // ghbar-gh-g
1019 return GhostGluonVertex(e, eta_s, eta_G);
1020 } else if(pi1==g && pi2==g && pi3==g) {
1021 // g^3
1022 return Gluon3Vertex(e, eta_s);
1023 } else if((pi3==g) && is_qbar_q(pi1, pi2)) {
1024 // qbar-q-g
1025 return QuarkGluonVertex(e, eta_s);
1026 } else if((pi3==A) && is_qbar_q(pi1, pi2)) {
1027 // qbar-q-A
1028 return QuarkAVertex(e, EL*Q[si], eta_e);
1029 } else if((pi3==A) && is_lbar_l(pi1, pi2)) {
1030 // lbar-l-A
1031 return LeptonAVertex(e, EL*Q[si], eta_e);
1032 } else if((pi3==Z) && is_qbar_q(pi1, pi2)) {
1033 // qbar-q-Z
1034 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));
1035 } else if((pi3==Z) && (is_lbar_l(pi1, pi2) || is_nbar_n(pi1, pi2))) {
1036 // lbar-l-Z & nbar-n-Z
1037 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)));
1038 } else if(pi1==Wp && pi2==Wm && pi3==A) {
1039 // Wp-Wm-A
1040 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)));
1041 } else if(pi1==Wp && pi2==Wm && pi3==Z) {
1042 // Wp-Wm-Z
1043 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)));
1044 } else if(is_uqbar(pi1) && is_dq(pi2) && pi3==Wp) {
1045 // uqbar-dq-Wp
1046 auto ckm = CKM(si1, si2);
1047 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));
1048 } else if(is_dqbar(pi1) && is_uq(pi2) && pi3==Wm) {
1049 // dqbar-uq-Wm
1050 auto ckm = CKM(si1, si2);
1051 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));
1052 } else if( (is_nbar_l(pi1, pi2) && pi3==Wp) || (is_lbar_n(pi1, pi2) && pi3==Wm) ) {
1053 // nbar-l-Wp & lbar-n-Wm
1054 return -I*eta*GEW/sqrt2*GMat(GAS(LI(fi3))-GAS(LI(fi3))*GAS(5),DI(fi1),DI(fi2))/2;
1055 } else if(is_qbar_q(pi1, pi2) && (pi3==H)) {
1056 // qbar-q-H
1057 return -I*GEW/2*M(si)/MW*GMat(GAS(1),DI(fi1),DI(fi2))*SP(TI(fi1),TI(fi2));
1058 } else if( (is_lbar_l(pi1, pi2) || is_nbar_n(pi1, pi2)) && (pi3==H)) {
1059 // lbar-l-H & nbar-n-H
1060 return -I*GEW/2*M(si)/MW*GMat(GAS(1),DI(fi1),DI(fi2));
1061 } else if(is_qbar_q(pi1, pi2) && (pi3==chi)) {
1062 // qbar-q-chi
1063 return -GEW*T3[si]*M(si)/MW*GMat(GAS(5),DI(fi1),DI(fi2))*SP(TI(fi1),TI(fi2));
1064 } else if( (is_lbar_l(pi1, pi2) || is_nbar_n(pi1, pi2)) && (pi3==chi)) {
1065 // lbar-l-chi & nbar-n-chi
1066 return -GEW*T3[si]*M(si)/MW*GMat(GAS(5),DI(fi1),DI(fi2));
1067 } else if(is_uqbar_dq(pi1, pi2) && pi3==phip) {
1068 // uqbar-dq-phip
1069 auto ckm = CKM(si1,si2);
1070 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));
1071 } else if(is_dqbar_uq(pi1, pi2) && pi3==phim) {
1072 // dqbar-uq-phim
1073 auto ckm = CKM(si1,si2);
1074 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));
1075 } else if(is_nbar_l(pi1, pi2) && pi3==phip) {
1076 // nbar-l-phip
1077 return -I*GEW/sqrt2*M(si2)/MW*GMat(GAS(1)+GAS(5),DI(fi1),DI(fi2))/2; // Mnl = 0
1078 } else if(is_lbar_n(pi1, pi2) && pi3==phim) {
1079 // lbar-n-phim
1080 return -I*GEW/sqrt2*M(si1)/MW*GMat(GAS(1)-GAS(5),DI(fi1),DI(fi2))/2; // Mnl = 0
1081 } else if(pi1==A && pi2==phip && pi3==phim) {
1082 // A-phip-phim
1083 return -I*eta_e*EL*SP(mom2-mom3,LI(fi1));
1084 } else if(pi1==Z && pi2==phip && pi3==phim) {
1085 // Z-phip-phim
1086 return -I*eta*eta_Z*GEW*C2W/(2*CW)*SP(mom2-mom3,LI(fi1));
1087 } else if(pi1==Wp && pi2==phim && pi3==H) {
1088 // Wp-phim-H
1089 return I/2*eta*GEW*SP(mom2-mom3,LI(fi1));
1090 } else if(pi1==Wm && pi2==phip && pi3==H) {
1091 // Wm-phip-H
1092 return -I/2*eta*GEW*SP(mom2-mom3,LI(fi1));
1093 } else if( ((pi1==Wp && pi2==phim) || (pi1==Wm && pi2==phip)) && pi3==chi) {
1094 // Wp-phim-chi & Wm-phip-chi
1095 return -eta*GEW/2*SP(mom2-mom3,LI(fi1));
1096 } else if(pi1==Z && pi2==chi && pi3==H) {
1097 // Z-chi-H
1098 return -eta*eta_Z*GEW/(2*CW)*SP(mom2-mom3,LI(fi1));
1099 } else if( ((pi1==phim && pi2==Wp) || (pi1==phip && pi2==Wm)) && pi3==A) {
1100 // phim-Wp-A & phip-Wm-A
1101 return I*eta_e*eta*EL*MW*SP(LI(fi2),LI(fi3));
1102 } else if( ((pi1==phim && pi2==Wp) || (pi1==phip && pi2==Wm)) && pi3==Z) {
1103 // phim-Wp-Z & phip-Wm-Z
1104 return -I*eta_Z*GEW*MZ*SW2*SP(LI(fi2),LI(fi3));
1105 } else if(pi1==H && pi2==Wp && pi3==Wm) {
1106 // H-Wp-Wm
1107 return I*GEW*MW*SP(LI(fi2),LI(fi3));
1108 } else if(pi1==H && pi2==Z && pi3==Z) {
1109 // H-Z-Z
1110 return I*GEW/CW*MZ*SP(LI(fi2),LI(fi3));
1111 } else if(pi1==H && pi2==phim && pi3==phip) {
1112 // H-phim-phip
1113 return -I/2*GEW*MH2/MW;
1114 } else if(pi1==H && pi2==H && pi3==H) {
1115 // H-H-H
1116 return -3/ex(2)*I*GEW*MH2/MW;
1117 } else if(pi1==H && pi2==chi && pi3==chi) {
1118 // H-chi-chi
1119 return -I/2*GEW*MH2/MW;
1120 } else if(pi1==ghWpbar && pi2==ghWp && pi3==A) {
1121 // ghWpbar-ghWp-A
1122 return I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1123 } else if(pi1==ghWmbar && pi2==ghWm && pi3==A) {
1124 // ghWmbar-ghWm-A
1125 return -I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1126 } else if(pi1==ghWpbar && pi2==ghWp && pi3==Z) {
1127 // ghWpbar-ghWp-Z
1128 return I*eta_G*eta*eta_Z*GEW*CW*SP(mom1,LI(fi3));
1129 } else if(pi1==ghWmbar && pi2==ghWm && pi3==Z) {
1130 // ghWmbar-ghWm-Z
1131 return -I*eta_G*eta*eta_Z*GEW*CW*SP(mom1,LI(fi3));
1132 } else if(pi1==ghWpbar && pi2==ghZ && pi3==Wp) {
1133 // ghWpbar-ghZ-Wp
1134 return -I*eta_G*eta*eta_Z*GEW*CW*SP(mom1,LI(fi3));
1135 } else if(pi1==ghWmbar && pi2==ghZ && pi3==Wm) {
1136 // ghWmbar-ghZ-Wm
1137 return I*eta_G*eta*eta_Z*GEW*CW*SP(mom1,LI(fi3));
1138 } else if(pi1==ghWpbar && pi2==ghA && pi3==Wp) {
1139 // ghWpbar-ghA-Wp
1140 return -I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1141 } else if(pi1==ghWmbar && pi2==ghA && pi3==Wm) {
1142 // ghWmbar-ghA-Wm
1143 return I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1144 } else if(pi1==ghZbar && pi2==ghWp && pi3==Wm) {
1145 // ghZbar-ghWp-Wm
1146 return -I*eta_G*eta*GEW*CW*SP(mom1,LI(fi3));
1147 } else if(pi1==ghZbar && pi2==ghWm && pi3==Wp) {
1148 // ghZbar-ghWm-Wp
1149 return I*eta_G*eta*GEW*CW*SP(mom1,LI(fi3));
1150 } else if(pi1==ghAbar && pi2==ghWp && pi3==Wm) {
1151 // ghAbar-ghWp-Wm
1152 return -I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1153 } else if(pi1==ghAbar && pi2==ghWm && pi3==Wp) {
1154 // ghAbar-ghWm-Wp
1155 return I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1156 } else if(pi1==ghWpbar && pi2==ghWp && pi3==chi) {
1157 // ghWpbar-ghWp-chi
1158 return eta_G*GEW/2*xi*MW;
1159 } else if(pi1==ghWmbar && pi2==ghWm && pi3==chi) {
1160 // ghWmbar-ghWm-chi
1161 return -eta_G*GEW/2*xi*MW;
1162 } else if( ((pi1==ghWpbar && pi2==ghWp) ||(pi1==ghWmbar && pi2==ghWm)) && pi3==H) {
1163 // ghWpbar-ghWp-H & ghWmbar-ghWm-H
1164 return -I/2*eta_G*GEW*xi*MW;
1165 } else if(pi1==ghZbar && pi2==ghZ && pi3==H) {
1166 // ghZbar-ghZ-H
1167 return -eta_G*I*GEW/(2*CW)*xi*MZ;
1168 } else if( pi1==ghZbar && ((pi2==ghWp && pi3==phim) || (pi2==ghWm && pi3==phip)) ) {
1169 // ghZbar-ghWp-phim & ghZbar-ghWm-phip
1170 return I/2*eta_G*eta_Z*GEW*xi*MZ;
1171 } else if( ((pi1==ghWpbar && pi3==phip) || (pi1==ghWmbar && pi3==phim)) && pi2==ghZ) {
1172 // ghWpbar-ghZ-phip & ghWmbar-ghZ-phim
1173 return -I*eta_G*eta_Z*GEW*C2W/(2*CW)*xi*MW;
1174 } else if( ((pi1==ghWpbar && pi3==phip) || (pi1==ghWmbar && pi3==phim)) && pi2==ghA) {
1175 // ghWpbar-ghA-phip & ghWmbar-ghA-phim
1176 return -I*eta_G*eta_e*eta*EL*xi*MW;
1177 } else {
1178 cout << endl << e << endl;
1179 throw Error("Vertex Not defined!");
1180 }
1181 } else if(nn==4) {
1182 auto pi4 = e.op(3).op(0);
1183 auto fi4 = e.op(3).op(1);
1184 auto mom4 = e.op(3).op(2);
1185
1186 if(pi1==g && pi2==g && pi3==g && pi4==g) {
1187 // g^4
1188 return Gluon4Vertex(e, eta_s);
1189 } else if(pi1==Wp && pi2==Wm && pi3==A && pi4==A) {
1190 // Wp-Wm-A-A
1191 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)));
1192 } else if(pi1==Wp && pi2==Wm && pi3==Z && pi4==Z) {
1193 // Wp-Wm-Z-Z
1194 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)));
1195 } else if(pi1==Wp && pi2==Wm && pi3==A && pi4==Z) {
1196 // Wp-Wm-A-Z
1197 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)));
1198 } else if(pi1==Wp && pi2==Wp && pi3==Wm && pi4==Wm) {
1199 // Wp-Wp-Wm-Wm
1200 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)));
1201 } else if(pi1==Wp && pi2==Wm && pi3==H && pi4==H) {
1202 // Wp-Wm-H-H
1203 return I/2*GEW2*SP(LI(fi1),LI(fi2));
1204 } else if(pi1==Wp && pi2==Wm && pi3==chi && pi4==chi) {
1205 // Wp-Wm-chi-chi
1206 return I/2*GEW2*SP(LI(fi1),LI(fi2));
1207 } else if(pi1==Z && pi2==Z && pi3==H && pi4==H) {
1208 // Z-Z-H-H
1209 return I/2*GEW2/CW2*SP(LI(fi1),LI(fi2));
1210 } else if(pi1==Z && pi2==Z && pi3==chi && pi4==chi) {
1211 // Z-Z-chi-chi
1212 return I/2*GEW2/CW2*SP(LI(fi1),LI(fi2));
1213 } else if(pi1==A && pi2==A && pi3==phip && pi4==phim) {
1214 // A-A-phip-phim
1215 return 2*I*EL2*SP(LI(fi1),LI(fi2));
1216 } else if(pi1==Z && pi2==Z && pi3==phip && pi4==phim) {
1217 // Z-Z-phip-phim
1218 return I/2*pow(GEW*C2W/CW,2)*SP(LI(fi1),LI(fi2));
1219 } else if(pi1==Wp && pi2==Wm && pi3==phip && pi4==phim) {
1220 // Wp-Wm-phip-phim
1221 return I/2*GEW2*SP(LI(fi1),LI(fi2));
1222 } else if( ((pi1==Wp && pi3==phim) || (pi1==Wm && pi3==phip)) && pi2==Z && pi4==H) {
1223 // Wp-Z-phim-H & Wm-Z-phip-H
1224 return -I*eta_Z*GEW2*SW2/(2*CW)*SP(LI(fi1),LI(fi2));
1225 } else if(pi1==Wm && pi2==Z && pi3==phip && pi4==chi) {
1226 // Wm-Z-phip-chi
1227 return -eta_Z*GEW2*SW2/(2*CW)*SP(LI(fi1),LI(fi2));
1228 } else if(pi1==Wp && pi2==Z && pi3==phim && pi4==chi) {
1229 // Wp-Z-phim-chi
1230 return eta_Z*GEW2*SW2/(2*CW)*SP(LI(fi1),LI(fi2));
1231 } else if( ((pi1==Wm && pi3==phip) || (pi1==Wp && pi3==phim)) && pi2==A && pi4==H) {
1232 // Wm-A-phip-H & Wp-A-phim-H
1233 return I/2*eta_e*eta*EL*GEW*SP(LI(fi1),LI(fi2));
1234 } else if(pi1==Wp && pi2==A && pi3==phim && pi4==chi) {
1235 // Wp-A-phim-chi
1236 return -1/ex(2)*eta_e*eta*EL*GEW*SP(LI(fi1),LI(fi2));
1237 } else if(pi1==Wm && pi2==A && pi3==phip && pi4==chi) {
1238 // Wm-A-phip-chi
1239 return 1/ex(2)*eta_e*eta*EL*GEW*SP(LI(fi1),LI(fi2));
1240 } else if(pi1==Z && pi2==A && pi3==phip && pi4==phim) {
1241 // Z-A-phip-phim
1242 return I*eta_e*eta*eta_Z*EL*GEW*C2W/CW*SP(LI(fi1),LI(fi2));
1243 } else if(pi1==phim && pi2==phip && pi3==phim && pi4==phip) {
1244 // phim-phip-phim-phip
1245 return -I/2*GEW2*MH2/MW2;
1246 } else if(pi1==H && pi2==H && pi3==phim && pi4==phip) {
1247 // H-H-phim-phip
1248 return -I/4*GEW2*MH2/MW2;
1249 } else if(pi1==chi && pi2==chi && pi3==phim && pi4==phip) {
1250 // chi-chi-phim-phip
1251 return -I/4*GEW2*MH2/MW2;
1252 } else if(pi1==H && pi2==H && pi3==H && pi4==H) {
1253 // H-H-H-H
1254 return -3/ex(4)*I*GEW2*MH2/MW2;
1255 } else if(pi1==H && pi2==H && pi3==chi && pi4==chi) {
1256 // H-H-chi-chi
1257 return -I/4*GEW2*MH2/MW2;
1258 } else if(pi1==chi && pi2==chi && pi3==chi && pi4==chi) {
1259 // chi-chi-chi-chi
1260 return -3/ex(4)*I*GEW2*MH2/MW2;
1261 } else {
1262 cout << endl << e << endl;
1263 throw Error("Vertex Not defined!");
1264 }
1265 } else {
1266 cout << endl << e << endl;
1267 throw Error("Vertex Not defined!");
1268 }
1269 } else return e.map(self);
1270 return e;
1271 });
1272 return fr(amp);
1273 }
1274
1275
1276 lst Models::FeynRulesQCD(const lst & amps, const ex & xi) {
1277 lst ret;
1278 for(auto item : amps) ret.append(FeynRulesQCD(item,xi));
1279 return ret;
1280 }
1281 ex Models::FeynRulesQCD(const ex & amp, const ex & xi) {
1282 if(is_a<lst>(amp)) return FeynRulesQCD(ex_to<lst>(amp), xi);
1283
1284 static Symbol A("A"), Q("Q"), Qbar("Qbar"), q("q"), qbar("qbar"), g("g"), gh("gh"), ghbar("ghbar");
1285 static Symbol m("m"), eq("eq"), eQ("eQ");
1286
1287 auto fr = MapFunction([&](const ex &e, MapFunction &self)->ex {
1288 if(isFunction(e,"OutField") || isFunction(e,"InField")) return 1;
1289 else if(isFunction(e, "Propagator")) {
1290 if(e.op(0).op(0)==q) {
1291 return QuarkPropagator(e, 0);
1292 } else if(e.op(0).op(0)==Q) {
1293 return QuarkPropagator(e, m);
1294 } else if(e.op(0).op(0)==g) {
1295 return GluonPropagator(e, xi);
1296 } else if(e.op(0).op(0)==gh) {
1297 return GluonGhostPropagator(e, xi);
1298 } else {
1299 cout << "expr: " << e << endl;
1300 throw Error("Propagator Not Found!");
1301 }
1302 } else if(isFunction(e, "Vertex")) {
1303 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) ) {
1304 // qbar-q-g or Qbar-Q-g or g -> A
1305 if(e.op(2).op(0)==g) return QuarkGluonVertex(e);
1306 else {
1307 auto fi1 = e.op(0).op(1);
1308 auto fi2 = e.op(1).op(1);
1309 auto fi3 = e.op(2).op(1);
1310 if(e.op(1).op(0)==q) QuarkAVertex(e, eq);
1311 else if(e.op(1).op(0)==Q) QuarkAVertex(e, eQ);
1312 else throw Error("Vertex Error.");
1313 }
1314 } else if(e.nops()==3 && e.op(0).op(0)==ghbar && e.op(1).op(0)==gh) {
1315 // ghbar-gh-g
1316 return GhostGluonVertex(e);
1317 } else if(e.nops()==3 && e.op(0).op(0)==g && e.op(1).op(0)==g && e.op(2).op(0)==g) {
1318 // g^3
1319 return Gluon3Vertex(e);
1320 } 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) {
1321 // g^4
1322 return Gluon4Vertex(e);
1323 } else {
1324 cout << "expr: " << e << endl;
1325 throw Error("Vertex Not Found!");
1326 }
1327 } else return e.map(self);
1328 return e;
1329 });
1330 return fr(amp);
1331 }
1332
1333}
1334
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:730
HepLib::QGRAF::HasLoop
bool HasLoop(ex amp, lst prop)
check a amplitude has a loop w.r.t. propapagtor
Definition QGRAF.cpp:427
HepLib::QGRAF::AntiGhostSumL
ex AntiGhostSumL(int qi)
polarization sum for anti-ghost, note that we will add extra - here
Definition QGRAF.cpp:771
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:652
HepLib::QGRAF::GluonPropagator
ex GluonPropagator(const ex &e, const ex &xi)
Propagator for gluon.
Definition QGRAF.cpp:490
HepLib::QGRAF::LeptonAVertex
ex LeptonAVertex(const ex &e, const ex &eq, const ex &eta_e)
l-lbar-A vertex
Definition QGRAF.cpp:639
HepLib::QGRAF::Gluon3Vertex
ex Gluon3Vertex(const ex &e, const ex &eta_s)
g-g-g vertex
Definition QGRAF.cpp:571
HepLib::QGRAF::GhostSumL
ex GhostSumL(int qi)
polarization sum for ghost
Definition QGRAF.cpp:753
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:780
HepLib::QGRAF::J1SumTr
ex J1SumTr(int si, const ex &p, const ex &n)
Definition QGRAF.cpp:795
HepLib::QGRAF::QuarkGluonVertex
ex QuarkGluonVertex(const ex &e, const ex &eta_s)
q-qbar-g vertex
Definition QGRAF.cpp:558
HepLib::QGRAF::ScalarPropagator
ex ScalarPropagator(const ex &e, const ex &m, const ex &xi)
Propagator for scalar.
Definition QGRAF.cpp:545
HepLib::QGRAF::J1SumL
ex J1SumL(int qi, ex p)
polarization sum for total angular momentum
Definition QGRAF.cpp:790
HepLib::QGRAF::LineTeX
map< ex, string, ex_is_less > LineTeX
Definition Init.cpp:185
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:477
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:339
HepLib::QGRAF::AntiQuarkSumR
ex AntiQuarkSumR(int qi, ex p, ex m, bool color)
polarization sum for anti-quark
Definition QGRAF.cpp:743
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:806
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:531
HepLib::QGRAF::VerTeX
map< ex, string, ex_is_less > VerTeX
Definition Init.cpp:186
HepLib::QGRAF::Gluon4Vertex
ex Gluon4Vertex(const ex &e, const ex &eta_s)
g-g-g-g vertex
Definition QGRAF.cpp:591
HepLib::QGRAF::IndexCC
ex IndexCC(ex e, bool all)
Definition QGRAF.cpp:670
HepLib::QGRAF::GhostGluonVertex
ex GhostGluonVertex(const ex &e, const ex &eta_s, const ex &eta_G)
ghbar-gh-g vertex
Definition QGRAF.cpp:610
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:464
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:717
HepLib::QGRAF::QuarkSumL
ex QuarkSumL(int qi, ex p, ex m, bool color)
polarization sum for quark
Definition QGRAF.cpp:704
HepLib::QGRAF::GhostSumR
ex GhostSumR(int qi)
polarization sum for ghost
Definition QGRAF.cpp:762
HepLib::QGRAF::LI
Index LI(ex fn)
Definition QGRAF.cpp:46
HepLib::QGRAF::InOutTeX
map< ex, string, ex_is_less > InOutTeX
Definition Init.cpp:187
HepLib::QGRAF::AZWPropagator
ex AZWPropagator(const ex &e, const ex &m, const ex &xi)
Propagator for A/Z/W boson.
Definition QGRAF.cpp:517
HepLib::QGRAF::QuarkAVertex
ex QuarkAVertex(const ex &e, const ex &eq, const ex &eta_e)
qbar-q-A vertex
Definition QGRAF.cpp:625
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:503
HepLib::QGRAF::FI
Index FI(ex fn)
Definition QGRAF.cpp:48
HepLib::QGRAF::DrawPDF
void DrawPDF(const lst &amps, string fn, int nr, bool single_page)
generate Feynman diagrams for the amplitudes, in PDF format
Definition QGRAF.cpp:169
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:123
HepLib::QGRAF::GluonSumR
ex GluonSumR(int qi, bool color)
polarization sum for gluon
Definition QGRAF.cpp:691
HepLib::QGRAF::GluonSumL
ex GluonSumL(int qi, bool color)
polarization sum for gluon
Definition QGRAF.cpp:680
HepLib::NA
const Symbol NA
HepLib::gs
const Symbol gs
HepLib::Combinations
void Combinations(int n, int m, std::function< void(const int *)> f)
Definition Functions.cpp:184
HepLib::NF
const Symbol NF
HepLib::ep
const Symbol ep
HepLib::GAS
ex GAS(const ex &expr, unsigned rl)
function similar to GAD/GSD in FeynClac
Definition DGamma.cpp:280
HepLib::file_exists
bool file_exists(string fn)
Definition BASIC.h:289
HepLib::Debug
bool Debug
Definition Init.cpp:144
HepLib::w
ex w
Definition Init.cpp:93
HepLib::vs
const Symbol vs
HepLib::d
const Symbol d
HepLib::isFunction
bool isFunction(const ex &e, string func_name)
Definition BASIC.h:658
HepLib::mu
const Symbol mu
HepLib::InstallPrefix
string InstallPrefix
Definition Init.cpp:164
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:1159
HepLib::SP
ex SP(const ex &a, bool use_map=false)
Definition Pair.cpp:166
Generated by doxygen 1.9.8