/* * Copyright (C) 2003-2009, CompHEP Collaboration * Copyright (C) 2003, Slava Bunichev * ------------------------------------------------------ */ #include #include #include "service2/include/chep_limits.h" #include "service2/include/tptcmac.h" #include "service2/include/unix_utils.h" #include "service2/include/read_func.h" #include "service2/include/syst.h" #include "service2/include/s_files.h" #include "service2/include/lbl.h" #include "service2/include/parser.h" #include "service2/include/getmem.h" #include "chep_crt/include/chep_crt.h" #include "physics.h" #include "ghosts.h" #include "cweight.h" #include "prepdiag.h" #include "diaprins.h" #include "chess.h" #include "out_service.h" #include "rfactor.h" #include "denominators.h" #include "process.h" #include "process_core.h" #include "r_code.h" #include "procvar.h" #include "out_c.h" #include "form_code.h" char momsubst[9] = { 0 }; char indsubst[9] = { 0 }; static int r_reading0 = FALSE; static int vertmap[3 * maxvert]; static int gammaflag; static char *parsedverts[2 * maxvert]; static int numberg5; static int indicatorFl[20]; static int FermLoopIndicator; void run_formprograms (void) { int i; DIR *dirp; struct direct *directp; char form_filename[25]; char cmd_string[50]; dirp = opendir ("results/"); while ((directp = readdir (dirp)) != NULL) { sscanf (directp->d_name, "%s", form_filename); for (i = 0; i < 25; i++) if (form_filename[i] == '.') { if ((form_filename[i + 1] == 'f') && (form_filename[i + 2] == 'r') && (form_filename[i + 3] == 'm')) { sprintf (cmd_string, "./form ./results/%s\n", form_filename); system (cmd_string); } } } (void) closedir (dirp); } static void * bact0 (char ch, pointer mm1, pointer mm2) { char *m1, *m2, *ans; int sgn; if (r_reading0 && (ch == '*')) { m1 = (char *) mm2; m2 = (char *) mm1; } else { m1 = (char *) mm1; m2 = (char *) mm2; } if (ch == '+' || ch == '-') { lShift (m1, 2); lShift (m2, 2); } else { if (m1[0] == 'P' || ch == '^') { lShift (m1, 1); m1[0] = '('; strcat (m1, ")"); } else lShift (m1, 2); if (m2[0] == 'P' || ch == '^') { lShift (m2, 1); m2[0] = '('; strcat (m2, ")"); } else lShift (m2, 2); } ans = (char *) m_alloc (strlen (m1) + strlen (m2) + 8); switch (ch) { case '+': if (m2[0] == '-') sprintf (ans, "P|%s%s", m1, m2); else sprintf (ans, "P|%s+%s", m1, m2); break; case '-': if (m2[0] == '-') sprintf (ans, "P|%s+%s", m1, m2 + 1); else sprintf (ans, "P|%s-%s", m1, m2); break; case '*': sgn = 1; if (m1[0] == '-') { lShift (m1, 1); sgn = -sgn; } if (m2[0] == '-') { lShift (m2, 1); sgn = -sgn; } if (sgn == 1) sprintf (ans, "M|%s*%s", m1, m2); else sprintf (ans, "M|-%s*%s", m1, m2); break; case '.': if (m2[0] != '-') sprintf (ans, "M|F(%s,%s)", m1, m2); else if (m1[0] != '-') sprintf (ans, "M|F(%s,%s)", m2, m1); else { lShift (m1, 1); lShift (m2, 1); sprintf (ans, "M|F(%s,%s)", m1, m2); } break; case '^': sprintf (ans, "M|%s^%s", m1, m2); } return (pointer) ans; } static void * uact0 (char *ch, pointer mm) { char *m, *ans; m = (char *) mm; ans = (char *) m_alloc (strlen (m) + 10); /* if (strcmp(ch,"-") == 0)*/ if (m[0] == 'M') { if (m[2] == '-') sprintf (ans, "M|%s", m + 3); else sprintf (ans, "M|-%s", m + 2); } else sprintf (ans, "M|-(%s)", m + 2); if (strcmp (ch, "G") == 0) { if (r_reading0) sprintf (ans, "M|-g_(ln,%s)", m + 2); else sprintf (ans, "M|g_(ln,%s)", m + 2); } return (pointer) ans; } static void * act_rcode (char *ch, int n, void **args) { if (n == 1) return uact0 (ch, args[0]); if (n == 2) return bact0 (ch[0], args[0], args[1]); if (n == 4 && !strcmp (ch, "eps")) { int l = 15 + strlen (args[0]) + strlen (args[1]) + strlen (args[2]) + strlen (args[3]); char *ans = (char *) m_alloc (l); sprintf (ans, "M|eps(%s,%s,%s,%s)", (char *) args[0] + 2, (char *) args[1] + 2, /* eps(a,b,c,d) */ (char *) args[2] + 2, (char *) args[3] + 2); return ans; } return 0; } static void * rd_rcode (char *s) { char *p; /*char * GammaString; */ int num; /*GammaString = (char *) m_alloc(12); */ p = (char *) m_alloc (12); p[0] = 0; if (strlen (s) == 2 && s[1] > '0' && s[1] <= '9') { switch (s[0]) { case 'p': case 'P': num = s[1] - '0'; num = momsubst[num - 1]; if (num > 0) sprintf (p, "M|p%d", num); else sprintf (p, "M|-p%d", -num); break; case 'm': num = s[1] - '0'; num = indsubst[num - 1]; sprintf (p, "M|m%d", num); break; case 'M': num = s[1] - '0'; num = indsubst[num - 1] - 1; sprintf (p, "M|m%d", num); } if (strcmp (s, "G5") == 0) strcpy (p, "M|g_(ln,5_)"); } if (!strlen (p)) sprintf (p, "M|%s", s); return (void *) p; } static char * pexpr (char p) { static char snum[8]; if (p >= 0) sprintf (snum, "p%d", p); else sprintf (snum, "(-p%d)", -p); return snum; } static char * iexpr (int i) { static char snum[6]; sbld (snum, "m%d", i); return snum; } static void var_def (void) { int k, i = 1; writeF ("#define dp(n) DP[(n)]\n"); for (k = 1; k <= nmodelvar; k++) { varlist modl = modelvars + k; if (!modl->func) { writeF ("#define %s va[%d]\n", modl->varname, i); i++; } } for (k = 1; k <= nmodelvar; k++) { varlist modl = modelvars + k; if (modl->func) { writeF ("#define %s va[%d]\n", modl->varname, i); i++; } } } static void head (void) { writeF ("*----------- VARIABLES ------------------ \n"); writeF (" S x,xx,xxx,y,mas,MZ,MW,Mt,Mb,Mc,Mu,Md,Ms,Mtop,Me,MH,"); writeF (" Mm,SW,EE,GG,Cw,CW,"); writeF (" Vus,Vcd,Vtb,Vcb,Vud,Vub,Vts,Vtd,Vcs,wtop,wZ,wW,wH,"); writeF (" s12,s23,s13,c12,c23,c13,Sqrt2,width,label;\n"); writeF (" S d0,...,d27;\n"); writeF (" Auto S cc;\n"); writeF (" V p1,...,p16,ZERO_,z,mom,moment;\n"); writeF (" I m1,...,m16,s,t,ind,j,ln;\n"); writeF (" CFunction dp,ANT,H;\n"); writeF (" Function F,propDen;\n"); writeF ("*\n"); writeF ("Off statistics;\n\n"); writeF ("#call Optimal()\n"); } static void writesubst (void) { int i, j, l; int c; int ntot = getntot (); writeF ("*--- Mass declarations and Moment substitutions ---\n\n"); writeF ("#procedure Substitution()\n"); /* writeF(" .sort\n"); */ writeF (" repeat;\n"); for (l = 1; l <= ntot; l++) writeF (" id p%d.p%d = %s*%s;\n", l, l, inoutmasses[l - 1], inoutmasses[l - 1]); writeF ("\n"); for (i = 1; i <= 3 * maxvert; i++) { l = momdep[i - 1][0]; /* For momdep this is a length V.E. */ if (l > 1) { writeF (" id p%d = ", i); for (j = 1; j <= l; j++) { c = /* (char) */ momdep[i - 1][j]; if (c > 0) writeF ("+p%d", c); else writeF ("-p%d", -c); } writeF (";\n"); } } writeF (" endrepeat;\n"); writeF ("#endprocedure\n\n\n"); } static int MaxDotProd (int i) { if (i == 2) return 0; if (i == 3) return 1; if (i == 4) return 4; if (i == 5) return 8; if (i == 6) return 13; if (i == 7) return 19; if (i == 8) return 26; return 0; } static void optimal (void) { int i, j, l, max, sign; int nin = getnin (); writeF ("#procedure Optimal()\n"); writeF (" .sort\n"); for (i = 1; i <= 3 * maxvert; i++) { l = momdep[i - 1][0]; if (l > 1) { max = i; i = 3 * maxvert + 1; } } writeF (" L PmaxPmax="); sign = 1; if ((max - 2) > nin) sign *= -1; if ((max - 1) > nin) sign *= -1; if (sign < 0) writeF ("-1*("); else writeF ("("); for (i = 1; i < max; i++) for (j = i; j < max; j++) if ((i != max - 2) || (j != max - 1)) { sign = 1; if (i > nin) sign *= -1; if (j > nin) sign *= -1; if (sign < 0) writeF ("+"); else writeF ("-"); writeF ("p%d.p%d", i, j); if (i == j) writeF ("/2"); } writeF ("+p%d.p%d/2);\n", max, max); writeF (" .sort\n"); writeF (" skip; nskip PmaxPmax;\n"); writeF (" #call Substitution()\n"); writeF (" .sort\n"); /* writeF(" print PmaxPmax;\n"); */ writeF (" #$MaxDpMom= p%d.p%d;\n", max - 2, max - 1); writeF (" #$MaxDpSub= PmaxPmax;\n"); writeF (" #$MaxDpInd= %d;\n", MaxDotProd (max - 1)); writeF (" .sort\n"); writeF (" Drop PmaxPmax;\n"); writeF (" .sort\n"); writeF ("#endprocedure\n\n\n"); } static void emitfactors (void) { char num[STRSIZ]; int c; int i, j, s, vln, pnum; writeF ("*---------- Factors ---------------\n"); /* ------- Symmetry Factor ------- */ writeF (" L SymmFact=%d/%d; * Diagram symmerty factor\n", vcs.symnum, vcs.symdenum); /* ----- average factor --------- */ c = 1; for (i = 1; i <= vcs.sizel; i++) { vln = vcs.valence[i - 1]; for (j = 1; j <= vln; j++) if (IN_PRTCL & vcs.vertlist[i - 1][j - 1].prop) { pnum = vcs.vertlist[i - 1][j - 1].partcl; s = prtclbase[pnum - 1].spin; switch (s) { case 1: if (strchr ("LR", prtclbase[pnum - 1].hlp) == NULL) c *= 2; break; case 2: if (zeromass (pnum)) c *= 2; else c *= 3; } c *= abs (prtclbase[pnum - 1].cdim); } } writeF (" L AverFact=1/%d; * Normalization factor of polarization average\n", c); /* ----- Fermion factor --------- */ c = 0; for (i = 0; i < vcs.sizel; i++) for (j = 0; j < vcs.valence[i]; j++) { pnum = vcs.vertlist[i][j].partcl; if (prtclbase[pnum - 1].spin == 1 && (IN_PRTCL & vcs.vertlist[i][j].prop)) ++(c); } if ((c & 1) == 1) strcpy (num, "-1"); else strcpy (num, "1"); writeF (" L FermFact=%s; * (-1)**(number of in-fermion particles)\n", num); /* ----- Vector factor --------- */ /* c = 0; for (i = 1; i <= vcs.sizet; i++) { vln = vcs.valence[i-1]; for (j = 1; j <= vln; j++) if (i < vcs.vertlist[i-1][j-1].nextvert.vno) { pnum = vcs.vertlist[i-1][j-1].partcl; if (prtclbase[pnum-1].spin == 2) ++(c); } } if ((c & 1) == 1) strcpy(num,"-1"); else strcpy(num,"1"); writeF(" Local VectFact=%s; * (-1)**(number of vector propagators)\n",num); */ /* ----- Color factor --------- */ writeF (" L ColorFact=%d/%d; * QCD color weight of diagram \n", vcs.clrnum, vcs.clrdenum); writeF ("*\n"); } /*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | Program for emitting denominators of propagators in CrossSection | | diagramm. | | September 19, 1989 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ static void emitdenoms (void) { int j; /* 1..2 * maxvert */ int k; /* 1..4 */ int denrno; /* EmitDenoms -- main */ for (j = 1; j <= vcs.sizet; j++) for (k = 1; k <= vcs.valence[j - 1]; k++) { edgeinvert *ln = &vcs.vertlist[j - 1][k - 1]; if (ln->moment > 0 && !(ln->prop & (IN_PRTCL | OUT_PRTCL)) && pseudop (ln->partcl)) writeF ("L totFactor=totFactor/%s**2;\n", prtclbase[ln->partcl - 1].massidnt); } writeF (" L denominator="); denrno = 0; for (j = 1; j <= vcs.sizet; j++) for (k = 1; k <= vcs.valence[j - 1]; k++) { edgeinvert *ln = &vcs.vertlist[j - 1][k - 1]; if (ln->moment > 0 && !(ln->prop & (IN_PRTCL | OUT_PRTCL))) { if (!pseudop (ln->partcl)) { char width[10] = "0"; if (denrno) writeF ("*"); denrno++; if (!ttypepropag (j, k)) strcpy (width, prtclbase[ln->partcl - 1].imassidnt); writeF ("propDen(p%d,%s,%s)", ln->moment, prtclbase[ln->partcl - 1].massidnt, width); } } } if (!denrno) writeF ("1"); writeF (";\n"); } /* EmitDenoms */ /* Programm for constructing vertex-oriented representation for ! ! amplitude diagram. ! ! A.Taranov 20.07.89 */ /* static void emitindex(indvertset index) { int m=1, n=0, l=1; if (setofb_eq0(index)) return; writeF(" Index "); do { if (insetb(m,index)) { if(n) writeF(","); n++; writeF("m%d",l); setofb_cpy(index,setofb_aun(index,setofb(m,_E))); } l++; if (!setofb_eq0(index)) m++; } while (!(n == 1 || setofb_eq0(index))); writeF("$\n"); } */ static void reducemult (char *nameres, char *name1, char *name2, indvertset index) { indvertset index1, index2; lvcpy (index1, index); lvcpy (index2, index); /* emitindex(index1); */ writeF (" L %s=%s*%s;\n", nameres, name1, name2); /* while (!setofb_eq0(index1)) { emitindex(index1); writeF(" %s:=%s$\n",nameres,nameres); } if (!setofb_eq0(index2)) { writeF(" RemInd "); l = 1; n = 1; while (!setofb_eq0(index2)) { if (insetb(l,index2)) { if (n != 1) writeF(","); ++(n); writeF("m%d",l); } setofb_cpy(index2,setofb_aun(index2,setofb(l,_E))); ++(l); } writeF("$\n"); } */ } static int fermvrt (int v) { int i; for (i = 0; i < vcs.valence[v - 1]; i++) if (fermionp (vcs.vertlist[v - 1][i].partcl)) return 1; return 0; } static char * subst (int ind, int arg) { int v, l, i; static shortstr ans; v = massindpos[ind - 1].vrt1; if (vertmap[v - 1] == arg) l = massindpos[ind - 1].ln1; else { v = massindpos[ind - 1].vrt2; l = massindpos[ind - 1].ln2; } sbld (ans, "%s,(", iexpr (ind)); if (gammaflag && fermvrt (v)) { for (i = 1; i <= vcs.valence[v - 1]; i++) if (i != l) sbld (ans, "%s+%s", ans, pexpr (vcs.vertlist[v - 1][i - 1].moment)); } else sbld (ans, "%s%s", ans, pexpr (-vcs.vertlist[v - 1][l - 1].moment)); sbld (ans, "%s)/%s", ans, prtclbase[vcs.vertlist[v - 1][l - 1].partcl - 1].massidnt); return ans; } static void r_mult (int arg1, int arg2, indvertset index1) { char name0[8], name1[8], name2[8]; indvertset mind; unsigned m, mm, maxmm; int pstn[15]; int n, l; indvertset index; /* R_mult */ setofb_cpy (index, index1); sbld (name1, "Vrt%d", arg1); sbld (name2, "Vrt%d", arg2); if (arg1 < arg2) strcpy (name0, name1); else strcpy (name0, name2); setofb_cpy (mind, setofb_its (index, setmassindex)); if (setofb_eq0 (mind)) { reducemult (name0, name1, name2, index); writeF (" .sort\n"); } else { reducemult ("Vrt0", name1, name2, index); l = 1; n = 0; while (!setofb_eq0 (mind)) { if (insetb (l, mind)) pstn[++n - 1] = l; setofb_cpy (mind, setofb_aun (mind, setofb (l, _E))); ++(l); } maxmm = (1 << n) - 1; for (mm = 1; mm <= maxmm; mm++) { writeF (" L VrtL=%s; L VrtR=%s;\n", name1, name2); setofb_zero (mind); m = mm; n = 1; while (m != 0) { if ((m & 1) != 0 && (l = pstn[n - 1]) != 0) { /* writeF(" VrtL=(VrtL where %s);\n",subst(l,arg1)); writeF(" VrtR=(VrtR where %s);\n",subst(l,arg2)); */ writeF (" #call RepIND(L,%s)\n", subst (l, arg1)); writeF (" #call RepIND(R,%s)\n", subst (l, arg2)); setofb_cpy (mind, setofb_uni (mind, setofb (l, _E))); } ++(n); m >>= 1; } reducemult ("Vrt0", "Vrt0 + VrtL", "VrtR", setofb_aun (index, mind)); writeF (" .sort\n"); /* writeF( " Clear Vrt_L $ Clear Vrt_R $\n"); */ } writeF (" L %s=Vrt0;\n", name0); writeF (" .sort\n"); } if (arg1 < arg2) { /* writeF(" Clear %s$\n",name2); */ for (n = 1; n <= vcs.sizet; n++) if (vertmap[n - 1] == arg2) vertmap[n - 1] = arg1; } if (arg2 < arg1) { /* writeF(" Clear %s$\n",name1); */ for (n = 1; n <= vcs.sizet; n++) if (vertmap[n - 1] == arg1) vertmap[n - 1] = arg2; } } /* R_mult */ /* Emit Reduce prg for ferm loop calculation */ static char * fermpropag (int v, int l) { static shortstr ans, mass; /*char * GammaString; */ edgeinvert *ln = &vcs.vertlist[v - 1][l - 1]; /*GammaString = (char *) m_alloc(18); */ strcpy (mass, prtclbase[ln->partcl - 1].massidnt); if (prtclbase[ln->partcl - 1].hlp == '*') return mass; sprintf (ans, "g_(ln,%s)", pexpr (ln->moment)); if (strcmp (mass, "0") != 0) sbld (ans, "(%s+%s)", ans, mass); else { if (fermionp (ln->partcl)) { if (prtclbase[ln->partcl - 1].hlp == 'L') strcat (ans, "*(1-g_(ln,5_))/2"); /*{ sprintf(GammaString,"*(1-g_(%d,5_))/2",FermLoopIndicator); strcat(ans,GammaString); } */ else if (prtclbase[ln->partcl - 1].hlp == 'R') strcat (ans, "*(1+g_(ln,5_))/2"); /*{ sprintf(GammaString,"*(1+g_(%d,5_))/2",FermLoopIndicator); strcat(ans,GammaString); } */ } else { if (prtclbase[ln->partcl - 1].hlp == 'R') strcat (ans, "*(1-g_(ln,5_))/2"); /*{ sprintf(GammaString,"*(1-g_(%d,5_))/2",FermLoopIndicator); strcat(ans,GammaString); } */ else if (prtclbase[ln->partcl - 1].hlp == 'L') strcat (ans, "*(1+g_(ln,5_))/2"); /*{ sprintf(GammaString,"*(1+g_(%d,5_))/2",FermLoopIndicator); strcat(ans,GammaString); } */ } } return ans; } static void fermprgemit (void) { int v, v1, l; int i, k, lpcount, nn; indvertset indexs; char *vertStr, *propStr; char *GammaString; GammaString = (char *) m_alloc (12); /* Nested function: fermpropag */ /* Main Procedure -- FermPrgEmit */ numberg5 = 0; if (nloop == 0) { gammaflag = FALSE; return; } gammaflag = TRUE; for (v = 1; v <= nloop; v++) { l = 1; fermloops[v - 1].g5 = fermloops[v - 1].lprtcl; while (!fermloops[v - 1].g5 && l <= fermloops[v - 1].len) { sprintf (GammaString, "g_(ln,5_)"); fermloops[v - 1].g5 = spos (GammaString, parsedverts[fermloops[v - 1].vv[l - 1] - 1]) != 0; ++(l); } if (fermloops[v - 1].g5) ++(numberg5); } for (v = 1; v <= 2 * maxvert; v++) vertmap[v - 1] = 0; for (lpcount = 1; lpcount <= nloop; lpcount++) { writeF ("*------- Fermion loop calculation ------- \n"); writeF (" \n"); FermLoopIndicator++; /* FermLoopIndicator */ for (nn = 1; nn <= fermloops[lpcount - 1].len; nn++) { if (vertexes[fermloops[lpcount - 1].vv[nn - 1] - 1].r_vert) writeF ("* reversed vertex\n"); vertStr = parsedverts[fermloops[lpcount - 1].vv[nn - 1] - 1]; propStr = fermpropag (fermloops[lpcount - 1].vv[nn - 1], fermloops[lpcount - 1].ll[nn - 1]); writeF (" L Vrt%d=%s*%s;\n", nn + lpcount - 1, vertStr, propStr); vertmap[fermloops[lpcount - 1].vv[nn - 1] - 1] = nn + lpcount - 1; } /*====================== if (fermloops[lpcount-1].intln[fermloops[lpcount-1].len-1] == 0) setofb_zero(indexs); else setofb_cpy (indexs,setofb(vcs.vertlist[fermloops[lpcount-1].vv[fermloops[lpcount-1].len-1]-1] [fermloops[lpcount-1]. intln[fermloops[lpcount-1].len-1]-1].lorentz ,_E)); ========================*/ writeF (" .sort\n"); for (v1 = 2; v1 <= fermloops[lpcount - 1].len; v1++) { int newInd = fermloops[lpcount - 1].intln[v1 - 1]; /* vector line */ if (newInd != 0) newInd = vcs.vertlist[fermloops[lpcount - 1].vv[v1 - 1] - 1][newInd - 1].lorentz; if (newInd == 0) setofb_zero (indexs); else setofb_cpy (indexs, setofb (newInd, _E)); newInd = fermloops[lpcount - 1].intln2[v1 - 1]; /* second vector line */ if (newInd != 0) newInd = vcs.vertlist[fermloops[lpcount - 1].vv[v1 - 1] - 1][newInd - 1].lorentz; if (newInd != 0) setofb_cpy (indexs, setofb_uni (indexs, setofb (newInd, _E))); r_mult (lpcount, lpcount + v1 - 1, indexs); } if (numberg5 == 1 && fermloops[lpcount - 1].g5) { writeF (" #call G5ravno0(%d)\n", lpcount); /* writeF( " Vrt_%d:=Sub(A=0*A,Vrt_%d)$\n" " Vrt_%d:=(Vrt_%d where A=>0*ZERO_)$\n" ,lpcount,lpcount); */ numberg5 = 0; fermloops[lpcount - 1].g5 = FALSE; } for (k = 1; k <= strlen (fermloops[lpcount - 1].invrt); k++) { v1 = fermloops[lpcount - 1].invrt[k - 1]; setofb_zero (indexs); for (l = 1; l <= vcs.valence[v1 - 1]; l++) { i = vcs.vertlist[v1 - 1][l - 1].lorentz; if (i != 0) setofb_cpy (indexs, setofb_uni (indexs, setofb (i, _E))); } vertmap[v1 - 1] = lpcount + 1; writeF (" Local Vrt%d=%s;\n", lpcount + 1, parsedverts[v1 - 1]); writeF (" #call DefineF()\n"); r_mult (lpcount, lpcount + 1, indexs); } /* writeF(" trace4,1;\n"); writeF(" contract;\n"); writeF(" .sort\n"); */ writeF (" L Fl%d=-Vrt%d;\n", lpcount, lpcount); writeF (" .sort\n"); indicatorFl[lpcount] = 1; writeF (" Drop;"); writeF (" Ndrop"); for (k = 0; k < 20; k++) if (indicatorFl[k] == 1) writeF (" Fl%d,", k); writeF (" numerator, GhostFact, totFactor;\n"); writeF (" .sort\n"); writeF (" #call Substitution()\n"); writeF (" #call RepG(%d,%d)\n", lpcount, FermLoopIndicator); /* writeF(" .sort\n"); */ writeF ("*\n"); } for (v1 = 1; v1 <= 2 * maxvert; v1++) fermmap[v1 - 1] = vertmap[v1 - 1]; writeF ("*\n"); gammaflag = FALSE; } static void formblocks (int *vrtsize) { int v; int l; int lori; int v1, v2, vv; int count; char vnum[11]; /* FormBlocks */ for (v = 1; v <= 2 * maxvert; v++) vertmap[v - 1] = fermmap[v - 1]; /* copy fermMap */ for (count = 1; count <= nloop; count++) /* ferm loop vertex */ { sbld (vnum, "%d", count); writeF (" L Vrt%s=Fl%s;\n", vnum, vnum); indicatorFl[count] = 0; } count = nloop; for (v = 1; v <= vcs.sizet; v++) /* simple vertex */ if (vertmap[v - 1] == 0) { ++(count); vertmap[v - 1] = count; writeF (" L Vrt%d:=%s;\n", count, parsedverts[v - 1]); } writeF (" #call DefineF()\n"); /* end of vertex filling */ *vrtsize = count; for (v = 1; v <= count; v++) { vertinfo[v - 1].vlnc = 0; vertinfo[v - 1].weit = 1; setofb_zero (vertinfo[v - 1].ind); vertinfo[v - 1].g5 = 0; } for (v = 1; v <= vcs.sizet; v++) { v1 = vertmap[v - 1]; if (v1 <= nloop) { vertinfo[v1 - 1].weit += 2; if (fermloops[v1 - 1].g5) { ++(vertinfo[v1 - 1].weit); vertinfo[v1 - 1].g5 = 1; } } for (l = 1; l <= vcs.valence[v - 1]; l++) { vv = vcs.vertlist[v - 1][l - 1].nextvert.vno; v2 = vertmap[vv - 1]; lori = vcs.vertlist[v - 1][l - 1].lorentz; if (lori != 0) { if (v1 != v2) { int np; ++(vertinfo[v1 - 1].vlnc); vertinfo[v1 - 1].link[vertinfo[v1 - 1].vlnc - 1] = v2; setofb_cpy (vertinfo[v1 - 1].ind, setofb_uni (vertinfo[v1 - 1].ind, setofb (lori, _E))); np = vcs.vertlist[v - 1][l - 1].partcl; if (prtclbase[np - 1].hlp == 't') { ++(vertinfo[v1 - 1].vlnc); vertinfo[v1 - 1].link[vertinfo[v1 - 1].vlnc - 1] = v2; setofb_cpy (vertinfo[v1 - 1].ind, setofb_uni (vertinfo[v1 - 1].ind, setofb (lori - 1, _E))); } } vertinfo[v1 - 1].weit += 2; if (insetb (lori, setmassindex)) ++(vertinfo[v1 - 1].weit); } } } /* End of Blk filling */ } /* FormBlocks */ static void emitreducecode (void) { int i, j, v0, v1, v2; indvertset ind1, ind2; int g5_exist; char *pstr; int sgn; for (i = 1; i <= vcs.sizet; i++) { for (j = 1; j <= MAX (1, vcs.valence[i - 1]); j++) { momsubst[j - 1] = vcs.vertlist[i - 1][vertexes[i - 1].subst[j - 1] - 1].moment; indsubst[j - 1] = vcs.vertlist[i - 1][vertexes[i - 1].subst[j - 1] - 1].lorentz; } r_reading0 = vertexes[i - 1].r_vert; pstr = (char *) readExpression (vertexes[i - 1].lgrnptr->description, rd_rcode, act_rcode, free); if (rderrcode != 0) { outFileClose (); finish ("*** Error! CompHEP is restarted"); exit (0); } sgn = 1; if (vcs.valence[i - 1] == 3 && prtclbase[vcs.vertlist[i - 1][0].partcl - 1].cdim == 8 && prtclbase[vcs.vertlist[i - 1][1].partcl - 1].cdim == 8 && prtclbase[vcs.vertlist[i - 1][2].partcl - 1].cdim == 8) { if (vertexes[i - 1].subst[0] > vertexes[i - 1].subst[1]) sgn = -sgn; if (vertexes[i - 1].subst[1] > vertexes[i - 1].subst[2]) sgn = -sgn; if (vertexes[i - 1].subst[0] > vertexes[i - 1].subst[2]) sgn = -sgn; } parsedverts[i - 1] = m_alloc (8 + strlen (pstr)); if (sgn == 1) sprintf (parsedverts[i - 1], "(%s)", pstr + 2); else sprintf (parsedverts[i - 1], "(-1)*(%s)", pstr + 2); free (pstr); } fermprgemit (); formblocks (&n_vrt); makeprgcode (); for (i = 0; i < vcs.sizet; ++i) free (parsedverts[i]); g5_exist = numberg5 != 0; for (i = n_vrt - 1; i >= 1; i--) { v1 = prgcode[i - 1][0]; v2 = prgcode[i - 1][1]; setofb_cpy (ind1, vertinfo[v1 - 1].ind); setofb_cpy (ind2, vertinfo[v2 - 1].ind); if (g5_exist && numberg5 == vertinfo[v1 - 1].g5 + vertinfo[v2 - 1].g5) { writeF (" L Vrt1r=Vrt%d;\n", v1); writeF (" L Vrt2r=Vrt%d;\n", v2); writeF (" #call Iravno0(%d,%d)\n", v1, v2); r_mult (v1, v2, setofb_its (ind1, ind2)); v0 = MIN (v1, v2); writeF (" L VrtTmp=Vrt%d; \n", v0); writeF (" .sort\n"); writeF (" L Vrt%d=Vrt1r; \n", v1); writeF (" .sort\n"); writeF (" L Vrt%d=Vrt2r; \n", v2); writeF (" .sort\n"); r_mult (v1, v2, setofb_its (ind1, ind2)); writeF (" L Vrt%d=Vrt%d+VrtTmp; \n", v0, v0); writeF (" .sort\n"); g5_exist = FALSE; } else r_mult (v1, v2, setofb_its (ind1, ind2)); vertinfo[MIN (v1, v2) - 1].g5 = vertinfo[v1 - 1].g5 + vertinfo[v2 - 1].g5; setofb_cpy (vertinfo[MIN (v1, v2) - 1].ind, setofb_aun (setofb_uni (ind1, ind2), setofb_its (ind1, ind2))); } for (i = 0; i < vcs.sizet; i++) for (j = 0; j < vcs.valence[i]; j++) { int mom, np; mom = vcs.vertlist[i][j].moment; np = vcs.vertlist[i][j].partcl; if (mom > 0 && prtclbase[np - 1].hlp == 't') { if (prtclbase[ghostmother (np) - 1].hlp == '*') writeF ("L Vrt1=Vrt1*(%s^2); \n", prtclbase[np - 1].massidnt); else writeF ("L Vrt1=Vrt1*(%s^2-p%d.p%d); \n", prtclbase[np - 1].massidnt, mom, mom); } } writeF ("*\n"); } void mk_formprograms (void) { int k; int ndel, ncalc, nrest; long nrecord, naxu; csdiagram csd; unsigned ncalctot; shortstr txt; hlpcsptr gstlist, c; vcsect vcs_copy; int i; s_listptr d_facts, df; rmptr t_fact; int NumberDiagrams = 0; int *info_for_sqme; int ServiceSqme; goto_xy (1, 21); scrcolor (Yellow, Blue); print (" FORM code generation \n"); scrcolor (Red, BGmain); print (" Generated........\n"); print (" current diagram :\n"); scrcolor (Yellow, Blue); print (" Press Esc to halt FORM codes generation "); scrcolor (FGmain, BGmain); diagrq = fopen (DIAGRQ_NAME, "rb"); ncalctot = 0; menuq = fopen (MENUQ_NAME, "rb"); /*-------------- service files --------------------------------*/ outFileOpen (scat ("%sresults%cvar_def.h", pathtouser, f_slash)); var_def (); outFileClose (); /*-------------------------------------------------------------*/ info_for_sqme = (int *) malloc ((subproc_sq + 1) * sizeof (int)); for (nsub = 1; nsub <= subproc_sq; nsub++) { rd_menu (menuq, 2, nsub, txt, &ndel, &ncalc, &nrest, &nrecord); fseek (diagrq, nrecord * sizeof (csdiagram), SEEK_SET); naxu = ndel + ncalc + nrest; info_for_sqme[nsub] = naxu; for (ndiagr = 1; ndiagr <= naxu; ndiagr++) { goto_xy (20, 22); print ("%u", ncalctot); goto_xy (20, 23); print ("%u", ndiagr); clr_eol (); FREAD1 (csd, diagrq); if (csd.status != -1) { outFileOpen (scat ("%sresults%cp%d_%d.frm", pathtouser, f_slash, nsub, ndiagr)); NumberDiagrams++; for (k = 0; k < 20; k++) indicatorFl[k] = 0; writeF ("#-\n"); writeF ("#: MaxTermSize 30000\n"); writeF ("#: ContinuationLines 30000\n"); writeLabel ('*'); writeF ("*\n"); transfdiagr (&csd, &vcs); writeF ("\n\n#include procedur.prc\n\n"); cwtarg (&vcs); if (vcs.clrnum == 0) { writeF ("*------- Zero color factor --------\n"); writeF ("L totFactor=0;\n"); writeF ("L numerator=0;\n"); writeF ("L denominator=1;\n"); writeF ("L NumDenominator=1;\n"); writeF (".sort\n"); } else { generateghosts (&vcs, &gstlist); if (gstlist == NULL) { writeF ("*------- non-existent diagram --------\n"); writeF ("L totFactor=0;\n"); writeF ("L numerator=0;\n"); writeF ("L denominator=1;\n"); writeF ("L NumDenominator=1;\n"); writeF (".sort\n"); } else { goto_xy (40, 23); print ("(%% %4d subdiagrams)", gstlist->maxnum); writeF ("* The total number of diagrams %d\n", gstlist->maxnum); preperdiagram (); writesubst (); optimal (); head (); emitfactors (); diagramsrfactors (gstlist, &d_facts, &t_fact); writeF (" L totFactor=%s;\n", rmonomtxt (*t_fact)); writeF (" .sort\n"); writeF (" L totFactor=" "totFactor*SymmFact*AverFact*FermFact*ColorFact;\n"); writeF (" .sort\n\n\n"); clrvm (t_fact->n.v); clrvm (t_fact->d.v); free (t_fact); writeF ("L numerator=0;\n"); c = gstlist; df = d_facts; vcs_copy = vcs; while (c != NULL) { FermLoopIndicator = 0; /* FermLoopIndicator */ coloringvcs (c); writeF ("* diagram number = %d\n", c->num); DiagramToOutFile (&vcs, TRUE, '*'); { int k; int sgn = c->sgn; for (k = 0; k < vcs.sizet; k++) sgn *= vertexes[k].lgrnptr->factor; writeF ("L GhostFact=%d;\n", sgn); } findReversVert (); attachvertexes (); emitreducecode (); writeF (" L numerator=numerator +(%s)*GhostFact*Vrt1;\n", smonomtxt (df->monom)); writeF (" .sort\n"); writeF (" Drop; Ndrop numerator, totFactor;\n"); writeF (" .sort\n"); for (i = 1; i <= FermLoopIndicator; i++) writeF (" trace4,%d;\n", i); writeF (" contract;\n"); writeF (" .sort\n"); writeF (" #call Substitution()\n"); writeF (" .sort\n"); writeF ("*\n"); vcs = vcs_copy; c = c->next; df = df->next; } eraseslist (d_facts); eraseghosts (gstlist); vcs = vcs_copy; writeF (" #call SubstOptimal()\n"); writeF (" #call OtherVar()\n"); writeF (" #call OutputTest(%d,%d)\n", nsub, ndiagr); writeF (" #call DotProduct()\n"); writeF ("*\n"); writeF ("*\n"); writeF (" .sort\n"); emitdenoms (); writeF (" #call DefineDenominator()\n"); } } writeF (" #call OutputC(%d)\n", NumberDiagrams); writeF (".end;\n"); outFileClose (); --nrest; ++ncalctot; if (escpressed ()) goto exi; } } } exi: fclose (diagrq); fclose (menuq); clrbox (1, 21, 70, 24); ServiceSqme = service_sqme (NumberDiagrams, info_for_sqme); }