FORM v5.0.0-35-g6318119
diawrap.cc
Go to the documentation of this file.
1
5/* #[ License : */
6/*
7 * Copyright (C) 1984-2026 J.A.M. Vermaseren
8 * When using this file you are requested to refer to the publication
9 * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
10 * This is considered a matter of courtesy as the development was paid
11 * for by FOM the Dutch physics granting agency and we would like to
12 * be able to track its scientific use to convince FOM of its value
13 * for the community.
14 *
15 * This file is part of FORM.
16 *
17 * FORM is free software: you can redistribute it and/or modify it under the
18 * terms of the GNU General Public License as published by the Free Software
19 * Foundation, either version 3 of the License, or (at your option) any later
20 * version.
21 *
22 * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
23 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
24 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
25 * details.
26 *
27 * You should have received a copy of the GNU General Public License along
28 * with FORM. If not, see <http://www.gnu.org/licenses/>.
29 */
30/* #] License : */
31// #[ Includes : diawrap.cc
32
33extern "C" {
34#include "form3.h"
35}
36
37#include "grccparam.h"
38#include "grcc.h"
39#include <map>
40
41#define MAXPOINTS 120
42
43typedef struct ToPoTyPe {
44 WORD *vert;
45 WORD *vertmax;
46 Options *opt;
47 int cldeg[MAXPOINTS], clnum[MAXPOINTS], clext[MAXPOINTS];
48 int cmind[MAXLEGS+1],cmaxd[MAXLEGS+1];
49 int ncl, nvert;
50} TOPOTYPE;
51
52static void ProcessDiagram(EGraph *eg, void *ti);
53static int processVertex(TOPOTYPE *TopoInf, int pointsremaining, int level);
54
55// #] Includes :
56// #[ LoadModel :
57
58int LoadModel(MODEL *m)
59{
60//
61// First check whether there was a model already.
62// In the Kaneko setup there can be only one at the same time.
63// Hence if there was one we need to remove it first.
64// Unless of course, it was already the model we want.
65//
66 if ( m->grccmodel != NULL ) return(0);
67
68 int i,j,k;
69//
70// First the information that goes into the Model struct.
71// Note that new Model takes over (does not copy) minp.cnamlist.
72//
73 MInput minp;
74 PInput pinp;
75 IInput iinp;
76 if ( m->ncouplings > GRCC_MAXNCPLG ) {
77 MesPrint("Too many coupling constants in model. Current limit is %d.",(WORD)GRCC_MAXNCPLG);
78 MesPrint("Suggestion: recompile Form with a larger value for GRCC_MAXNCPLG");
79 return(-1);
80 }
81 minp.defpart = GRCC_DEFBYCODE;
82 minp.name = (char *)(m->name);
83 minp.ncouple = m->ncouplings;
84 for ( i = 0; i < GRCC_MAXNCPLG; i++ ) minp.cnamlist[i] = NULL;
85 for ( i = 0; i < minp.ncouple; i++ )
86 minp.cnamlist[i] = (char *)(VARNAME(symbols,m->couplings[i]));
87 Model *mdl = new Model(&minp);
88//
89// Now the particles
90//
91 for ( i = 0; i < m->nparticles; i++ ) {
92 if ( minp.defpart == GRCC_DEFBYCODE ) {
93 pinp.name = NULL;
94 pinp.aname = NULL;
95 pinp.pcode = m->vertices[i]->particles[0].number;
96 pinp.acode = m->vertices[i]->particles[1].number;
97 switch ( m->vertices[i]->particles[0].spin ) {
98 case 1:
99 pinp.ptypec = GRCC_PT_Scalar;
100 break;
101 case -1:
102 pinp.ptypec = GRCC_PT_Ghost;
103 break;
104 case 2:
105 if ( m->vertices[i]->particles[0].type == 0 )
106 pinp.ptypec = GRCC_PT_Majorana;
107 else
108 pinp.ptypec = GRCC_PT_Undef;
109 break;
110 case -2:
111 if ( m->vertices[i]->particles[0].type == 0 )
112 pinp.ptypec = GRCC_PT_Majorana;
113 else
114 pinp.ptypec = GRCC_PT_Dirac;
115 break;
116 case 3:
117 pinp.ptypec = GRCC_PT_Vector;
118 break;
119 default:
120 pinp.ptypec = GRCC_PT_Undef;
121 break;
122 }
123 }
124 else {
125 pinp.name = (char *)(VARNAME(functions,m->vertices[i]->particles[0].number));
126 pinp.aname = (char *)(VARNAME(functions,m->vertices[i]->particles[1].number));
127 switch ( m->vertices[i]->particles[0].spin ) {
128 case 1:
129 pinp.ptypen = "scalar";
130 break;
131 case -1:
132 pinp.ptypen = "ghost";
133 break;
134 case 2:
135 if ( m->vertices[i]->particles[0].type == 0 )
136 pinp.ptypen = "majorana";
137 else
138 pinp.ptypen = "undef";
139 break;
140 case -2:
141 if ( m->vertices[i]->particles[0].type == 0 )
142 pinp.ptypen = "majorana";
143 else
144 pinp.ptypen = "dirac";
145 break;
146 case 3:
147 pinp.ptypen = "vector";
148 break;
149 default:
150 pinp.ptypen = "undef";
151 break;
152 }
153 }
154 pinp.extonly = m->vertices[i]->externonly;
155 mdl->addParticle(&pinp);
156 }
157 mdl->addParticleEnd();
158//
159// Now the vertices
160//
161 for ( i = m->nparticles; i < m->invertices; i++ ) {
162 VERTEX *v = m->vertices[i];
163 if ( minp.defpart == GRCC_DEFBYCODE ) {
164 iinp.icode = NODEFUNCTION+i;
165 iinp.name = NULL;
166 }
167 else {
168 iinp.name = "node_";
169 }
170 iinp.nplistn = v->nparticles;
171 for ( j = 0; j < iinp.nplistn; j++ ) {
172 if ( minp.defpart == GRCC_DEFBYCODE ) {
173 iinp.plistc[j] = v->particles[j].number;
174 }
175 else {
176 iinp.plistn[j] = (char *)(VARNAME(functions,v->particles[j].number));
177 }
178 }
179//
180// We need a properly ordered list of coupling constants
181// The ordered list is in m->couplings.
182// For each vertex they are in v->couplings.
183//
184// iinp.cvallist = (int *)Malloc1(m->ncouplings*sizeof(int),"couplings");
185
186 for ( j = 0; j < m->ncouplings; j++ ) {
187 iinp.cvallist[j] = 0;
188 for ( k = 0; k < v->ncouplings; k += 2 ) {
189 if ( v->couplings[k] == m->couplings[j] ) {
190 iinp.cvallist[j] = v->couplings[k+1];
191 break;
192 }
193 }
194 }
195 mdl->addInteraction(&iinp);
196 }
197 mdl->addInteractionEnd();
198 m->grccmodel = (void *)mdl;
199 return(0);
200}
201
202// #] LoadModel :
203// #[ ConvertParticle :
204
205int ConvertParticle(Model *model,int formnum)
206{
207//
208// Returns the grcc number of the particle, because grcc does not convert
209//
210 int i;
211 for ( i = 0; i < model->nParticles; i++ ) {
212 if ( model->particles[i]->pcode == formnum ) { return(i); }
213 else if ( model->particles[i]->acode == formnum ) { return(-i); }
214 }
215 MesPrint("Particle %d not found in model %s",formnum,model->name);
216 Terminate(-1);
217 return(0);
218}
219
220// #] ConvertParticle :
221// #[ ReConvertParticle :
222
223int ReConvertParticle(Model *model,int grccnum)
224{
225//
226// Returns the grcc number of the particle, because grcc does not convert
227//
228 if ( grccnum < 0 ) { return(model->particles[-grccnum]->acode); }
229 else { return(model->particles[grccnum]->pcode); }
230}
231
232// #] ReConvertParticle :
233// #[ numParticle :
234
235int numParticle(MODEL *m,WORD n)
236{
237 int i;
238 for ( i = 0; i < m->nparticles; i++ ) {
239 if ( m->vertices[i]->particles[0].number == n ) return(i);
240 if ( m->vertices[i]->particles[1].number == n ) return(i);
241 }
242 MesPrint("numParticle: particle %d not found in model",n);
243 Terminate(-1);
244 return(-1);
245}
246
247// #] numParticle :
248// #[ ProcessDiagram :
249
250void ProcessDiagram(EGraph *eg, void *ti)
251{
252//
253// This is the routine that gets a complete diagram and passes it on
254// to Form (Generator) for further algebraic manipulations.
255// The term is picked up from AT.diaterm and a new term is constructed
256// in the workspace.
257//
258 GETIDENTITY
259 TERMINFO *info = (TERMINFO *)ti;
260
261 if ( ( info->flags & TOPOLOGIESONLY ) == TOPOLOGIESONLY ) return;
262
263 WORD *term = info->term, *newterm, *oldworkpointer = AT.WorkPointer;
264 WORD *tdia = term + info->diaoffset;
265 WORD *tail = tdia + tdia[1];
266 WORD *tend = term + *term;
267 WORD *fill, *startfill, *cfill, *afill;
268 int i, j, intr;
269 Model *model = (Model *)info->currentModel;
270 MODEL *m = (MODEL *)info->currentMODEL;
271 int numlegs, vect, edge, maxmom = 0;
272
273 newterm = term + *term;
274 for ( i = 1; i < info->diaoffset; i++ ) newterm[i] = term[i];
275 fill = newterm + info->diaoffset;
276//
277// Now get the nodes
278//
279 if ( ( info->flags & WITHOUTNODES ) == 0 ) {
280 for ( i = 0; i < eg->nNodes; i++ ) {
281//
282// node_(number,coupling,particle_1(momentum_1),...,particle_n(momentum_n))
283//
284 numlegs = eg->nodes[i]->deg;
285 startfill = fill;
286 *fill++ = NODEFUNCTION;
287 *fill++ = 0;
288 FILLFUN(fill)
289 *fill++ = -SNUMBER; *fill++ = i+1;
290//
291// Now we put the coupling constants. This is done inside the
292// function for when we want to work with counterterms.
293//
294 if ( !eg->isExternal(i) ) {
295 afill = fill; *fill++ = 0; *fill++ = 1; FILLARG(fill)
296 cfill = fill; *fill++ = 0;
297 intr = eg->nodes[i]->intrct;
298 for ( j = 0; j < model->interacts[intr]->nclist; j++ ) {
299 if ( model->interacts[intr]->clist[j] != 0 ) {
300 *fill++ = SYMBOL; *fill++ = 4;
301 *fill++ = m->couplings[j];
302 *fill++ = model->interacts[intr]->clist[j];
303 }
304 }
305 *fill++ = 1; *fill++ = 1; *fill++ = 3;
306 *cfill = fill - cfill;
307 *afill = fill - afill;
308 }
309 else {
310 *fill++ = -SNUMBER;
311 *fill++ = 1;
312 }
313//
314// Now the particles and their momenta.
315//
316 for ( j = 0; j < numlegs; j++ ) {
317 *fill++ = ARGHEAD+FUNHEAD+6;
318 *fill++ = 0;
319 FILLARG(fill)
320 edge = eg->nodes[i]->edges[j];
321 vect = ABS(edge)-1;
322 *fill++ = 6+FUNHEAD;
323 int a;
324 if ( edge < 0 ) { a = ReConvertParticle(model,-eg->edges[vect]->ptcl); }
325 else { a = ReConvertParticle(model,eg->edges[vect]->ptcl); }
326 *fill++ = a;
327 *fill++ = FUNHEAD+2;
328 FILLFUN(fill)
329 *fill++ = edge < 0 ? -MINVECTOR: -VECTOR;
330 if ( numlegs == 1 || vect < info->numextern ) { // Look up in set of external momenta
331 *fill++ = SetElements[Sets[info->externalset].first+vect];
332 }
333 else { // Look up in set of internal momenta set
334 *fill++ = SetElements[Sets[info->internalset].first+(vect-eg->nExtern)];
335 // determine the number of momenta required from internalset:
336 maxmom = MaX(maxmom, vect-eg->nExtern);
337 }
338 *fill++ = 1; *fill++ = 1; *fill++ = 3;
339 }
340 startfill[1] = fill-startfill;
341 }
342 }
343 if ( ( info->flags & WITHEDGES ) == WITHEDGES ) {
344 for ( i = 0; i < eg->nEdges; i++ ) {
345 int n1 = eg->edges[i]->nodes[0];
346 int n2 = eg->edges[i]->nodes[1];
347// int l1 = eg->edges[i]->nlegs[0];
348// int l2 = eg->edges[i]->nlegs[1];
349
350 startfill = fill;
351 *fill++ = EDGE;
352 *fill++ = 0;
353 FILLFUN(fill)
354 *fill++ = -SNUMBER; *fill++ = i+1; // number of the edge
355//
356 *fill++ = ARGHEAD+FUNHEAD+6;
357 *fill++ = 0;
358 FILLARG(fill)
359 *fill++ = 6+FUNHEAD;
360 int a = ReConvertParticle(model,eg->edges[i]->ptcl);
361 *fill++ = a;
362 *fill++ = FUNHEAD+2;
363 FILLFUN(fill)
364 *fill++ = -VECTOR;
365 // Look up in set of internal momenta set
366 if ( i < info->numextern ) { // Look up in set of external momenta
367 *fill++ = SetElements[Sets[info->externalset].first+i];
368 }
369 else { // Look up in set of internal momenta set
370 *fill++ = SetElements[Sets[info->internalset].first+(i-eg->nExtern)];
371 maxmom = MaX(maxmom, i-eg->nExtern);
372 }
373 *fill++ = 1; *fill++ = 1; *fill++ = 3;
374//
375 *fill++ = -SNUMBER; *fill++ = n1+1; // number of the node from
376 *fill++ = -SNUMBER; *fill++ = n2+1; // number of the node to
377 startfill[1] = fill - startfill;
378 }
379 }
380 if ( ( info->flags & WITHBLOCKS ) == WITHBLOCKS ) {
381 for ( i = 0; i < eg->econn->nblocks; i++ ) {
382 startfill = fill;
383 *fill++ = BLOCK;
384 *fill++ = 0;
385 FILLFUN(fill);
386 *fill++ = -SNUMBER;
387 *fill++ = i+1;
388 *fill++ = -SNUMBER;
389 *fill++ = eg->econn->blocks[i].loop;
390//
391// Now we have to make a list of all nodes inside this block
392//
393 int bnodes[GRCC_MAXNODES], k;
394 WORD *argfill = fill, *funfill;
395 *fill++ = 0; *fill++ = 0; FILLARG(fill)
396 *fill++ = 0;
397 for ( k = 0; k < GRCC_MAXNODES; k++ ) bnodes[k] = 0;
398 for ( k = 0; k < eg->econn->blocks[i].nmedges; k++ ) {
399 bnodes[eg->econn->blocks[i].edges[k][0]] = 1;
400 bnodes[eg->econn->blocks[i].edges[k][1]] = 1;
401 }
402 for ( k = 0; k < GRCC_MAXNODES; k++ ) {
403 if ( bnodes[k] == 0 ) continue;
404//
405// Now we put the node inside this argument.
406//
407 funfill = fill;
408 *fill++ = NODEFUNCTION;
409 *fill++ = 0;
410 FILLFUN(fill)
411 *fill++ = -SNUMBER; *fill++ = k+1;
412 numlegs = eg->nodes[k]->deg;
413 for ( j = 0; j < numlegs; j++ ) {
414 edge = eg->nodes[k]->edges[j];
415 vect = ABS(edge)-1;
416 *fill++ = -VECTOR;
417 if ( numlegs == 1 || vect < info->numextern ) { // Look up in set of external momenta
418 *fill++ = SetElements[Sets[info->externalset].first+vect];
419 }
420 else { // Look up in set of internal momenta set
421 *fill++ = SetElements[Sets[info->internalset].first+(vect-info->numextern)];
422 maxmom = MaX(maxmom, vect-info->numextern);
423 }
424 }
425 funfill[1] = fill-funfill;
426 }
427 *fill++ = 1; *fill++ = 1; *fill++ = 3;
428 *argfill = fill - argfill;
429 argfill[ARGHEAD] = argfill[0] - ARGHEAD;
430 startfill[1] = fill-startfill;
431 }
432 }
433 if ( ( info->flags & WITHONEPISETS ) == WITHONEPISETS ) {
434 for ( i = 0; i < eg->econn->nopic; i++ ) {
435 startfill = fill;
436 *fill++ = ONEPI;
437 *fill++ = 0;
438 FILLFUN(fill);
439 *fill++ = -SNUMBER;
440 *fill++ = i+1;
441 *fill++ = -ONEPI;
442 for ( j = 0; j < eg->econn->opics[i].nnodes; j++ ) {
443 *fill++ = -SNUMBER;
444 *fill++ = eg->econn->opics[i].nodes[j]+1;
445 }
446 startfill[1] = fill-startfill;
447 }
448 }
449//
450// Topology counter. We have exaggerated a bit with the eye on the far future.
451//
452 if ( info->numtopo-1 < MAXPOSITIVE ) {
453 *fill++ = TOPO; *fill++ = FUNHEAD+2; FILLFUN(fill)
454 *fill++ = -SNUMBER; *fill++ = (WORD)(info->numtopo-1);
455 }
456 else if ( info->numtopo-1 < FULLMAX-1 ) {
457 *fill++ = TOPO; *fill++ = FUNHEAD+ARGHEAD+4; FILLFUN(fill)
458 *fill++ = ARGHEAD+4; *fill++ = 0; FILLARG(fill)
459 *fill++ = 4;
460 *fill++ = (WORD)((info->numtopo-1) & WORDMASK);
461 *fill++ = 1; *fill++ = 3;
462 }
463 else { // for now: science fiction
464 *fill++ = TOPO; *fill++ = FUNHEAD+ARGHEAD+6; FILLFUN(fill)
465 *fill++ = ARGHEAD+6; *fill++ = 0; FILLARG(fill)
466 *fill++ = 6; *fill++ = (WORD)((info->numtopo-1) >> BITSINWORD);
467 *fill++ = (WORD)((info->numtopo-1) & WORDMASK);
468 *fill++ = 0; *fill++ = 1; *fill++ = 5;
469 }
470//
471// Symmetry factors. We let Normalize do the multiplication.
472//
473 if ( eg->nsym != 1 ) {
474 *fill++ = SNUMBER; *fill++ = 4; *fill++ = (WORD)eg->nsym; *fill++ = -1;
475 }
476 if ( eg->esym != 1 ) {
477 *fill++ = SNUMBER; *fill++ = 4; *fill++ = (WORD)eg->esym; *fill++ = -1;
478 }
479 if ( eg->extperm != 1 ) {
480 *fill++ = SNUMBER; *fill++ = 4; *fill++ = (WORD)eg->extperm; *fill++ = 1;
481 }
482//
483// verify internalset has sufficient momenta:
484//
485 if ( maxmom >= Sets[info->internalset].last - Sets[info->internalset].first ) {
486 MLOCK(ErrorMessageLock);
487 MesPrint("&Insufficient internal momenta in diagrams_");
488 MUNLOCK(ErrorMessageLock);
489 Terminate(-1);
490 }
491//
492// finish it off
493//
494 while ( tail < tend ) *fill++ = *tail++;
495 if ( eg->fsign < 0 ) fill[-1] = -fill[-1];
496 *newterm = fill - newterm;
497 AT.WorkPointer = fill;
498
499 Generator(BHEAD newterm,info->level);
500 AT.WorkPointer = oldworkpointer;
501}
502
503// #] ProcessDiagram :
504// #[ fendMG :
505
506Bool fendMG(EGraph *eg, void *ti)
507{
508 DUMMYUSE(eg);
509 DUMMYUSE(ti);
510 return True;
511}
512
513// #] fendMG :
514// #[ ProcessTopology :
515
516Bool ProcessTopology(EGraph *eg, void *ti)
517{
518//
519// This routine is called for each new topology.
520// New convention: return True; generate the corresponding diagrams if needed
521// return False; skip diagram generation (when asked for).
522//
523 TERMINFO *info = (TERMINFO *)ti;
524
525// This seems to work properly. It was disabled before.
526#define WITHEARLYVETO
527#ifdef WITHEARLYVETO
528 if ( ( ( info->flags & CHECKEXTERN ) == CHECKEXTERN ) && info->currentMODEL != NULL ) {
529 int i, j;
530 int numlegs, vect, edge;
531 for ( i = 0; i < eg->nNodes; i++ ) {
532 if ( eg->isExternal(i) ) continue;
533 numlegs = eg->nodes[i]->deg;
534 for ( j = 0; j < numlegs; j++ ) {
535 edge = eg->nodes[i]->edges[j];
536 vect = ABS(edge)-1;
537 if ( vect < info->numextern && info->legcouple[vect][numlegs] == 0 ) {
538
539// This cannot be.
540
541 info->numtopo++;
542 return False;
543 }
544 }
545 }
546 }
547#endif
548
549 if ( ( info->flags & TOPOLOGIESONLY ) == 0 ) {
550 info->numtopo++;
551 return True;
552 }
553//
554// Now we are just generating topologies.
555//
556 GETIDENTITY
557 WORD *term = info->term, *newterm, *oldworkpointer = AT.WorkPointer;
558 WORD *tdia = term + info->diaoffset;
559 WORD *tail = tdia + tdia[1];
560 WORD *tend = term + *term;
561 WORD *fill, *startfill;
562 Model *model = (Model *)info->currentModel;
563 MODEL *m = (MODEL *)info->currentMODEL;
564 int i, j;
565 int numlegs, vect, edge, maxmom = 0;
566
567 newterm = term + *term;
568 for ( i = 1; i < info->diaoffset; i++ ) newterm[i] = term[i];
569 fill = newterm + info->diaoffset;
570//
571// Now get the nodes
572//
573 for ( i = 0; i < eg->nNodes; i++ ) {
574//
575// node_(number,momentum_1,...,momentum_n)
576//
577 numlegs = eg->nodes[i]->deg;
578 startfill = fill;
579 *fill++ = NODEFUNCTION;
580 *fill++ = 0;
581 FILLFUN(fill)
582 *fill++ = -SNUMBER; *fill++ = i+1;
583 if ( model != NULL && m != NULL ) {
584 if ( !eg->isExternal(i) ) {
585 WORD *afill = fill; *fill++ = 0; *fill++ = 1; FILLARG(fill)
586 WORD *cfill = fill; *fill++ = 0;
587 int cpl = 2*eg->nodes[i]->extloop+eg->nodes[i]->deg-2;
588 *fill++ = SYMBOL; *fill++ = 4;
589 *fill++ = m->couplings[0];
590 *fill++ = cpl;
591 *fill++ = 1; *fill++ = 1; *fill++ = 3;
592 *cfill = fill - cfill;
593 *afill = fill - afill;
594 if ( *afill == ARGHEAD+8 && afill[ARGHEAD+4] == 1 ) {
595 fill = afill; *fill++ = -SYMBOL; *fill++ = afill[ARGHEAD+3];
596 }
597 }
598 else {
599 *fill++ = -SNUMBER;
600 *fill++ = 1;
601 }
602 }
603//
604// Now the momenta.
605//
606 for ( j = 0; j < numlegs; j++ ) {
607 edge = eg->nodes[i]->edges[j];
608 vect = ABS(edge)-1;
609 *fill++ = edge < 0 ? -MINVECTOR: -VECTOR;
610 if ( numlegs == 1 || vect < info->numextern ) { // Look up in set of external momenta
611 *fill++ = SetElements[Sets[info->externalset].first+vect];
612 }
613 else { // Look up in set of internal momenta set
614 *fill++ = SetElements[Sets[info->internalset].first+(vect-info->numextern)];
615 // determine the number of momenta required from internalset:
616 maxmom = MaX(maxmom, vect-info->numextern);
617 }
618 }
619 startfill[1] = fill-startfill;
620 }
621 if ( ( info->flags & WITHEDGES ) == WITHEDGES ) {
622 for ( i = 0; i < eg->nEdges; i++ ) {
623 int n1 = eg->edges[i]->nodes[0];
624 int n2 = eg->edges[i]->nodes[1];
625// int l1 = eg->edges[i]->nlegs[0];
626// int l2 = eg->edges[i]->nlegs[1];
627 startfill = fill;
628 *fill++ = EDGE;
629 *fill++ = 0;
630 FILLFUN(fill)
631 *fill++ = -SNUMBER; *fill++ = i+1; // number of the edge
632//
633 *fill++ = -VECTOR;
634 if ( i < info->numextern ) { // Look up in set of external momenta
635 *fill++ = SetElements[Sets[info->externalset].first+i];
636 }
637 else { // Look up in set of internal momenta set
638 *fill++ = SetElements[Sets[info->internalset].first+(i-eg->nExtern)];
639 maxmom = MaX(maxmom, i-eg->nExtern);
640 }
641//
642 *fill++ = -SNUMBER; *fill++ = n1+1; // number of the node from
643 *fill++ = -SNUMBER; *fill++ = n2+1; // number of the node to
644 startfill[1] = fill - startfill;
645 }
646 }
647//
648// Block information
649//
650 if ( ( info->flags & WITHBLOCKS ) == WITHBLOCKS ) {
651 for ( i = 0; i < eg->econn->nblocks; i++ ) {
652 startfill = fill;
653 *fill++ = BLOCK;
654 *fill++ = 0;
655 FILLFUN(fill);
656 *fill++ = -SNUMBER;
657 *fill++ = i+1;
658 *fill++ = -SNUMBER;
659 *fill++ = eg->econn->blocks[i].loop;
660//
661// Now we have to make a list of all nodes inside this block
662//
663 int bnodes[GRCC_MAXNODES], k;
664 WORD *argfill = fill, *funfill;
665 *fill++ = 0; *fill++ = 0; FILLARG(fill)
666 *fill++ = 0;
667 for ( k = 0; k < GRCC_MAXNODES; k++ ) bnodes[k] = 0;
668 for ( k = 0; k < eg->econn->blocks[i].nmedges; k++ ) {
669 bnodes[eg->econn->blocks[i].edges[k][0]] = 1;
670 bnodes[eg->econn->blocks[i].edges[k][1]] = 1;
671 }
672 for ( k = 0; k < GRCC_MAXNODES; k++ ) {
673 if ( bnodes[k] == 0 ) continue;
674//
675// Now we put the node inside this argument.
676//
677 funfill = fill;
678 *fill++ = NODEFUNCTION;
679 *fill++ = 0;
680 FILLFUN(fill)
681 *fill++ = -SNUMBER; *fill++ = k+1;
682 numlegs = eg->nodes[k]->deg;
683 for ( j = 0; j < numlegs; j++ ) {
684 edge = eg->nodes[k]->edges[j];
685 vect = ABS(edge)-1;
686 *fill++ = -VECTOR;
687 if ( numlegs == 1 || vect < info->numextern ) { // Look up in set of external momenta
688 *fill++ = SetElements[Sets[info->externalset].first+vect];
689 }
690 else { // Look up in set of internal momenta set
691 *fill++ = SetElements[Sets[info->internalset].first+(vect-info->numextern)];
692 maxmom = MaX(maxmom, vect-info->numextern);
693 }
694 }
695 funfill[1] = fill-funfill;
696 }
697 *fill++ = 1; *fill++ = 1; *fill++ = 3;
698 *argfill = fill - argfill;
699 argfill[ARGHEAD] = argfill[0] - ARGHEAD;
700 startfill[1] = fill-startfill;
701 }
702
703// if ( eg->econn->narticuls > 0 ) {
704// startfill = fill;
705// *fill++ = BLOCK;
706// *fill++ = 0;
707// FILLFUN(fill);
708// for ( i = 0; i < eg->econn->snodes; i++ ) {
709// if ( eg->econn->articuls[i] != 0 ) {
710// *fill++ = -SNUMBER;
711// *fill++ = i+1;
712// }
713// }
714// startfill[1] = fill-startfill;
715// }
716 }
717 if ( ( info->flags & WITHONEPISETS ) == WITHONEPISETS ) {
718 for ( i = 0; i < eg->econn->nopic; i++ ) {
719 startfill = fill;
720 *fill++ = ONEPI;
721 *fill++ = 0;
722 FILLFUN(fill);
723 *fill++ = -SNUMBER;
724 *fill++ = i+1;
725 *fill++ = -ONEPI;
726 for ( j = 0; j < eg->econn->opics[i].nnodes; j++ ) {
727 *fill++ = -SNUMBER;
728 *fill++ = eg->econn->opics[i].nodes[j]+1;
729 }
730 startfill[1] = fill-startfill;
731 }
732 }
733//
734// Topology counter. We have exaggerated a bit with the eye on the far future.
735//
736 if ( info->numtopo < MAXPOSITIVE ) {
737 *fill++ = TOPO; *fill++ = FUNHEAD+2; FILLFUN(fill)
738 *fill++ = -SNUMBER; *fill++ = (WORD)(info->numtopo);
739 }
740 else if ( info->numtopo < FULLMAX-1 ) {
741 *fill++ = TOPO; *fill++ = FUNHEAD+ARGHEAD+4; FILLFUN(fill)
742 *fill++ = ARGHEAD+4; *fill++ = 0; FILLARG(fill)
743 *fill++ = 4;
744 *fill++ = (WORD)(info->numtopo & WORDMASK);
745 *fill++ = 1; *fill++ = 3;
746 }
747 else { // for now: science fiction
748 *fill++ = TOPO; *fill++ = FUNHEAD+ARGHEAD+6; FILLFUN(fill)
749 *fill++ = ARGHEAD+6; *fill++ = 0; FILLARG(fill)
750 *fill++ = 6; *fill++ = (WORD)(info->numtopo >> BITSINWORD);
751 *fill++ = (WORD)(info->numtopo & WORDMASK);
752 *fill++ = 0; *fill++ = 1; *fill++ = 5;
753 }
754//
755// verify internalset has sufficient momenta:
756//
757 if ( maxmom >= Sets[info->internalset].last - Sets[info->internalset].first ) {
758 MLOCK(ErrorMessageLock);
759 MesPrint("&Insufficient internal momenta in diagrams_");
760 MUNLOCK(ErrorMessageLock);
761 Terminate(-1);
762 }
763//
764// finish it off
765//
766 while ( tail < tend ) *fill++ = *tail++;
767 if ( eg->fsign < 0 ) fill[-1] = -fill[-1];
768 *newterm = fill - newterm;
769 AT.WorkPointer = fill;
770
771 Generator(BHEAD newterm,info->level);
772 AT.WorkPointer = oldworkpointer;
773 info->numtopo++;
774 return False;
775}
776
777// #] ProcessTopology :
778// #[ SetDualOpts :
779void SetDualOpts(int *opt, const WORD num, const int key, const char* key_name,
780 const int dual, const char* dual_name, const int val, const int dval) {
781
782 if ( ( num & key ) == key ) {
783 if ( ( num & dual ) == dual ) {
784 MLOCK(ErrorMessageLock);
785 MesPrint("&Conflicting diagram filters: %s and %s.", key_name, dual_name);
786 MUNLOCK(ErrorMessageLock);
787 Terminate(-1);
788 }
789 else {
790 *opt = val;
791 }
792 }
793 else {
794 if ( ( num & dual ) == dual ) {
795 *opt = dval;
796 }
797 else {
798 // The default value is always 0.
799 *opt = 0;
800 }
801 }
802}
803// #] SetDualOpts :
804// #[ GenDiagrams :
805
806int GenDiagrams(PHEAD WORD *term, WORD level)
807{
808 Model *model;
809 MODEL *m;
810 Options *opt;
811 Process *proc;
812 int pid = 1, x;
813 int babble = AC.GrccVerbose ? 2 : 0;
814 TERMINFO info;
815 WORD inset,outset,*coupl,setnum,optionnumber = 0;
816 int i, j, cpl[GRCC_MAXNCPLG];
817 int ninitl, initlPart[GRCC_MAXLEGS], nfinal, finalPart[GRCC_MAXLEGS];
818 for ( i = 0; i < GRCC_MAXNCPLG; i++ ) cpl[i] = 0;
819 std::map<int,int> momlist;
820//
821// Here we create an object of type Option and load it up.
822// Next we run the diagram generation on it.
823//
824 info.term = term;
825 info.level = level;
826 info.diaoffset = AR.funoffset;
827 info.externalset = term[info.diaoffset+FUNHEAD+7];
828 info.internalset = term[info.diaoffset+FUNHEAD+9];
829 info.flags = 0;
830 inset = term[info.diaoffset+FUNHEAD+3];
831 outset = term[info.diaoffset+FUNHEAD+5];
832 coupl = term + info.diaoffset + FUNHEAD + 10;
833 if ( *coupl < 0 ) {
834 if ( term[info.diaoffset+1] > FUNHEAD + 12 ) {
835 optionnumber = term[info.diaoffset+FUNHEAD+13];
836 }
837 }
838 else {
839 if ( term[info.diaoffset+1] > *coupl+FUNHEAD+10 )
840 optionnumber = term[info.diaoffset+*coupl+FUNHEAD+11];
841 }
842 setnum = term[info.diaoffset+FUNHEAD+1];
843
844 m = AC.models[SetElements[Sets[setnum].first]];
845 LoadModel(m);
846 model = (Model *)m->grccmodel;
847
848 info.currentModel = (void *)model;
849 info.currentMODEL = (void *)m;
850 info.numdia = 0;
851 info.numtopo = 1;
852 info.flags = optionnumber;
853
854 opt = new Options();
855
856 opt->setOutAG(ProcessDiagram, &info);
857 opt->setOutMG(ProcessTopology, &info);
858
859 opt->values[GRCC_OPT_SymmInitial] = ( optionnumber & WITHSYMMETRIZEI ) == WITHSYMMETRIZEI;
860 opt->values[GRCC_OPT_SymmFinal] = ( optionnumber & WITHSYMMETRIZEF ) == WITHSYMMETRIZEF;
861
862 // WITHBLOCKS controls output formatting. We could introduce an extra filtering option
863 // corresponding to GRCC_OPT_Block, which is somewhat like Qgraf "onevi" but not quite
864 // the same currently.
865 //opt->values[GRCC_OPT_Block] = ;
866
867 // Now the "qgraf-compatible filtering options":
868 int qgopt[GRCC_QGRAF_OPT_Size];
869 SetDualOpts(&qgopt[GRCC_QGRAF_OPT_ONEPI], optionnumber,ONEPARTI, "ONEPI_", ONEPARTR, "ONEPR_", 1,-1);
870 SetDualOpts(&qgopt[GRCC_QGRAF_OPT_ONSHELL], optionnumber,ONSHELL, "ONSHELL_", OFFSHELL, "OFFSHELL_", 1,-1);
871 SetDualOpts(&qgopt[GRCC_QGRAF_OPT_NOSIGMA], optionnumber,NOSIGMA, "NOSIGMA_", SIGMA, "SIGMA_", 1,-1);
872 SetDualOpts(&qgopt[GRCC_QGRAF_OPT_NOSNAIL], optionnumber,NOSNAIL, "NOSNAIL_", SNAIL, "SNAIL_", 1,-1);
873 SetDualOpts(&qgopt[GRCC_QGRAF_OPT_NOTADPOLE],optionnumber,NOTADPOLE,"NOTADPOLE_",TADPOLE , "TADPOLE_", 1,-1);
874 SetDualOpts(&qgopt[GRCC_QGRAF_OPT_SIMPLE], optionnumber,SIMPLE, "SIMPLE_", NOTSIMPLE,"NOTSIMPLE_",1,-1);
875 SetDualOpts(&qgopt[GRCC_QGRAF_OPT_BIPART], optionnumber,BIPART, "BIPART_", NONBIPART,"NONBIPART_",1,-1);
876 SetDualOpts(&qgopt[GRCC_QGRAF_OPT_CYCLI], optionnumber,CYCLI, "CYCLI_", CYCLR, "CYCLR_", 1,-1);
877 SetDualOpts(&qgopt[GRCC_QGRAF_OPT_FLOOP], optionnumber,FLOOP, "FLOOP_", NOTFLOOP, "NOTFLOOP_", 1,-1);
878 // Now set the options internally:
879 opt->setQGrafOpt(qgopt);
880
881 opt->setOutputF(False,"");
882 opt->setOutputP(False,"");
883 opt->printLevel(babble);
884
885// Load the various arrays.
886 ninitl = Sets[inset].last - Sets[inset].first;
887 for ( i = 0; i < ninitl; i++ ) {
888 x = SetElements[Sets[inset].first+i];
889 initlPart[i] = ConvertParticle(model,x);
890 info.legcouple[i] = m->vertices[numParticle(m,x)]->couplings;
891 }
892 nfinal = Sets[outset].last - Sets[outset].first;
893 for ( i = 0; i < nfinal; i++ ) {
894 x = SetElements[Sets[outset].first+i];
895 finalPart[i] = ConvertParticle(model,x);
896 info.legcouple[i+ninitl] = m->vertices[numParticle(m,x)]->couplings;
897 }
898 info.numextern = ninitl + nfinal;
899 for ( i = 2; i <= MAXLEGS; i++ ) {
900 if ( m->legcouple[i] == 1 ) {
901 for ( j = 0; j < info.numextern; j++ ) {
902 if ( info.legcouple[j][i] == 0 ) { info.flags |= CHECKEXTERN; goto Go_on; }
903 }
904 }
905 }
906
907 // Check that we have sufficient external momenta in the set:
908 if ( info.numextern > Sets[info.externalset].last - Sets[info.externalset].first ) {
909 MLOCK(ErrorMessageLock);
910 MesPrint("&Insufficient external momenta in diagrams_");
911 MUNLOCK(ErrorMessageLock);
912 Terminate(-1);
913 }
914
915 // Check that none of the supplied momenta are negative or repeated:
916 for ( i = 0; i < Sets[info.externalset].last - Sets[info.externalset].first; i++ ) {
917 const int momcode = SetElements[Sets[info.externalset].first + i];
918 if ( momcode < AM.OffsetVector ) {
919 MLOCK(ErrorMessageLock);
920 MesPrint("&Invalid negative external momentum in diagrams_: -%s",
921 VARNAME(vectors, momcode+WILDMASK-AM.OffsetVector));
922 MUNLOCK(ErrorMessageLock);
923 Terminate(-1);
924 }
925 momlist[momcode]++;
926 if ( momlist[momcode] != 1 ) {
927 MLOCK(ErrorMessageLock);
928 MesPrint("&Invalid repeated momentum in diagrams_: %s",
929 VARNAME(vectors, momcode-AM.OffsetVector));
930 MUNLOCK(ErrorMessageLock);
931 Terminate(-1);
932 }
933 }
934 for ( i = 0; i < Sets[info.internalset].last - Sets[info.internalset].first; i++ ) {
935 const int momcode = SetElements[Sets[info.internalset].first + i];
936 if ( momcode < AM.OffsetVector ) {
937 MLOCK(ErrorMessageLock);
938 MesPrint("&Invalid negative internal momentum in diagrams_: -%s",
939 VARNAME(vectors, momcode+WILDMASK-AM.OffsetVector));
940 MUNLOCK(ErrorMessageLock);
941 Terminate(-1);
942 }
943 momlist[momcode]++;
944 if ( momlist[momcode] != 1 ) {
945 MLOCK(ErrorMessageLock);
946 MesPrint("&Invalid repeated momentum in diagrams_: %s",
947 VARNAME(vectors, momcode-AM.OffsetVector));
948 MUNLOCK(ErrorMessageLock);
949 Terminate(-1);
950 }
951 }
952
953Go_on:;
954//
955// Now we have to sort out the coupling constants.
956// The argument at coupl can be of type -SNUMBER, -SYMBOL or generic
957// It has however already be tested for syntax.
958// Note that one cannot have 1 for the coupling constants.
959// In that case one should select 0 loops or something equivalent.
960//
961 if ( *coupl == -SNUMBER ) { // Number of loops
962//
963// This is the complicated case.
964// We have to compute the number of coupling constants and then
965// generate diagrams for all combinations with the proper power.
966//
967 int nc = coupl[1]*2 + ninitl + nfinal - 2;
968 int *scratch = (int *)Malloc1(nc*sizeof(int),"DistrN");
969 scratch[0] = -2; // indicating startup cq first call.
970
971 if ( ( info.flags & TOPOLOGIESONLY ) == 0 ) {
972 while ( DistrN(nc,cpl,m->ncouplings,scratch) ) {
973 proc = new Process(pid, model, opt, ninitl, initlPart, nfinal, finalPart, cpl);
974 delete proc;
975 info.numtopo = 1;
976 }
977 }
978 else {
979 cpl[0] = nc;
980 proc = new Process(pid, model, opt, ninitl, initlPart, nfinal, finalPart, cpl);
981 delete proc;
982 }
983 M_free(scratch,"DistrN");
984 opt->end();
985 delete opt;
986 return(0);
987 }
988 else if ( *coupl == -SYMBOL ) { // Just a single power of one constant
989 for ( i = 0; i < m->ncouplings; i++ ) {
990 if ( m->couplings[i] == coupl[1] ) {
991 cpl[i] = 1;
992 break;
993 }
994 }
995 }
996 else { // One term with powers of coupling constants
997 WORD *t, *tstop;
998 t = coupl + ARGHEAD+3;
999 tstop = coupl+*coupl; tstop -= ABS(tstop[-1]);
1000 while ( t < tstop ) {
1001 for ( i = 0; i < m->ncouplings; i++ ) {
1002 if ( m->couplings[i] == *t ) {
1003 cpl[i] = t[1];
1004 break;
1005 }
1006 }
1007 t += 2;
1008 }
1009 }
1010/*
1011 And now the generation:
1012*/
1013 proc = new Process(pid, model, opt, ninitl, initlPart, nfinal, finalPart, cpl);
1014 opt->end();
1015 delete proc;
1016 delete opt;
1017 return(0);
1018}
1019
1020// #] GenDiagrams :
1021// #[ processVertex :
1022
1023// Routine is to be used recursively to work its way through a list
1024// of possible vertices. The array of vertices is in TopoInf->vert
1025// with TopoInf->nvert the number of possible vertices.
1026// Currently we allow in TopoInf->vert only vertices with 3 or more edges.
1027//
1028// We work with a point system. Each n-point vertex contributes n-2 points.
1029// When all points are assigned, we can call mgraph->generate().
1030//
1031// The number of vertices of a given number of edges is stored in
1032// TopoInf->clnum[..] but the loop that determines how many there are
1033// may be limited by the corresponding element in TopoInf->vertmax[level]
1034
1035int processVertex(TOPOTYPE *TopoInf, int pointsremaining, int level)
1036{
1037 int i, j;
1038
1039 for ( i = pointsremaining, j = 0; i >= 0; i -= TopoInf->vert[level]-2, j++ ) {
1040 if ( TopoInf->vertmax && TopoInf->vertmax[level] >= 0
1041 && j > TopoInf->vertmax[level] ) break;
1042 if ( i == 0 ) { // We got one!
1043 TopoInf->cldeg[TopoInf->ncl] = TopoInf->vert[level];
1044 TopoInf->clnum[TopoInf->ncl] = j;
1045 TopoInf->clext[TopoInf->ncl] = 0;
1046 TopoInf->ncl++;
1047
1048 MGraph *mgraph = new MGraph(1, TopoInf->ncl, TopoInf->cldeg,
1049 TopoInf->clnum, TopoInf->clext,
1050 TopoInf->cmind, TopoInf->cmaxd, TopoInf->opt);
1051
1052 mgraph->generate();
1053
1054 delete mgraph;
1055
1056 TopoInf->ncl--;
1057
1058 break;
1059 }
1060 if ( level < TopoInf->nvert-1 ) {
1061 if ( j > 0 ) {
1062 TopoInf->cldeg[TopoInf->ncl] = TopoInf->vert[level];
1063 TopoInf->clnum[TopoInf->ncl] = j;
1064 TopoInf->clext[TopoInf->ncl] = 0;
1065 TopoInf->ncl++;
1066 }
1067 if ( processVertex(TopoInf,i,level+1) < 0 ) return(-1);
1068 if ( j > 0 ) { TopoInf->ncl--; }
1069 }
1070 }
1071 return(0);
1072}
1073
1074// #] processVertex :
1075
EGraph
Definition grcc.h:624
MGraph
Definition grcc.h:793
Model
Definition grcc.h:284
Generator
int Generator(PHEAD WORD *, WORD)
Definition proces.c:3249
Generated by doxygen 1.9.8