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) {
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 exmap bend_map;
187 std::map<ex,int,ex_is_less> vtex_map; // vertex option, only once
188 for(auto l : lines) {
189 lst ll = lst{l.op(0), l.op(1)};
190 bool isExt = (is_a<numeric>(ll.op(0)) && ll.op(0)<0) || (is_a<numeric>(ll.op(1)) && ll.op(1)<0);
191 ll.sort();
192 bend_map[ll] = bend_map[ll] + 1;
193
194 auto fidL = (is_a<numeric>(l.op(1)) ? l.op(1) : l.op(1).op(0));
195 out << "\"" << fidL << "\"";
196 if(fidL<0) {
197 out << "[particle=";
198 if(InOutTeX[fidL].length()>0) out << InOutTeX[fidL];
199 else out << fidL;
200 out << "]";
201 } else if(vtex_map[l.op(1).op(1)]==0) {
202 out << VerTeX[l.op(1).op(1)];
203 vtex_map[l.op(1).op(1)]=1;
204 }
205
206 out << " --[";
207 auto f = l.op(2).op(0);
208 if(LineTeX[f].length()>0) {
209 if(!isExt) out << LineTeX[f];
210 else {
211 auto cpos = LineTeX[f].find(", edge");
212 if(cpos>0) out << LineTeX[f].substr(0,cpos);
213 else out << LineTeX[f];
214 }
215 }
216 if(bend_map[ll]>2) out << ",half right";
217 else if(bend_map[ll]>1) out << ",half left";
218 if(is_zero(l.op(0)-l.op(1))) out << ",loop,distance=2cm";
219 out << "]";
220
221 auto fidR = (is_a<numeric>(l.op(0)) ? l.op(0) : l.op(0).op(0));
222 out << " \"" << fidR << "\"";
223 if(fidR<0) {
224 out << "[particle=";
225 if(InOutTeX[fidR].length()>0) out << InOutTeX[fidR];
226 else out<< fidR;
227 out << "]";
228 } else if(vtex_map[l.op(0).op(1)]==0) {
229 out << VerTeX[l.op(0).op(1)];
230 vtex_map[l.op(0).op(1)]=1;
231 }
232 out << ";" << endl;
233 }
234 out << "};" << endl;
235 out << "\\end{document}" << endl;
236 out.close();
237 auto rc = system(("cd "+tex_path+" && echo X | lualatex " + to_string(idx) + " 1>/dev/null").c_str());
238 return 0;
239 }, "TeX");
240
241 ofstream out(tex_path+"diagram.tex");
242 if(single_page) {
243 out << "\\let\\mypdfximage\\pdfximage" << endl;
244 out << "\\def\\pdfximage{\\immediate\\mypdfximage}" << endl;
245 out << "\\documentclass{standalone}" << endl;
246 } else {
247 out << "\\documentclass{article}" << endl;
248 }
249 out << "\\usepackage{graphicx}" << endl;
250 out << "\\usepackage{adjustbox}" << endl;
251 out << "\\usepackage{longtable}" << endl;
252 out << "\\usepackage{scalefnt}" << endl;
253 out << "\\begin{document}" << endl;
254 //out << "\\scalefont{" << .5/nr << "}" << endl;
255 if(single_page) {
256 out << "\\begin{adjustbox}{valign=T,width=\\textwidth}" << endl;
257 out << "\\begin{tabular}{|"; for(int i=0; i<nr; i++) out << "cc|"; out << "}" << endl;
258 } else {
259 out << "\\setlength\\LTleft{-2.5cm} \\setlength\\LTright{0pt}" << endl;
260 out << "\\begin{longtable}{|"; for(int i=0; i<nr; i++) out << "cc|"; out << "}" << endl;
261 }
262 out << "\\hline" << endl;
263 int total = amps.nops();
264 int namps = total;
265 if((total%nr)!=0) total = (total/nr+1)*nr;
266 for(int i=0; i<total; i++) {
267
268 //int limit = 300;
269 //if((i!=0) && (i+1!=total)) {
270 // out << "\\end{tabular}" << endl << endl;
271 // out << "\\begin{tabular}{|"; for(int i=0; i<nr; i++) out << "cc|"; out << "}" << endl;
272 // out << "\\hline" << endl;
273 //}
274
275 out << "{\\tiny " << i+1 << "}&" << endl;
276 if(i<namps) {
277 out << "\\includegraphics[keepaspectratio,";
278 out << "height=" << 1.0/nr << "\\textwidth,";
279 out << "width=" << 1.0/nr << "\\textwidth]";
280 out << "{"<<i<<".pdf}" << endl;
281 }
282 if((i+1)%nr==0) out << "\\\\ \\hline";
283 else out << "&";
284 }
285 if(single_page) {
286 out << "\\end{tabular}" << endl;
287 out << "\\end{adjustbox}" << endl;
288 } else {
289 out << "\\end{longtable}" << endl;
290 }
291 out << "\\end{document}" << endl;
292 out.close();
293 if(Debug) rc = system(("ulimit -n unlimited > /dev/null; cd "+tex_path+" && pdflatex diagram && mv diagram.pdf ../"+fn).c_str());
294 else rc = system(("ulimit -n unlimited > /dev/null; cd "+tex_path+" && echo X | pdflatex diagram 1>/dev/null && mv diagram.pdf ../"+fn).c_str());
295
296 if(!file_exists(fn)) {
297 cout << "DrawPDF failed, files not removed, please check the TeX files!" << endl;
298 cout << "latex command: pdflatex diagram.tex" << endl;
299 cout << "if there are many pdf files, one needs ulimint -n <large number>." << endl;
300 exit(1);
301 } else if(!Debug) {
302 rc = system(("rm -r "+tex_path).c_str());
303 }
304 }
305
313 vector<lst> ShrinkCut(ex amp, lst prop, int n) {
314 vector<lst> ret;
315 if(prop.nops()<1) throw Error("ShrinkCut: no cut provided!");
316 auto tls = TopoLines(amp);
317 vector<int> cls_vec;
318 for(int i=0; i<tls.nops(); i++) {
319 auto pi = tls.op(i).op(2);
320 if(pi.nops()<2) continue;
321 if(!is_a<lst>(prop.op(0))) { // prop : lst{ g, g }
322 if(is_zero(pi.op(0)-prop.op(0)) && is_zero(pi.op(1)-prop.op(1))) cls_vec.push_back(i);
323 } else {
324 for(auto iprop : prop) {
325 if(is_zero(pi.op(0)-iprop.op(0)) && is_zero(pi.op(1)-iprop.op(1))) cls_vec.push_back(i);
326 }
327 }
328 }
329 if(cls_vec.size()<n) return ret;
330
331 Combinations(cls_vec.size(), n, [n,&ret,cls_vec,tls](const int * is)->void {
332 int cls[n];
333 for(int i=0; i<n; i++) cls[i] = cls_vec[is[i]];
334
335 // cut each line into 2 half-lines, labeled with 0
336 auto tls2 = tls;
337 for(auto ci : cls) {
338 auto ol = tls2.op(ci);
339 tls2[ci] = lst{ol.op(0), 0, lst{ol.op(2).op(0)}, ol.op(3)};
340 tls2.append(lst{0, ol.op(1), lst{ol.op(2).op(1)}, ol.op(3)} );
341 }
342
343 // shrink internal lines
344 int last = 0;
345 int ntls2 = tls2.nops();
346 while(true) {
347 ex lp = 0;
348 for(int i=last; i<ntls2; i++) {
349 auto li = tls2.op(i);
350 if(is_zero(li) || li.op(2).nops()<2) continue;
351 last = i;
352 lp = li;
353 tls2[last] = 0;
354 break;
355 }
356 if(is_zero(lp)) break;
357
358 for(int i=0; i<ntls2; i++) {
359 if(is_zero(tls2.op(i))) continue;
360 if(is_zero(tls2.op(i).op(0)-lp.op(0))) tls2[i][0] = lp.op(1);
361 if(is_zero(tls2.op(i).op(1)-lp.op(0))) tls2[i][1] = lp.op(1);
362 }
363 }
364
365 // final connected parts
366 map<int, lst> con_map;
367 for(auto li : tls2) {
368 if(is_zero(li)) continue;
369 ex key, val;
370 ex fiL = li.op(0);
371 if(is_a<lst>(fiL)) fiL = fiL.op(0);
372 ex fiR = li.op(1);
373 if(is_a<lst>(fiR)) fiR = fiR.op(0);
374 if(fiL>0 && fiR<0) {
375 con_map[ex_to<numeric>(fiL).to_int()].append(fiR);
376 } else if(fiR>0 && fiL<0) {
377 con_map[ex_to<numeric>(fiR).to_int()].append(fiL);
378 } else {
379 if(fiL>0 && is_zero(fiR)) key = fiL;
380 else if(fiR>0 && is_zero(fiL)) key = fiR;
381 else throw Error("ShrinkCut: unexpcected point reached.");
382 val = li.op(2).op(0);
383 con_map[ex_to<numeric>(key).to_int()].append(val);
384 }
385 }
386
387 lst item;
388 for(auto kv : con_map) item.append(kv.second.sort());
389 ret.push_back(item);
390 });
391
392 return ret;
393 }
394
401 bool HasLoop(ex amp, lst prop) {
402 auto tls = TopoLines(amp);
403 int ntls = tls.nops();
404 // shrink internal lines of prop
405 int last = 0;
406 while(true) {
407 ex lp = 0;
408 for(int i=last; i<ntls; i++) {
409 auto li = tls.op(i);
410 if(is_zero(li) || li.op(2).nops()<2) continue; // external line
411 auto cpi = li.op(2);
412 bool m1 = !(is_zero(cpi.op(0)-prop.op(0)) && is_zero(cpi.op(1)-prop.op(1)));
413 bool m2 = !(is_zero(cpi.op(0)-prop.op(1)) && is_zero(cpi.op(1)-prop.op(0)));
414 if(m1 && m2) continue; // other propagators
415 if(is_zero(li.op(0)-li.op(1))) return true; // a loop found
416 last = i;
417 lp = li;
418 tls[last] = 0;
419 break;
420 }
421 if(is_zero(lp)) break;
422
423 for(int i=0; i<ntls; i++) {
424 if(is_zero(tls.op(i))) continue;
425 if(is_zero(tls.op(i).op(0)-lp.op(0))) tls[i][0] = lp.op(1);
426 if(is_zero(tls.op(i).op(1)-lp.op(0))) tls[i][1] = lp.op(1);
427 }
428 }
429 return false;
430 }
431
438 ex QuarkPropagator(const ex & e, const ex & m) {
439 auto fi1 = e.op(0).op(1);
440 auto fi2 = e.op(1).op(1);
441 auto mom = e.op(2);
442 return I * SP(TI(fi1),TI(fi2)) * GMat(GAS(mom)+GAS(1)*m, DI(fi1),DI(fi2)) / (SP(mom)-m*m);
443 }
444
451 ex LeptonPropagator(const ex & e, const ex & m) {
452 auto fi1 = e.op(0).op(1);
453 auto fi2 = e.op(1).op(1);
454 auto mom = e.op(2);
455 return I * GMat(GAS(mom)+GAS(1)*m, DI(fi1),DI(fi2)) / (SP(mom)-m*m);
456 }
457
464 ex GluonPropagator(const ex & e, const ex & xi) {
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(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) );
469 }
470
477 ex GluonGhostPropagator(const ex & e, 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 I * eta_G * SP(CI(fi1),CI(fi2)) / SP(mom);
482 }
483
491 ex AZWPropagator(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)-m*m) * (SP(LI(fi1),LI(fi2))-(1-xi)*SP(LI(fi1))*SP(LI(fi2))/(SP(mom)-xi*m*m));
496 }
497
505 ex AZWGhostPropagator(const ex & e, const ex & m, const ex & xi, const ex & eta_G) {
506 auto fi1 = e.op(0).op(1);
507 auto fi2 = e.op(1).op(1);
508 auto mom = e.op(2);
509 return eta_G * I / (SP(mom)-xi*m*m);
510 }
511
519 ex ScalarPropagator(const ex & e, const ex & m, const ex & xi) {
520 auto fi1 = e.op(0).op(1);
521 auto fi2 = e.op(1).op(1);
522 auto mom = e.op(2);
523 return I / (SP(mom)-xi*m*m);
524 }
525
532 ex QuarkGluonVertex(const ex & e, const ex & eta_s) {
533 auto fi1 = e.op(0).op(1);
534 auto fi2 = e.op(1).op(1);
535 auto fi3 = e.op(2).op(1);
536 return -I*eta_s*gs*GMat(GAS(LI(fi3)),DI(fi1),DI(fi2))*SUNT(CI(fi3),TI(fi1),TI(fi2));
537 }
538
545 ex Gluon3Vertex(const ex & e, const ex & eta_s) {
546 auto fi1 = e.op(0).op(1);
547 auto fi2 = e.op(1).op(1);
548 auto fi3 = e.op(2).op(1);
549 auto mom1 = e.op(0).op(2);
550 auto mom2 = e.op(1).op(2);
551 auto mom3 = e.op(2).op(2);
552 return -eta_s*gs*SUNF(CI(fi1),CI(fi2),CI(fi3))*(
553 SP(mom1-mom2,LI(fi3))*SP(LI(fi1),LI(fi2)) +
554 SP(mom2-mom3,LI(fi1))*SP(LI(fi2),LI(fi3)) +
555 SP(mom3-mom1,LI(fi2))*SP(LI(fi3),LI(fi1))
556 );
557 }
558
565 ex Gluon4Vertex(const ex & e, const ex & eta_s) {
566 auto fi1 = e.op(0).op(1);
567 auto fi2 = e.op(1).op(1);
568 auto fi3 = e.op(2).op(1);
569 auto fi4 = e.op(3).op(1);
570 return -I*gs*gs*(
571 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))) +
572 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))) +
573 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)))
574 );
575 }
576
584 ex GhostGluonVertex(const ex & e, const ex & eta_s, const ex & eta_G) {
585 auto fi1 = e.op(0).op(1);
586 auto fi2 = e.op(1).op(1);
587 auto fi3 = e.op(2).op(1);
588 auto mom1 = e.op(0).op(2);
589 return eta_s*eta_G*gs*SUNF(CI(fi1),CI(fi2),CI(fi3))*SP(mom1,LI(fi3));
590 }
591
599 ex QuarkAVertex(const ex & e, const ex & eq, const ex & eta_e) {
600 auto fi1 = e.op(0).op(1);
601 auto fi2 = e.op(1).op(1);
602 auto fi3 = e.op(2).op(1);
603 return -I*eta_e*eq*GMat(GAS(LI(fi3)),DI(fi1),DI(fi2))*SP(TI(fi1),TI(fi2));
604 }
605
613 ex LeptonAVertex(const ex & e, const ex & eq, const ex & eta_e) {
614 auto fi1 = e.op(0).op(1);
615 auto fi2 = e.op(1).op(1);
616 auto fi3 = e.op(2).op(1);
617 return -I*eta_e*eq*GMat(GAS(LI(fi3)),DI(fi1),DI(fi2));
618 }
619
626 ex IndexL2R(ex e, bool all) {
627 static MapFunction map([all](const ex &e, MapFunction &self)->ex {
628 if(!Index::has(e)) return e;
629 else if(is_a<Index>(e)) {
630 auto idx = ex_to<Index>(e);
631 auto nstr = idx.name.get_name();
632 if(!all && nstr.rfind("lim",0)==0) return e;
633 else if(!all && nstr.rfind("dim",0)==0) return e;
634 else if(!all && nstr.rfind("cim",0)==0) return e;
635 else if(!all && nstr.rfind("tim",0)==0) return e;
636 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);
637 else return e;
638 }
639 else return e.map(self);
640 });
641 return map(e);
642 }
643
644 ex IndexCC(ex e, bool all) {
645 return IndexL2R(e,all);
646 }
647
654 ex GluonSumL(int qi, bool color) {
655 if(color) return -SP(LI(qi), RLI(qi)) * SP(CI(qi), RCI(qi));
656 else return -SP(LI(qi), RLI(qi));
657 }
658
665 ex GluonSumR(int qi, bool color) {
666 if(color) return -SP(RLI(qi),LI(qi)) * SP(RCI(qi),CI(qi));
667 else return -SP(RLI(qi),LI(qi));
668 }
669
678 ex QuarkSumL(int qi, ex p, ex m, bool color) {
679 if(color) return GMat(GAS(p)+m*GAS(1), RDI(qi), DI(qi)) * SP(RTI(qi), TI(qi));
680 else return GMat(GAS(p)+m*GAS(1), RDI(qi), DI(qi));
681 }
682
691 ex QuarkSumR(int qi, ex p, ex m, bool color) {
692 if(color) return GMat(GAS(p)+m*GAS(1), DI(qi), RDI(qi)) * SP(TI(qi),RTI(qi));
693 else return GMat(GAS(p)+m*GAS(1), DI(qi), RDI(qi));
694 }
695
704 ex AntiQuarkSumL(int qi, ex p, ex m, bool color) {
705 if(color) return GMat(GAS(p)-m*GAS(1), DI(qi), RDI(qi)) * SP(TI(qi), RTI(qi));
706 else return GMat(GAS(p)-m*GAS(1), DI(qi), RDI(qi));
707 }
708
717 ex AntiQuarkSumR(int qi, ex p, ex m, bool color) {
718 if(color) return GMat(GAS(p)-m*GAS(1), RDI(qi), DI(qi)) * SP(RTI(qi), TI(qi));
719 else return GMat(GAS(p)-m*GAS(1), RDI(qi), DI(qi));
720 }
721
727 ex GhostSumL(int qi) {
728 return SP(CI(qi), RCI(qi));
729 }
730
736 ex GhostSumR(int qi) {
737 return SP(RCI(qi),CI(qi));
738 }
739
745 ex AntiGhostSumL(int qi) {
746 return -SP(CI(qi), RCI(qi));
747 }
748
754 ex AntiGhostSumR(int qi) {
755 return -SP(RCI(qi),CI(qi));
756 }
757
764 ex J1SumL(int qi, ex p) {
765 if(is_zero(p)) return -SP(LI(qi),RLI(qi));
766 else return -SP(LI(qi),RLI(qi)) + SP(p,LI(qi)) * SP(p,RLI(qi)) / SP(p);
767 }
768
775 ex J1SumR(int qi, ex p) {
776 if(is_zero(p)) return -SP(RLI(qi),LI(qi));
777 else return -SP(RLI(qi),LI(qi)) + SP(p,RLI(qi)) * SP(p,LI(qi)) / SP(p);
778 }
779
780 // https://arxiv.org/abs/1209.6213v2
781 lst Models::FeynRulesSM(const lst & amps, const ex & xi) {
782 lst ret;
783 for(auto item : amps) ret.append(FeynRulesSM(item,xi));
784 return ret;
785 }
786 ex Models::FeynRulesSM(const ex & amp, const ex & xi) {
787 if(is_a<lst>(amp)) return FeynRulesSM(ex_to<lst>(amp), xi);
788
789 static Symbol CW("CW"); // cos(theta)
790 static Symbol SW("SW"); // sin(theta)
791 //static Symbol C2W("C2W"); // cos(2theta)
792 static ex CW2 = CW*CW;
793 static ex SW2 = SW*SW;
794 static ex C2W = CW2-SW2;
795
796 static Symbol U("U"); // U-quark
797 static Symbol Ubar("Ubar"); // anti U-quark
798 static Symbol D("D"); // D-quark
799 static Symbol Dbar("Dbar"); // anti D-quark
800 static Symbol C("C"); // C-quark
801 static Symbol Cbar("Cbar"); // anti C-quark
802 static Symbol S("S"); // S-quark
803 static Symbol Sbar("Sbar"); // anti S-quark
804 static Symbol T("T"); // T-quark
805 static Symbol Tbar("Tbar"); // anti T-quark
806 static Symbol B("B"); // B-quark
807 static Symbol Bbar("Bbar"); // anti B-quark
808
809 static Symbol g("g"); // gluon
810 static Symbol gh("gh"); // gluon ghost
811 static Symbol ghbar("ghbar"); // anti gluon ghost
812
813 static Symbol A("A"); // photon
814 static Symbol Wm("Wm"); // W-
815 static Symbol Wp("Wp"); // W+
816 static Symbol Z("Z"); // Z
817
818 static Symbol ghA("ghA"); // photon ghost
819 static Symbol ghAbar("ghAbar"); // anti photon ghost
820 static Symbol ghWm("ghWm"); // W- ghost
821 static Symbol ghWmbar("ghWmbar"); // anti W- ghost
822 static Symbol ghWp("ghWp"); // W+ ghost
823 static Symbol ghWpbar("ghWpbar"); // anti W+ ghost
824 static Symbol ghZ("ghZ"); // Z ghost
825 static Symbol ghZbar("ghZbar"); // anti Z ghost
826
827 static Symbol em("em"); // e-
828 static Symbol ep("ep"); // e+
829 static Symbol ne("ne"); // e-neutrino
830 static Symbol nebar("nebar"); // anti e-neutrino
831 static Symbol mum("mum"); // mu-
832 static Symbol mup("mup"); // mu+
833 static Symbol nmu("nmu"); // mu-neutrino
834 static Symbol nmubar("nmubar"); // anti mu-neutrino
835 static Symbol taum("taum"); // tau-
836 static Symbol taup("taup"); // tau+
837 static Symbol ntau("ntau"); // tau-neutrino
838 static Symbol ntaubar("ntaubar"); // anti tau-neutrino
839
840 static Symbol chi("chi"); // Z goldstone
841 static Symbol phim("phim"); // W- goldstone
842 static Symbol phip("phip"); // W+ goldstone
843 static Symbol H("H"); // higgs
844
845 static auto M = [](const string & si)->ex {
846 return Symbol("M"+si);
847 };
848 static ex MW = M("W");
849 static ex MZ = M("Z");
850 static ex MH = M("H");
851
852 static Symbol EL("e"); // charge e
853
854 static ex eta_s = -1;
855 // Peskin's book convention
856 static ex eta = -1;
857 static ex eta_prime = -1;
858 static ex eta_Z = 1;
859 static ex eta_theta = 1;
860 static ex eta_Y = 1;
861 static ex eta_e = -1;
862 static ex eta_G = 1;
863
864 static ex GEW = eta_e*eta*eta_theta * EL/SW;
865 static ex EL2 = EL*EL;
866 static ex MW2 = MW*MW;
867 static ex MZ2 = MZ*MZ;
868 static ex MH2 = MH*MH;
869 static ex GEW2 = GEW*GEW;
870 static ex sqrt2 = sqrt(ex(2));
871
872 map<string,ex> T3;
873 T3["U"] = T3["C"] = T3["T"] = ex(1)/2;
874 T3["D"] = T3["S"] = T3["B"] = -ex(1)/2;
875 T3["ne"] = T3["nmu"] = T3["ntau"] = ex(1)/2;
876 T3["em"] = T3["mum"] = T3["taum"] = -ex(1)/2;
877 T3["ep"] = T3["mup"] = T3["taup"] = -ex(1)/2; // for sure
878
879 map<string,ex> Q;
880 Q["U"] = Q["C"] = Q["T"] = ex(2)/3;
881 Q["D"] = Q["S"] = Q["B"] = -ex(1)/3;
882 Q["ne"] = Q["nmu"] = Q["ntau"] = 0;
883 Q["em"] = Q["mum"] = Q["taum"] = -1;
884 Q["ep"] = Q["mup"] = Q["taup"] = -1; // for sure
885
886 auto is_qbar_q = [&](const ex pi1, const ex pi2)->bool {
887 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);
888 };
889 auto is_lbar_l = [&](const ex pi1, const ex pi2)->bool {
890 return (pi1==ep && pi2==em) || (pi1==mup && pi2==mum) || (pi1==taup && pi2==taum);
891 };
892 auto is_lbar_n = [&](const ex pi1, const ex pi2)->bool {
893 return (pi1==ep && pi2==ne) || (pi1==mup && pi2==nmu) || (pi1==taup && pi2==ntau);
894 };
895 auto is_nbar_n = [&](const ex pi1, const ex pi2)->bool {
896 return (pi1==nebar && pi2==ne) || (pi1==nmubar && pi2==nmu) || (pi1==ntaubar && pi2==ntau);
897 };
898 auto is_nbar_l = [&](const ex pi1, const ex pi2)->bool {
899 return (pi1==nebar && pi2==em) || (pi1==nmubar && pi2==mum) || (pi1==ntaubar && pi2==taum);
900 };
901 auto is_uqbar = [&](const ex pi)->bool {
902 return (pi==Ubar) || (pi==Cbar) || (pi==Tbar);
903 };
904 auto is_uq = [&](const ex pi)->bool {
905 return (pi==U) || (pi==C) || (pi==T);
906 };
907 auto is_dqbar = [&](const ex pi)->bool {
908 return (pi==Dbar) || (pi==Sbar) || (pi==Bbar);
909 };
910 auto is_dq = [&](const ex pi)->bool {
911 return (pi==D) || (pi==S) || (pi==B);
912 };
913 auto is_uqbar_dq = [&](const ex pi1, const ex pi2)->bool {
914 return is_uqbar(pi1) && is_dq(pi2);
915 };
916 auto is_dqbar_uq = [&](const ex pi1, const ex pi2)->bool {
917 return is_dqbar(pi1) && is_uq(pi2);
918 };
919 auto CKM = [&](const string & si1, const string & si2)->ex {
920 return Symbol("V"+si1+si2);
921 };
922
923 auto gfV = [&](const string & si)->ex {
924 return T3[si]/2-Q[si]*SW2;
925 };
926 auto gfA = [&](const string & si)->ex {
927 return T3[si]/2;
928 };
929
930 auto fr = MapFunction([&](const ex &e, MapFunction &self)->ex {
931 if(isFunction(e,"OutField") || isFunction(e,"InField")) return 1;
932 else if(isFunction(e, "Propagator")) {
933 auto pi = e.op(0).op(0);
934 auto si = ex2str(pi);
935 if(pi==U || pi==D || pi==C || pi==S || pi==T || pi==B) {
936 return QuarkPropagator(e, M(si));
937 } else if(pi==g) {
938 return GluonPropagator(e, xi);
939 } else if(pi==gh) {
940 return GluonGhostPropagator(e, xi, eta_G);
941 } else if(pi==A) {
942 return AZWPropagator(e, 0, xi);
943 } else if(pi==Z) {
944 return AZWPropagator(e, MZ, xi);
945 } else if(pi==Wm) {
946 return AZWPropagator(e, MW, xi);
947 } else if(pi==ghA) {
948 return AZWGhostPropagator(e, 0, xi, eta_G);
949 } else if(pi==ghZ) {
950 return AZWGhostPropagator(e, MZ, xi, eta_G);
951 } else if(pi==ghWm || pi==ghWp) {
952 return AZWGhostPropagator(e, MW, xi, eta_G);
953 } else if(pi==em || pi==mum || pi==taum) {
954 return LeptonPropagator(e, M(si));
955 } else if(pi==ne || pi==nmu || pi==ntau) {
956 return LeptonPropagator(e, 0);
957 } else if(pi==H) {
958 return ScalarPropagator(e, MH, 1);
959 } else if(pi==chi) {
960 return ScalarPropagator(e, MZ, xi);
961 } else if(pi==phim) {
962 return ScalarPropagator(e, MW, xi);
963 } else {
964 cout << endl << e << endl;
965 throw Error("Propagator Not defined!");
966 }
967 } else if(isFunction(e, "Vertex")) {
968 auto pi1 = e.op(0).op(0);
969 auto pi2 = e.op(1).op(0);
970 auto pi3 = e.op(2).op(0);
971 auto fi1 = e.op(0).op(1);
972 auto fi2 = e.op(1).op(1);
973 auto fi3 = e.op(2).op(1);
974 auto mom1 = e.op(0).op(2);
975 auto mom2 = e.op(1).op(2);
976 auto mom3 = e.op(2).op(2);
977 int nn = e.nops();
978
979 auto si1 = ex2str(pi1);
980 auto si2 = ex2str(pi2);
981 string_replace_all(si1, "bar", "");
982 string_replace_all(si2, "bar", "");
983 auto si = si2;
984
985 if(nn==3) {
986 if(pi1==ghbar && pi2==gh && pi3==g) {
987 // ghbar-gh-g
988 return GhostGluonVertex(e, eta_s, eta_G);
989 } else if(pi1==g && pi2==g && pi3==g) {
990 // g^3
991 return Gluon3Vertex(e, eta_s);
992 } else if((pi3==g) && is_qbar_q(pi1, pi2)) {
993 // qbar-q-g
994 return QuarkGluonVertex(e, eta_s);
995 } else if((pi3==A) && is_qbar_q(pi1, pi2)) {
996 // qbar-q-A
997 return QuarkAVertex(e, EL*Q[si], eta_e);
998 } else if((pi3==A) && is_lbar_l(pi1, pi2)) {
999 // lbar-l-A
1000 return LeptonAVertex(e, EL*Q[si], eta_e);
1001 } else if((pi3==Z) && is_qbar_q(pi1, pi2)) {
1002 // qbar-q-Z
1003 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));
1004 } else if((pi3==Z) && (is_lbar_l(pi1, pi2) || is_nbar_n(pi1, pi2))) {
1005 // lbar-l-Z & nbar-n-Z
1006 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)));
1007 } else if(pi1==Wp && pi2==Wm && pi3==A) {
1008 // Wp-Wm-A
1009 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)));
1010 } else if(pi1==Wp && pi2==Wm && pi3==Z) {
1011 // Wp-Wm-Z
1012 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)));
1013 } else if(is_uqbar(pi1) && is_dq(pi2) && pi3==Wp) {
1014 // uqbar-dq-Wp
1015 auto ckm = CKM(si1, si2);
1016 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));
1017 } else if(is_dqbar(pi1) && is_uq(pi2) && pi3==Wm) {
1018 // dqbar-uq-Wm
1019 auto ckm = CKM(si1, si2);
1020 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));
1021 } else if( (is_nbar_l(pi1, pi2) && pi3==Wp) || (is_lbar_n(pi1, pi2) && pi3==Wm) ) {
1022 // nbar-l-Wp & lbar-n-Wm
1023 return -I*eta*GEW/sqrt2*GMat(GAS(LI(fi3))-GAS(LI(fi3))*GAS(5),DI(fi1),DI(fi2))/2;
1024 } else if(is_qbar_q(pi1, pi2) && (pi3==H)) {
1025 // qbar-q-H
1026 return -I*GEW/2*M(si)/MW*GMat(GAS(1),DI(fi1),DI(fi2))*SP(TI(fi1),TI(fi2));
1027 } else if( (is_lbar_l(pi1, pi2) || is_nbar_n(pi1, pi2)) && (pi3==H)) {
1028 // lbar-l-H & nbar-n-H
1029 return -I*GEW/2*M(si)/MW*GMat(GAS(1),DI(fi1),DI(fi2));
1030 } else if(is_qbar_q(pi1, pi2) && (pi3==chi)) {
1031 // qbar-q-chi
1032 return -GEW*T3[si]*M(si)/MW*GMat(GAS(5),DI(fi1),DI(fi2))*SP(TI(fi1),TI(fi2));
1033 } else if( (is_lbar_l(pi1, pi2) || is_nbar_n(pi1, pi2)) && (pi3==chi)) {
1034 // lbar-l-chi & nbar-n-chi
1035 return -GEW*T3[si]*M(si)/MW*GMat(GAS(5),DI(fi1),DI(fi2));
1036 } else if(is_uqbar_dq(pi1, pi2) && pi3==phip) {
1037 // uqbar-dq-phip
1038 auto ckm = CKM(si1,si2);
1039 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));
1040 } else if(is_dqbar_uq(pi1, pi2) && pi3==phim) {
1041 // dqbar-uq-phim
1042 auto ckm = CKM(si1,si2);
1043 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));
1044 } else if(is_nbar_l(pi1, pi2) && pi3==phip) {
1045 // nbar-l-phip
1046 return -I*GEW/sqrt2*M(si2)/MW*GMat(GAS(1)+GAS(5),DI(fi1),DI(fi2))/2; // Mnl = 0
1047 } else if(is_lbar_n(pi1, pi2) && pi3==phim) {
1048 // lbar-n-phim
1049 return -I*GEW/sqrt2*M(si1)/MW*GMat(GAS(1)-GAS(5),DI(fi1),DI(fi2))/2; // Mnl = 0
1050 } else if(pi1==A && pi2==phip && pi3==phim) {
1051 // A-phip-phim
1052 return -I*eta_e*EL*SP(mom2-mom3,LI(fi1));
1053 } else if(pi1==Z && pi2==phip && pi3==phim) {
1054 // Z-phip-phim
1055 return -I*eta*eta_Z*GEW*C2W/(2*CW)*SP(mom2-mom3,LI(fi1));
1056 } else if(pi1==Wp && pi2==phim && pi3==H) {
1057 // Wp-phim-H
1058 return I/2*eta*GEW*SP(mom2-mom3,LI(fi1));
1059 } else if(pi1==Wm && pi2==phip && pi3==H) {
1060 // Wm-phip-H
1061 return -I/2*eta*GEW*SP(mom2-mom3,LI(fi1));
1062 } else if( ((pi1==Wp && pi2==phim) || (pi1==Wm && pi2==phip)) && pi3==chi) {
1063 // Wp-phim-chi & Wm-phip-chi
1064 return -eta*GEW/2*SP(mom2-mom3,LI(fi1));
1065 } else if(pi1==Z && pi2==chi && pi3==H) {
1066 // Z-chi-H
1067 return -eta*eta_Z*GEW/(2*CW)*SP(mom2-mom3,LI(fi1));
1068 } else if( ((pi1==phim && pi2==Wp) || (pi1==phip && pi2==Wm)) && pi3==A) {
1069 // phim-Wp-A & phip-Wm-A
1070 return I*eta_e*eta*EL*MW*SP(LI(fi2),LI(fi3));
1071 } else if( ((pi1==phim && pi2==Wp) || (pi1==phip && pi2==Wm)) && pi3==Z) {
1072 // phim-Wp-Z & phip-Wm-Z
1073 return -I*eta_Z*GEW*MZ*SW2*SP(LI(fi2),LI(fi3));
1074 } else if(pi1==H && pi2==Wp && pi3==Wm) {
1075 // H-Wp-Wm
1076 return I*GEW*MW*SP(LI(fi2),LI(fi3));
1077 } else if(pi1==H && pi2==Z && pi3==Z) {
1078 // H-Z-Z
1079 return I*GEW/CW*MZ*SP(LI(fi2),LI(fi3));
1080 } else if(pi1==H && pi2==phim && pi3==phip) {
1081 // H-phim-phip
1082 return -I/2*GEW*MH2/MW;
1083 } else if(pi1==H && pi2==H && pi3==H) {
1084 // H-H-H
1085 return -3/ex(2)*I*GEW*MH2/MW;
1086 } else if(pi1==H && pi2==chi && pi3==chi) {
1087 // H-chi-chi
1088 return -I/2*GEW*MH2/MW;
1089 } else if(pi1==ghWpbar && pi2==ghWp && pi3==A) {
1090 // ghWpbar-ghWp-A
1091 return I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1092 } else if(pi1==ghWmbar && pi2==ghWm && pi3==A) {
1093 // ghWmbar-ghWm-A
1094 return -I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1095 } else if(pi1==ghWpbar && pi2==ghWp && pi3==Z) {
1096 // ghWpbar-ghWp-Z
1097 return I*eta_G*eta*eta_Z*GEW*CW*SP(mom1,LI(fi3));
1098 } else if(pi1==ghWmbar && pi2==ghWm && pi3==Z) {
1099 // ghWmbar-ghWm-Z
1100 return -I*eta_G*eta*eta_Z*GEW*CW*SP(mom1,LI(fi3));
1101 } else if(pi1==ghWpbar && pi2==ghZ && pi3==Wp) {
1102 // ghWpbar-ghZ-Wp
1103 return -I*eta_G*eta*eta_Z*GEW*CW*SP(mom1,LI(fi3));
1104 } else if(pi1==ghWmbar && pi2==ghZ && pi3==Wm) {
1105 // ghWmbar-ghZ-Wm
1106 return I*eta_G*eta*eta_Z*GEW*CW*SP(mom1,LI(fi3));
1107 } else if(pi1==ghWpbar && pi2==ghA && pi3==Wp) {
1108 // ghWpbar-ghA-Wp
1109 return -I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1110 } else if(pi1==ghWmbar && pi2==ghA && pi3==Wm) {
1111 // ghWmbar-ghA-Wm
1112 return I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1113 } else if(pi1==ghZbar && pi2==ghWp && pi3==Wm) {
1114 // ghZbar-ghWp-Wm
1115 return -I*eta_G*eta*GEW*CW*SP(mom1,LI(fi3));
1116 } else if(pi1==ghZbar && pi2==ghWm && pi3==Wp) {
1117 // ghZbar-ghWm-Wp
1118 return I*eta_G*eta*GEW*CW*SP(mom1,LI(fi3));
1119 } else if(pi1==ghAbar && pi2==ghWp && pi3==Wm) {
1120 // ghAbar-ghWp-Wm
1121 return -I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1122 } else if(pi1==ghAbar && pi2==ghWm && pi3==Wp) {
1123 // ghAbar-ghWm-Wp
1124 return I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1125 } else if(pi1==ghWpbar && pi2==ghWp && pi3==chi) {
1126 // ghWpbar-ghWp-chi
1127 return eta_G*GEW/2*xi*MW;
1128 } else if(pi1==ghWmbar && pi2==ghWm && pi3==chi) {
1129 // ghWmbar-ghWm-chi
1130 return -eta_G*GEW/2*xi*MW;
1131 } else if( ((pi1==ghWpbar && pi2==ghWp) ||(pi1==ghWmbar && pi2==ghWm)) && pi3==H) {
1132 // ghWpbar-ghWp-H & ghWmbar-ghWm-H
1133 return -I/2*eta_G*GEW*xi*MW;
1134 } else if(pi1==ghZbar && pi2==ghZ && pi3==H) {
1135 // ghZbar-ghZ-H
1136 return -eta_G*I*GEW/(2*CW)*xi*MZ;
1137 } else if( pi1==ghZbar && ((pi2==ghWp && pi3==phim) || (pi2==ghWm && pi3==phip)) ) {
1138 // ghZbar-ghWp-phim & ghZbar-ghWm-phip
1139 return I/2*eta_G*eta_Z*GEW*xi*MZ;
1140 } else if( ((pi1==ghWpbar && pi3==phip) || (pi1==ghWmbar && pi3==phim)) && pi2==ghZ) {
1141 // ghWpbar-ghZ-phip & ghWmbar-ghZ-phim
1142 return -I*eta_G*eta_Z*GEW*C2W/(2*CW)*xi*MW;
1143 } else if( ((pi1==ghWpbar && pi3==phip) || (pi1==ghWmbar && pi3==phim)) && pi2==ghA) {
1144 // ghWpbar-ghA-phip & ghWmbar-ghA-phim
1145 return -I*eta_G*eta_e*eta*EL*xi*MW;
1146 } else {
1147 cout << endl << e << endl;
1148 throw Error("Vertex Not defined!");
1149 }
1150 } else if(nn==4) {
1151 auto pi4 = e.op(3).op(0);
1152 auto fi4 = e.op(3).op(1);
1153 auto mom4 = e.op(3).op(2);
1154
1155 if(pi1==g && pi2==g && pi3==g && pi4==g) {
1156 // g^4
1157 return Gluon4Vertex(e, eta_s);
1158 } else if(pi1==Wp && pi2==Wm && pi3==A && pi4==A) {
1159 // Wp-Wm-A-A
1160 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)));
1161 } else if(pi1==Wp && pi2==Wm && pi3==Z && pi4==Z) {
1162 // Wp-Wm-Z-Z
1163 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)));
1164 } else if(pi1==Wp && pi2==Wm && pi3==A && pi4==Z) {
1165 // Wp-Wm-A-Z
1166 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)));
1167 } else if(pi1==Wp && pi2==Wp && pi3==Wm && pi4==Wm) {
1168 // Wp-Wp-Wm-Wm
1169 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)));
1170 } else if(pi1==Wp && pi2==Wm && pi3==H && pi4==H) {
1171 // Wp-Wm-H-H
1172 return I/2*GEW2*SP(LI(fi1),LI(fi2));
1173 } else if(pi1==Wp && pi2==Wm && pi3==chi && pi4==chi) {
1174 // Wp-Wm-chi-chi
1175 return I/2*GEW2*SP(LI(fi1),LI(fi2));
1176 } else if(pi1==Z && pi2==Z && pi3==H && pi4==H) {
1177 // Z-Z-H-H
1178 return I/2*GEW2/CW2*SP(LI(fi1),LI(fi2));
1179 } else if(pi1==Z && pi2==Z && pi3==chi && pi4==chi) {
1180 // Z-Z-chi-chi
1181 return I/2*GEW2/CW2*SP(LI(fi1),LI(fi2));
1182 } else if(pi1==A && pi2==A && pi3==phip && pi4==phim) {
1183 // A-A-phip-phim
1184 return 2*I*EL2*SP(LI(fi1),LI(fi2));
1185 } else if(pi1==Z && pi2==Z && pi3==phip && pi4==phim) {
1186 // Z-Z-phip-phim
1187 return I/2*pow(GEW*C2W/CW,2)*SP(LI(fi1),LI(fi2));
1188 } else if(pi1==Wp && pi2==Wm && pi3==phip && pi4==phim) {
1189 // Wp-Wm-phip-phim
1190 return I/2*GEW2*SP(LI(fi1),LI(fi2));
1191 } else if( ((pi1==Wp && pi3==phim) || (pi1==Wm && pi3==phip)) && pi2==Z && pi4==H) {
1192 // Wp-Z-phim-H & Wm-Z-phip-H
1193 return -I*eta_Z*GEW2*SW2/(2*CW)*SP(LI(fi1),LI(fi2));
1194 } else if(pi1==Wm && pi2==Z && pi3==phip && pi4==chi) {
1195 // Wm-Z-phip-chi
1196 return -eta_Z*GEW2*SW2/(2*CW)*SP(LI(fi1),LI(fi2));
1197 } else if(pi1==Wp && pi2==Z && pi3==phim && pi4==chi) {
1198 // Wp-Z-phim-chi
1199 return eta_Z*GEW2*SW2/(2*CW)*SP(LI(fi1),LI(fi2));
1200 } else if( ((pi1==Wm && pi3==phip) || (pi1==Wp && pi3==phim)) && pi2==A && pi4==H) {
1201 // Wm-A-phip-H & Wp-A-phim-H
1202 return I/2*eta_e*eta*EL*GEW*SP(LI(fi1),LI(fi2));
1203 } else if(pi1==Wp && pi2==A && pi3==phim && pi4==chi) {
1204 // Wp-A-phim-chi
1205 return -1/ex(2)*eta_e*eta*EL*GEW*SP(LI(fi1),LI(fi2));
1206 } else if(pi1==Wm && pi2==A && pi3==phip && pi4==chi) {
1207 // Wm-A-phip-chi
1208 return 1/ex(2)*eta_e*eta*EL*GEW*SP(LI(fi1),LI(fi2));
1209 } else if(pi1==Z && pi2==A && pi3==phip && pi4==phim) {
1210 // Z-A-phip-phim
1211 return I*eta_e*eta*eta_Z*EL*GEW*C2W/CW*SP(LI(fi1),LI(fi2));
1212 } else if(pi1==phim && pi2==phip && pi3==phim && pi4==phip) {
1213 // phim-phip-phim-phip
1214 return -I/2*GEW2*MH2/MW2;
1215 } else if(pi1==H && pi2==H && pi3==phim && pi4==phip) {
1216 // H-H-phim-phip
1217 return -I/4*GEW2*MH2/MW2;
1218 } else if(pi1==chi && pi2==chi && pi3==phim && pi4==phip) {
1219 // chi-chi-phim-phip
1220 return -I/4*GEW2*MH2/MW2;
1221 } else if(pi1==H && pi2==H && pi3==H && pi4==H) {
1222 // H-H-H-H
1223 return -3/ex(4)*I*GEW2*MH2/MW2;
1224 } else if(pi1==H && pi2==H && pi3==chi && pi4==chi) {
1225 // H-H-chi-chi
1226 return -I/4*GEW2*MH2/MW2;
1227 } else if(pi1==chi && pi2==chi && pi3==chi && pi4==chi) {
1228 // chi-chi-chi-chi
1229 return -3/ex(4)*I*GEW2*MH2/MW2;
1230 } else {
1231 cout << endl << e << endl;
1232 throw Error("Vertex Not defined!");
1233 }
1234 } else {
1235 cout << endl << e << endl;
1236 throw Error("Vertex Not defined!");
1237 }
1238 } else return e.map(self);
1239 return e;
1240 });
1241 return fr(amp);
1242 }
1243
1244
1245 lst Models::FeynRulesQCD(const lst & amps, const ex & xi) {
1246 lst ret;
1247 for(auto item : amps) ret.append(FeynRulesQCD(item,xi));
1248 return ret;
1249 }
1250 ex Models::FeynRulesQCD(const ex & amp, const ex & xi) {
1251 if(is_a<lst>(amp)) return FeynRulesQCD(ex_to<lst>(amp), xi);
1252
1253 static Symbol A("A"), Q("Q"), Qbar("Qbar"), q("q"), qbar("qbar"), g("g"), gh("gh"), ghbar("ghbar");
1254 static Symbol m("m"), eq("eq"), eQ("eQ");
1255
1256 auto fr = MapFunction([&](const ex &e, MapFunction &self)->ex {
1257 if(isFunction(e,"OutField") || isFunction(e,"InField")) return 1;
1258 else if(isFunction(e, "Propagator")) {
1259 if(e.op(0).op(0)==q) {
1260 return QuarkPropagator(e, 0);
1261 } else if(e.op(0).op(0)==Q) {
1262 return QuarkPropagator(e, m);
1263 } else if(e.op(0).op(0)==g) {
1264 return GluonPropagator(e, xi);
1265 } else if(e.op(0).op(0)==gh) {
1266 return GluonGhostPropagator(e, xi);
1267 } else {
1268 cout << "expr: " << e << endl;
1269 throw Error("Propagator Not Found!");
1270 }
1271 } else if(isFunction(e, "Vertex")) {
1272 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) ) {
1273 // qbar-q-g or Qbar-Q-g or g -> A
1274 if(e.op(2).op(0)==g) return QuarkGluonVertex(e);
1275 else {
1276 auto fi1 = e.op(0).op(1);
1277 auto fi2 = e.op(1).op(1);
1278 auto fi3 = e.op(2).op(1);
1279 if(e.op(1).op(0)==q) QuarkAVertex(e, eq);
1280 else if(e.op(1).op(0)==Q) QuarkAVertex(e, eQ);
1281 else throw Error("Vertex Error.");
1282 }
1283 } else if(e.nops()==3 && e.op(0).op(0)==ghbar && e.op(1).op(0)==gh) {
1284 // ghbar-gh-g
1285 return GhostGluonVertex(e);
1286 } else if(e.nops()==3 && e.op(0).op(0)==g && e.op(1).op(0)==g && e.op(2).op(0)==g) {
1287 // g^3
1288 return Gluon3Vertex(e);
1289 } 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) {
1290 // g^4
1291 return Gluon4Vertex(e);
1292 } else {
1293 cout << "expr: " << e << endl;
1294 throw Error("Vertex Not Found!");
1295 }
1296 } else return e.map(self);
1297 return e;
1298 });
1299 return fr(amp);
1300 }
1301
1302}
1303
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:704
HepLib::QGRAF::HasLoop
bool HasLoop(ex amp, lst prop)
check a amplitude has a loop w.r.t. propapagtor
Definition QGRAF.cpp:401
HepLib::QGRAF::AntiGhostSumL
ex AntiGhostSumL(int qi)
polarization sum for anti-ghost, note that we will add extra - here
Definition QGRAF.cpp:745
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:626
HepLib::QGRAF::GluonPropagator
ex GluonPropagator(const ex &e, const ex &xi)
Propagator for gluon.
Definition QGRAF.cpp:464
HepLib::QGRAF::LeptonAVertex
ex LeptonAVertex(const ex &e, const ex &eq, const ex &eta_e)
l-lbar-A vertex
Definition QGRAF.cpp:613
HepLib::QGRAF::Gluon3Vertex
ex Gluon3Vertex(const ex &e, const ex &eta_s)
g-g-g vertex
Definition QGRAF.cpp:545
HepLib::QGRAF::GhostSumL
ex GhostSumL(int qi)
polarization sum for ghost
Definition QGRAF.cpp:727
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:754
HepLib::QGRAF::QuarkGluonVertex
ex QuarkGluonVertex(const ex &e, const ex &eta_s)
q-qbar-g vertex
Definition QGRAF.cpp:532
HepLib::QGRAF::ScalarPropagator
ex ScalarPropagator(const ex &e, const ex &m, const ex &xi)
Propagator for scalar.
Definition QGRAF.cpp:519
HepLib::QGRAF::J1SumL
ex J1SumL(int qi, ex p)
polarization sum for total angular momentum
Definition QGRAF.cpp:764
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:451
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:313
HepLib::QGRAF::AntiQuarkSumR
ex AntiQuarkSumR(int qi, ex p, ex m, bool color)
polarization sum for anti-quark
Definition QGRAF.cpp:717
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:775
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:505
HepLib::QGRAF::Gluon4Vertex
ex Gluon4Vertex(const ex &e, const ex &eta_s)
g-g-g-g vertex
Definition QGRAF.cpp:565
HepLib::QGRAF::IndexCC
ex IndexCC(ex e, bool all)
Definition QGRAF.cpp:644
HepLib::QGRAF::GhostGluonVertex
ex GhostGluonVertex(const ex &e, const ex &eta_s, const ex &eta_G)
ghbar-gh-g vertex
Definition QGRAF.cpp:584
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:438
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:691
HepLib::QGRAF::QuarkSumL
ex QuarkSumL(int qi, ex p, ex m, bool color)
polarization sum for quark
Definition QGRAF.cpp:678
HepLib::QGRAF::GhostSumR
ex GhostSumR(int qi)
polarization sum for ghost
Definition QGRAF.cpp:736
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:491
HepLib::QGRAF::QuarkAVertex
ex QuarkAVertex(const ex &e, const ex &eq, const ex &eta_e)
qbar-q-A vertex
Definition QGRAF.cpp:599
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:477
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:665
HepLib::QGRAF::GluonSumL
ex GluonSumL(int qi, bool color)
polarization sum for gluon
Definition QGRAF.cpp:654
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:281
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::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:1125
HepLib::SP
ex SP(const ex &a, bool use_map=false)
Definition Pair.cpp:166
Generated by doxygen 1.9.8