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
168 void DrawPDF(const lst & amps, string fn, int nr) {
169 int id=0, rc;
170 exvector amp_vec;
171 for(auto item : amps) amp_vec.push_back(item);
172 string tex_path = to_string(getpid()) + "_TeX/";
173 if(!dir_exists(tex_path)) rc = system(("mkdir -p "+tex_path).c_str());
174
175 GiNaC_Parallel(amp_vec.size(), [&amp_vec,tex_path](int idx)->ex {
176 auto amp = amp_vec[idx];
177 ofstream out(tex_path+to_string(idx)+".tex");
178 out << "\\documentclass[tikz]{standalone}" << endl;
179 out << "\\usepackage{tikz-feynman}" << endl;
180 out << "\\tikzfeynmanset{compat=1.1.0}" << endl;
181 out << "\\begin{document}" << endl;
182 out << "\\feynmandiagram{" << endl;
183 auto lines = TopoLines(amp);
184
185 exmap bend_map;
186 std::map<ex,int,ex_is_less> vtex_map; // vertex option, only once
187 for(auto l : lines) {
188 lst ll = lst{l.op(0), l.op(1)};
189 bool isExt = (is_a<numeric>(ll.op(0)) && ll.op(0)<0) || (is_a<numeric>(ll.op(1)) && ll.op(1)<0);
190 ll.sort();
191 bend_map[ll] = bend_map[ll] + 1;
192
193 auto fidL = (is_a<numeric>(l.op(1)) ? l.op(1) : l.op(1).op(0));
194 out << "\"" << fidL << "\"";
195 if(fidL<0) {
196 out << "[particle=";
197 if(InOutTeX[fidL].length()>0) out << InOutTeX[fidL];
198 else out << fidL;
199 out << "]";
200 } else if(vtex_map[l.op(1).op(1)]==0) {
201 out << VerTeX[l.op(1).op(1)];
202 vtex_map[l.op(1).op(1)]=1;
203 }
204
205 out << " --[";
206 auto f = l.op(2).op(0);
207 if(LineTeX[f].length()>0) {
208 if(!isExt) out << LineTeX[f];
209 else {
210 auto cpos = LineTeX[f].find(", edge");
211 if(cpos>0) out << LineTeX[f].substr(0,cpos);
212 else out << LineTeX[f];
213 }
214 }
215 if(bend_map[ll]>2) out << ",half right";
216 else if(bend_map[ll]>1) out << ",half left";
217 if(is_zero(l.op(0)-l.op(1))) out << ",loop,distance=2cm";
218 out << "]";
219
220 auto fidR = (is_a<numeric>(l.op(0)) ? l.op(0) : l.op(0).op(0));
221 out << " \"" << fidR << "\"";
222 if(fidR<0) {
223 out << "[particle=";
224 if(InOutTeX[fidR].length()>0) out << InOutTeX[fidR];
225 else out<< fidR;
226 out << "]";
227 } else if(vtex_map[l.op(0).op(1)]==0) {
228 out << VerTeX[l.op(0).op(1)];
229 vtex_map[l.op(0).op(1)]=1;
230 }
231 out << ";" << endl;
232 }
233 out << "};" << endl;
234 out << "\\end{document}" << endl;
235 out.close();
236 auto rc = system(("cd "+tex_path+" && echo X | lualatex " + to_string(idx) + " 1>/dev/null").c_str());
237 return 0;
238 }, "TeX");
239
240 ofstream out(tex_path+"diagram.tex");
241 out << "\\let\\mypdfximage\\pdfximage" << endl;
242 out << "\\def\\pdfximage{\\immediate\\mypdfximage}" << endl;
243 out << "\\documentclass{standalone}" << endl;
244 out << "\\usepackage{graphicx}" << endl;
245 out << "\\usepackage{adjustbox}" << endl;
246 out << "\\usepackage{scalefnt}" << endl;
247 out << "\\begin{document}" << endl;
248 //out << "\\scalefont{" << .5/nr << "}" << endl;
249 out << "\\begin{adjustbox}{valign=T,width=\\textwidth}" << endl;
250 out << "\\begin{tabular}{|"; for(int i=0; i<nr; i++) out << "cc|"; out << "}" << endl;
251 out << "\\hline" << endl;
252 int total = amps.nops();
253 int namps = total;
254 if((total%nr)!=0) total = (total/nr+1)*nr;
255 for(int i=0; i<total; i++) {
256
257 //int limit = 300;
258 //if((i!=0) && (i+1!=total)) {
259 // out << "\\end{tabular}" << endl << endl;
260 // out << "\\begin{tabular}{|"; for(int i=0; i<nr; i++) out << "cc|"; out << "}" << endl;
261 // out << "\\hline" << endl;
262 //}
263
264 out << "{\\tiny " << i+1 << "}&" << endl;
265 if(i<namps) {
266 out << "\\includegraphics[keepaspectratio,";
267 out << "height=" << 1.0/nr << "\\textwidth,";
268 out << "width=" << 1.0/nr << "\\textwidth]";
269 out << "{"<<i<<".pdf}" << endl;
270 }
271 if((i+1)%nr==0) out << "\\\\ \\hline";
272 else out << "&";
273 }
274 out << "\\end{tabular}" << endl;
275 out << "\\end{adjustbox}" << endl;
276 out << "\\end{document}" << endl;
277 out.close();
278 if(Debug) rc = system(("cd "+tex_path+" && pdflatex diagram && mv diagram.pdf ../"+fn).c_str());
279 else rc = system(("cd "+tex_path+" && echo X | pdflatex diagram 1>/dev/null && mv diagram.pdf ../"+fn).c_str());
280 if(!Debug) rc = system(("rm -r "+tex_path).c_str());
281 }
282
290 vector<lst> ShrinkCut(ex amp, lst prop, int n) {
291 vector<lst> ret;
292 if(prop.nops()<1) throw Error("ShrinkCut: no cut provided!");
293 auto tls = TopoLines(amp);
294 vector<int> cls_vec;
295 for(int i=0; i<tls.nops(); i++) {
296 auto pi = tls.op(i).op(2);
297 if(pi.nops()<2) continue;
298 if(!is_a<lst>(prop.op(0))) { // prop : lst{ g, g }
299 if(is_zero(pi.op(0)-prop.op(0)) && is_zero(pi.op(1)-prop.op(1))) cls_vec.push_back(i);
300 } else {
301 for(auto iprop : prop) {
302 if(is_zero(pi.op(0)-iprop.op(0)) && is_zero(pi.op(1)-iprop.op(1))) cls_vec.push_back(i);
303 }
304 }
305 }
306 if(cls_vec.size()<n) return ret;
307
308 Combinations(cls_vec.size(), n, [n,&ret,cls_vec,tls](const int * is)->void {
309 int cls[n];
310 for(int i=0; i<n; i++) cls[i] = cls_vec[is[i]];
311
312 // cut each line into 2 half-lines, labeled with 0
313 auto tls2 = tls;
314 for(auto ci : cls) {
315 auto ol = tls2.op(ci);
316 tls2.let_op(ci) = lst{ol.op(0), 0, lst{ol.op(2).op(0)}, ol.op(3)};
317 tls2.append(lst{0, ol.op(1), lst{ol.op(2).op(1)}, ol.op(3)} );
318 }
319
320 // shrink internal lines
321 int last = 0;
322 int ntls2 = tls2.nops();
323 while(true) {
324 ex lp = 0;
325 for(int i=last; i<ntls2; i++) {
326 auto li = tls2.op(i);
327 if(is_zero(li) || li.op(2).nops()<2) continue;
328 last = i;
329 lp = li;
330 tls2.let_op(last) = 0;
331 break;
332 }
333 if(is_zero(lp)) break;
334
335 for(int i=0; i<ntls2; i++) {
336 if(is_zero(tls2.op(i))) continue;
337 if(is_zero(tls2.op(i).op(0)-lp.op(0))) tls2.let_op(i).let_op(0) = lp.op(1);
338 if(is_zero(tls2.op(i).op(1)-lp.op(0))) tls2.let_op(i).let_op(1) = lp.op(1);
339 }
340 }
341
342 // final connected parts
343 map<int, lst> con_map;
344 for(auto li : tls2) {
345 if(is_zero(li)) continue;
346 ex key, val;
347 ex fiL = li.op(0);
348 if(is_a<lst>(fiL)) fiL = fiL.op(0);
349 ex fiR = li.op(1);
350 if(is_a<lst>(fiR)) fiR = fiR.op(0);
351 if(fiL>0 && fiR<0) {
352 con_map[ex_to<numeric>(fiL).to_int()].append(fiR);
353 } else if(fiR>0 && fiL<0) {
354 con_map[ex_to<numeric>(fiR).to_int()].append(fiL);
355 } else {
356 if(fiL>0 && is_zero(fiR)) key = fiL;
357 else if(fiR>0 && is_zero(fiL)) key = fiR;
358 else throw Error("ShrinkCut: unexpcected point reached.");
359 val = li.op(2).op(0);
360 con_map[ex_to<numeric>(key).to_int()].append(val);
361 }
362 }
363
364 lst item;
365 for(auto kv : con_map) item.append(kv.second.sort());
366 ret.push_back(item);
367 });
368
369 return ret;
370 }
371
378 bool HasLoop(ex amp, lst prop) {
379 auto tls = TopoLines(amp);
380 int ntls = tls.nops();
381 // shrink internal lines of prop
382 int last = 0;
383 while(true) {
384 ex lp = 0;
385 for(int i=last; i<ntls; i++) {
386 auto li = tls.op(i);
387 if(is_zero(li) || li.op(2).nops()<2) continue; // external line
388 auto cpi = li.op(2);
389 bool m1 = !(is_zero(cpi.op(0)-prop.op(0)) && is_zero(cpi.op(1)-prop.op(1)));
390 bool m2 = !(is_zero(cpi.op(0)-prop.op(1)) && is_zero(cpi.op(1)-prop.op(0)));
391 if(m1 && m2) continue; // other propagators
392 if(is_zero(li.op(0)-li.op(1))) return true; // a loop found
393 last = i;
394 lp = li;
395 tls.let_op(last) = 0;
396 break;
397 }
398 if(is_zero(lp)) break;
399
400 for(int i=0; i<ntls; i++) {
401 if(is_zero(tls.op(i))) continue;
402 if(is_zero(tls.op(i).op(0)-lp.op(0))) tls.let_op(i).let_op(0) = lp.op(1);
403 if(is_zero(tls.op(i).op(1)-lp.op(0))) tls.let_op(i).let_op(1) = lp.op(1);
404 }
405 }
406 return false;
407 }
408
415 ex QuarkPropagator(const ex & e, const ex & m) {
416 auto fi1 = e.op(0).op(1);
417 auto fi2 = e.op(1).op(1);
418 auto mom = e.op(2);
419 return I * SP(TI(fi1),TI(fi2)) * GMat(GAS(mom)+GAS(1)*m, DI(fi1),DI(fi2)) / (SP(mom)-m*m);
420 }
421
428 ex LeptonPropagator(const ex & e, const ex & m) {
429 auto fi1 = e.op(0).op(1);
430 auto fi2 = e.op(1).op(1);
431 auto mom = e.op(2);
432 return I * GMat(GAS(mom)+GAS(1)*m, DI(fi1),DI(fi2)) / (SP(mom)-m*m);
433 }
434
441 ex GluonPropagator(const ex & e, const ex & xi) {
442 auto fi1 = e.op(0).op(1);
443 auto fi2 = e.op(1).op(1);
444 auto mom = e.op(2);
445 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) );
446 }
447
454 ex GluonGhostPropagator(const ex & e, const ex & xi, const ex & eta_G) {
455 auto fi1 = e.op(0).op(1);
456 auto fi2 = e.op(1).op(1);
457 auto mom = e.op(2);
458 return I * eta_G * SP(CI(fi1),CI(fi2)) / SP(mom);
459 }
460
468 ex AZWPropagator(const ex & e, const ex & m, const ex & xi) {
469 auto fi1 = e.op(0).op(1);
470 auto fi2 = e.op(1).op(1);
471 auto mom = e.op(2);
472 return -I/(SP(mom)-m*m) * (SP(LI(fi1),LI(fi2))-(1-xi)*SP(LI(fi1))*SP(LI(fi2))/(SP(mom)-xi*m*m));
473 }
474
482 ex AZWGhostPropagator(const ex & e, const ex & m, const ex & xi, const ex & eta_G) {
483 auto fi1 = e.op(0).op(1);
484 auto fi2 = e.op(1).op(1);
485 auto mom = e.op(2);
486 return eta_G * I / (SP(mom)-xi*m*m);
487 }
488
496 ex ScalarPropagator(const ex & e, const ex & m, const ex & xi) {
497 auto fi1 = e.op(0).op(1);
498 auto fi2 = e.op(1).op(1);
499 auto mom = e.op(2);
500 return I / (SP(mom)-xi*m*m);
501 }
502
509 ex QuarkGluonVertex(const ex & e, const ex & eta_s) {
510 auto fi1 = e.op(0).op(1);
511 auto fi2 = e.op(1).op(1);
512 auto fi3 = e.op(2).op(1);
513 return -I*eta_s*gs*GMat(GAS(LI(fi3)),DI(fi1),DI(fi2))*SUNT(CI(fi3),TI(fi1),TI(fi2));
514 }
515
522 ex Gluon3Vertex(const ex & e, const ex & eta_s) {
523 auto fi1 = e.op(0).op(1);
524 auto fi2 = e.op(1).op(1);
525 auto fi3 = e.op(2).op(1);
526 auto mom1 = e.op(0).op(2);
527 auto mom2 = e.op(1).op(2);
528 auto mom3 = e.op(2).op(2);
529 return -eta_s*gs*SUNF(CI(fi1),CI(fi2),CI(fi3))*(
530 SP(mom1-mom2,LI(fi3))*SP(LI(fi1),LI(fi2)) +
531 SP(mom2-mom3,LI(fi1))*SP(LI(fi2),LI(fi3)) +
532 SP(mom3-mom1,LI(fi2))*SP(LI(fi3),LI(fi1))
533 );
534 }
535
542 ex Gluon4Vertex(const ex & e, const ex & eta_s) {
543 auto fi1 = e.op(0).op(1);
544 auto fi2 = e.op(1).op(1);
545 auto fi3 = e.op(2).op(1);
546 auto fi4 = e.op(3).op(1);
547 return -I*gs*gs*(
548 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))) +
549 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))) +
550 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)))
551 );
552 }
553
561 ex GhostGluonVertex(const ex & e, const ex & eta_s, const ex & eta_G) {
562 auto fi1 = e.op(0).op(1);
563 auto fi2 = e.op(1).op(1);
564 auto fi3 = e.op(2).op(1);
565 auto mom1 = e.op(0).op(2);
566 return eta_s*eta_G*gs*SUNF(CI(fi1),CI(fi2),CI(fi3))*SP(mom1,LI(fi3));
567 }
568
576 ex QuarkAVertex(const ex & e, const ex & eq, const ex & eta_e) {
577 auto fi1 = e.op(0).op(1);
578 auto fi2 = e.op(1).op(1);
579 auto fi3 = e.op(2).op(1);
580 return -I*eta_e*eq*GMat(GAS(LI(fi3)),DI(fi1),DI(fi2))*SP(TI(fi1),TI(fi2));
581 }
582
590 ex LeptonAVertex(const ex & e, const ex & eq, const ex & eta_e) {
591 auto fi1 = e.op(0).op(1);
592 auto fi2 = e.op(1).op(1);
593 auto fi3 = e.op(2).op(1);
594 return -I*eta_e*eq*GMat(GAS(LI(fi3)),DI(fi1),DI(fi2));
595 }
596
603 ex IndexL2R(ex e, bool all) {
604 static MapFunction map([all](const ex &e, MapFunction &self)->ex {
605 if(!Index::has(e)) return e;
606 else if(is_a<Index>(e)) {
607 auto idx = ex_to<Index>(e);
608 auto nstr = idx.name.get_name();
609 if(!all && nstr.rfind("lim",0)==0) return e;
610 else if(!all && nstr.rfind("dim",0)==0) return e;
611 else if(!all && nstr.rfind("cim",0)==0) return e;
612 else if(!all && nstr.rfind("tim",0)==0) return e;
613 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);
614 else return e;
615 }
616 else return e.map(self);
617 });
618 return map(e);
619 }
620
621 ex IndexCC(ex e, bool all) {
622 return IndexL2R(e,all);
623 }
624
631 ex GluonSumL(int qi, bool color) {
632 if(color) return -SP(LI(qi), RLI(qi)) * SP(CI(qi), RCI(qi));
633 else return -SP(LI(qi), RLI(qi));
634 }
635
642 ex GluonSumR(int qi, bool color) {
643 if(color) return -SP(RLI(qi),LI(qi)) * SP(RCI(qi),CI(qi));
644 else return -SP(RLI(qi),LI(qi));
645 }
646
655 ex QuarkSumL(int qi, ex p, ex m, bool color) {
656 if(color) return GMat(GAS(p)+m*GAS(1), RDI(qi), DI(qi)) * SP(RTI(qi), TI(qi));
657 else return GMat(GAS(p)+m*GAS(1), RDI(qi), DI(qi));
658 }
659
668 ex QuarkSumR(int qi, ex p, ex m, bool color) {
669 if(color) return GMat(GAS(p)+m*GAS(1), DI(qi), RDI(qi)) * SP(TI(qi),RTI(qi));
670 else return GMat(GAS(p)+m*GAS(1), DI(qi), RDI(qi));
671 }
672
681 ex AntiQuarkSumL(int qi, ex p, ex m, bool color) {
682 if(color) return GMat(GAS(p)-m*GAS(1), DI(qi), RDI(qi)) * SP(TI(qi), RTI(qi));
683 else return GMat(GAS(p)-m*GAS(1), DI(qi), RDI(qi));
684 }
685
694 ex AntiQuarkSumR(int qi, ex p, ex m, bool color) {
695 if(color) return GMat(GAS(p)-m*GAS(1), RDI(qi), DI(qi)) * SP(RTI(qi), TI(qi));
696 else return GMat(GAS(p)-m*GAS(1), RDI(qi), DI(qi));
697 }
698
704 ex GhostSumL(int qi) {
705 return SP(CI(qi), RCI(qi));
706 }
707
713 ex GhostSumR(int qi) {
714 return SP(RCI(qi),CI(qi));
715 }
716
722 ex AntiGhostSumL(int qi) {
723 return -SP(CI(qi), RCI(qi));
724 }
725
731 ex AntiGhostSumR(int qi) {
732 return -SP(RCI(qi),CI(qi));
733 }
734
741 ex J1SumL(int qi, ex p) {
742 if(is_zero(p)) return -SP(LI(qi),RLI(qi));
743 else return -SP(LI(qi),RLI(qi)) + SP(p,LI(qi)) * SP(p,RLI(qi)) / SP(p);
744 }
745
752 ex J1SumR(int qi, ex p) {
753 if(is_zero(p)) return -SP(RLI(qi),LI(qi));
754 else return -SP(RLI(qi),LI(qi)) + SP(p,RLI(qi)) * SP(p,LI(qi)) / SP(p);
755 }
756
757 // https://arxiv.org/abs/1209.6213v2
758 lst Models::FeynRulesSM(const lst & amps, const ex & xi) {
759 lst ret;
760 for(auto item : amps) ret.append(FeynRulesSM(item,xi));
761 return ret;
762 }
763 ex Models::FeynRulesSM(const ex & amp, const ex & xi) {
764 if(is_a<lst>(amp)) return FeynRulesSM(ex_to<lst>(amp), xi);
765
766 static Symbol CW("CW"); // cos(theta)
767 static Symbol SW("SW"); // sin(theta)
768 //static Symbol C2W("C2W"); // cos(2theta)
769 static ex CW2 = CW*CW;
770 static ex SW2 = SW*SW;
771 static ex C2W = CW2-SW2;
772
773 static Symbol U("U"); // U-quark
774 static Symbol Ubar("Ubar"); // anti U-quark
775 static Symbol D("D"); // D-quark
776 static Symbol Dbar("Dbar"); // anti D-quark
777 static Symbol C("C"); // C-quark
778 static Symbol Cbar("Cbar"); // anti C-quark
779 static Symbol S("S"); // S-quark
780 static Symbol Sbar("Sbar"); // anti S-quark
781 static Symbol T("T"); // T-quark
782 static Symbol Tbar("Tbar"); // anti T-quark
783 static Symbol B("B"); // B-quark
784 static Symbol Bbar("Bbar"); // anti B-quark
785
786 static Symbol g("g"); // gluon
787 static Symbol gh("gh"); // gluon ghost
788 static Symbol ghbar("ghbar"); // anti gluon ghost
789
790 static Symbol A("A"); // photon
791 static Symbol Wm("Wm"); // W-
792 static Symbol Wp("Wp"); // W+
793 static Symbol Z("Z"); // Z
794
795 static Symbol ghA("ghA"); // photon ghost
796 static Symbol ghAbar("ghAbar"); // anti photon ghost
797 static Symbol ghWm("ghWm"); // W- ghost
798 static Symbol ghWmbar("ghWmbar"); // anti W- ghost
799 static Symbol ghWp("ghWp"); // W+ ghost
800 static Symbol ghWpbar("ghWpbar"); // anti W+ ghost
801 static Symbol ghZ("ghZ"); // Z ghost
802 static Symbol ghZbar("ghZbar"); // anti Z ghost
803
804 static Symbol em("em"); // e-
805 static Symbol ep("ep"); // e+
806 static Symbol ne("ne"); // e-neutrino
807 static Symbol nebar("nebar"); // anti e-neutrino
808 static Symbol mum("mum"); // mu-
809 static Symbol mup("mup"); // mu+
810 static Symbol nmu("nmu"); // mu-neutrino
811 static Symbol nmubar("nmubar"); // anti mu-neutrino
812 static Symbol taum("taum"); // tau-
813 static Symbol taup("taup"); // tau+
814 static Symbol ntau("ntau"); // tau-neutrino
815 static Symbol ntaubar("ntaubar"); // anti tau-neutrino
816
817 static Symbol chi("chi"); // Z goldstone
818 static Symbol phim("phim"); // W- goldstone
819 static Symbol phip("phip"); // W+ goldstone
820 static Symbol H("H"); // higgs
821
822 static auto M = [](const string & si)->ex {
823 return Symbol("M"+si);
824 };
825 static ex MW = M("W");
826 static ex MZ = M("Z");
827 static ex MH = M("H");
828
829 static Symbol EL("e"); // charge e
830
831 static ex eta_s = -1;
832 // Peskin's book convention
833 static ex eta = -1;
834 static ex eta_prime = -1;
835 static ex eta_Z = 1;
836 static ex eta_theta = 1;
837 static ex eta_Y = 1;
838 static ex eta_e = -1;
839 static ex eta_G = 1;
840
841 static ex GEW = eta_e*eta*eta_theta * EL/SW;
842 static ex EL2 = EL*EL;
843 static ex MW2 = MW*MW;
844 static ex MZ2 = MZ*MZ;
845 static ex MH2 = MH*MH;
846 static ex GEW2 = GEW*GEW;
847 static ex sqrt2 = sqrt(ex(2));
848
849 map<string,ex> T3;
850 T3["U"] = T3["C"] = T3["T"] = ex(1)/2;
851 T3["D"] = T3["S"] = T3["B"] = -ex(1)/2;
852 T3["ne"] = T3["nmu"] = T3["ntau"] = ex(1)/2;
853 T3["em"] = T3["mum"] = T3["taum"] = -ex(1)/2;
854 T3["ep"] = T3["mup"] = T3["taup"] = -ex(1)/2; // for sure
855
856 map<string,ex> Q;
857 Q["U"] = Q["C"] = Q["T"] = ex(2)/3;
858 Q["D"] = Q["S"] = Q["B"] = -ex(1)/3;
859 Q["ne"] = Q["nmu"] = Q["ntau"] = 0;
860 Q["em"] = Q["mum"] = Q["taum"] = -1;
861 Q["ep"] = Q["mup"] = Q["taup"] = -1; // for sure
862
863 auto is_qbar_q = [&](const ex pi1, const ex pi2)->bool {
864 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);
865 };
866 auto is_lbar_l = [&](const ex pi1, const ex pi2)->bool {
867 return (pi1==ep && pi2==em) || (pi1==mup && pi2==mum) || (pi1==taup && pi2==taum);
868 };
869 auto is_lbar_n = [&](const ex pi1, const ex pi2)->bool {
870 return (pi1==ep && pi2==ne) || (pi1==mup && pi2==nmu) || (pi1==taup && pi2==ntau);
871 };
872 auto is_nbar_n = [&](const ex pi1, const ex pi2)->bool {
873 return (pi1==nebar && pi2==ne) || (pi1==nmubar && pi2==nmu) || (pi1==ntaubar && pi2==ntau);
874 };
875 auto is_nbar_l = [&](const ex pi1, const ex pi2)->bool {
876 return (pi1==nebar && pi2==em) || (pi1==nmubar && pi2==mum) || (pi1==ntaubar && pi2==taum);
877 };
878 auto is_uqbar = [&](const ex pi)->bool {
879 return (pi==Ubar) || (pi==Cbar) || (pi==Tbar);
880 };
881 auto is_uq = [&](const ex pi)->bool {
882 return (pi==U) || (pi==C) || (pi==T);
883 };
884 auto is_dqbar = [&](const ex pi)->bool {
885 return (pi==Dbar) || (pi==Sbar) || (pi==Bbar);
886 };
887 auto is_dq = [&](const ex pi)->bool {
888 return (pi==D) || (pi==S) || (pi==B);
889 };
890 auto is_uqbar_dq = [&](const ex pi1, const ex pi2)->bool {
891 return is_uqbar(pi1) && is_dq(pi2);
892 };
893 auto is_dqbar_uq = [&](const ex pi1, const ex pi2)->bool {
894 return is_dqbar(pi1) && is_uq(pi2);
895 };
896 auto CKM = [&](const string & si1, const string & si2)->ex {
897 return Symbol("V"+si1+si2);
898 };
899
900 auto gfV = [&](const string & si)->ex {
901 return T3[si]/2-Q[si]*SW2;
902 };
903 auto gfA = [&](const string & si)->ex {
904 return T3[si]/2;
905 };
906
907 auto fr = MapFunction([&](const ex &e, MapFunction &self)->ex {
908 if(isFunction(e,"OutField") || isFunction(e,"InField")) return 1;
909 else if(isFunction(e, "Propagator")) {
910 auto pi = e.op(0).op(0);
911 auto si = ex2str(pi);
912 if(pi==U || pi==D || pi==C || pi==S || pi==T || pi==B) {
913 return QuarkPropagator(e, M(si));
914 } else if(pi==g) {
915 return GluonPropagator(e, xi);
916 } else if(pi==gh) {
917 return GluonGhostPropagator(e, xi, eta_G);
918 } else if(pi==A) {
919 return AZWPropagator(e, 0, xi);
920 } else if(pi==Z) {
921 return AZWPropagator(e, MZ, xi);
922 } else if(pi==Wm) {
923 return AZWPropagator(e, MW, xi);
924 } else if(pi==ghA) {
925 return AZWGhostPropagator(e, 0, xi, eta_G);
926 } else if(pi==ghZ) {
927 return AZWGhostPropagator(e, MZ, xi, eta_G);
928 } else if(pi==ghWm || pi==ghWp) {
929 return AZWGhostPropagator(e, MW, xi, eta_G);
930 } else if(pi==em || pi==mum || pi==taum) {
931 return LeptonPropagator(e, M(si));
932 } else if(pi==ne || pi==nmu || pi==ntau) {
933 return LeptonPropagator(e, 0);
934 } else if(pi==H) {
935 return ScalarPropagator(e, MH, 1);
936 } else if(pi==chi) {
937 return ScalarPropagator(e, MZ, xi);
938 } else if(pi==phim) {
939 return ScalarPropagator(e, MW, xi);
940 } else {
941 cout << endl << e << endl;
942 throw Error("Propagator Not defined!");
943 }
944 } else if(isFunction(e, "Vertex")) {
945 auto pi1 = e.op(0).op(0);
946 auto pi2 = e.op(1).op(0);
947 auto pi3 = e.op(2).op(0);
948 auto fi1 = e.op(0).op(1);
949 auto fi2 = e.op(1).op(1);
950 auto fi3 = e.op(2).op(1);
951 auto mom1 = e.op(0).op(2);
952 auto mom2 = e.op(1).op(2);
953 auto mom3 = e.op(2).op(2);
954 int nn = e.nops();
955
956 auto si1 = ex2str(pi1);
957 auto si2 = ex2str(pi2);
958 string_replace_all(si1, "bar", "");
959 string_replace_all(si2, "bar", "");
960 auto si = si2;
961
962 if(nn==3) {
963 if(pi1==ghbar && pi2==gh && pi3==g) {
964 // ghbar-gh-g
965 return GhostGluonVertex(e, eta_s, eta_G);
966 } else if(pi1==g && pi2==g && pi3==g) {
967 // g^3
968 return Gluon3Vertex(e, eta_s);
969 } else if((pi3==g) && is_qbar_q(pi1, pi2)) {
970 // qbar-q-g
971 return QuarkGluonVertex(e, eta_s);
972 } else if((pi3==A) && is_qbar_q(pi1, pi2)) {
973 // qbar-q-A
974 return QuarkAVertex(e, EL*Q[si], eta_e);
975 } else if((pi3==A) && is_lbar_l(pi1, pi2)) {
976 // lbar-l-A
977 return LeptonAVertex(e, EL*Q[si], eta_e);
978 } else if((pi3==Z) && is_qbar_q(pi1, pi2)) {
979 // qbar-q-Z
980 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));
981 } else if((pi3==Z) && (is_lbar_l(pi1, pi2) || is_nbar_n(pi1, pi2))) {
982 // lbar-l-Z & nbar-n-Z
983 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)));
984 } else if(pi1==Wp && pi2==Wm && pi3==A) {
985 // Wp-Wm-A
986 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)));
987 } else if(pi1==Wp && pi2==Wm && pi3==Z) {
988 // Wp-Wm-Z
989 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)));
990 } else if(is_uqbar(pi1) && is_dq(pi2) && pi3==Wp) {
991 // uqbar-dq-Wp
992 auto ckm = CKM(si1, si2);
993 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));
994 } else if(is_dqbar(pi1) && is_uq(pi2) && pi3==Wm) {
995 // dqbar-uq-Wm
996 auto ckm = CKM(si1, si2);
997 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));
998 } else if( (is_nbar_l(pi1, pi2) && pi3==Wp) || (is_lbar_n(pi1, pi2) && pi3==Wm) ) {
999 // nbar-l-Wp & lbar-n-Wm
1000 return -I*eta*GEW/sqrt2*GMat(GAS(LI(fi3))-GAS(LI(fi3))*GAS(5),DI(fi1),DI(fi2))/2;
1001 } else if(is_qbar_q(pi1, pi2) && (pi3==H)) {
1002 // qbar-q-H
1003 return -I*GEW/2*M(si)/MW*GMat(GAS(1),DI(fi1),DI(fi2))*SP(TI(fi1),TI(fi2));
1004 } else if( (is_lbar_l(pi1, pi2) || is_nbar_n(pi1, pi2)) && (pi3==H)) {
1005 // lbar-l-H & nbar-n-H
1006 return -I*GEW/2*M(si)/MW*GMat(GAS(1),DI(fi1),DI(fi2));
1007 } else if(is_qbar_q(pi1, pi2) && (pi3==chi)) {
1008 // qbar-q-chi
1009 return -GEW*T3[si]*M(si)/MW*GMat(GAS(5),DI(fi1),DI(fi2))*SP(TI(fi1),TI(fi2));
1010 } else if( (is_lbar_l(pi1, pi2) || is_nbar_n(pi1, pi2)) && (pi3==chi)) {
1011 // lbar-l-chi & nbar-n-chi
1012 return -GEW*T3[si]*M(si)/MW*GMat(GAS(5),DI(fi1),DI(fi2));
1013 } else if(is_uqbar_dq(pi1, pi2) && pi3==phip) {
1014 // uqbar-dq-phip
1015 auto ckm = CKM(si1,si2);
1016 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));
1017 } else if(is_dqbar_uq(pi1, pi2) && pi3==phim) {
1018 // dqbar-uq-phim
1019 auto ckm = CKM(si1,si2);
1020 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));
1021 } else if(is_nbar_l(pi1, pi2) && pi3==phip) {
1022 // nbar-l-phip
1023 return -I*GEW/sqrt2*M(si2)/MW*GMat(GAS(1)+GAS(5),DI(fi1),DI(fi2))/2; // Mnl = 0
1024 } else if(is_lbar_n(pi1, pi2) && pi3==phim) {
1025 // lbar-n-phim
1026 return -I*GEW/sqrt2*M(si1)/MW*GMat(GAS(1)-GAS(5),DI(fi1),DI(fi2))/2; // Mnl = 0
1027 } else if(pi1==A && pi2==phip && pi3==phim) {
1028 // A-phip-phim
1029 return -I*eta_e*EL*SP(mom2-mom3,LI(fi1));
1030 } else if(pi1==Z && pi2==phip && pi3==phim) {
1031 // Z-phip-phim
1032 return -I*eta*eta_Z*GEW*C2W/(2*CW)*SP(mom2-mom3,LI(fi1));
1033 } else if(pi1==Wp && pi2==phim && pi3==H) {
1034 // Wp-phim-H
1035 return I/2*eta*GEW*SP(mom2-mom3,LI(fi1));
1036 } else if(pi1==Wm && pi2==phip && pi3==H) {
1037 // Wm-phip-H
1038 return -I/2*eta*GEW*SP(mom2-mom3,LI(fi1));
1039 } else if( ((pi1==Wp && pi2==phim) || (pi1==Wm && pi2==phip)) && pi3==chi) {
1040 // Wp-phim-chi & Wm-phip-chi
1041 return -eta*GEW/2*SP(mom2-mom3,LI(fi1));
1042 } else if(pi1==Z && pi2==chi && pi3==H) {
1043 // Z-chi-H
1044 return -eta*eta_Z*GEW/(2*CW)*SP(mom2-mom3,LI(fi1));
1045 } else if( ((pi1==phim && pi2==Wp) || (pi1==phip && pi2==Wm)) && pi3==A) {
1046 // phim-Wp-A & phip-Wm-A
1047 return I*eta_e*eta*EL*MW*SP(LI(fi2),LI(fi3));
1048 } else if( ((pi1==phim && pi2==Wp) || (pi1==phip && pi2==Wm)) && pi3==Z) {
1049 // phim-Wp-Z & phip-Wm-Z
1050 return -I*eta_Z*GEW*MZ*SW2*SP(LI(fi2),LI(fi3));
1051 } else if(pi1==H && pi2==Wp && pi3==Wm) {
1052 // H-Wp-Wm
1053 return I*GEW*MW*SP(LI(fi2),LI(fi3));
1054 } else if(pi1==H && pi2==Z && pi3==Z) {
1055 // H-Z-Z
1056 return I*GEW/CW*MZ*SP(LI(fi2),LI(fi3));
1057 } else if(pi1==H && pi2==phim && pi3==phip) {
1058 // H-phim-phip
1059 return -I/2*GEW*MH2/MW;
1060 } else if(pi1==H && pi2==H && pi3==H) {
1061 // H-H-H
1062 return -3/ex(2)*I*GEW*MH2/MW;
1063 } else if(pi1==H && pi2==chi && pi3==chi) {
1064 // H-chi-chi
1065 return -I/2*GEW*MH2/MW;
1066 } else if(pi1==ghWpbar && pi2==ghWp && pi3==A) {
1067 // ghWpbar-ghWp-A
1068 return I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1069 } else if(pi1==ghWmbar && pi2==ghWm && pi3==A) {
1070 // ghWmbar-ghWm-A
1071 return -I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1072 } else if(pi1==ghWpbar && pi2==ghWp && pi3==Z) {
1073 // ghWpbar-ghWp-Z
1074 return I*eta_G*eta*eta_Z*GEW*CW*SP(mom1,LI(fi3));
1075 } else if(pi1==ghWmbar && pi2==ghWm && pi3==Z) {
1076 // ghWmbar-ghWm-Z
1077 return -I*eta_G*eta*eta_Z*GEW*CW*SP(mom1,LI(fi3));
1078 } else if(pi1==ghWpbar && pi2==ghZ && pi3==Wp) {
1079 // ghWpbar-ghZ-Wp
1080 return -I*eta_G*eta*eta_Z*GEW*CW*SP(mom1,LI(fi3));
1081 } else if(pi1==ghWmbar && pi2==ghZ && pi3==Wm) {
1082 // ghWmbar-ghZ-Wm
1083 return I*eta_G*eta*eta_Z*GEW*CW*SP(mom1,LI(fi3));
1084 } else if(pi1==ghWpbar && pi2==ghA && pi3==Wp) {
1085 // ghWpbar-ghA-Wp
1086 return -I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1087 } else if(pi1==ghWmbar && pi2==ghA && pi3==Wm) {
1088 // ghWmbar-ghA-Wm
1089 return I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1090 } else if(pi1==ghZbar && pi2==ghWp && pi3==Wm) {
1091 // ghZbar-ghWp-Wm
1092 return -I*eta_G*eta*GEW*CW*SP(mom1,LI(fi3));
1093 } else if(pi1==ghZbar && pi2==ghWm && pi3==Wp) {
1094 // ghZbar-ghWm-Wp
1095 return I*eta_G*eta*GEW*CW*SP(mom1,LI(fi3));
1096 } else if(pi1==ghAbar && pi2==ghWp && pi3==Wm) {
1097 // ghAbar-ghWp-Wm
1098 return -I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1099 } else if(pi1==ghAbar && pi2==ghWm && pi3==Wp) {
1100 // ghAbar-ghWm-Wp
1101 return I*eta_G*eta_e*EL*SP(mom1,LI(fi3));
1102 } else if(pi1==ghWpbar && pi2==ghWp && pi3==chi) {
1103 // ghWpbar-ghWp-chi
1104 return eta_G*GEW/2*xi*MW;
1105 } else if(pi1==ghWmbar && pi2==ghWm && pi3==chi) {
1106 // ghWmbar-ghWm-chi
1107 return -eta_G*GEW/2*xi*MW;
1108 } else if( ((pi1==ghWpbar && pi2==ghWp) ||(pi1==ghWmbar && pi2==ghWm)) && pi3==H) {
1109 // ghWpbar-ghWp-H & ghWmbar-ghWm-H
1110 return -I/2*eta_G*GEW*xi*MW;
1111 } else if(pi1==ghZbar && pi2==ghZ && pi3==H) {
1112 // ghZbar-ghZ-H
1113 return -eta_G*I*GEW/(2*CW)*xi*MZ;
1114 } else if( pi1==ghZbar && ((pi2==ghWp && pi3==phim) || (pi2==ghWm && pi3==phip)) ) {
1115 // ghZbar-ghWp-phim & ghZbar-ghWm-phip
1116 return I/2*eta_G*eta_Z*GEW*xi*MZ;
1117 } else if( ((pi1==ghWpbar && pi3==phip) || (pi1==ghWmbar && pi3==phim)) && pi2==ghZ) {
1118 // ghWpbar-ghZ-phip & ghWmbar-ghZ-phim
1119 return -I*eta_G*eta_Z*GEW*C2W/(2*CW)*xi*MW;
1120 } else if( ((pi1==ghWpbar && pi3==phip) || (pi1==ghWmbar && pi3==phim)) && pi2==ghA) {
1121 // ghWpbar-ghA-phip & ghWmbar-ghA-phim
1122 return -I*eta_G*eta_e*eta*EL*xi*MW;
1123 } else {
1124 cout << endl << e << endl;
1125 throw Error("Vertex Not defined!");
1126 }
1127 } else if(nn==4) {
1128 auto pi4 = e.op(3).op(0);
1129 auto fi4 = e.op(3).op(1);
1130 auto mom4 = e.op(3).op(2);
1131
1132 if(pi1==g && pi2==g && pi3==g && pi4==g) {
1133 // g^4
1134 return Gluon4Vertex(e, eta_s);
1135 } else if(pi1==Wp && pi2==Wm && pi3==A && pi4==A) {
1136 // Wp-Wm-A-A
1137 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)));
1138 } else if(pi1==Wp && pi2==Wm && pi3==Z && pi4==Z) {
1139 // Wp-Wm-Z-Z
1140 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)));
1141 } else if(pi1==Wp && pi2==Wm && pi3==A && pi4==Z) {
1142 // Wp-Wm-A-Z
1143 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)));
1144 } else if(pi1==Wp && pi2==Wp && pi3==Wm && pi4==Wm) {
1145 // Wp-Wp-Wm-Wm
1146 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)));
1147 } else if(pi1==Wp && pi2==Wm && pi3==H && pi4==H) {
1148 // Wp-Wm-H-H
1149 return I/2*GEW2*SP(LI(fi1),LI(fi2));
1150 } else if(pi1==Wp && pi2==Wm && pi3==chi && pi4==chi) {
1151 // Wp-Wm-chi-chi
1152 return I/2*GEW2*SP(LI(fi1),LI(fi2));
1153 } else if(pi1==Z && pi2==Z && pi3==H && pi4==H) {
1154 // Z-Z-H-H
1155 return I/2*GEW2/CW2*SP(LI(fi1),LI(fi2));
1156 } else if(pi1==Z && pi2==Z && pi3==chi && pi4==chi) {
1157 // Z-Z-chi-chi
1158 return I/2*GEW2/CW2*SP(LI(fi1),LI(fi2));
1159 } else if(pi1==A && pi2==A && pi3==phip && pi4==phim) {
1160 // A-A-phip-phim
1161 return 2*I*EL2*SP(LI(fi1),LI(fi2));
1162 } else if(pi1==Z && pi2==Z && pi3==phip && pi4==phim) {
1163 // Z-Z-phip-phim
1164 return I/2*pow(GEW*C2W/CW,2)*SP(LI(fi1),LI(fi2));
1165 } else if(pi1==Wp && pi2==Wm && pi3==phip && pi4==phim) {
1166 // Wp-Wm-phip-phim
1167 return I/2*GEW2*SP(LI(fi1),LI(fi2));
1168 } else if( ((pi1==Wp && pi3==phim) || (pi1==Wm && pi3==phip)) && pi2==Z && pi4==H) {
1169 // Wp-Z-phim-H & Wm-Z-phip-H
1170 return -I*eta_Z*GEW2*SW2/(2*CW)*SP(LI(fi1),LI(fi2));
1171 } else if(pi1==Wm && pi2==Z && pi3==phip && pi4==chi) {
1172 // Wm-Z-phip-chi
1173 return -eta_Z*GEW2*SW2/(2*CW)*SP(LI(fi1),LI(fi2));
1174 } else if(pi1==Wp && pi2==Z && pi3==phim && pi4==chi) {
1175 // Wp-Z-phim-chi
1176 return eta_Z*GEW2*SW2/(2*CW)*SP(LI(fi1),LI(fi2));
1177 } else if( ((pi1==Wm && pi3==phip) || (pi1==Wp && pi3==phim)) && pi2==A && pi4==H) {
1178 // Wm-A-phip-H & Wp-A-phim-H
1179 return I/2*eta_e*eta*EL*GEW*SP(LI(fi1),LI(fi2));
1180 } else if(pi1==Wp && pi2==A && pi3==phim && pi4==chi) {
1181 // Wp-A-phim-chi
1182 return -1/ex(2)*eta_e*eta*EL*GEW*SP(LI(fi1),LI(fi2));
1183 } else if(pi1==Wm && pi2==A && pi3==phip && pi4==chi) {
1184 // Wm-A-phip-chi
1185 return 1/ex(2)*eta_e*eta*EL*GEW*SP(LI(fi1),LI(fi2));
1186 } else if(pi1==Z && pi2==A && pi3==phip && pi4==phim) {
1187 // Z-A-phip-phim
1188 return I*eta_e*eta*eta_Z*EL*GEW*C2W/CW*SP(LI(fi1),LI(fi2));
1189 } else if(pi1==phim && pi2==phip && pi3==phim && pi4==phip) {
1190 // phim-phip-phim-phip
1191 return -I/2*GEW2*MH2/MW2;
1192 } else if(pi1==H && pi2==H && pi3==phim && pi4==phip) {
1193 // H-H-phim-phip
1194 return -I/4*GEW2*MH2/MW2;
1195 } else if(pi1==chi && pi2==chi && pi3==phim && pi4==phip) {
1196 // chi-chi-phim-phip
1197 return -I/4*GEW2*MH2/MW2;
1198 } else if(pi1==H && pi2==H && pi3==H && pi4==H) {
1199 // H-H-H-H
1200 return -3/ex(4)*I*GEW2*MH2/MW2;
1201 } else if(pi1==H && pi2==H && pi3==chi && pi4==chi) {
1202 // H-H-chi-chi
1203 return -I/4*GEW2*MH2/MW2;
1204 } else if(pi1==chi && pi2==chi && pi3==chi && pi4==chi) {
1205 // chi-chi-chi-chi
1206 return -3/ex(4)*I*GEW2*MH2/MW2;
1207 } else {
1208 cout << endl << e << endl;
1209 throw Error("Vertex Not defined!");
1210 }
1211 } else {
1212 cout << endl << e << endl;
1213 throw Error("Vertex Not defined!");
1214 }
1215 } else return e.map(self);
1216 return e;
1217 });
1218 return fr(amp);
1219 }
1220
1221
1222 lst Models::FeynRulesQCD(const lst & amps, const ex & xi) {
1223 lst ret;
1224 for(auto item : amps) ret.append(FeynRulesQCD(item,xi));
1225 return ret;
1226 }
1227 ex Models::FeynRulesQCD(const ex & amp, const ex & xi) {
1228 if(is_a<lst>(amp)) return FeynRulesQCD(ex_to<lst>(amp), xi);
1229
1230 static Symbol A("A"), Q("Q"), Qbar("Qbar"), q("q"), qbar("qbar"), g("g"), gh("gh"), ghbar("ghbar");
1231 static Symbol m("m"), eq("eq"), eQ("eQ");
1232
1233 auto fr = MapFunction([&](const ex &e, MapFunction &self)->ex {
1234 if(isFunction(e,"OutField") || isFunction(e,"InField")) return 1;
1235 else if(isFunction(e, "Propagator")) {
1236 if(e.op(0).op(0)==q) {
1237 return QuarkPropagator(e, 0);
1238 } else if(e.op(0).op(0)==Q) {
1239 return QuarkPropagator(e, m);
1240 } else if(e.op(0).op(0)==g) {
1241 return GluonPropagator(e, xi);
1242 } else if(e.op(0).op(0)==gh) {
1243 return GluonGhostPropagator(e, xi);
1244 } else {
1245 cout << "expr: " << e << endl;
1246 throw Error("Propagator Not Found!");
1247 }
1248 } else if(isFunction(e, "Vertex")) {
1249 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) ) {
1250 // qbar-q-g or Qbar-Q-g or g -> A
1251 if(e.op(2).op(0)==g) return QuarkGluonVertex(e);
1252 else {
1253 auto fi1 = e.op(0).op(1);
1254 auto fi2 = e.op(1).op(1);
1255 auto fi3 = e.op(2).op(1);
1256 if(e.op(1).op(0)==q) QuarkAVertex(e, eq);
1257 else if(e.op(1).op(0)==Q) QuarkAVertex(e, eQ);
1258 else throw Error("Vertex Error.");
1259 }
1260 } else if(e.nops()==3 && e.op(0).op(0)==ghbar && e.op(1).op(0)==gh) {
1261 // ghbar-gh-g
1262 return GhostGluonVertex(e);
1263 } else if(e.nops()==3 && e.op(0).op(0)==g && e.op(1).op(0)==g && e.op(2).op(0)==g) {
1264 // g^3
1265 return Gluon3Vertex(e);
1266 } 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) {
1267 // g^4
1268 return Gluon4Vertex(e);
1269 } else {
1270 cout << "expr: " << e << endl;
1271 throw Error("Vertex Not Found!");
1272 }
1273 } else return e.map(self);
1274 return e;
1275 });
1276 return fr(amp);
1277 }
1278
1279}
1280
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:681
HepLib::QGRAF::HasLoop
bool HasLoop(ex amp, lst prop)
check a amplitude has a loop w.r.t. propapagtor
Definition QGRAF.cpp:378
HepLib::QGRAF::AntiGhostSumL
ex AntiGhostSumL(int qi)
polarization sum for anti-ghost, note that we will add extra - here
Definition QGRAF.cpp:722
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:603
HepLib::QGRAF::GluonPropagator
ex GluonPropagator(const ex &e, const ex &xi)
Propagator for gluon.
Definition QGRAF.cpp:441
HepLib::QGRAF::LeptonAVertex
ex LeptonAVertex(const ex &e, const ex &eq, const ex &eta_e)
l-lbar-A vertex
Definition QGRAF.cpp:590
HepLib::QGRAF::Gluon3Vertex
ex Gluon3Vertex(const ex &e, const ex &eta_s)
g-g-g vertex
Definition QGRAF.cpp:522
HepLib::QGRAF::GhostSumL
ex GhostSumL(int qi)
polarization sum for ghost
Definition QGRAF.cpp:704
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:731
HepLib::QGRAF::QuarkGluonVertex
ex QuarkGluonVertex(const ex &e, const ex &eta_s)
q-qbar-g vertex
Definition QGRAF.cpp:509
HepLib::QGRAF::ScalarPropagator
ex ScalarPropagator(const ex &e, const ex &m, const ex &xi)
Propagator for scalar.
Definition QGRAF.cpp:496
HepLib::QGRAF::J1SumL
ex J1SumL(int qi, ex p)
polarization sum for total angular momentum
Definition QGRAF.cpp:741
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:428
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:290
HepLib::QGRAF::AntiQuarkSumR
ex AntiQuarkSumR(int qi, ex p, ex m, bool color)
polarization sum for anti-quark
Definition QGRAF.cpp:694
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:752
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:482
HepLib::QGRAF::Gluon4Vertex
ex Gluon4Vertex(const ex &e, const ex &eta_s)
g-g-g-g vertex
Definition QGRAF.cpp:542
HepLib::QGRAF::DrawPDF
void DrawPDF(const lst &amps, string fn, int nr)
generate Feynman diagrams for the amplitudes, in PDF format
Definition QGRAF.cpp:168
HepLib::QGRAF::IndexCC
ex IndexCC(ex e, bool all)
Definition QGRAF.cpp:621
HepLib::QGRAF::GhostGluonVertex
ex GhostGluonVertex(const ex &e, const ex &eta_s, const ex &eta_G)
ghbar-gh-g vertex
Definition QGRAF.cpp:561
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:415
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:668
HepLib::QGRAF::QuarkSumL
ex QuarkSumL(int qi, ex p, ex m, bool color)
polarization sum for quark
Definition QGRAF.cpp:655
HepLib::QGRAF::GhostSumR
ex GhostSumR(int qi)
polarization sum for ghost
Definition QGRAF.cpp:713
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:468
HepLib::QGRAF::QuarkAVertex
ex QuarkAVertex(const ex &e, const ex &eq, const ex &eta_e)
qbar-q-A vertex
Definition QGRAF.cpp:576
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:454
HepLib::QGRAF::FI
Index FI(ex fn)
Definition QGRAF.cpp:48
HepLib::QGRAF::RFI
Index RFI(ex fn)
Definition QGRAF.cpp:54
HepLib::QGRAF::TopoLines
lst TopoLines(const ex &amp)
generate the topological lines for the amplitude
Definition QGRAF.cpp:123
HepLib::QGRAF::GluonSumR
ex GluonSumR(int qi, bool color)
polarization sum for gluon
Definition QGRAF.cpp:642
HepLib::QGRAF::GluonSumL
ex GluonSumL(int qi, bool color)
polarization sum for gluon
Definition QGRAF.cpp:631
HepLib::gs
const Symbol gs
HepLib::Combinations
void Combinations(int n, int m, std::function< void(const int *)> f)
Definition Functions.cpp:184
HepLib::ep
const Symbol ep
HepLib::GAS
ex GAS(const ex &expr, unsigned rl)
function similar to GAD/GSD in FeynClac
Definition DGamma.cpp:257
HepLib::Debug
bool Debug
Definition Init.cpp:144
HepLib::w
ex w
Definition Init.cpp:93
HepLib::vs
const Symbol vs
HepLib::isFunction
bool isFunction(const ex &e, string func_name)
Definition BASIC.h:658
HepLib::InstallPrefix
string InstallPrefix
Definition Init.cpp:162
HepLib::string_replace_all
void string_replace_all(string &str, const string &from, const string &to)
Definition Functions.cpp:148
HepLib::ex2str
string ex2str(const ex &expr)
convert ex to output string, the defalut printer format will be used
Definition BASIC.cpp:715
HepLib::dir_exists
bool dir_exists(string dir)
Definition BASIC.h:297
HepLib::GiNaC_Parallel
exvector GiNaC_Parallel(int ntotal, std::function< ex(int)> f, const string &key)
GiNaC Parallel Evaluation using fork.
Definition BASIC.cpp:260
HepLib::REGISTER_FUNCTION
REGISTER_FUNCTION(GMat, do_not_evalf_params().print_func< FCFormat >(&GMat_fc_print).conjugate_func(mat_conj).set_return_type(return_types::commutative)) bool IsZero(const ex &e)
Definition Basic.cpp:1026
HepLib::SP
ex SP(const ex &a, bool use_map=false)
Definition Pair.cpp:166
Generated by doxygen 1.9.8