!2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! ! ---------------------------- ! beowulfsubs.f Version 4.0 ! ---------------------------- ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine version ! use beowulf_parameters implicit none write(nout,*)'beowulf 4.0 subroutines 4 September 2005' write(nout,*)'These subroutines are unchanged from' write(nout,*)'beowulf 3.3 subroutines 23 December 2004' write(nout,*)'Coulomb gauge and Feynman gauge plus showers' RETURN END subroutine version ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! ! Subroutines for numerical integration of jet cross sections in ! electron-positron annihilation. -- D. E. Soper ! ! First code: 29 November 1992. ! Latest revision: see Version subroutine above. ! ! The main program and subroutines that a user might want to modify ! are contained in the companion package, beowulf.f. In particular, a ! user may very well want to modify parameter settings in the main ! program and to change the observables calculated in the subroutine ! HROTHGAR and in the functions CALStype(nparticles,kf,index). Subroutines ! that can be modified only at extreme peril to the reliability of ! the results are in this package, beowulfsubs.f. ! ! There are two parallel calculations. Program beowulf calculates a ! sample integral, which by default is the average value of ! (1 - thrust)^2. These are summed in the variable INTEGRAL and ! reported upon completion of the program. The program also computes ! a simple check integral in order to check on the jacobians etc. ! In the meantime, for each point in loop space and each final ! state cut, the program reports the corresponding point in the space ! of final state momenta along with the corresponding weight (Feynman ! diagram times jacobian factors) to the subroutine HROTHGAR, which ! multiplies by the measurement functions CALS corresponding to the ! measurements desired and accumulates the results. ! ! There are four function subroutines representing the Feynman graphs ! that are generated by Mathematica routines. These are ! ! function feynman(graphnumber,flavorsetnumber,kc,cut,mumsbar,flag) ! function feynman0(graphnumber,flavorsetnumber,kin,cut) ! function feynmanSH0(graphnumber,flavorsetnumber,kin,cut, & ! vquark,vqbar,tglue) ! function softsubtraction(k,absk, & ! graphnumber,flavorsetnumber,cutnumber) ! ! In order to control roundoff errors, a point in loop space is rejected ! if the point is too near a singularity or if there is too much ! cancellation in the contribution from that point to INTEGRAL. ! ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! ! PROGRAM STRUCTURE !main program beowulf ! graphcountinit ! getcounts and getcountsf ! hrothgar... (interface) ! randominit ! daytime ! version ! reno... (main calculation) ! diagnostic ! !subroutine reno ! timing ! hrothgar ! getnewgraph ! findtypes (maptypes for chosing points) ! getnewcut ! finda (matrix to get propagator momenta from loop momenta) ! newpoint ! checkpoint ! calculate... (calculate at chosen point) ! !subroutine calculate ! function density ! getnewcut ! deform ! finda (matrix propagator energies from final state energies) ! checkcalc (the known integral) ! hrothgar ! functions feynmanf, feynman, feynman0f, feynman0 ! ! Common blocks (used for data that does not change): ! ! Set in main program beowulf: ! common /sizes/ nloops1,nprops1,nverts1,cutmax1, & ! nloops2,nprops2,nverts2,cutmax2 ! common /limits/ badnesslimit,cancellimit,thrustcut ! common /maxtime/ timelimit ! common /programmode/ mode ! common /whichgraphs/ usegraph ! common /montecarlo/groupsize,groupsizegraph,groupsizetotal ! common /calculook/ report,details ! common /colorfactors/ nc,nf,cf ! common /gaugechoice/ gauge ! common /physicsdata/ alphasofmz,sinsqthetaw,mz,widthz,externalrts ! common /deformscales/deformalpha,deformbeta,deformgamma ! common /smearparms/ smearfctr,lowpwr,highpwr ! common /renormalize/ muoverrts ! ! Set in graphcountinit (called from main program beowulf): ! common /graphcounts/ numberofgraphs,numberofcuts,numberofmaps ! ! Here is a listing of the subroutines and functions included. ! ! Subroutines and functions in beowulf file: ! module beowulf_parameters ! program beowulf ! subroutine hrothgar ! function calsthrust ! function calstmoments ! function thrust ! function thrustdist ! function kn0 ! function kn ! function calsymoments ! subroutine combinejets ! function cals0 ! subroutine daytime ! subroutine timing ! subroutine getcounts ! subroutine getcountsf ! ! Subroutines and functions in beowulfsubs file: ! subroutine version ! subroutine reno ! subroutine getnewgraph ! subroutine getnewflavorset ! subroutine finda ! function propsign ! subroutine findtypes ! subroutine newpoint ! subroutine checkpoint ! subroutine axes ! subroutine calculate ! subroutine checkcalc ! function density ! subroutine choose3 ! function rho3 ! subroutine choose2to2s ! function rho2to2s ! subroutine choose2to2t ! function rho2to2t ! subroutine choose2to3d ! function rho2to3d ! subroutine choose2to3e ! function rho2to3e ! subroutine choose2to1 ! function rho2to1 ! subroutine graphcountinit ! subroutine getnewcut ! subroutine deform ! function delta ! function smear ! function feynmanf ! function feynman0f ! subroutine twopointgf ! subroutine twopointqf ! subroutine vertexf ! subroutine twopt2f ! function feynman ! function feynman0 ! subroutine twopointg ! subroutine twopointq ! subroutine vertex ! subroutine twopt2 ! subroutine epsilon4 ! subroutine epsilon3 ! subroutine epsilon2 ! subroutine epsilon2n ! subroutine epsilon1n ! subroutine epsilont2 ! subroutine epsilont1 ! subroutine diagnostic ! function alpi ! function xxreal ! function xximag ! function complexsqrt ! function factorial ! function sinhinv ! function expm1 ! function sqrtm1 ! function logsqint ! function invlogsqint ! subroutine randominit ! subroutine newran ! function random ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! ! A brief introduction to the variables used: ! ! Size of the calculation: ! NLOOPS = number of loops (in cut photon self energy graph). ! NPROPS = number of propagators in graph, = 3 * NLOOPS - 1. ! NVERTS = number of vertices in graph, = 2 * NLOOPS. ! CUTMAX = NLOOPS + 1 ! = maximum number of cut propagators; ! = number of independent loop momenta needed to determine the ! propagator momenta, counting the virtual photon momentum. ! The current program is restricted to 0 and 1 virtual loops. ! ! Labels: ! L = index of loop momenta, L = 0,1,...,NLOOPS. ! L = 0 normally denontes the virtual photon momentum. ! P = index of propagator, P = 0,1,...,NPROPS. ! P = 0 denotes the virtual photon momentum. ! Q(L) = index P of propagator carrying the Lth loop momentum. ! V = index of vertices, V = 1,...,NVERTS ! ! Momentum variables (MU = 0,1,2,3): ! K(P,MU) = Momentum of Pth propagator. ! For P = 0, this is the virtual photon momentum: ! K(0,MU) = 0 for MU = 1, 2, 3 while K(0,0) = RTS. ! ABSK(P) = Square of the three momentum of Pth propagator. ! KINLOOP(J,MU) = K(LOOPINDEX(J),MU) = momenta of loop propagators. ! KCUT(I,MU) = K(CUTINDEX(I),MU) = momenta of cut propagators. ! K(Q(L),MU) = Lth loop momentum, L = 0,...,NLOOPS; ! KC(P,MU) = complex propagator momenta. ! A(P,L) = Matrix relating propagator momenta to loop momenta. ! K(P,MU) = SUM_{L=0}^{NLOOPS} A(P,L) K(Q(L),MU) ! kf(i,mu) = momenta of final state particles, i = 1,nparticles ! equal to the kcut(i,mu) if no showering. ! ! Variables from NEWGRAPH: ! VRTX(P,I) = Index of vertex at beginning (i= 1) and end (I = 2) of ! of propagator P. Specifies the supergraph. ! PROP(V,I) = Index of Ith propagator attached to vertex V, I = 1,2,3. ! Also specifies the supergraph. ! ! Variables associated with NEWPOINT and FINDTYPES: ! NMAPS = Number of different maps from random x's to momenta. ! MAPNMUMBER = Number labelling a certain map. ! QS(MAPNUMBER,II) = Label of the IIth propagator that is special ! in map number MAPNUMBER. ! QSIGNS(MAPNUMBER,II) = sign needed to relate the conventional ! direction of the propagator to that in an elementary scattering ! MAPTYPES(MAPNUMBER) = T2TO3, T2TO2T, T2TO2S, T2TO1. ! ! JACNEWPOINT =1/DENSITY(GRAPHNUMBER,K,QS,QSIGNS,MAPTYPES,NMAPS,ORDER) ! = Jacobian for loop momenta L. ! ! Variables from NEWCUT: ! NEWCUTINIT: .TRUE. tells NEWCUT to initialize itself. ! NCUT = Number of cut propagators. ! ISIGN(P) = +1 if propagator P is left of cut, -1 if right, 0 if cut. ! CUTINDEX(I) = Index P of cut propagator I, I = 1,...,CUTMAX. ! CUTSIGN(I) = Sign of cut propagator I I = 1,...,CUTMAX. ! (+1 if K(P,0) >0 for cut propagator.) ! LEFTLOOP = True iff there is a virtual loop to the left of the cut. ! RIGHTLOOP = True iff there is a virtual loop to the right of the cut. ! NINLOOP = Number of propagators in loop. ! LOOPINDEX(NP) = Index P of NPth propagator around the loop. ! LOOPSIGN(NP) = 1 if propagator direction is same as loop direction. ! -1 if direction is opposite to loop direction. ! NP = JCUT: Propagator cut by loopcut. ! CUTFOUND: .TRUE. if NEWCUT found a new cut. ! ! In RENO we use CUTINDEX to define CUT(P) = True if propagator ! P is cut. ! ! Solving for the propagator energies: ! For NCUT = CUTMAX, cut propagators are P = CUTINDEX(I). ! with direction of positive energy given by CUTSIGN(I). ! For NCUT = CUTMAX - 1, we define a "loopcut" on the propagator ! numbered JCUT in order around the loop, 1.LE.JCUT.LE.NINLOOP: ! CUTINDEX(CUTMAX) = LOOPINDEX(JCUT) and ! CUTSIGN(CUTMAX) = LOOPSIGN(JCUT). ! Energies of cut propagators are ! E(I-1) = K(CUTINDEX(I),0) for I = 1,...,CUTMAX. ! and are determined from ! E(I-1) = CUTSIGN(I) * SQRT( Sum_J [ K(CUTINDEX(I),J)**2 ] ). ! This gives energies E(L) for L = 0,...,NLOOPS. We consider the ! propagators designated by QE(L) = CUTINDEX(L+1) as independent ! and generate the matrix AE(P,L) that gives the propagator energies ! in terms of these independent momenta. This gives the propagator ! energies. ! ! Contour deformation: ! NEWKINLOOP(MU) = addition to the momentum going around the loop ! caused by deforming the contour. We have ! Im[ KC(LOOPINDEX(J,MU)) ] = LOOPSIGN(LOOPINDEX(J)) ! * Im[ NEWKINLOOP(J,MU) ] for MU = 1,2,3. ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine renoinit(sums) use beowulf_parameters use beowulf_structures implicit none ! Out: type(renoresults) :: sums ! integer :: i ! sums%i = 0.0d0 sums%r = 0.0d0 sums%bis = 0.0d0 sums%chkr = 0.0d0 sums%chki = 0.0d0 sums%rsq = 0.0d0 sums%isq = 0.0d0 sums%bissq = 0.0d0 sums%chkrsq = 0.0d0 sums%chkisq = 0.0d0 sums%included = 0 sums%extracount = 0 sums%omitted = 0.0d0 sums%nvalpt = (/(0,i = -9,6)/) sums%valptmax = 0.0d0 sums%fluct = reshape((/(0,i=1,maxgraphs*maxmaps)/),(/maxgraphs,maxmaps/)) ! sums%badpointinfo does not need to be initialized ! end subroutine renoinit ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine reno(sums) ! use beowulf_parameters use beowulf_structures implicit none ! In and Out: type(renoresults) :: sums ! ! When called many times, computes the cross section integral ! by Monte Carlo integration. ! ! Latest revision 16 November 2003 ! integer :: nloops,nprops integer :: nloops1,nprops1,nverts1,cutmax1 integer :: nloops2,nprops2,nverts2,cutmax2 common /sizes/ nloops1,nprops1,nverts1,cutmax1, & nloops2,nprops2,nverts2,cutmax2 real(kind=dbl) :: badnesslimit,cancellimit,thrustcut common /limits/ badnesslimit,cancellimit,thrustcut real(kind=dbl) :: timelimit common /maxtime/ timelimit ! What the program should do character(len=8) :: mode common /programmode/ mode ! Graphs to include logical :: usegraph(maxgraphs) common /whichgraphs/ usegraph ! How many graphs and how many maps for each: integer :: numberofgraphs integer :: numberofmaps(maxgraphs) common /graphcounts/ numberofgraphs,numberofmaps ! Order of perturbation theory integer :: order ! Momenta: real(kind=dbl) :: k(0:3*size-1,0:3),absk(0:3*size-1) ! Matrices: integer :: a(0:3*size-1,0:size) ! A blank shower for hrothgar: type(showerlist) :: blankshower ! NEWGRAPH variables: type(graphstructure) :: graph logical :: graphfound integer :: graphnumber ! FINDA variable: logical :: qok ! MAP variables: integer :: nmaps,mapnumber integer :: qs(maxmaps,0:size),qsigns(maxmaps,0:size) integer :: q(0:size),qsign(0:size) character(len=6) :: maptypes(maxmaps) character(len=6) :: maptype ! Variable from CHECKPOINT: real(kind=dbl) :: badness ! Problem report from NEWPOINT logical :: badnewpoint ! Logical variables to tell how to treat point: logical :: xtrapointq, badpointq ! Functions: real(kind=dbl) :: xxreal,xximag ! Index variables: integer :: l,p,mu ! Reno size and counting variables: integer :: groupsize(maxgraphs,maxmaps) integer :: groupsizegraph(maxgraphs) integer :: groupsizetotal common /montecarlo/groupsize,groupsizegraph,groupsizetotal integer :: point ! The standard integrals to be calculated complex(kind=dbl) :: integral,integralchk real(kind=dbl) :: integralbis ! Calculate variables: complex(kind=dbl) :: value,valuechk real(kind=dbl) :: maxpart real(kind=dbl) :: valpt,logvalpt ! This gives us access to the random number generator's internal workings. ! We do not change these numbers in subroutine reno! But we save the ! information needed to set the random number generator properly in ! subroutine diagnostic. integer :: jrandom,irandom(250) real(kind=dbl) :: rrandom(250) common/rando/ rrandom,jrandom,irandom real(kind=dbl) :: temp integer :: tempjr,tempir(250) real(kind=dbl) :: temprr(250) real(kind=dbl) :: random ! save ! Save everything ! !------------------------------ Begin ---------------------------------- ! ! Initialize integrals for this reno group. ! integral = (0.0d0,0.0d0) integralbis = 0.0d0 integralchk = (0.0d0,0.0d0) ! ! Call Hrothgar to tell him to that we are starting a new group. ! call hrothgar(blankshower,1.0d0,1,'startgroup') ! ! Get a new graph. The starting value for GRAPHNUMBER depends ! on the order we want. Graphs for ORDER = 2 are numbered 1,...,10 ! and those for ORDER = 1 are numbered 11,12. For MODE = nlo, ! we will do the ORDER = 2 graphs first, then continue with the ! first order graphs. Thus we wait for NEWGRAPH to return ! GRAPHFOUND = false, then reset it (see ELSE part of IF(GRAPHFOUND)). ! graphfound = .true. ! IF (mode.EQ.'born ') THEN order = 1 nloops = nloops1 nprops = nprops1 graphnumber = 10 ELSE IF (mode.EQ.'hocoef ') THEN order = 2 nloops = nloops2 nprops = nprops2 graphnumber = 0 ELSE IF (mode.EQ.'nlo ') THEN order = 2 nloops = nloops2 nprops = nprops2 graphnumber = 0 ELSE IF ((mode.EQ.'showerI ').OR.(mode.EQ.'showerII')) THEN order = 2 nloops = nloops2 nprops = nprops2 graphnumber = 0 ELSE write(nout,*)'Not programmed for this mode. ',mode stop END IF ! DO call getnewgraph(order,graph,graphfound) IF (.NOT.graphfound) THEN IF (((mode.EQ.'nlo ').OR.(mode.EQ.'showerI ').OR.(mode.EQ.'showerII')) & .AND.(order.EQ.2)) THEN ! We have finished with the order 2 graphs and should do order 1. graphfound = .true. order = 1 nloops = nloops1 nprops = nprops1 call getnewgraph(order,graph,graphfound) IF (.NOT.graphfound) THEN write(nout,*)'What, graph 11 disappeared?' STOP END IF ELSE EXIT ! We have finished with all of the graphs. END IF END IF ! We have a graph to work with. graphnumber = graphnumber + 1 IF (.NOT.(graphnumber.EQ.graph%graphnumber)) THEN write(nout,*)'Oops, we lost track of graphnumber.' STOP END IF ! ! Check if we were supposed to use this graph (USEGRAPH is ! set in the main program.) ! IF (usegraph(graphnumber)) THEN ! ! Calculate number of maps NMAPS, index arrays QS, ! types MAPTYPES, and signs QSIGNS associated with the maps. ! call findtypes(graph,nmaps,qs,qsigns,maptypes) ! ! Loop over choices of maps from x's to loop momenta. ! DO mapnumber = 1,nmaps ! maptype = maptypes(mapnumber) DO l = 0,nloops q(l) = qs(mapnumber,l) qsign(l) = qsigns(mapnumber,l) END DO ! call finda(graph%vrtx,q,nloops,order,a,qok) ! ! Loop over Reno points within a group. ! DO point = 1,groupsize(graphnumber,mapnumber) ! ! Call Hrothgar to tell him that we are starting a new point. ! call hrothgar(blankshower,1.0d0,1,'startpoint') ! ! Get a new point. Check on its badness. If it is too bad, ! or if NEWPOINT reported a problem, we omit the point after ! notifying Hrothgar. ! badpointq = .false. xtrapointq = .false. call newpoint(a,qsign,maptype,order,k,absk,badnewpoint) IF (badnewpoint) THEN call hrothgar(blankshower,1.0d0,1,'badpoint ') badpointq = .true. END IF call checkpoint(k,absk,graph%prop,order,badness) IF (badness.GT.100*badnesslimit) THEN call hrothgar(blankshower,1.0d0,1,'badpoint ') badpointq = .true. ELSE IF (badness.GT.badnesslimit) THEN call hrothgar(blankshower,1.0d0,1,'xtrapoint ') xtrapointq = .true. END IF ! ! If the point is not too bad, we can call CALCULATE. ! The final state momenta found, along with the corresponding ! weights, are reported to Hrothgar by CACULATE. ! Then call Hrothgar to tell him that we are done with this point. ! IF (.NOT.badpointq) THEN temp = random(1) tempjr = jrandom ! This saves the state of the random temprr = rrandom ! number generator just before tempir = irandom ! calculate was called. call calculate(graph,k,absk, & qs,qsigns,maptypes,nmaps,value,maxpart,valuechk) END IF ! ! Add contribution from this point to integral. ! We count the point if Maxvalue/|Value| < Cancellimit. ! IF (.NOT.badpointq) THEN IF ( maxpart .GT. 100*cancellimit*abs(xxreal(value)) ) THEN call hrothgar(blankshower,1.0d0,1,'badpoint ') badpointq = .true. ELSE IF ( maxpart .GT. cancellimit*abs(xxreal(value)) ) THEN call hrothgar(blankshower,1.0d0,1,'xtrapoint ') xtrapointq = .true. END IF END IF ! IF ( (.NOT.badpointq).AND.(.NOT.xtrapointq) ) THEN integral = integral + value sums%fluct(graphnumber,mapnumber) = sums%fluct(graphnumber,mapnumber) & + xxreal(value)**2/groupsize(graphnumber,mapnumber) integralchk = integralchk + valuechk sums%included = sums%included + 1 ! ! For diagnostic purposes, we need VALPT, the contribution to ! the integral being calculated from this point, normalized such ! that the integral is the sum over all points chosen of VALPT ! divided by the total number of points, NRENO * GROUPSIZETOTAL. ! valpt = abs(xxreal(value))*groupsizetotal IF (valpt .GT. 0.0d0) then logvalpt = log10(valpt) ELSE logvalpt = -100.0d0 END IF DO l = -9,6 IF((logvalpt.GE.l).AND.(logvalpt.LT.(l+1))) THEN sums%nvalpt(l) = sums%nvalpt(l) + 1 END IF END DO IF (valpt.GT.sums%valptmax) THEN ! ! We save the information about the worst point for later analysis. ! sums%valptmax = valpt DO mu = 0,3 sums%badpointinfo%k(0,mu) = 0.0d0 END DO DO p = 1,nprops sums%badpointinfo%k(p,0) = 0.0d0 DO mu = 1,3 sums%badpointinfo%k(p,mu) = k(p,mu) END DO END DO sums%badpointinfo%graphnumber = graphnumber sums%badpointinfo%mapnumber = mapnumber sums%badpointinfo%jr = tempjr ! This saves the state of the random sums%badpointinfo%rr = temprr ! number generator just before calculate sums%badpointinfo%ir = tempir ! was called for the bad point. ! END IF ELSE IF ((.NOT.badpointq).AND.(xtrapointq) ) THEN ! ! For points that are 'extra', we include the value of ! the integrand in the INTEGRALBIS, which will provide an estimate ! or the effect of the cutoffs. ! integralbis = integralbis + xxreal(value) sums%extracount = sums%extracount + 1 ! ELSE sums%omitted = sums%omitted + 1 END IF ! ! End of loop over POINT. ! call hrothgar(blankshower,1.0d0,1,'pointdone ') END DO ! ! End of loop over MAPNUMBER. ! END DO ! ! End for IF (USEGRAPH(GRAPHNUMBER)) THEN ! END IF ! ! End of loop DO ... call getnewgraph() with EXIT if we ran out of graphs. ! END DO ! ! Call Hrothgar to tell him that we are done with this reno group. ! call hrothgar(blankshower,1.0d0,1,'groupdone ') ! ! Add results from this group to the SUM variables. ! sums%r = sums%r + xxreal(integral) sums%i = sums%i + xximag(integral) sums%bis = sums%bis + integralbis sums%chkr = sums%chkr + xxreal(integralchk) sums%chki = sums%chki + xximag(integralchk) ! sums%rsq = sums%rsq + xxreal(integral)**2 sums%isq = sums%isq + xximag(integral)**2 sums%bissq = sums%bissq + integralbis**2 sums%chkrsq = sums%chkrsq + xxreal(integralchk)**2 sums%chkisq = sums%chkisq + xximag(integralchk)**2 ! RETURN END subroutine reno ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! Subroutine GETNEWGRAPH !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine getnewgraph(orderwanted,graph,graphfound) ! use beowulf_parameters use beowulf_structures implicit none ! In: integer :: orderwanted ! Out: type(graphstructure) :: graph logical :: graphfound ! ! VRTX(P,I) = Index of vertex at beginning (i= 1) and end (I = 2) of ! of propagator P. Specifies the supergraph. ! PROP(V,I) = Index of Ith propagator attached to vertex V, I = 1,2,3. ! Also specifies the supergraph. ! Use ! C(V,I) = Index of Ith vertex connected to vertex V. ! V = 1,...,NVERTS; I =1,2,3; C(V,I) = 1,...,NVERTS and -1,-2. ! Here C(V,1).LE.C(V,2).LE.C(V,3). ! This is the fundamental specification of the supergraph. ! Vertex 1 is automatically connected to the photon "-1":C(1,1) = -1. ! Vertex 2 is automatically connected to the photon "-2":C(2,1) = -2. ! Connections to the external boson: ! In C(V,I) we assign the first connection of vertex 1 to be vertex "-1" ! while the first connection of vertex 2 is vertex "-2." This numbering ! is convenient for working out C(V,I). In reporting the results, ! however, we label the external boson with propagator 0, so that ! PROP(1,1) = PROP(2,1) = 0. Then propagator 0 attaches to vertices ! 1 and 2: VERT(0,1) = 2, VERT(0,2) = 1. ! ! 13 June 2002 ! integer, dimension(2*size,3) :: c ! integer :: nprops,nverts integer :: nloops1,nprops1,nverts1,cutmax1 integer :: nloops2,nprops2,nverts2,cutmax2 common /sizes/ nloops1,nprops1,nverts1,cutmax1, & nloops2,nprops2,nverts2,cutmax2 integer, save :: graphnumber = 0 integer :: vv,p,i,va,vb integer, dimension(2*size) :: nused ! ! Initializations. ! IF (orderwanted.EQ.2) THEN IF (graphnumber.GE.10) THEN graphfound = .false. graphnumber = 0 RETURN END IF nprops = nprops2 nverts = nverts2 ELSE IF (orderwanted.EQ.1) THEN IF (graphnumber.GE.12) THEN graphfound = .false. graphnumber = 0 RETURN ELSE IF (graphnumber.LT.10) THEN graphnumber = 10 ! Starting value for order = 1. END IF nprops = nprops1 nverts = nverts1 ELSE write(nout,*)'Order must be 1 or 2.',orderwanted END IF graphnumber = graphnumber + 1 ! ! We know which graph, and everything is OK, so we proceed. ! graphfound = .true. graph%graphnumber = graphnumber graph%order = orderwanted ! ! Data are from getnewgraph.nb ! IF (graphnumber.EQ.1) THEN c(1,:) = (/-1,2,3/) c(2,:) = (/-2,1,4/) c(3,:) = (/1,4,5/) c(4,:) = (/2,3,6/) c(5,:) = (/3,6,6/) c(6,:) = (/4,5,5/) ELSE IF (graphnumber.EQ.2) THEN c(1,:) = (/-1,2,3/) c(2,:) = (/-2,1,4/) c(3,:) = (/1,5,5/) c(4,:) = (/2,6,6/) c(5,:) = (/3,3,6/) c(6,:) = (/4,4,5/) ELSE IF (graphnumber.EQ.3) THEN c(1,:) = (/-1,2,3/) c(2,:) = (/-2,1,4/) c(3,:) = (/1,5,6/) c(4,:) = (/2,5,6/) c(5,:) = (/3,4,6/) c(6,:) = (/3,4,5/) ELSE IF (graphnumber.EQ.4) THEN c(1,:) = (/-1,3,4/) c(2,:) = (/-2,3,4/) c(3,:) = (/1,2,5/) c(4,:) = (/1,2,6/) c(5,:) = (/3,6,6/) c(6,:) = (/4,5,5/) ELSE IF (graphnumber.EQ.5) THEN c(1,:) = (/-1,3,4/) c(2,:) = (/-2,3,5/) c(3,:) = (/1,2,4/) c(4,:) = (/1,3,6/) c(5,:) = (/2,6,6/) c(6,:) = (/4,5,5/) ELSE IF (graphnumber.EQ.6) THEN c(1,:) = (/-1,3,4/) c(2,:) = (/-2,3,5/) c(3,:) = (/1,2,5/) c(4,:) = (/1,6,6/) c(5,:) = (/2,3,6/) c(6,:) = (/4,4,5/) ELSE IF (graphnumber.EQ.7) THEN c(1,:) = (/-1,3,4/) c(2,:) = (/-2,3,5/) c(3,:) = (/1,2,6/) c(4,:) = (/1,5,6/) c(5,:) = (/2,4,6/) c(6,:) = (/3,4,5/) ELSE IF (graphnumber.EQ.8) THEN c(1,:) = (/-1,3,4/) c(2,:) = (/-2,5,6/) c(3,:) = (/1,4,5/) c(4,:) = (/1,3,6/) c(5,:) = (/2,3,6/) c(6,:) = (/2,4,5/) ELSE IF (graphnumber.EQ.9) THEN c(1,:) = (/-1,3,4/) c(2,:) = (/-2,5,6/) c(3,:) = (/1,5,5/) c(4,:) = (/1,6,6/) c(5,:) = (/2,3,3/) c(6,:) = (/2,4,4/) ELSE IF (graphnumber.EQ.10) THEN c(1,:) = (/-1,3,4/) c(2,:) = (/-2,5,6/) c(3,:) = (/1,5,6/) c(4,:) = (/1,5,6/) c(5,:) = (/2,3,4/) c(6,:) = (/2,3,4/) ELSE IF (graphnumber.EQ.11) THEN c(1,:) = (/-1,2,3/) c(2,:) = (/-2,1,4/) c(3,:) = (/1,4,4/) c(4,:) = (/2,3,3/) ELSE IF (graphnumber.EQ.12) THEN c(1,:) = (/-1,3,4/) c(2,:) = (/-2,3,4/) c(3,:) = (/1,2,4/) c(4,:) = (/1,2,3/) END IF ! ! Exit. We translate the results for C(V,I) into VRTX(P,I), I = 1,2, ! and PROP(P,I), I = 1,2,3. Here NUSED(V) denotes how many propagators ! we have so far assigned connecting to vertex V. ! DO vv = 1,nverts nused(vv) = 0 END DO graph%vrtx(0,1) = 2 graph%vrtx(0,2) = 1 graph%prop(1,1) = 0 nused(1) = 1 graph%prop(2,1) = 0 nused(2) = 1 p = 1 DO vv = 1,nverts DO i = 1,3 IF (c(vv,i).GT.vv) THEN va = vv vb = c(vv,i) graph%vrtx(p,1) = va nused(va) = nused(va) + 1 graph%prop(va,nused(va)) = p graph%vrtx(p,2) = vb nused(vb) = nused(vb) + 1 graph%prop(vb,nused(vb)) = p p = p+1 END IF END DO END DO IF (p.NE.nprops+1) THEN write(nout,*)'Snafu in getnewgraph.',p-1,nprops write(nout,*)'Graphnumber',graph%graphnumber stop END IF DO vv = 1,nverts IF (nused(vv).NE.3) THEN write(nout,*)'Problem in getnewwgraph.',vv,nused(vv) write(nout,*)'Graphnumber',graph%graphnumber stop END IF END DO ! ! OK. We are ready to return. ! RETURN END subroutine getnewgraph ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine getnewflavorset(graphnumber,flavorset,flavorsetfound) ! use beowulf_parameters use beowulf_structures implicit none ! In: integer :: graphnumber ! Out: type (flavorstructure) :: flavorset logical :: flavorsetfound ! ! Latest revisions: ! 28 June 2002 ! 8 November 2003 ! 27 August 2004 ! real(kind=dbl) :: nc,nf,cf common /colorfactors/ nc,nf,cf ! Physics data: real(kind=dbl) :: alphasofmz,sinsqthetaw,mz,widthz,externalrts common /physicsdata/ alphasofmz,sinsqthetaw,mz,widthz,externalrts ! integer, save :: flavorsetnumber = 0 ! character(len=6), parameter :: gluon = 'gluon' character(len=6), parameter :: quark = 'quark' character(len=6), parameter :: qbar = 'qbar ' character(len=6), parameter :: none = 'none ' ! real(kind=dbl), dimension(6), parameter :: & charge = (/0.666666666666667d0,-0.333333333333333d0, & -0.333333333333333d0,0.666666666666667d0, & -0.333333333333333d0,0.666666666666667d0/) real(kind=dbl), dimension(6), parameter :: & t3 = (/ 0.5d0,-0.5d0, & -0.5d0, 0.5d0, & -0.5d0, 0.5d0/) real(kind=dbl) :: s,cossqthetaw,denomz,kappaz,ve,ae real(kind=dbl) :: effchrgsqtot,a,v,q real(kind=dbl), dimension(6) :: effchrgsq ! logical, save :: initialize = .true. integer, save :: nflavors integer :: p,i,flavorindex real(kind=dbl) :: temp,r,random real(kind=dbl), dimension(6), save :: probabilitysum integer :: quarkflavor character(len=5) :: thetype ! flavorsetfound = .false. flavorsetnumber = flavorsetnumber + 1 ! IF (graphnumber.EQ.1) THEN IF (flavorsetnumber.EQ.1) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,quark,qbar,gluon,gluon,qbar,quark/) ELSE IF(flavorsetnumber.EQ.2) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,quark,qbar,gluon,gluon,gluon,gluon/) ELSE IF(flavorsetnumber.EQ.3) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,quark,gluon,qbar,quark,gluon,qbar/) ELSE IF(flavorsetnumber.EQ.4) THEN flavorsetfound = .true. flavorset%type = & (/qbar,quark,qbar,quark,gluon,gluon,qbar,quark/) ELSE IF(flavorsetnumber.EQ.5) THEN flavorsetfound = .true. flavorset%type = & (/qbar,quark,qbar,quark,gluon,gluon,gluon,gluon/) ELSE IF(flavorsetnumber.EQ.6) THEN flavorsetfound = .true. flavorset%type = & (/qbar,quark,qbar,gluon,quark,qbar,gluon,quark/) END IF ELSE IF(graphnumber.EQ.2) THEN IF (flavorsetnumber.EQ.1) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,quark,gluon,qbar,gluon,quark,qbar/) ELSE IF(flavorsetnumber.EQ.2) THEN flavorsetfound = .true. flavorset%type = & (/qbar,quark,qbar,gluon,quark,gluon,qbar,quark/) END IF ELSE IF(graphnumber.EQ.3) THEN IF (flavorsetnumber.EQ.1) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,quark,gluon,qbar,quark,gluon,quark/) ELSE IF(flavorsetnumber.EQ.2) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,quark,gluon,qbar,gluon,quark,gluon/) ELSE IF(flavorsetnumber.EQ.3) THEN flavorsetfound = .true. flavorset%type = & (/qbar,quark,qbar,gluon,quark,qbar,gluon,qbar/) ELSE IF(flavorsetnumber.EQ.4) THEN flavorsetfound = .true. flavorset%type = & (/qbar,quark,qbar,gluon,quark,gluon,qbar,gluon/) END IF ELSE IF(graphnumber.EQ.4) THEN IF (flavorsetnumber.EQ.1) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,qbar,quark,gluon,gluon,qbar,quark/) ELSE IF(flavorsetnumber.EQ.2) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,qbar,quark,gluon,gluon,gluon,gluon/) END IF ELSE IF(graphnumber.EQ.5) THEN IF (flavorsetnumber.EQ.1) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,qbar,quark,gluon,qbar,gluon,quark/) ELSE IF(flavorsetnumber.EQ.2) THEN flavorsetfound = .true. flavorset%type = & (/qbar,quark,quark,qbar,gluon,quark,gluon,qbar/) END IF ELSE IF(graphnumber.EQ.6) THEN IF (flavorsetnumber.EQ.1) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,qbar,quark,gluon,gluon,qbar,quark/) ELSE IF(flavorsetnumber.EQ.2) THEN flavorsetfound = .true. flavorset%type = & (/qbar,quark,quark,qbar,gluon,gluon,quark,qbar/) END IF ELSE IF(graphnumber.EQ.7) THEN IF (flavorsetnumber.EQ.1) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,qbar,quark,gluon,qbar,gluon,gluon/) ELSE IF(flavorsetnumber.EQ.2) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,qbar,quark,gluon,gluon,qbar,quark/) ELSE IF(flavorsetnumber.EQ.3) THEN flavorsetfound = .true. flavorset%type = & (/qbar,quark,quark,qbar,gluon,quark,gluon,gluon/) ELSE IF(flavorsetnumber.EQ.4) THEN flavorsetfound = .true. flavorset%type = & (/qbar,quark,quark,qbar,gluon,gluon,quark,qbar/) END IF ELSE IF(graphnumber.EQ.8) THEN IF (flavorsetnumber.EQ.1) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,quark,qbar,quark,gluon,gluon,quark/) ELSE IF(flavorsetnumber.EQ.2) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,qbar,quark,quark,gluon,gluon,qbar/) ELSE IF(flavorsetnumber.EQ.3) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,qbar,quark,gluon,quark,qbar,gluon/) END IF ELSE IF(graphnumber.EQ.9) THEN IF (flavorsetnumber.EQ.1) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,qbar,quark,gluon,quark,gluon,qbar/) END IF ELSE IF(graphnumber.EQ.10) THEN IF (flavorsetnumber.EQ.1) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,quark,qbar,gluon,quark,qbar,gluon/) END IF ELSE IF(graphnumber.EQ.11) THEN IF (flavorsetnumber.EQ.1) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,quark,gluon,qbar,none,none,none/) ELSE IF(flavorsetnumber.EQ.2) THEN flavorsetfound = .true. flavorset%type = & (/qbar,quark,qbar,gluon,quark,none,none,none/) END IF ELSE IF(graphnumber.EQ.12) THEN IF (flavorsetnumber.EQ.1) THEN flavorsetfound = .true. flavorset%type = & (/quark,qbar,qbar,quark,gluon,none,none,none/) END IF ELSE write(nout,*)'Bad graphnumber in getnewflavorset.' STOP END IF ! IF (flavorsetfound) THEN flavorset%setnumber = flavorsetnumber ELSE flavorsetnumber = 0 END IF ! IF (flavorsetfound) THEN ! ! Choose flavors for the quarks. ! IF (initialize) THEN initialize = .false. s = externalrts**2 cossqthetaw = 1.0d0 - sinsqthetaw denomz = (s - mz**2)**2 + mz**2 * widthz**2 kappaz = 1/4.0d0/sinsqthetaw/cossqthetaw ve = -0.5d0 + 2.0d0 * sinsqthetaw ae = -0.5d0 nflavors = nint(nf) effchrgsqtot = 0.0d0 DO i = 1,nflavors a = t3(i) v = t3(i) - 2.0d0*charge(i)*sinsqthetaw q = charge(i) effchrgsq(i) = & denomz*q**2 & + kappaz**2 * s**2 * (ve**2 + ae**2)*(v**2 + a**2) & - 2.0d0*kappaz*s*(s-mz**2)*ve*q*v effchrgsqtot = effchrgsqtot + effchrgsq(i) END DO temp = 0.0d0 DO i = 1,nflavors temp = temp + effchrgsq(i)/effchrgsqtot probabilitysum(i) = temp ! we save this END DO END IF ! r = random(1) quarkflavor = 1 DO i = 1,nflavors - 1 IF (r.GT.probabilitysum(i)) THEN quarkflavor = quarkflavor + 1 END IF END DO DO p = 1,8 thetype = flavorset%type(p) IF (thetype.EQ.gluon) THEN flavorset%flavor(p) = 0 ELSE IF (thetype.EQ.quark) THEN flavorset%flavor(p) = quarkflavor ELSE IF (thetype.EQ.qbar) THEN flavorset%flavor(p) = - quarkflavor ELSE IF (thetype.EQ.none) THEN flavorset%flavor(p) = - 777 ! For p = 6,7,8 for Born graphs. ELSE write(nout,*)'Problem with parton type in getnewflavorset. ',thetype STOP END IF END DO ! ! For graphs 1 and 4, we have a separate quark loop comprising propagators ! 7 and 8, for which the flavor should be chosen at random. ! IF ((graphnumber.EQ.1).OR.(graphnumber.EQ.4)) THEN r = random(1) flavorindex = ceiling(random(1)*nf) ! choose the flavor at random DO p = 7,8 thetype = flavorset%type(p) IF (thetype.EQ.gluon) THEN flavorset%flavor(p) = 0 ELSE IF (thetype.EQ.quark) THEN flavorset%flavor(p) = flavorindex ELSE IF (thetype.EQ.qbar) THEN flavorset%flavor(p) = - flavorindex ELSE write(nout,*)'Problem with parton type in getnewflavorset. ',thetype STOP END IF END DO END IF END IF ! (flavorsetfound) ! RETURN END subroutine getnewflavorset ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 subroutine getcolors(graphnumber,flavorsetnumber,cutnumber,strings) use beowulf_parameters implicit none ! In: integer :: graphnumber,flavorsetnumber,cutnumber ! Out: integer, dimension(8) :: strings ! ! For a given graphnumber,flavorsetnumber,cutnumber with the cut such that there ! are four particles in the final state, finds a color structure for the final ! state. The output is an eight component array of integers, strings, with ! index of string connected to quark index of parton 1 = strings(1) ! index of string connected to qbar index of parton 1 = strings(2) ! index of string connected to quark index of parton 2 = strings(3) ! index of string connected to qbar index of parton 2 = strings(4) ! index of string connected to quark index of parton 3 = strings(5) ! index of string connected to qbar index of parton 3 = strings(6) ! index of string connected to quark index of parton 4 = strings(7) ! index of string connected to qbar index of parton 4 = strings(8) ! ! 21 November 2003 ! type colorstructure integer :: n integer, dimension(8) :: cnct1 integer, dimension(8) :: cnct2 end type colorstructure type(colorstructure), dimension(10,6,3:9), save :: color ! indices denote color(graphnumber,flavorsetnumber,cutnumber) logical, save :: init = .true. real(kind=dbl) :: random integer :: x = -500 ! IF (init) THEN color(1,1,3)%n = 1 color(1,1,3)%cnct1 = (/1,x,x,2,x,1,2,x/) color(1,2,3)%n = 2 color(1,2,3)%cnct1 = (/2,1,3,2,x,3,1,x/) color(1,2,3)%cnct2 = (/3,2,2,1,x,3,1,x/) color(1,3,3)%n = 1 color(1,3,3)%cnct1 = (/x,3,3,2,2,1,1,x/) color(1,4,3)%n = 1 color(1,4,3)%cnct1 = (/1,x,x,2,2,x,x,1/) color(1,5,3)%n = 2 color(1,5,3)%cnct1 = (/2,1,3,2,1,x,x,3/) color(1,5,3)%cnct2 = (/3,2,2,1,1,x,x,3/) color(1,6,3)%n = 1 color(1,6,3)%cnct1 = (/1,x,2,1,3,2,x,3/) color(3,1,3)%n = 2 color(3,1,3)%cnct1 = (/x,3,2,1,3,2,1,x/) color(3,1,3)%cnct2 = (/x,3,3,2,2,1,1,x/) color(3,1,4)%n = 2 color(3,1,4)%cnct1 = (/1,x,x,1,x,2,2,x/) color(3,1,4)%cnct2 = (/1,x,x,2,x,1,2,x/) color(3,2,3)%n = 1 color(3,2,3)%cnct1 = (/3,2,x,3,2,1,1,x/) color(3,2,4)%n = 1 color(3,2,4)%cnct1 = (/3,2,2,1,x,3,1,x/) color(3,3,3)%n = 2 color(3,3,3)%cnct1 = (/1,x,2,1,3,2,x,3/) color(3,3,3)%cnct2 = (/1,x,3,2,2,1,x,3/) color(3,3,4)%n = 2 color(3,3,4)%cnct1 = (/x,1,1,x,2,x,x,2/) color(3,3,4)%cnct2 = (/x,2,1,x,2,x,x,1/) color(3,4,3)%n = 1 color(3,4,3)%cnct1 = (/2,1,1,x,3,2,x,3/) color(3,4,4)%n = 1 color(3,4,4)%cnct1 = (/2,1,3,2,1,x,x,3/) color(4,1,5)%n = 1 color(4,1,5)%cnct1 = (/x,2,1,x,x,1,2,x/) color(4,1,6)%n = 1 color(4,1,6)%cnct1 = (/1,x,x,2,2,x,x,1/) color(4,2,5)%n = 2 color(4,2,5)%cnct1 = (/2,1,3,2,x,3,1,x/) color(4,2,5)%cnct2 = (/3,2,2,1,x,3,1,x/) color(4,2,6)%n = 2 color(4,2,6)%cnct1 = (/2,1,3,2,1,x,x,3/) color(4,2,6)%cnct2 = (/3,2,2,1,1,x,x,3/) color(5,1,5)%n = 1 color(5,1,5)%cnct1 = (/x,3,3,2,2,1,1,x/) color(5,2,5)%n = 1 color(5,2,5)%cnct1 = (/1,x,2,1,3,2,x,3/) color(6,1,5)%n = 1 color(6,1,5)%cnct1 = (/x,3,3,2,2,1,1,x/) color(6,2,5)%n = 1 color(6,2,5)%cnct1 = (/1,x,2,1,3,2,x,3/) color(7,1,5)%n = 1 color(7,1,5)%cnct1 = (/3,2,2,1,x,3,1,x/) color(7,1,6)%n = 1 color(7,1,6)%cnct1 = (/3,2,x,3,2,1,1,x/) color(7,1,7)%n = 1 color(7,1,7)%cnct1 = (/2,1,3,2,1,x,x,3/) color(7,1,8)%n = 1 color(7,1,8)%cnct1 = (/3,2,x,3,2,1,1,x/) color(7,2,5)%n = 2 color(7,2,5)%cnct1 = (/1,x,x,1,x,2,2,x/) color(7,2,5)%cnct2 = (/1,x,x,2,x,1,2,x/) color(7,2,6)%n = 2 color(7,2,6)%cnct1 = (/x,3,2,1,3,2,1,x/) color(7,2,6)%cnct2 = (/x,3,3,2,2,1,1,x/) color(7,2,7)%n = 2 color(7,2,7)%cnct1 = (/x,1,1,x,2,x,x,2/) color(7,2,7)%cnct2 = (/x,2,1,x,2,x,x,1/) color(7,2,8)%n = 2 color(7,2,8)%cnct1 = (/x,3,2,1,3,2,1,x/) color(7,2,8)%cnct2 = (/x,3,3,2,2,1,1,x/) color(7,3,5)%n = 1 color(7,3,5)%cnct1 = (/2,1,3,2,1,x,x,3/) color(7,3,6)%n = 1 color(7,3,6)%cnct1 = (/2,1,1,x,3,2,x,3/) color(7,3,7)%n = 1 color(7,3,7)%cnct1 = (/3,2,2,1,x,3,1,x/) color(7,3,8)%n = 1 color(7,3,8)%cnct1 = (/2,1,1,x,3,2,x,3/) color(7,4,5)%n = 2 color(7,4,5)%cnct1 = (/x,1,1,x,2,x,x,2/) color(7,4,5)%cnct2 = (/x,2,1,x,2,x,x,1/) color(7,4,6)%n = 2 color(7,4,6)%cnct1 = (/1,x,2,1,3,2,x,3/) color(7,4,6)%cnct2 = (/1,x,3,2,2,1,x,3/) color(7,4,7)%n = 2 color(7,4,7)%cnct1 = (/1,x,x,1,x,2,2,x/) color(7,4,7)%cnct2 = (/1,x,x,2,x,1,2,x/) color(7,4,8)%n = 2 color(7,4,8)%cnct1 = (/1,x,2,1,3,2,x,3/) color(7,4,8)%cnct2 = (/1,x,3,2,2,1,x,3/) color(8,1,5)%n = 1 color(8,1,5)%cnct1 = (/x,2,x,1,1,x,2,x/) color(8,1,6)%n = 1 color(8,1,6)%cnct1 = (/1,x,2,x,x,2,x,1/) color(8,2,5)%n = 1 color(8,2,5)%cnct1 = (/1,x,x,1,x,2,2,x/) color(8,2,6)%n = 1 color(8,2,6)%cnct1 = (/x,1,1,x,2,x,x,2/) color(8,3,5)%n = 1 color(8,3,5)%cnct1 = (/3,2,2,1,x,3,1,x/) color(8,3,6)%n = 1 color(8,3,6)%cnct1 = (/2,1,3,2,1,x,x,3/) color(9,1,5)%n = 1 color(9,1,5)%cnct1 = (/x,3,3,2,1,x,2,1/) color(10,1,5)%n = 2 color(10,1,5)%cnct1 = (/2,1,3,2,x,3,1,x/) color(10,1,5)%cnct2 = (/3,2,2,1,x,3,1,x/) color(10,1,6)%n = 2 color(10,1,6)%cnct1 = (/x,1,x,2,1,x,2,x/) color(10,1,6)%cnct2 = (/x,2,x,1,1,x,2,x/) color(10,1,7)%n = 2 color(10,1,7)%cnct1 = (/1,x,2,x,x,1,x,2/) color(10,1,7)%cnct2 = (/1,x,2,x,x,2,x,1/) color(10,1,8)%n = 2 color(10,1,8)%cnct1 = (/2,1,3,2,1,x,x,3/) color(10,1,8)%cnct2 = (/3,2,2,1,1,x,x,3/) color(10,1,9)%n = 2 color(10,1,9)%cnct1 = (/2,1,x,3,1,x,3,2/) color(10,1,9)%cnct2 = (/3,2,x,3,1,x,2,1/) END IF ! (init) ! IF (color(graphnumber,flavorsetnumber,cutnumber)%n.EQ.2) THEN IF (random(1).GT.0.5d0) THEN strings = color(graphnumber,flavorsetnumber,cutnumber)%cnct1 + 500 ELSE strings = color(graphnumber,flavorsetnumber,cutnumber)%cnct2 + 500 END IF ELSE strings = color(graphnumber,flavorsetnumber,cutnumber)%cnct1 + 500 END IF RETURN END ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine finda(vrtx,q,nq,order,a,qok) ! use beowulf_parameters implicit none ! In: integer :: vrtx(0:3*size-1,2),q(0:size),nq,order ! Out: integer :: a(0:3*size-1,0:size) logical :: qok ! ! Finds matrix A relating propagator momenta to loop momenta. ! ! VRTX(P,N) specifies the graph considered ! Q(L) specifies the propagators to be considered independent ! NQ specifies how many entries of Q should be considered ! NQ = NLOOPS all the entries in Q should be considered. ! If Q(0),Q(1),...,Q(NLOOPS) are independent then ! FINDA generates the matrix A and sets QOK = .TRUE. ! Otherwise the generation of A fails and QOK = .FALSE. ! NQ < NLOOPS only first NQ entries in Q should be considered. ! If Q(0),Q(1),...,Q(NQ) are independent then ! FINDA sets QOK = .TRUE. ! Otherwise QOK = .FALSE. ! In either case, a complete A is not generated. ! ! L index of loop momenta, L = 0,1,...,NLOOPS. ! L = 0 normally denontes the virtual photon momentum. ! P index of propagator, P = 0,1,...,NPROPS. ! P = 0 denotes the virtual photon momentum. ! V index of vertices, V = 1,...,NVERTS ! A(P,L) matrix relating propagator momenta to loop momenta. ! K(P) = Sum_L A(P,L) K(Q(L)). ! VRTX(P,1) = V means that the vertex connected to the tail of ! propagator P is V. ! VRTX(P,2) = V means that the vertex connected to the head of ! propagator P is V. ! Q(L) = P means that we consider the Lth loop momentum to ! be that carried by propagator P. ! CONNECTED(V,J) = P means that the Jth propagator connected to ! vertex V is P. ! FIXED(P) = True means that we have determined the momentum carried ! by propagator P. ! FINISHED(V) = True means that we have determined the momenta carried ! by all the propagators connected to vertex V. ! PROPSIGN(VRTX,P,V) is a function that returns +1 if the head of ! propagator P is at V, -1 if the tail is at V. ! COUNT is the number of propagators connected to the vertex ! under consideration such that FIXED(P) = True. If ! COUNT = 2, then we can fix another propagator momentum. ! 3 July 1994 ! 19 December 1995 ! integer :: nloops,nprops,nverts integer :: nloops1,nprops1,nverts1,cutmax1 integer :: nloops2,nprops2,nverts2,cutmax2 common /sizes/ nloops1,nprops1,nverts1,cutmax1, & nloops2,nprops2,nverts2,cutmax2 ! integer :: l,p,v,j,l1,l2 integer :: connected(2*size,3) logical :: fixed(0:3*size-1),finished(2*size) logical :: change integer :: propsign,sign integer :: sum(0:size) integer :: count integer :: ptofix ! IF (order.EQ.1) THEN nloops = nloops1 nprops = nprops1 nverts = nverts1 ELSE IF (order.EQ.2) THEN nloops = nloops2 nprops = nprops2 nverts = nverts2 ELSE write(nout,*)'Order must be 1 or 2.',order END IF ! IF((nq.LT.1).OR.(nq.GT.nloops)) THEN write(nout,*)'Nq out of range in finda.' END IF ! ! First check to see that the same propagator hasn't been ! assigned to two loop variables. ! DO l1 = 0,nq-1 DO l2 = l1+1,nq IF (q(l1).EQ.q(l2)) THEN qok = .false. RETURN END IF END DO END DO ! ! Initialization. ! qok = .false. ! DO v = 1,nverts j = 1 DO p = 0,nprops IF( (vrtx(p,1).EQ.v).OR.(vrtx(p,2).EQ.v) ) THEN connected(v,j) = p j = j+1 END IF END DO END DO ! DO p = 0,nprops DO l = 0,nloops a(p,l) = 0 END DO END DO DO l = 0,nq a(q(l),l) = 1 END DO ! DO p = 0,nprops fixed(p) = .false. END DO DO l = 0,nq fixed(q(l)) = .true. END DO ! DO v = 1,nverts finished(v) = .false. END DO ! change = .true. ! ! Start. ! DO WHILE (change) change = .false. ! DO v = 1,nverts IF (.NOT.finished(v)) THEN count = 0 DO j = 1,3 p = connected(v,j) IF ( fixed(p) ) THEN count = count + 1 END IF END DO ! ! There are 3 already fixed propagators conencted to this vertex, so ! we must check to see if the momenta coming into the vertex sum to ! zero. ! IF (count.EQ.3) THEN DO l = 0,nq sum(l) = 0 END DO DO j = 1,3 p = connected(v,j) sign = propsign(vrtx,p,v) DO l = 0,nq sum(l) = sum(l) + sign * a(p,l) END DO END DO DO l = 0,nq ! ! Dependent propagators given to FINDA. ! IF (sum(l).NE.0) THEN qok = .false. RETURN ! END IF END DO finished(v) = .true. change = .true. ! ! There are two already fixed propagators connected to this vertex, ! so we should determine the momentum carried by the remaining, ! unfixed, propagator. ! ELSE IF (count.EQ.2) THEN DO l = 0,nq sum(l) = 0 END DO DO j = 1,3 p = connected(v,j) IF ( fixed(p) ) THEN sign = propsign(vrtx,p,v) DO l = 0,nq sum(l) = sum(l) + sign * a(p,l) END DO ELSE ptofix = p END IF END DO sign = propsign(vrtx,ptofix,v) DO l = 0,nq a(ptofix,l) = - sign * sum(l) END DO fixed(ptofix) = .true. finished(v) = .true. change = .true. END IF ! ! Close loop DO V = 1,NVERTS ; IF (.NOT.FINISHED(V)) THEN. ! END IF END DO ! ! Close loop DO WHILE (CHANGE) ! END DO ! ! At this point, we have not found a contradiction with momentum ! conservation, so the Q's must have been OK: ! qok = .true. ! ! If we had been given a complete set of Q's, then we should have ! fixed each propagator at each vertex. Check just to make sure. ! IF (nq.EQ.nloops) THEN DO v = 1,nverts IF (.NOT.finished(v) ) THEN write(nout,*)'Snafu in finda.' write(nout,*)'v = ',v,' nq =',nq write(nout,*)'q =',q(0),q(1),q(2),q(3) stop END IF END DO END IF ! RETURN ! END subroutine finda ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! function propsign(vrtx,p,v) ! use beowulf_parameters implicit none ! In: integer :: vrtx(0:3*size-1,2) integer :: p,v ! Out: integer :: propsign ! ! IF ( vrtx(p,1).EQ.v ) THEN propsign = -1 RETURN ELSE IF ( vrtx(p,2).EQ.v ) THEN propsign = 1 RETURN ELSE write(nout,*)'Propsign called for p not connected to v.' stop END IF END function propsign ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine findtypes(graph,nmaps,qs,qsigns,maptypes) ! use beowulf_parameters use beowulf_structures implicit none ! In: type(graphstructure) :: graph ! Out: integer :: nmaps,qs(maxmaps,0:size),qsigns(maxmaps,0:size) character(len=6) :: maptypes(maxmaps) ! ! Given a graph specified by VRTX and PROP, this subroutine finds the ! characteristics of each map, labelled by an index MAPNUMBER. For a given map, ! it finds labels Q of 'special' propagators and the corresponding signs QSIGN ! and the MAPTYPE. The subroutine does finds the total number of maps, NMAPS, ! and fills the corresponding arrays QS, QSIGNS, and MAPTYPES, each of which ! carries a MAPNUMBER index. The definition is that the loop momenta are ! ! ELL(L,mu) = K(Q(L),mu) for L = 1,2,3. ! ! The other K(P,mu) are related to these by K(P,mu) = Sum A_{P,L} ELL(L,mu). ! Following the paper "Choosing integration points for QCD calculations by ! numerical integration" (Phys. Rev. D 64, 034018 (2001)), we also define ! vectors ELL1(mu), ELL2(mu) = P1(mu), and ELL3(mu) = P2(mu) which equal the ! respective ELL(L,mu) up to signs: ! ! ELL(1,mu) = QSIGN(1)*ELL1(mu) ! ELL(2,mu) = QSIGN(2)*ELL2(mu) ! ELL(3,mu) = QSIGN(3)*ELL3(mu) ! ! The paper does not, unfortunately, take note of these signs. The momentum ! flow directions shown in the paper are those of {ELL1(mu),ELL2(mu),ELL3(mu)}. ! The signs arise because the momentum flow directions for the propagator ! momenta K(P,mu) have independent definitions. ! ! The possibilities for maptypes are as follows: ! ! 1) T2TO2T used for k1 + k2 -> p1 + p2 with a virtual parton ! with momentum q exchanged, q = k1 - p1. Then P(Q(1)) = q, ! P(Q(2)) = p1, P(Q(3)) = p2. We will generate points with ! CHOOSET2TO2T(p1,p2,ell1,ok). ! ! 2) T2TO2S used for k1 + k2 -> p1 + p2 with a *no* virtual ! parton with momentum q = k1 - p1 exchanged. Then P(Q(1)) = k1, ! P(Q(2)) = p1, P(Q(3)) = p2. We will generate points with ! CHOOSE2TO2S(p1,p2,ell1,ok). ! ! 3) T2TO3 used for k1 + k2 -> p1 + p2 + p3 with k2 = - k1. ! Then P(Q(1)) = k1, P(Q(2)) = p1, P(Q(3)) = p2. We will generate ! points with CHOOSET2TO3(p1,p2,ell1,ok). ! ! 4) T2TO1A used for k1 + k2 -> p1 on shell for a virtual self-energy ! diagram. We will choose points with CHOOSEST2TO1A(p1,p2,ell1,ok). ! ! 4) T2TO1B used for k1 + k2 -> p1 on shell for a virtual three or four ! point function. We will choose points with CHOOSEST2TO1B(p1,p2,ell1,ok). ! This case is not discussed in the paper "Choosing integration points ! for QCD calculations by numerical integration." It is designed for ! the weak 1/\theta collinear singularities that occur in diagrams other ! than self energy diagrams in Coulomb gauge before real-virtual cancellations. ! ! We also have the possibility of Born graphs, for which the maptype ! is BORN and Q(1) and Q(2) are chosen as two of the cut propagators. ! ! 20 December 2000 ! 20 March 2001 ! 1 February 2002 ! 13 June 2002 ! 7 December 2002 ! integer :: cutmax integer :: nloops1,nprops1,nverts1,cutmax1 integer :: nloops2,nprops2,nverts2,cutmax2 common /sizes/ nloops1,nprops1,nverts1,cutmax1, & nloops2,nprops2,nverts2,cutmax2 ! integer :: mapnumber ! Newcut variables type (cutstructure) :: cut logical :: cutfound ! integer :: l,p,kj,k1,k2,kdirect,ptest,pleaving,pp1,pp2 integer :: l1,signl1,ltesta,ltestb integer :: i,j,jfound1,jfound2,jfound integer :: v1,v2,v3,vother,vv1,vv2 integer :: sign0,sign1,sign2 integer :: timesfound1,timesfound2,timesfound integer :: pl,signl,vl,kl,signkl logical :: notinloop ! !---------------------------------- ! IF (graph%order.EQ.1) THEN cutmax = cutmax1 ELSE IF (graph%order.EQ.2) THEN cutmax = cutmax2 ELSE write(nout,*)'Order must be 1 or 2.',graph%order END IF ! mapnumber = 0 DO call getnewcut(graph%graphnumber,cut,cutfound) ! IF (cutfound) THEN ! IF (graph%order.EQ.1) THEN ! ! First, we have the code for what to do for Born graphs. There ! is no Q(3) or QSIGN(3) in this case. ! mapnumber = mapnumber + 1 qs(mapnumber,0) = 0 qsigns(mapnumber,0) = +1 qs(mapnumber,1) = cut%cutindex(1) qsigns(mapnumber,1) = cut%cutsign(1) qs(mapnumber,2) = cut%cutindex(2) qsigns(mapnumber,2) = cut%cutsign(2) qs(mapnumber,3) = 137 qsigns(mapnumber,3) = 137 maptypes(mapnumber) = 'born ' ! ! Alternative for IF (ORDER.EQ.1) THEN ! ELSE IF (graph%order.EQ.2) THEN ! ! We want to do something only if there is a virtual loop: ! IF (cut%ncut.EQ.(cutmax-1)) THEN !--- ! Case of 4 propagators in the loop !--- IF (cut%ninloop.EQ.4) THEN ! ! For a 4 propagator loop there are three 2 to 1 maps. ! DO l = 1,3 mapnumber = mapnumber + 1 ! ! We let pl be the lth propagator around the loop and vl be the lth ! vertex around the loop. Also, signl = +1 if propagator pl is directed ! in the positive loop direction. ! pl = cut%loopindex(l) signl = cut%loopsign(l) IF (signl.EQ.1) THEN vl = graph%vrtx(pl,2) ELSE vl = graph%vrtx(pl,1) END IF ! ! We find the cut propagator kl connected to vl along ! with signkl = +1 if the cut propagator kl is leaving vertex vl ! and signkl = -1 if the cut propagator kl is entering vertex vl. Just ! as a check, we define TIMESFOUND to see if we find kl exactly ! once. ! timesfound = 0 DO j = 1,3 kj = cut%cutindex(j) IF (graph%vrtx(kj,1).EQ.vl) THEN kl = kj signkl = +1 jfound = j timesfound = timesfound + 1 ELSE IF (graph%vrtx(kj,2).EQ.vl) THEN kl = kj signkl = -1 jfound = j timesfound = timesfound + 1 END IF END DO IF (timesfound.NE.1) THEN write(nout,*) 'Failure in findtypes.' stop END IF ! ! Now we record this information. The propagator Q(3) is one of the ! propagators in the final state other than that connected to vl. The ! corresponding sign is +1 if this propagator crosses the final ! state cut in the same direction as the propagator connected to vl. ! qs(mapnumber,0) = 0 qsigns(mapnumber,0) = +1 qs(mapnumber,1) = pl qsigns(mapnumber,1) = signl qs(mapnumber,2) = kl qsigns(mapnumber,2) = signkl IF (cut%cutindex(1).NE.kl) THEN qs(mapnumber,3) = cut%cutindex(1) qsigns(mapnumber,3) = cut%cutsign(1)*cut%cutsign(jfound)*signkl ELSE qs(mapnumber,3) = cut%cutindex(2) qsigns(mapnumber,3) = cut%cutsign(2)*cut%cutsign(jfound)*signkl END IF maptypes(mapnumber) = 't2to1b' ! ! Close DO l = 1,3 END DO ! ! For a 4 propagator loop there are two ellipse maps (T2T02T) and ! one circle map (T2TO3). We do the two ellipse maps first. ! DO l = 2,3 mapnumber = mapnumber + 1 p = cut%loopindex(l) v1 = graph%vrtx(p,1) v2 = graph%vrtx(p,2) ! ! We find the cut propagators K1 and K2 connected to V1 and V2 along ! with the sign = +1 if the cut propagator Kj is leaving vertex Vj ! and sign = -1 if the cut propagator Kj is entering vertex Vj. Just ! as a check, we define FOUNDJ to see if we find K1 and K2 exactly ! once. ! timesfound1 = 0 timesfound2 = 0 DO j = 1,3 kj = cut%cutindex(j) IF (graph%vrtx(kj,1).EQ.v1) THEN k1 = kj sign1 = +1 timesfound1 = timesfound1+1 ELSE IF (graph%vrtx(kj,2).EQ.v1) THEN k1 = kj sign1 = -1 timesfound1 = timesfound1+1 ELSE IF (graph%vrtx(kj,1).EQ.v2) THEN k2 = kj sign2 = +1 timesfound2 = timesfound2+1 ELSE IF (graph%vrtx(kj,2).EQ.v2) THEN k2 = kj sign2 = -1 timesfound2 = timesfound2+1 END IF END DO IF ((timesfound1.NE.1).OR.(timesfound2.NE.1)) THEN write(nout,*) 'Failure in findtypes.' stop END IF ! ! Now we record this information. ! qs(mapnumber,0) = 0 qsigns(mapnumber,0) = +1 qs(mapnumber,1) = p qsigns(mapnumber,1) = +1 qs(mapnumber,2) = k1 qsigns(mapnumber,2) = sign1 qs(mapnumber,3) = k2 qsigns(mapnumber,3) = sign2 maptypes(mapnumber) = 't2to2t' ! ! End DO L = 2,3 for the choice of two ellipse maps. ! END DO ! ! Now we do the circle map. ! Our definition for the circle map T2TO3E is that Q(1) is ! LOOPINDEX(1) the first propagator in the loop starting from the ! current vertex. Then Q(2) is the label of the propagator that ! enters the final state and connects to the vertex at the head ! of propagator Q(1). Then Q(3) is the label of the propagator ! that enters the final state and connects to the propagator with ! label LOOPINDEX(4), the last propagator in the loop. The sign ! QSIGN(1) = +1 since this propagator always points *from* the ! current vertex. For QSIGN(2) and QSIGN(3) a plus sign indicates ! that the propagator points toward the final state, a minus sign ! indicates the opposite. ! IF(cut%loopsign(1).NE.1) THEN write(nout,*)'loopsign(1) not 1 in findtypes.' stop ELSE IF(cut%loopsign(4).NE.-1) THEN write(nout,*)'loopsign(4) not -1 in findtypes.' stop END IF ! v1 = graph%vrtx(cut%loopindex(1),2) v2 = graph%vrtx(cut%loopindex(4),2) timesfound1 = 0 timesfound2 = 0 DO j = 1,3 kj = cut%cutindex(j) IF (graph%vrtx(kj,1).EQ.v1) THEN k1 = kj sign1 = +1 timesfound1 = timesfound1+1 ELSE IF (graph%vrtx(kj,2).EQ.v1) THEN k1 = kj sign1 = -1 timesfound1 = timesfound1+1 ELSE IF (graph%vrtx(kj,1).EQ.v2) THEN k2 = kj sign2 = +1 timesfound2 = timesfound2+1 ELSE IF (graph%vrtx(kj,2).EQ.v2) THEN k2 = kj sign2 = -1 timesfound2 = timesfound2+1 END IF END DO IF ((timesfound1.NE.1).OR.(timesfound2.NE.1)) THEN write(nout,*) 'Oops, failure in findtypes.', & timesfound1,timesfound2 stop END IF ! mapnumber = mapnumber + 1 qs(mapnumber,0) = 0 qsigns(mapnumber,0) = +1 qs(mapnumber,1) = cut%loopindex(1) qsigns(mapnumber,1) = +1 qs(mapnumber,2) = k1 qsigns(mapnumber,2) = sign1 qs(mapnumber,3) = k2 qsigns(mapnumber,3) = sign2 maptypes(mapnumber) = 't2to3e' ! !--- ! Case of 3 propagators in the loop !--- ELSE IF (cut%ninloop.EQ.3) THEN ! ! With three propagators in the loop there are one or two 2 to 1 maps. ! DO l = 1,2 ! ! We let pl be the lth propagator around the loop and vl be the lth ! vertex around the loop. Also, signl = +1 if propagator pl is directed ! in the positive loop direction. ! pl = cut%loopindex(l) signl = cut%loopsign(l) IF (signl.EQ.1) THEN vl = graph%vrtx(pl,2) ELSE vl = graph%vrtx(pl,1) END IF ! ! We find the cut propagator kl connected to vl along ! with signkl = +1 if the cut propagator kl is leaving vertex vl ! and signkl = -1 if the cut propagator kl is entering vertex vl. ! We define TIMESFOUND to see if we find kl exactly once. It can ! happen that we will not find final state propagator that connects ! to vl. Then we do not record a 2 to 1 map. ! timesfound = 0 DO j = 1,3 kj = cut%cutindex(j) IF (graph%vrtx(kj,1).EQ.vl) THEN kl = kj signkl = +1 jfound = j timesfound = timesfound + 1 ELSE IF (graph%vrtx(kj,2).EQ.vl) THEN kl = kj signkl = -1 jfound = j timesfound = timesfound + 1 END IF END DO IF (timesfound.GT.1) THEN write(nout,*) 'Failure in findtypes.' stop ELSE IF (timesfound.EQ.1) THEN ! ! We have a 2 to 1 map and we need to record this information. ! The propagator Q(3) is one of the propagators in the final state other ! than that connected to vl. The corresponding sign is +1 if this propagator ! crosses the final state cut in the same direction as the propagator ! connected to vl. ! mapnumber = mapnumber + 1 qs(mapnumber,0) = 0 qsigns(mapnumber,0) = +1 qs(mapnumber,1) = pl qsigns(mapnumber,1) = signl qs(mapnumber,2) = kl qsigns(mapnumber,2) = signkl IF (cut%cutindex(1).NE.kl) THEN qs(mapnumber,3) = cut%cutindex(1) qsigns(mapnumber,3) = cut%cutsign(1)*cut%cutsign(jfound)*signkl ELSE qs(mapnumber,3) = cut%cutindex(2) qsigns(mapnumber,3) = cut%cutsign(2)*cut%cutsign(jfound)*signkl END IF maptypes(mapnumber) = 't2to1b' ! ! Close ELSE IF (timesfound.EQ.1) THEN END IF ! Close DO l = 1,3 END DO ! ! Now continuing with the virtual three point function, we will need ! either a 2to2(t) map or a 2to2(s) map and a 2to3 map. ! We are not sure which of two possibilities we have, but we proceed ! as if we had the case of a virtual loop that connects to two ! propagators that go into the final state. ! mapnumber = mapnumber + 1 p = cut%loopindex(2) v1 = graph%vrtx(p,1) v2 = graph%vrtx(p,2) ! ! We find the cut propagators K1 and K2 connected to V1 and V2 along ! with the sign = +1 if the cut propagator Kj is leaving vertex Vj ! and sign = -1 if the cut propagator Kj is entering vertex Vj. We ! check using FOUNDJ to see if we find K1 and K2 exactly once. ! timesfound1 = 0 timesfound2 = 0 DO j = 1,3 kj = cut%cutindex(j) IF (graph%vrtx(kj,1).EQ.v1) THEN k1 = kj sign1 = +1 timesfound1 = timesfound1+1 ELSE IF (graph%vrtx(kj,2).EQ.v1) THEN k1 = kj sign1 = -1 timesfound1 = timesfound1+1 ELSE IF (graph%vrtx(kj,1).EQ.v2) THEN k2 = kj sign2 = +1 timesfound2 = timesfound2+1 ELSE IF (graph%vrtx(kj,2).EQ.v2) THEN k2 = kj sign2 = -1 timesfound2 = timesfound2+1 END IF END DO ! ! Now we figure out what to do based on what we found. ! IF ((timesfound1.GT.1).OR.(timesfound2.GT.1)) THEN write(nout,*) 'Failure in findtypes.' stop ELSE IF ((timesfound1.LT.1).AND.(timesfound2.LT.1)) THEN write(nout,*) 'Failure in findtypes.' stop ! ELSE IF ((timesfound1.EQ.1).AND.(timesfound2.EQ.1)) THEN ! ! This is the case we were looking for. Now we record the information. ! qs(mapnumber,0) = 0 qsigns(mapnumber,0) = +1 qs(mapnumber,1) = p qsigns(mapnumber,1) = +1 qs(mapnumber,2) = k1 qsigns(mapnumber,2) = sign1 qs(mapnumber,3) = k2 qsigns(mapnumber,3) = sign2 maptypes(mapnumber) = 't2to2t' ! ELSE ! ! Either Found1 = 1 and Found2 = 0 or Found2 = 1 and Found1 = 0. ! In these cases our loop does *not* connect to two propagators ! that go to the final state. The label of the propagator ! that enters the final state will be called Kdirect and the ! vertex that does not connect to this propagator will be called ! Vother. We take sign0 = +1 if our loop propagator points from ! Kdirect to the s-channel propagator that splits into two ! propagators that go to the final state. Otherwise sign0 = -1. ! IF ((timesfound1.EQ.1).AND.(timesfound2.EQ.0)) THEN kdirect = k1 vother = v2 sign0 = +1 ELSE IF ((timesfound1.EQ.0).AND.(timesfound2.EQ.1)) THEN kdirect = k2 vother = v1 sign0 = -1 END IF ! ! Now we deal with this case. ! IF (cut%cutindex(1).EQ.kdirect) THEN k1 = cut%cutindex(2) k2 = cut%cutindex(3) ELSE IF (cut%cutindex(2).EQ.kdirect) THEN k1 = cut%cutindex(3) k2 = cut%cutindex(1) ELSE IF (cut%cutindex(3).EQ.kdirect) THEN k1 = cut%cutindex(1) k2 = cut%cutindex(2) ELSE write(nout,*)'We are in real trouble here.' stop END IF ! ! We have K1 and K2, but we need the corresponding signs. ! Find the index Pleaving of the propagator leaving the loop toward ! the final state. ! timesfound = 0 DO j = 1,3 ptest = graph%prop(vother,j) notinloop = .true. DO i = 1,3 IF (ptest.EQ.cut%loopindex(i)) THEN notinloop = .false. END IF END DO IF (notinloop) THEN pleaving = ptest timesfound = timesfound + 1 END IF END DO IF (timesfound.NE.1) THEN write(nout,*)'pleaving not found or found twice.' stop END IF ! ! Let V3 be the vertex not in the loop at the end of propagator ! Pleaving. Two propagators in the final state must connect to this ! vertex. ! v3 = graph%vrtx(pleaving,1) IF (v3.EQ.vother) THEN v3 = graph%vrtx(pleaving,2) END IF ! ! We use V3 to get the proper signs. ! IF (graph%vrtx(k1,1).EQ.v3) THEN sign1 = +1 ELSE IF (graph%vrtx(k1,2).EQ.v3) THEN sign1 = -1 ELSE write(nout,*)'Yikes, this is bad.' stop END IF IF (graph%vrtx(k2,1).EQ.v3) THEN sign2 = +1 ELSE IF (graph%vrtx(k2,2).EQ.v3) THEN sign2 = -1 ELSE write(nout,*)'Yikes, this is also bad.' stop END IF ! ! Now we record the information. ! Recall that P = LOOPINDEX(2) and that SIGN0 = +1 if propagator P ! points toward propagator Pleaving. ! qs(mapnumber,0) = 0 qsigns(mapnumber,0) = +1 qs(mapnumber,1) = p qsigns(mapnumber,1) = sign0 qs(mapnumber,2) = k1 qsigns(mapnumber,2) = sign1 qs(mapnumber,3) = k2 qsigns(mapnumber,3) = sign2 maptypes(mapnumber) = 't2to2s' ! ! But we are not done, because in this case we need a circle map too. ! Our definition for the circle map T2TO3D is that Q(1) is ! LOOPINDEX(1) or LOOPINDEX(3), one of the two propagators that ! connects to a propagator that connects to the current vertex. ! We take the one that connects to vertex Vother that connects to ! a propagator Pleaving that connects vertex V3 that, finally, ! connects to to two propagators that enter the final state. Then ! Q(3) and Q(3) are these two propagators that enter the final ! state from vertex V3. For QSIGN(2) and QSIGN(3) a plus sign ! indicates that the propagator points toward the final state, a ! minus sign indicates the opposite. The sign QSIGN(1) is + 1 if ! this propagator points toward the final state, -1 in the ! opposite circumstance. ! IF (cut%loopsign(1).NE.1) THEN write(nout,*)'loopsign not 1 in findtypes.' stop END IF ! ! The loop momentum with label L1 is the one that ! attaches to VOTHER (the vertex that connects to a propagator ! that splits before going to the final state.) We take ! SIGNL1 = +1 if this propagator points towards VOTHER. ! ltesta = cut%loopindex(1) ltestb = cut%loopindex(3) IF (graph%vrtx(ltesta,2).EQ.vother) THEN l1 = ltesta signl1 = +1 ELSE IF (graph%vrtx(ltesta,1).EQ.vother) THEN l1 = ltesta signl1 = -1 ELSE IF (graph%vrtx(ltestb,2).EQ.vother) THEN l1 = ltestb signl1 = +1 ELSE IF (graph%vrtx(ltestb,1).EQ.vother) THEN l1 = ltestb signl1 = -1 ELSE write(nout,*)'Cannot seem to find l1.' stop END IF ! mapnumber = mapnumber + 1 qs(mapnumber,0) = 0 qsigns(mapnumber,0) = +1 qs(mapnumber,1) = l1 qsigns(mapnumber,1) = signl1 qs(mapnumber,2) = k1 qsigns(mapnumber,2) = sign1 qs(mapnumber,3) = k2 qsigns(mapnumber,3) = sign2 maptypes(mapnumber) = 't2to3d' ! ! Close the IF structure for the case of three particles in the loop, ! IF (FOUND1.GT.1).OR.(FOUND2.GT.1) THEN ... ! END IF !--- ! Case of 2 propagators in the loop !--- ELSE IF (cut%ninloop.EQ.2) THEN ! ! We are not sure which of two possibilities we have, but we proceed ! as if we had the case of a virtual loop that connects to two ! propagators that go into the final state. ! mapnumber = mapnumber + 1 p = cut%loopindex(1) v1 = graph%vrtx(p,1) v2 = graph%vrtx(p,2) ! ! We find the cut propagators K1 or K2 connected to V1 or V2 along ! with the sign = +1 if the cut propagator Kj is leaving vertex Vj ! and sign = -1 if the cut propagator Kj is entering vertex Vj. We ! check using FOUNDJ to see if we find K1 or K2 exactly once. ! timesfound1 = 0 timesfound2 = 0 DO j = 1,3 kj = cut%cutindex(j) IF (graph%vrtx(kj,1).EQ.v1) THEN k1 = kj sign1 = +1 jfound1 = j timesfound1 = timesfound1+1 ELSE IF (graph%vrtx(kj,2).EQ.v1) THEN k1 = kj sign1 = -1 jfound1 = j timesfound1 = timesfound1+1 ELSE IF (graph%vrtx(kj,1).EQ.v2) THEN k2 = kj sign2 = +1 jfound2 = j timesfound2 = timesfound2+1 ELSE IF (graph%vrtx(kj,2).EQ.v2) THEN k2 = kj sign2 = -1 jfound2 = j timesfound2 = timesfound2+1 END IF END DO ! ! Now we figure out what to do based on what we found. ! IF ((timesfound1.GT.1).OR.(timesfound2.GT.1)) THEN write(nout,*) 'Failure in findtypes.' stop ! ELSE IF ((timesfound1.EQ.1).AND.(timesfound2.EQ.0)) THEN ! ! This is one of the cases we were looking for. Now we record the ! information. The propagator Q(3) is one of the propagators ! in the final state other than that connected to our loop. The ! corresponding sign is +1 if this propagator crosses the final ! state cut in the same direction as the propagator connected to ! our loop. ! qs(mapnumber,0) = 0 qsigns(mapnumber,0) = +1 qs(mapnumber,1) = p qsigns(mapnumber,1) = -1 qs(mapnumber,2) = k1 qsigns(mapnumber,2) = sign1 IF (cut%cutindex(1).NE.k1) THEN qs(mapnumber,3) = cut%cutindex(1) qsigns(mapnumber,3) = cut%cutsign(1)*cut%cutsign(jfound1)*sign1 ELSE qs(mapnumber,3) = cut%cutindex(2) qsigns(mapnumber,3) = cut%cutsign(2)*cut%cutsign(jfound1)*sign1 END IF maptypes(mapnumber) = 't2to1a' ! ELSE IF ((timesfound1.EQ.0).AND.(timesfound2.EQ.1)) THEN ! ! This is one of the cases we were looking for. Now we record the ! information. The propagator Q(3) is one of the propagators ! in the final state other than that connected to our loop. The ! corresponding sign is +1 if this propagator crosses the final ! state cut in the same direction as the propagator connected to ! our loop. ! qs(mapnumber,0) = 0 qsigns(mapnumber,0) = +1 qs(mapnumber,1) = p qsigns(mapnumber,1) = +1 qs(mapnumber,2) = k2 qsigns(mapnumber,2) = sign2 IF (cut%cutindex(1).NE.k2) THEN qs(mapnumber,3) = cut%cutindex(1) qsigns(mapnumber,3) = cut%cutsign(1)*cut%cutsign(jfound2)*sign2 ELSE qs(mapnumber,3) = cut%cutindex(2) qsigns(mapnumber,3) = cut%cutsign(2)*cut%cutsign(jfound2)*sign2 END IF maptypes(mapnumber) = 't2to1a' ! ELSE IF ((timesfound1.EQ.0).AND.(timesfound2.EQ.0)) THEN ! ! Here TimesFound1 = 0 and TimesFound2 = 0, so our loop does *not* ! connect to a propagator that goes to the final state. ! Find the indices PP1 and PP2 of the propagators connected to ! our loop. ! timesfound = 0 DO j = 1,3 ptest = graph%prop(v1,j) notinloop = .true. DO i = 1,2 IF (ptest.EQ.cut%loopindex(i)) THEN notinloop = .false. END IF END DO IF (notinloop) THEN pp1 = ptest timesfound = timesfound + 1 END IF END DO IF (timesfound.NE.1) THEN write(nout,*)'pp1 not found or found twice.' stop END IF timesfound = 0 ! DO j = 1,3 ptest = graph%prop(v2,j) notinloop = .true. DO i = 1,2 IF (ptest.EQ.cut%loopindex(i)) THEN notinloop = .false. END IF END DO IF (notinloop) THEN pp2 = ptest timesfound = timesfound + 1 END IF END DO IF (timesfound.NE.1) THEN write(nout,*)'pp2 not found or found twice.' stop END IF ! ! Let VV1 and VV2 be the vertices not in the loop at the end of ! propagators PP1 and PP2 respectively. Two propagators in the final ! state must connect to one of these vertices. ! vv1 = graph%vrtx(pp1,1) IF (vv1.EQ.v1) THEN vv1 = graph%vrtx(pp1,2) END IF vv2 = graph%vrtx(pp2,1) IF (vv2.EQ.v2) THEN vv2 = graph%vrtx(pp2,2) END IF ! ! We have VV1 and VV2. A slight hitch is that one of them might be ! the vertex 1 or 2 that connect to the photon. In this case, ! in the next step we do *not* want to find the final state ! propagator that connects to this vertex. A cure is to set the ! vertex number to something impossible. ! IF ((vv1.EQ.1).OR.(vv1.EQ.2)) THEN vv1 = -17 END IF IF ((vv2.EQ.1).OR.(vv2.EQ.2)) THEN vv2 = -17 END IF ! ! Now we find two final state propagators connected to VV1 or ! else two final state propagators connected to VV2. ! timesfound = 0 DO j = 1,3 kj = cut%cutindex(j) IF (graph%vrtx(kj,1).EQ.vv1) THEN IF(timesfound.EQ.0) THEN k1 = kj sign1 = +1 ELSE k2 = kj sign2 = +1 END IF sign0 = -1 timesfound = timesfound+1 ELSE IF (graph%vrtx(kj,1).EQ.vv2) THEN IF(timesfound.EQ.0) THEN k1 = kj sign1 = +1 ELSE k2 = kj sign2 = +1 END IF sign0 = +1 timesfound = timesfound+1 ELSE IF (graph%vrtx(kj,2).EQ.vv1) THEN IF(timesfound.EQ.0) THEN k1 = kj sign1 = -1 ELSE k2 = kj sign2 = -1 END IF sign0 = -1 timesfound = timesfound+1 ELSE IF (graph%vrtx(kj,2).EQ.vv2) THEN IF(timesfound.EQ.0) THEN k1 = kj sign1 = -1 ELSE k2 = kj sign2 = -1 END IF sign0 = +1 timesfound = timesfound+1 END IF END DO IF (timesfound.NE.2) THEN write(nout,*)'Where are those tricky propagators?',timesfound stop END IF ! ! Now we record the information. ! Recall that P = LOOPINDEX(1) and that SIGN0 = +1 if propagator P ! points toward propagators in the final state. ! qs(mapnumber,0) = 0 qsigns(mapnumber,0) = +1 qs(mapnumber,1) = p qsigns(mapnumber,1) = sign0 qs(mapnumber,2) = k1 qsigns(mapnumber,2) = sign1 qs(mapnumber,3) = k2 qsigns(mapnumber,3) = sign2 maptypes(mapnumber) = 't2to2s' ! ! Close IF ((TIMESFOUND1.GT.1).OR.(TIMESFOUND2.GT.1)) THEN ... ! END IF ! ! Case of less than 2 propagators in the loop ! ELSE write(nout,*) 'Looped the loop in finddqs.' stop ! ! End IF (NINLOOP.EQ. n ) series ! END IF ! ! End IF (there is a virtual loop) THEN ... ! END IF ! ! Close IF (ORDER.EQ.1) THEN ... ELSE IF (ORDER.EQ.2) THEN ... ! ELSE write(nout,*)'Order in findtypes needed to be 1 or 2.' stop END IF ! ! End IF (CUTFOUND) THEN ... If the cut was not found, then we are done. ! ELSE EXIT END IF ! ! End main loop: DO call getnewcut() ! END DO ! nmaps = mapnumber RETURN END subroutine findtypes ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890!2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine newpoint(a,qsign,maptype,order,k,absk,badpoint) ! use beowulf_parameters implicit none ! In: integer :: a(0:3*size-1,0:size),qsign(0:size) character(len=6) :: maptype integer :: order ! Out: real(kind=dbl) :: k(0:3*size-1,0:3),absk(0:3*size-1) logical :: badpoint ! ! Chooses a new Monte Carlo point in the space of loop 3-momenta. ! 4 March 1993 ! 12 July 1993 ! 17 July 1994 ! 2 May 1996 ! 5 February 1997 ! 4 February 1999 ! 10 March 1999 ! 9 April 1999 ! 20 August 1999 ! 21 December 2000 ! 20 March 2001 ! 8 February 2002 ! 7 December 2002 ! integer :: nprops integer :: nloops1,nprops1,nverts1,cutmax1 integer :: nloops2,nprops2,nverts2,cutmax2 common /sizes/ nloops1,nprops1,nverts1,cutmax1, & nloops2,nprops2,nverts2,cutmax2 ! real(kind=dbl) :: p1(3),p2(3),p3(3),ell1(3) integer :: p,mu real(kind=dbl) :: temp,ksq logical :: ok ! !------------ ! IF (order.EQ.1) THEN nprops = nprops1 ELSE IF (order.EQ.2) THEN nprops = nprops2 ELSE write(nout,*)'Order must be 1 or 2.',order END IF ! badpoint = .false. ! IF (order.EQ.1) THEN ! ! We deal with the case of a Born graph first. ! call choose3(p1,p2,p3,ok) IF(.NOT.ok) THEN DO p = 1,nprops DO mu = 0,3 k(p,mu) = 0.0d0 END DO END DO badpoint = .true. RETURN END IF DO p = 1,nprops ksq = 0.0d0 DO mu = 1,3 temp = a(p,1)*qsign(1)*p1(mu) & + a(p,2)*qsign(2)*p2(mu) k(p,mu) = temp ksq = ksq + temp**2 END DO absk(p) = sqrt(ksq) k(p,0) = 0.0d0 END DO DO mu = 0,3 k(0,mu) = 0.0d0 END DO absk(0) = 0.0d0 ! ! Alternative for IF (ORDER.EQ.1) THEN ! ELSE IF (order.EQ.2) THEN ! ! Here is what we do for order alpha_s^2 graphs. ! ! Our notation: The definition is that the loop momenta are ! ! ELL(L,mu) = K(Q(L),mu) for L = 1,2,3. ! ! The other K(P,mu) are related to these by K(P,mu) = Sum A_{P,L} ELL(L,mu). ! Following the paper "Choosing integration points for QCD calculations by ! numerical integration" (Phys. Rev. D 64, 034018 (2001)), we also define ! vectors ELL1(mu), ELL2(mu) = P1(mu), and ELL3(mu) = P2(mu) which equal the ! respective ELL(L,mu) up to signs: ! ! ELL(1,mu) = QSIGN(1)*ELL1(mu) ! ELL(2,mu) = QSIGN(2)*ELL2(mu) ! ELL(3,mu) = QSIGN(3)*ELL3(mu) ! ! The paper does not, unfortunately, take note of these signs. The momentum ! flow directions shown in the paper are those of {ELL1(mu),ELL2(mu),ELL3(mu)}. ! The signs arise because the momentum flow directions for the propagator ! momenta K(P,mu) have independent definitions. We use {ELL1,P1,P2} in this ! subroutine directly to generate the K(P,mu). ! ! First, we need to generate a three parton final state. ! Abort if we get a not OK signal. ! call choose3(p1,p2,p3,ok) IF(.NOT.ok) THEN DO p = 1,nprops DO mu = 0,3 k(p,mu) = 0.0d0 END DO END DO badpoint = .true. RETURN END IF ! ! Then we generate the loop momentum, ell1. ! Abort if we get a not OK signal. ! IF (maptype.EQ.'t2to3d') THEN call choose2to3d(p1,p2,ell1,ok) ELSE IF (maptype.EQ.'t2to3e') THEN call choose2to3e(p1,p2,ell1,ok) ELSE IF (maptype.EQ.'t2to2t') THEN call choose2to2t(p1,p2,ell1,ok) ELSE IF (maptype.EQ.'t2to2s') THEN call choose2to2s(p1,p2,ell1,ok) ELSE IF (maptype.EQ.'t2to1a') THEN call choose2to1a(p1,p2,ell1,ok) ELSE IF (maptype.EQ.'t2to1b') THEN call choose2to1b(p1,p2,ell1,ok) ELSE write(nout,*)'Bad maptype in newpoint.' stop END IF IF (.NOT.ok) THEN DO p = 1,nprops DO mu = 0,3 k(p,mu) = 0.0d0 END DO END DO badpoint = .true. RETURN END IF ! ! Now we have ELL1(mu), P1(mu), and P2(mu) and we need to translate to ! the propagator momenta K(P,MU). ! DO p = 1,nprops ksq = 0.0d0 DO mu = 1,3 temp = a(p,1)*qsign(1)*ell1(mu) & + a(p,2)*qsign(2)*p1(mu) & + a(p,3)*qsign(3)*p2(mu) k(p,mu) = temp ksq = ksq + temp**2 END DO absk(p) = sqrt(ksq) k(p,0) = 0.0d0 END DO DO mu = 0,3 k(0,mu) = 0.0d0 END DO absk(0) = 0.0d0 ! ! Close IF (ORDER.EQ.1) THEN ... ELSE IF (ORDER.EQ.2) THEN ... ! ELSE write(nout,*)'Order should have been 1 or 2 in newpoint.' stop END IF ! RETURN END subroutine newpoint ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine checkpoint(k,absk,prop,order,badness) ! use beowulf_parameters implicit none ! In: real(kind=dbl) :: k(0:3*size-1,0:3),absk(0:3*size-1) integer :: prop(2*size,3) integer :: order ! Out: real(kind=dbl) :: badness ! ! Calculates the BADNESS of a point chosen by NEWPOINT. If there ! are very collinear particles meeting at a vertex or of there is a ! very soft particle, then the badness is big. Specifically, for ! each vertex V we order the momenta entering the vertex Kmin, Kmid ! Kmax in order of |K|. Then ! ! Kmin (Kmin + Kmid - Kmax )/Kmax^2 ! ! is the 1/badness^2 for that vertex. The badness for the point is the ! largest of the badness values of all the vertices. ! ! Revised 13 may 1998 ! integer :: nverts integer :: nloops1,nprops1,nverts1,cutmax1 integer :: nloops2,nprops2,nverts2,cutmax2 common /sizes/ nloops1,nprops1,nverts1,cutmax1, & nloops2,nprops2,nverts2,cutmax2 ! real(kind=dbl) :: smallnessv,smallness integer :: v real(kind=dbl) :: kmin,kmid,kmax,k1,k2,k3 ! IF (order.EQ.1) THEN nverts = nverts1 ELSE IF (order.EQ.2) THEN nverts = nverts2 ELSE write(nout,*)'Order must be 1 or 2.',order END IF ! smallness = 1.0d0 DO v = 3,nverts k1 = absk(prop(v,1)) k2 = absk(prop(v,2)) k3 = absk(prop(v,3)) IF (k1.LT.k2) THEN kmin = k1 kmax = k2 ELSE kmin = k2 kmax = k1 END IF IF (k3.LT.kmin) THEN kmid = kmin kmin = k3 ELSE IF (k3.GT.kmax) THEN kmid = kmax kmax = k3 ELSE kmid = k3 END IF smallnessv = kmin * (kmin + kmid - kmax) /kmax**2 IF( smallnessv .LT. smallness ) THEN smallness = smallnessv END IF END DO IF (smallness.LT.1.0d-30) THEN badness = 1.0d15 ELSE badness = sqrt(1.0d0/smallness) END IF RETURN END subroutine checkpoint ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine axes(ea,eb,ec) ! use beowulf_parameters implicit none ! In: real(kind=dbl) :: ea(3) ! Out: real(kind=dbl) :: eb(3),ec(3) ! ! Given a unit vector E_a(mu), generates unit vectors E_b(mu) and ! E_c(mu) such that E_a*E_b = E_b*E_c = E_c*E_a = 0. ! ! The vector E_b will lie in the plane formed by the z-axis and ! E_a unless E_a itself is nearly aligned along the z-axis, in which ! case E_b will lie in the plane formed by the x-axis and E_a. ! ! 18 April 1996 ! ! real(kind=dbl) :: costhetasq,sinthetainv ! ! For check ! integer :: mu real(kind=dbl) :: dotaa,dotbb,dotcc,dotab,dotac,dotbc ! costhetasq = ea(3)**2 IF(costhetasq.LT.0.9d0) THEN sinthetainv = 1.0d0/sqrt(1.0d0 - costhetasq) ec(1) = - ea(2)*sinthetainv ec(2) = ea(1)*sinthetainv ec(3) = 0.0d0 ELSE costhetasq = ea(1)**2 sinthetainv = 1.0d0/sqrt(1.0d0 - costhetasq) ec(1) = 0.0d0 ec(2) = - ea(3)*sinthetainv ec(3) = ea(2)*sinthetainv END IF eb(1) = ec(2)*ea(3) - ec(3)*ea(2) eb(2) = ec(3)*ea(1) - ec(1)*ea(3) eb(3) = ec(1)*ea(2) - ec(2)*ea(1) ! ! Check: ! dotaa = 0.0d0 dotbb = 0.0d0 dotcc = 0.0d0 dotab = 0.0d0 dotac = 0.0d0 dotbc = 0.0d0 DO mu = 1,3 dotaa = dotaa + ea(mu)*ea(mu) dotbb = dotbb + eb(mu)*eb(mu) dotcc = dotcc + ec(mu)*ec(mu) dotab = dotab + ea(mu)*eb(mu) dotac = dotac + ea(mu)*ec(mu) dotbc = dotbc + eb(mu)*ec(mu) END DO IF (abs(dotaa - 1.0d0).GT.1.0d20) THEN write(nout,*)'dotaa messed up in axes.' stop ELSE IF (abs(dotbb - 1.0d0).GT.1.0d20) THEN write(nout,*)'dotbb messed up in axes.' stop ELSE IF (abs(dotcc - 1.0d0).GT.1.0d20) THEN write(nout,*)'dotcc messed up in axes.' stop ELSE IF (abs(dotab).GT.1.0d20) THEN write(nout,*)'dotab messed up in axes.' stop ELSE IF (abs(dotac).GT.1.0d20) THEN write(nout,*)'dotac messed up in axes.' stop ELSE IF (abs(dotbc).GT.1.0d20) THEN write(nout,*)'dotbc messed up in axes.' stop END IF ! RETURN END subroutine axes ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! Subroutine to calculate integrand !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine calculate(graph,kin,abskin, & qs,qsigns,maptypes,nmaps,value,maxpart,valuechk) ! use beowulf_parameters use beowulf_structures implicit none ! In: type(graphstructure) :: graph real(kind=dbl) :: kin(0:3*size-1,0:3),abskin(0:3*size-1) integer :: qs(maxmaps,0:size),qsigns(maxmaps,0:size) character(len=6) :: maptypes(maxmaps) integer :: nmaps ! Out: complex(kind=dbl) :: value,valuechk real(kind=dbl) :: maxpart ! ! Calculates the value of the graph specified by VRTX at the point K, ! returning result in VALUE, which includes the division by the density ! of points and the jacobian for deforming the contour. Also reports ! MAXPART, the biggest absolute value of the contributions to Re(VALUE). ! This helps us to keep track of cancellations and thus to abort the ! calculation if too much cancellation among terms will be required. ! ! Latest revision: 11 May 1996 ! 24 October 1996 (call to CHECKDEFORM) ! 15 November 1996 (remove finite 'i epsilon') ! 18 November 1996 (add CHECKVALUE) ! 22 November 1996 Bug fixed. ! 27 November 1996 (complex checkvalue) ! 29 November 1996 (branchcut check; better checkvalue) ! 27 February 1997 renormalization; reporting ! 25 July 1997 renormalization; self-energy graphs ! 17 September 1997 more renormalization & self-energy ! 21 September 1997 finish DEFORM ! 24 September 1997 fix bugs ! 19 October 1997 fix cutsign bug ! 22 October 1997 fix renormalization sign bug ! 6 November 1997 improvements for deformation ! 28 November 1997 more work on deformation ! 2 December 1997 more precision in "report" numbers ! 4 January 1998 revisions for self-energy graphs ! 11 January 1998 renormalizaion for self-energy graphs ! 27 February 1998 use Hrothgar for output ! 5 March 1998 integrate Hrothgar ! 14 March 1998 restore checks of deformation direction ! 24 July 1998 use countfactor(graphnumber) ! 4 August 1998 better CHECKDEFORM ! 5 August 1998 change to groupsize(graphnumber) ! 22 August 1998 add color factors ! 22 December 1998 precalculate cut structure in RENO ! 26 April 1999 omit REFLECT except as option ! 22 December 2000 omit REFLECT entirely ! 22 December 2000 change method of choosing points ! 19 December 2001 call FEYNMANF, new organization. ! 31 December 2001 add Coulomb gauge. ! 11 February 2002 Add Born calculation. ! 14 June 2002 Revisions for f95 ! 30 June 2002 Calculate separately for each flavor set. ! 30 October 2002 Add showers. ! 16 May 2003 Showers with soft gluon radiation. ! 22 November 2003 Work on interface to Pythia ! 20 February 2004 more on interface with Pythia ! 1 October 2004 correlation with beam direction etc !---------------------------------- ! integer :: nloops,nprops,cutmax integer :: nloops1,nprops1,nverts1,cutmax1 integer :: nloops2,nprops2,nverts2,cutmax2 common /sizes/ nloops1,nprops1,nverts1,cutmax1, & nloops2,nprops2,nverts2,cutmax2 real(kind=dbl) :: muoverrts common /renormalize/ muoverrts logical :: report,details logical :: pythia character(len=10):: filename common /pythiainfo/ pythia,filename ! For pythia information common /calculook/ report,details real(kind=dbl) :: nc,nf,cf common /colorfactors/ nc,nf,cf character(len=7) :: gauge common /gaugechoice/ gauge ! Information on the current cut. type (cutstructure) :: cut logical :: cutfound logical :: keepcut ! Information on the current flavor set. type (flavorstructure) flavorset logical :: flavorsetfound character(len=5) :: partontype integer :: partonflavor ! Showers variables ! ! Below, theshower%multfactor is the factor that will multiply the ! Born graph. From subroutine softradiate, we start with ! fijk/(4 Pi * softintegral)/density when we add the soft ! gluon to the shower; then we further divide by the densities ! of points {qbarsq,x,phi} for each splitting in the first ! level splittings of the hard partons from a Born graph. ! The Sudakov factor is not included, nor are 1/qbarsq/(2 pi). ! type(showerlist) :: theshower ! real(kind=dbl) :: onemthrust,msoftsq ! Factor in msoftsq =(msoftfactor*rts0*onemthrust)**2 real(kind=dbl) :: msoftfactor common /softcutoff/ msoftfactor ! Factor for relating kappasq to onemthrust for small onemthrust real(kind=dbl) :: onemthrustcrit common /thrustfactor/ onemthrustcrit ! ! What the program should do character(len=8) :: mode common /programmode/ mode ! Physics data: real(kind=dbl) :: alphasofmz,sinsqthetaw,mz,widthz,externalrts common /physicsdata/ alphasofmz,sinsqthetaw,mz,widthz,externalrts ! ! Labels: integer :: qe(0:size) ! Momenta: real(kind=dbl) :: k(0:3*size-1,0:3),absk(0:3*size-1) real(kind=dbl) :: koff(0:3*size-1,0:3) real(kind=dbl) :: kinloop(size+1,0:3) real(kind=dbl) :: kcut(size+1,0:3) complex(kind=dbl) :: newkinloop(0:3) complex(kind=dbl) :: kc(0:3*size-1,0:3) complex(kind=dbl) :: ellsq,ell real(kind=dbl) :: e(0:size),rts,rts0 ! Squared momenta for factors F_j for showers real(kind=dbl) :: k1sq,k2sq ! Factors F_j for graph 9. real(kind=dbl) :: graph9factor ! Number of final state particles and their momenta integer :: nparticles real(kind=dbl) :: kf(maxparticles,0:3) ! Renormalization: real(kind=dbl) :: mumsbar ! Matrices: integer :: ae(0:3*size-1,0:size) ! FINDA variable: logical :: qok ! DENSITY variables: real(kind=dbl) :: jacnewpoint,density ! Loopcut variables: logical :: calcmore integer :: jcut,index,loopcutsign ! DEFORM variables: complex(kind=dbl) :: jacdeform ! Color connection information for final state integer, dimension(8) :: strings ! Functions: real(kind=dbl) :: cals0,smear real(kind=dbl) :: xxreal,xximag complex(kind=dbl) :: complexsqrt real(kind=dbl) :: alpi ! Index variables: integer :: p,mu,i,j,n ! Propagator properties logical :: cutQ(3*size-1) ! Flag for feynman function character(len=16) :: flag ! Results variables: real(kind=dbl) :: calsval,thrust real(kind=dbl) :: weight,maxweight complex(kind=dbl) :: feynman,feynmanf,feynman0,feynman0f,feynmanval complex(kind=dbl) :: softsubtraction real(kind=dbl) :: prefactor,prefactor0 complex(kind=dbl) :: integrand complex(kind=dbl) :: check ! Useful constants: real(kind=dbl), parameter :: metric(0:3) = (/1.0d0,-1.0d0,-1.0d0,-1.0d0/) real(kind=dbl), parameter :: pi = 3.141592653589793239d0 real(kind=dbl) :: temp ! For calculating matrix element square for shower from Born graph. integer :: nfound real(kind=dbl), dimension(3) :: qvec,kplusvec,kminusvec real(kind=dbl), dimension(0:3) :: vquark,vqbar real(kind=dbl), dimension(0:3,0:3) :: tglue integer :: n0,n1,n2 character(len=5) :: f1,f2 character(len=9) :: kind2pt complex(kind=dbl) :: feynmanSH0 ! Reno size and counting variables: integer :: groupsize(maxgraphs,maxmaps) integer :: groupsizegraph(maxgraphs) integer :: groupsizetotal common /montecarlo/groupsize,groupsizegraph,groupsizetotal ! type(softinformation), dimension(3) :: softinfos ! For makeshowerI type(softinformation) :: softinfo ! For shower...prop real(kind=dbl) :: cosine(3,3) integer :: j1,j2 ! !------------------------------------------------------------------------------- ! IF (graph%order.EQ.1) THEN nloops = nloops1 nprops = nprops1 cutmax = cutmax1 ELSE IF (graph%order.EQ.2) THEN nloops = nloops2 nprops = nprops2 cutmax = cutmax2 ELSE write(nout,*)'Order must be 1 or 2.',graph%order END IF ! ! We do not want to change the value of KIN and ABSKIN. (Actually, ! all that we change in the current version in k(p,0), so this step ! is not really necessary). ! DO p = 1,nprops absk(p) = abskin(p) DO mu = 0,3 k(p,mu) = kin(p,mu) END DO END DO ! ! Initialize contribution to integral from this point. Also initialize ! BIGGEST, which will be the biggest absolute value of the contributions ! to VALUE. This helps us to keep track of cancellations and thus to ! abort the calculation if too much cancellation among terms will be ! required. ! maxpart = 0.0d0 value = (0.0d0,0.0d0) valuechk = (0.0d0,0.0d0) ! ! Calculate jacobian. ! jacnewpoint = & 1.0d0/density(graph%graphnumber,k,qs,qsigns,maptypes,nmaps,graph%order) ! ! Loop over flavorsets. ! DO IF (gauge.EQ.'feynman') THEN ! In the case of Feynman gauge, we exit at the end of the first pass. flavorset%setnumber = 1 flavorset%type = & (/'none ','none ','none ','none ','none ','none ','none ','none '/) flavorset%flavor = & (/-777,-777,-777,-777,-777,-777,-777,-777/) ELSE IF (gauge.EQ.'coulomb') THEN call getnewflavorset(graph%graphnumber,flavorset,flavorsetfound) IF (.not.flavorsetfound) EXIT ELSE write(nout,*)'That gauge does not exist.' stop END IF ! ! Loop over cuts. ! DO IF ((mode.NE.'showerI ').AND.(mode.NE.'showerII')) THEN call getnewcut(graph%graphnumber,cut,cutfound) ELSE ! ! In the case of a shower mode, for certain cuts of certain graphs ! we keep only the wide angle splitting part, while other cuts are ! omitted entirely. Subroutine checknewcut finds out for us, and ! also calculates the factor F_1 or F_2 needed for graph 9. ! DO call getnewcut(graph%graphnumber,cut,cutfound) IF (.not.cutfound) EXIT call checknewcut(graph%graphnumber,cut,k,absk,keepcut,graph9factor) IF (keepcut) EXIT END DO ! END IF ! (mode.NE.'showerI ').AND.(mode.NE.'showerII') IF (.not.cutfound) EXIT !.... IF (report) THEN IF (graph%order.EQ.2) THEN write(nout,301)cut%ncut,cut%cutindex(1),cut%cutindex(2), & cut%cutindex(3),cut%cutindex(4) ELSE write(nout,301)cut%ncut,cut%cutindex(1),cut%cutindex(2), & cut%cutindex(3) END IF 301 format('ncut =',i2,' cutindex =',4i2) END IF !'''' ! ! Calculate final state momenta. ! Then we can also calculate CALSVAL and the PREFACTOR. ! DO i = 1,cut%ncut kcut(i,0) = 0.0d0 !We calculate the energy later. DO mu = 1,3 kcut(i,mu) = cut%cutsign(i) * k(cut%cutindex(i),mu) END DO END DO ! ! Set final state energies kcut(i,0) to their on-shell values and define ! final state momenta kf(i,mu) to be kcut(i,mu). Define rts to be the sum ! of these final state energies. Also, nparticles is the number of these ! final state particles. We define rts0 to be this "before showers" version ! of rts. If we later generate showers, we will overide these values. ! Look for "! Here we redefine rts, kcut(i,0), kf(i,mu), and nparticles." ! DO i = 1,cut%ncut kcut(i,0) = absk(cut%cutindex(i)) END DO nparticles = cut%ncut rts = 0.0d0 DO i = 1,nparticles DO mu = 0,3 kf(i,mu) = kcut(i,mu) END DO rts = rts + kf(i,0) END DO rts0 = rts mumsbar = muoverrts * rts0 ! The renormalization scale. ! ! Record starting values for the shower list. If we are not doing showering, ! just a plain purturbative calculation, then hrothgar will get this ! shower list as the representation of the final state. ! ! If we do generate showers, the variable onemthrust ! controls the width of showers generated by makeshowerII from ! graph%order.EQ.2 graphs and also the limit on 'soft' gluons in makeshowerI. ! Also, we need the soft scale and (1 - thrust). Here is where msoftsq ! is set, with a standard parameter msoftfactor set above. ! onemthrust = (1.0d0 - thrust(kf,nparticles)) IF (report) then write(nout,*)"1 - thrust before showering is ",onemthrust END IF msoftsq = (msoftfactor*rts0*onemthrust)**2 ! theshower%length = cut%ncut theshower%nstrings = 7777 ! This will be changed. theshower%rts0 = rts0 theshower%onemthrust = onemthrust theshower%msoftsq = msoftsq theshower%multfactor = 1.0d0 ! This will be changed. theshower%pii = 0 ! This will be changed. theshower%pjj = 0 ! This will be changed. DO n = 1,theshower%length partontype = flavorset%type(cut%cutindex(n)) partonflavor = flavorset%flavor(cut%cutindex(n)) IF (cut%cutsign(n).EQ.-1) THEN partonflavor = - partonflavor IF (partontype.EQ.'quark') THEN partontype = 'qbar ' ELSE IF (partontype.EQ.'qbar ') THEN partontype = 'quark' END IF END IF theshower%ptn(n)%type = partontype theshower%ptn(n)%flavor = partonflavor theshower%ptn(n)%stringquark = 500 ! If we make showers, this will change. theshower%ptn(n)%stringqbar = 500 ! If we make showers, this will change. ! theshower%ptn(n)%self = n theshower%ptn(n)%parent = 0 theshower%ptn(n)%ancestor = n theshower%ptn(n)%child1 = -1 theshower%ptn(n)%child2 = -1 theshower%ptn(n)%childless = .true. theshower%ptn(n)%done = .false. DO mu = 1,3 theshower%ptn(n)%momentum(mu) = kcut(n,mu) END DO IF (onemthrust.LT.onemthrustcrit) THEN ! default is onemthrustcrit = 0.05 ! for graph%order.EQ.1, theshower%ptn(n)%kappasq, n=1,2,3, is not used. theshower%ptn(n)%kappasq = kcut(n,0)**2 * onemthrust**2/onemthrustcrit ELSE theshower%ptn(n)%kappasq = kcut(n,0)**2 * onemthrust ENDIF END DO ! IF ((mode.EQ.'showerI ').OR.(mode.EQ.'showerII')) THEN ! ! Firt, string analysis for three parton final states. ! IF (theshower%length.EQ.3) THEN theshower%nstrings = 2 + 500 ! The 500 is to give strings labels 500 + n. DO n = 1,3 IF (theshower%ptn(n)%type.EQ.'gluon') THEN theshower%ptn(n)%stringquark = 502 theshower%ptn(n)%stringqbar = 501 ELSE IF (theshower%ptn(n)%type.EQ.'quark') THEN theshower%ptn(n)%stringquark = 501 theshower%ptn(n)%stringqbar = 0 ELSE IF (theshower%ptn(n)%type.EQ.'qbar ') THEN theshower%ptn(n)%stringquark = 0 theshower%ptn(n)%stringqbar = 502 ELSE write(nout,*)'Confused about parton type in calculate. ',partontype STOP END IF END DO END IF ! (theshower%length.EQ.3) ! ! String analysis for four parton final states. ! IF (theshower%length.EQ.4) THEN theshower%nstrings = 3 + 500 ! The 500 is to give strings labels 500 + n. call getcolors(graph%graphnumber,flavorset%setnumber, & cut%cutnumber,strings) theshower%ptn(1)%stringquark = strings(1) theshower%ptn(1)%stringqbar = strings(2) theshower%ptn(2)%stringquark = strings(3) theshower%ptn(2)%stringqbar = strings(4) theshower%ptn(3)%stringquark = strings(5) theshower%ptn(3)%stringqbar = strings(6) theshower%ptn(4)%stringquark = strings(7) theshower%ptn(4)%stringqbar = strings(8) END IF ! (theshower%length.EQ.4) ! IF (graph%order.EQ.1) THEN ! ! We will create showering. ! ! For each of the final hard state partons i = 1,2,3, we need its parton type, ! which we call softinfos(i)%type0. We also need the types of ! the other two, which we call softinfos(i)%type1 and softinfos(i)%type2 ! and we need the cosine of the angle from the direction of parton i to the ! directions of these two other partons, which we call softinfos(i)%cos1 ! and softinfo(i)s%cos2. This information will be used for the soft gluon ! subtractions. We use the starting values of the final state momenta for ! this, with exact three body kinematics. ! DO i = 1,3 DO j = i+1,3 cosine(i,j) = 0.0d0 DO mu = 1,3 cosine(i,j) = cosine(i,j) + kcut(i,mu) * kcut(j,mu) END DO cosine(i,j) = cosine(i,j)/kcut(i,0)/kcut(j,0) cosine(j,i) = cosine(i,j) END DO END DO ! DO i = 1,3 j1 = mod(i,3) + 1 ! 1 -> 2, 2 -> 3, 3 -> 1 j2 = mod(i+1,3) + 1 ! 1 -> 3, 2 -> 1, 3 -> 2 softinfos(i)%cos1 = cosine(i,j1) softinfos(i)%cos2 = cosine(i,j2) softinfos(i)%type0 = theshower%ptn(i)%type softinfos(i)%type1 = theshower%ptn(j1)%type softinfos(i)%type2 = theshower%ptn(j2)%type softinfos(i)%msoftsq = theshower%msoftsq END DO ! call softradiate(theshower) ! Radiate a soft gluon ! Note that subroutine softradiate has shifted the parton momenta ! so we modify kcut accordingly. DO i = 1,4 DO mu = 1,3 kcut(i,mu) = theshower%ptn(i)%momentum(mu) END DO END DO ! call makeshowerI(theshower,softinfos) ! First splittings. IF (mode.EQ.'showerII') THEN call makeshowerII(theshower) ! Shower with small angle approx. END IF ELSE ! Order is 2. We make a shower from the final state particles, but ! only if we have enabled soft gluon subtractions, or else cancellations ! will fail. IF (mode.EQ.'showerII') THEN call makeshowerII(theshower) ! Shower with small angle approx. END IF END IF ! ! Now we have a primary shower if the graph was first order plus a ! secondary shower if desired. We need to find the final state partons. ! kf(i,mu) is the four-momentum of the ith final state parton. We ! compute rts as the sum of the energies of all the final state partons. ! Also, we compute the energy, kcut(j,0) of the jth original cut ! parton as sum of the energies of its descendants in the final state. ! Here we redefine rts, kcut(i,0), kf(i,mu), and nparticles. ! rts = 0.0d0 DO i = 1,cut%ncut kcut(i,0) = 0.0d0 END DO IF (graph%order.EQ.1) THEN kcut(4,0) = 0.0d0 ! initialization also for the extra soft gluon. END IF i = 0 DO n = 1,theshower%length IF (theshower%ptn(n)%childless) THEN i = i+1 temp = 0.0d0 DO mu = 1,3 kf(i,mu) = theshower%ptn(n)%momentum(mu) temp = temp + kf(i,mu)**2 END DO temp = sqrt(temp) kf(i,0) = temp rts = rts + temp j = theshower%ptn(n)%ancestor kcut(j,0) = kcut(j,0) + temp END IF END DO nparticles = i ! IF (report) then write(nout,*)"1 - thrust after showering is ", & (1.0d0 - thrust(kf,nparticles)) END IF ! ! We also calculate momenta koff(p,mu), which is like k(p,mu) except that ! the final state particles are given energies equal to the total energy ! of the shower that originates from that particle and the momentum of ! the soft gluon shower is appropriately routed. This is only for the ! Born graphs. ! IF (graph%order.EQ.1) THEN call findkoff(graph%graphnumber,cut,kcut,theshower%pii,theshower%pjj,koff) DO mu = 1,3 IF (abs(koff(0,mu)).GT.1.0d-12) THEN write(nout,*)'Something is rotten in findkoff.', & graph%graphnumber,cut%cutindex,theshower%pii,theshower%pjj stop END IF END DO IF ( abs(rts - koff(0,0)).GT.1.0d-8 ) THEN write(nout,*)'Oops, the calculation of rts did not work.',rts,koff(0,0) stop END IF END IF ! (graph%order.EQ.1) ! END IF ! ((mode.EQ.'showerI ').OR.(mode.EQ.'showerII')) ! calsval = cals0(nparticles,kf) prefactor = 1.0d0 / (nc * rts**2 * (2.0d0 * pi)**nloops ) ! IF (mode.EQ.'born ') THEN prefactor = prefactor * alpi(muoverrts*externalrts) ELSE IF (mode.EQ.'nlo ') THEN prefactor = prefactor * alpi(muoverrts*externalrts)**graph%order ELSE IF (mode.EQ.'hocoef ') THEN continue ELSE IF ((mode.EQ.'showerI ').OR.(mode.EQ.'showerII')) THEN prefactor = prefactor * alpi(muoverrts*externalrts)**graph%order IF (graph%graphnumber.EQ.9) THEN prefactor = prefactor * graph9factor END IF IF (graph%order.EQ.1) THEN prefactor = prefactor*theshower%multfactor prefactor0 = prefactor ! Use this to modify prefactor for graph 11. END IF ELSE write(nout,*)'calculate not programmed for this mode. ',mode stop END IF ! ! Calculate momenta around loop (if any). In case NINLOOP = 0, this ! DO loop is skipped. ! DO j = 1,cut%ninloop DO mu = 1,3 kinloop(j,mu) = cut%loopsign(j) * k(cut%loopindex(j),mu) END DO END DO ! ! Please note that at this point the energy in the loop, KINLOOP(J,0), ! is not calculated. We have to wait until we have a loop cut to ! do this. ! ! Now KINLOOP(J,MU) gets an imaginary part for MU = 1,2,3. ! DEFORM calculates NEWKINLOOP and the associated jacobian, JACDEFORM. ! In case NINLOOP = 0, this subroutine just returns NEWKINLOOP(MU) = 0 ! and JACDEFORM = 1. ! call deform(graph%vrtx,cut%loopindex,rts,cut%leftloop,cut%rightloop, & cut%ninloop,kinloop,newkinloop,jacdeform) ! ! If there is a loop, we need to go around the loop and generate ! a "loopcut." There are four cases. ! ! 1) NINLOOP = 0, with NCUT = CUTMAX. ! Then we are ready to proceed, and we should calculate only once ! before going back to NEWCUT. Therefore we set CALCMORE to .FALSE. ! so that we do not enter this code again. ! ! In the other three cases, there is a loop with NINLOOP = 2, 3, or 4. ! We generate a loopcut specified by the index JCUT = 1, 2, ... around ! the loop: CUTINDEX(CUTMAX) = LOOPINDEX(JCUT). ! ! 2) NINLOOP = 2, with NCUT = CUTMAX - 1. ! Then the loop is a self-energy subgraph and, with our dispersive ! treatment of these graphs, there is one term. We need to calculate ! energies, so we put JCUT = 1, but this is just a convention: the ! choice JCUT = 1 or 2 affects only the energy in the loop and the ! two point function depends only on the 3-momentum in the loop. ! ! 3) NINLOOP = 3, with NCUT = CUTMAX - 1. ! We will go around the loop twice with JCUT = 1,2,3 and then ! with JCUT = 1,2,3 again. We set CUTSIGN(CUTMAX) = LOOPSIGN(JCUT) ! the first time and CUTSIGN(CUTMAX) = - LOOPSIGN(JCUT) the second time. ! This corresponds to doing the energy integral with the contour closed in ! the upper half plane and also in the lower half plane. We *average* over ! the two sign choices, so we will need to multiply FEYNMAN by 1/2 for ! NINLOOP = 3. When we are done with this we set CALCMORE to .FALSE. . ! We will also need the renormaliztion counter term. This is done at the ! end. Before the IF(calcmore)... loop is closed, we test for calcmore ! and if it is false we check whether graph%order is 2 and cut%ninloop ! is 3. If so, we calculate and add the counter term. There is a flag ! for this purpose: if flag = ' no flag set', then the functions ! fenyman and feynmanf calculate the normal terms, but if ! flag = 'renormalize 3 pt', these functions calculate the counter terms. ! ! 3) NINLOOP = 4, with NCUT = CUTMAX - 1. ! We will go around the loop twice with JCUT = 1,2,3,4 and then with ! JCUT = 1,2,3,4 again. We set CUTSIGN(CUTMAX) = LOOPSIGN(JCUT) the first ! time and CUTSIGN(CUTMAX) = - LOOPSIGN(JCUT) the second time. This ! corresponds to doing the energy integral with the contour closed in the ! upper half plane and also in the lower half plane. We *average* over the ! two sign choices, so we will need to multiply FEYNMAN by 1/2 for ! NINLOOP = 4. When we are done with this we set CALCMORE to .FALSE. . ! ! We initialize the weight, then add to it the contributions from ! each pass through this loop. In the last pass, we also add the ! renormalization counter term to the weight, if appropriate. ! weight = 0.0d0 maxweight = 0.0d0 ! index = 1 ! to count where we are calcmore = .true. ! DO WHILE (calcmore) ! IF (cut%ninloop.EQ.0) THEN calcmore = .false. ELSE IF (cut%ninloop.EQ.2) THEN jcut = 1 cut%cutindex(cutmax) = cut%loopindex(jcut) cut%cutsign(cutmax) = cut%loopsign(jcut) loopcutsign = 1 calcmore = .false. ELSE IF (cut%ninloop.EQ.3) THEN IF (index.LE.3) THEN jcut = index cut%cutindex(cutmax) = cut%loopindex(jcut) cut%cutsign(cutmax) = cut%loopsign(jcut) loopcutsign = 1 ELSE jcut = index - 3 cut%cutindex(cutmax) = cut%loopindex(jcut) cut%cutsign(cutmax) = - cut%loopsign(jcut) loopcutsign = -1 END IF index = index + 1 IF (index.GT.6) THEN calcmore = .false. END IF ELSE IF (cut%ninloop.EQ.4) THEN IF (index.LE.4) THEN jcut = index cut%cutindex(cutmax) = cut%loopindex(jcut) cut%cutsign(cutmax) = cut%loopsign(jcut) loopcutsign = 1 ELSE jcut = index - 4 cut%cutindex(cutmax) = cut%loopindex(jcut) cut%cutsign(cutmax) = - cut%loopsign(jcut) loopcutsign = -1 END IF index = index + 1 IF (index.GT.8) THEN calcmore = .false. END IF ELSE write(nout,*)'Impossible case in calculate.' stop END IF ! ! Calculate matrix AE(P,I) relating propagator energies to energies of ! cut lines. NOTE that the index I here is displaced by 1. ! DO i = 0,nloops qe(i) = cut%cutindex(i+1) END DO call finda(graph%vrtx,qe,nloops,graph%order,ae,qok) IF (.NOT.qok) THEN write(nout,*)'ae not found.' write(nout,*)'Graph number',graph%graphnumber,' cut number',cut%cutnumber write(nout,*) qe stop END IF ! ! Define logical cut variables: ! cutQ(P) = .TRUE. if propagator P crosses the final state cut ! OR if it crosses the loopcut. ! DO p = 1,nprops cutQ(p) = .false. END DO DO i = 1,cutmax cutQ(cut%cutindex(i)) = .true. END DO ! ! Calculate part of the propagator energies corresponding to the ! real part of the loop three-momenta. Here we calculate for with the ! energy of each final state particle from the original graph set ! equal to the absolute value of its momentum. That is, the mass of ! the shower originating from this particle (if any) is neglected. ! ! Note that i is displaced by 1 in the e(i) array in order to work with ! the matrix AE(p,i), for which i should run from 0 to cutmax - 1. ! DO i = 1,cutmax e(i-1) = cut%cutsign(i) * absk(cut%cutindex(i)) END DO ! DO p = 0,nprops k(p,0) = 0.0d0 DO i = 0,nloops k(p,0) = k(p,0) + ae(p,i) * e(i) END DO END DO IF ( abs(rts0 - k(0,0)).GT.1.0d-8 ) THEN write(nout,*)'Oops, the calculation of rts0 did not work.',rts0,k(0,0) stop END IF ! ! Calculate the added complex loop energy. Check that we do not ! cross the cut of Sqrt(ELLSQ) by using COMPLEXSQRT(ELLSQ). ! IF (cut%ninloop.GT.0) THEN kinloop(jcut,0) = cut%loopsign(jcut) * k(cut%loopindex(jcut),0) ellsq = (0.0d0,0.0d0) DO mu = 1,3 ellsq = ellsq + ( kinloop(jcut,mu) + newkinloop(mu) )**2 END DO ell = complexsqrt(ellsq) newkinloop(0) = loopcutsign*ell - kinloop(jcut,0) ELSE newkinloop(0) = (0.0d0,0.0d0) END IF !.... IF (report) THEN IF( details .AND. (cut%ninloop.GT.0) ) THEN write(nout,340)newkinloop(0),xximag(newkinloop(1)), & xximag(newkinloop(2)),xximag(newkinloop(3)) 340 format('newkinloop =',2(1p g12.3),' and',3(1p g12.3)) END IF END IF !'''' ! Now we calculate the complex propagator momenta, kc(p,mu). These are ! based on k(p,mu) for which the energy of each final state particle ! from the original graph is set equal to the absolute value of its ! momentum. ! DO p = 0,nprops DO mu = 0,3 kc(p,mu) = k(p,mu) END DO END DO ! DO j = 1,cut%ninloop DO mu = 0,3 kc(cut%loopindex(j),mu) = kc(cut%loopindex(j),mu) & + cut%loopsign(j) * newkinloop(mu) END DO END DO ! IF ((mode.EQ.'showerI ').OR.(mode.EQ.'showerII')) THEN ! ! There are some preliminary calculations for the case of shower mode. ! ! First, we multiply by factor F_j for showers for the first Born ! graph, graphnumber 11. To compute F_j, we use momenta kcut(i,mu), ! which represent the momenta of the final state particles 1,2,3,4 ! including the showers that develop from them. The graph structure ! for graph 11 is that particles 1 and 2 form a self cut energy ! subdiagram, particle 3 is the third hard parton in the Born graph ! and particle 4 is the soft gluon. We use the momentum of particle 4 ! in the calculation of F_j if it forms a self-energy subdiagram. ! Then we add it as a soft parton and subtract the same diagram ! from the splitting, so we want to treat it kinematically just ! like a splitting. ! IF (graph%graphnumber.EQ.11) THEN IF ((cut%cutindex(1).NE.5).OR. & (cut%cutindex(2).NE.4).OR. & (cut%cutindex(3).NE.1)) THEN write(nout,*)'Bad labelling in F_j calculation' stop END IF k1sq = 0.0d0 k2sq = 0.0d0 IF ( (theshower%pii.EQ.1).AND.(theshower%pjj.EQ.1) ) THEN DO mu = 0,3 k1sq = k1sq + kcut(3,mu)**2 * metric(mu) k2sq = k2sq + (kcut(1,mu) + kcut(2,mu) + kcut(4,mu))**2 * metric(mu) END DO ELSE IF ( (theshower%pii.EQ.2).AND.(theshower%pjj.EQ.2) ) THEN DO mu = 0,3 k1sq = k1sq + kcut(3,mu)**2 * metric(mu) k2sq = k2sq + (kcut(1,mu) + kcut(2,mu) + kcut(4,mu))**2 * metric(mu) END DO ELSE IF ( (theshower%pii.EQ.3).AND.(theshower%pjj.EQ.3) ) THEN DO mu = 0,3 k1sq = k1sq + (kcut(3,mu) + kcut(4,mu))**2 * metric(mu) k2sq = k2sq + (kcut(1,mu) + kcut(2,mu))**2 * metric(mu) END DO ELSE DO mu = 0,3 k1sq = k1sq + kcut(3,mu)**2 * metric(mu) k2sq = k2sq + (kcut(1,mu) + kcut(2,mu))**2 * metric(mu) END DO END IF prefactor = prefactor0 * k2sq**2/(k1sq**2 + k2sq**2) END IF ! IF (graph%order.EQ.1) THEN ! ! Calculate tglue, vquark, vqbar, which are needed for the first ! level of showering from a Born graph. ! DO n0=1,3 ! We want the three hard partons. qvec = theshower%ptn(n0)%momentum n1 = theshower%ptn(n0)%child1 n2 = theshower%ptn(n0)%child2 f1 = theshower%ptn(n1)%type f2 = theshower%ptn(n2)%type kplusvec = theshower%ptn(n1)%momentum kminusvec = theshower%ptn(n2)%momentum ! We need the softinfo for parton n0. We can tell by its type. nfound = 0 DO j = 1,3 IF (softinfos(j)%type0.EQ.theshower%ptn(n0)%type) THEN nfound = nfound + 1 softinfo = softinfos(j) END IF END DO IF (nfound.NE.1) THEN write(nout,*) 'Oops, types messed up in the shower' STOP END IF IF (theshower%ptn(n0)%type.EQ.'gluon') THEN IF ((f1.EQ.'quark').AND.(f2.EQ.'qbar ')) THEN kind2pt = 'quarkloop' ELSE IF ((f1.EQ.'gluon').AND.(f2.EQ.'gluon')) THEN kind2pt = 'gluonloop' ELSE write(nout,*)'Oops, that parton type combination should not exist.' STOP END IF call showerglueprop(kind2pt,softinfo,qvec,kplusvec,kminusvec,tglue) ELSE IF (theshower%ptn(n0)%type.EQ.'quark') THEN IF (.NOT.(f1.EQ.'gluon').AND.(f2.EQ.'quark')) THEN write(nout,*)'Oops, that parton type combination should not exist.' STOP END IF call showerquarkprop(softinfo,qvec,kplusvec,kminusvec,vquark) ELSE IF (theshower%ptn(n0)%type.EQ.'qbar ') THEN IF (.NOT.(f1.EQ.'gluon').AND.(f2.EQ.'qbar ')) THEN write(nout,*)'Oops, that parton type combination should not exist.' STOP END IF call showerquarkprop(softinfo,qvec,kplusvec,kminusvec,vqbar) DO mu = 0,3 vqbar(mu) = - vqbar(mu) END DO ELSE write(nout,*)'The parton type must be gluon, quark, or qbar.' STOP END IF END DO ! Close DO n0=1,3 ! END IF ! (graph%order.EQ.1) END IF ! ((mode.EQ.'showerI ').OR.(mode.EQ.'showerII')) ! ! Calculate graph and add to contribution for this point. ! IF (graph%order.EQ.2) THEN flag = ' no flag set' ! This is not a counter-term. IF (gauge.EQ.'feynman') THEN feynmanval = feynmanf(graph%graphnumber,kc,cutQ,mumsbar,flag) ELSE IF (gauge.EQ.'coulomb') THEN feynmanval = & feynman(graph%graphnumber,flavorset%setnumber,kc,cutQ,mumsbar,flag) ELSE write(nout,*)'That gauge does not exist.' stop END IF ELSE IF (graph%order.EQ.1) THEN ! ! There are a lot of alternatives in the case of a Born graph. ! IF (gauge.EQ.'feynman') THEN ! ! This gets the Feynman gauge matrix element squared. ! feynmanval = feynman0f(graph%graphnumber,k,cutQ) ! ELSE IF (gauge.EQ.'coulomb') THEN ! IF ((mode.EQ.'showerI ').OR.(mode.EQ.'showerII')) THEN ! ! We have previously calculated vquark,vqbar, and tglue. ! Here we use koff(p,mu) which is like k(p,mu) except that ! the final state particles are given energies equal to the total ! energy of the shower that originates from that particle. ! This gets the appropriate matrix element squared with a one ! level shower. ! feynmanval = feynmanSH0(graph%graphnumber,flavorset%setnumber, & koff,cutQ,vquark,vqbar,tglue) ELSE ! ! We have Coulomb gauge but not mode .EQ. 'showerI' or 'showerII'. ! This gets the Coulomb gauge matrix element squared. ! feynmanval = & feynman0(graph%graphnumber,flavorset%setnumber,k,cutQ) ! ! Close IF ((mode.EQ.'showerI ').OR.(mode.EQ.'showerII')) THEN ... ELSE ... ! END IF ! ELSE write(nout,*)'That gauge does not exist.' stop ! ! Close ! IF (gauge.EQ.'feynman') THEN ... ! ELSE IF (gauge.EQ.'coulomb') THEN ... ! ELSE ... END IF ELSE write(nout,*)'Order should have been 1 or 2.' stop ! ! Close !IF (graph%order.EQ.2) THEN ... !ELSE IF (graph%order.EQ.1) THEN ... !ELSE END IF ! integrand = prefactor * jacnewpoint * jacdeform * feynmanval * smear(rts) ! ! If we have a 3 or 4 point virtual loop, then we are averaging over ! closing the energy integral contour in the upper and lower ! half planes and we supply a 1/2. ! IF (cut%ninloop.GT.2) THEN integrand = 0.5d0*integrand END IF ! maxweight = max(maxweight,abs(xxreal(integrand))) weight = weight + xxreal(integrand) ! integrand = integrand * calsval maxpart = max(maxpart,abs(xxreal(integrand))) value = value + integrand !.... IF (report) THEN IF (details) THEN write(nout,370) 370 format('prefactor * jacnewpoint * (jacdeform-r jacdeform-i) * ', & ' (feynman-r feynman-i) * calsval * smear(rts)') write(nout,371)prefactor,jacnewpoint,jacdeform, & feynmanval,calsval,smear(rts) 371 format(9(1p g12.3)) END IF IF (cut%ninloop.GT.0) THEN write(nout,372)cut%loopindex(jcut),integrand*groupsizetotal 372 format(i3,' contribution:',2(1p g18.10)) ELSE write(nout,373)integrand*groupsizetotal 373 format(' contribution:',2(1p g18.10)) END IF IF (details) THEN write(nout,*)' ' END IF END IF !'''' ! ! Compute a known integral to see if we have it right. ! Subroutine CHECKCALC calculates CHECK. ! IF (flavorset%setnumber.EQ.1) THEN IF (flag.NE.'renormalize 3 pt') THEN call & checkcalc(graph%graphnumber,cut%cutindex,kc,jacnewpoint,jacdeform,check) IF (cut%ninloop.GT.2) THEN check = 0.5d0*check END IF valuechk = valuechk + check END IF END IF ! ! IF(.NOT.calcmore) THEN ! ! We are almost done with this cut. However, there may be a renormalization ! counter term and/or a soft subtraction needed. ! !.... IF (report) THEN write(nout,341)weight*calsval*groupsizetotal 341 format(' cut total so far:',1p g18.10) IF (graph%order.EQ.1) THEN write(nout,*)'Point:' DO p=1,nprops write(nout,1703) p,koff(p,0),koff(p,1),koff(p,2),koff(p,3) 1703 format('p =',i2,' koff = ',4(1p g12.3)) END DO write(nout,*)'pii = ',theshower%pii,' pjj = ',theshower%pjj write(nout,1704) kcut(4,0),kcut(4,1),kcut(4,2),kcut(4,3) 1704 format('kcut(4,mu) = ',4(1p g12.3)) END IF END IF !.... ! IF (graph%order.EQ.2) THEN IF (cut%ninloop.EQ.3) THEN ! ! We need the renormalization counter term: ! flag = 'renormalize 3 pt' IF (gauge.EQ.'feynman') THEN feynmanval = feynmanf(graph%graphnumber,kc,cutQ,mumsbar,flag) ELSE IF (gauge.EQ.'coulomb') THEN feynmanval = & feynman(graph%graphnumber,flavorset%setnumber,kc,cutQ,mumsbar,flag) ELSE write(nout,*)'That gauge does not exist.' stop END IF integrand = prefactor * jacnewpoint * jacdeform * feynmanval * smear(rts) maxweight = max(maxweight,abs(xxreal(integrand))) weight = weight + xxreal(integrand) integrand = integrand * calsval maxpart = max(maxpart,abs(xxreal(integrand))) value = value + integrand !.... IF (report) THEN IF (details) THEN write(nout,370) write(nout,371)prefactor,jacnewpoint,jacdeform, & feynmanval,calsval,smear(rts) END IF write(nout,374)integrand*groupsizetotal 374 format(' uv counter term:',2(1p g18.10)) write(nout,341)weight*calsval*groupsizetotal END IF !'''' ! END IF ! (cut%ninloop.EQ.3) ! IF ((mode.EQ.'showerI ').OR.(mode.EQ.'showerII')) THEN ! ! We need the soft subtraction: ! feynmanval = softsubtraction(k,absk, & graph%graphnumber,flavorset%setnumber,cut%cutnumber) integrand = prefactor * jacnewpoint * feynmanval * smear(rts) maxweight = max(maxweight,abs(xxreal(integrand))) weight = weight + xxreal(integrand) integrand = integrand * calsval maxpart = max(maxpart,abs(xxreal(integrand))) value = value + integrand !.... IF (report) THEN IF (details) THEN write(nout,416) 416 format('prefactor * jacnewpoint *', & ' (feynman-r feynman-i) * calsval * smear(rts)') write(nout,417)prefactor,jacnewpoint,feynmanval,calsval,smear(rts) 417 format(7(1p g12.3)) END IF write(nout,375)integrand*groupsizetotal 375 format(' soft subtraction:',2(1p g18.10)) write(nout,341)weight*calsval*groupsizetotal END IF !.... END IF ! ((mode.EQ.'showerI ').OR.(mode.EQ.'showerII')) ! END IF ! (graph%order.EQ.2) ! END IF ! (.NOT.calcmore) ! ! End of loop DO WHILE (CALCMORE) that runs over loopcuts. ! END DO ! ! We are ready to call Hrothgar to process the result for this cut. ! call hrothgar(theshower,weight,1,'newresult ') ! ! Close loop DO ... call getnewcut(); IF (.NOT.cutfound) EXIT ! END DO ! ! Close loop DO with call getnewflavorset(); IF (.NOT.flavorsetfound) EXIT ! for Coulomb gauge. ! IF (gauge.EQ.'feynman') EXIT END DO ! RETURN END subroutine calculate ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! End of subroutine to calculate integrand !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine & checkcalc(graphnumber,cutindex,kc,jacnewpoint,jacdeform,check) ! use beowulf_parameters implicit none ! In: integer :: graphnumber,cutindex(size+1) complex(kind=dbl) :: kc(0:3*size-1,0:3) real(kind=dbl) :: jacnewpoint complex(kind=dbl) :: jacdeform ! Out: complex(kind=dbl) :: check ! ! Compute a known integral to see if we have it right. ! This subroutine calculates the integrand. ! The check is based on ! Int d^3 p [p^2 + M^2]^(-3) = Pi^2/ (4 M^3). ! Int d^3 p [p^2 (p^2 + M^2)]^(-1) = 2 Pi^2 /M ! Note that we look at just one term in the sum over cuts ! and loopcuts: ! For graph 10, we take Cutindex = (7,5,4,1); ! For graph 8, we take Cutindex = (8,6,4,1), etc. ! ! Latest modification: 11 February 2002. ! ! Reno size and counting variables: integer :: groupsize(maxgraphs,maxmaps) integer :: groupsizegraph(maxgraphs) integer :: groupsizetotal common /montecarlo/groupsize,groupsizegraph,groupsizetotal ! real(kind=dbl), parameter :: mm = 3.0d-1 real(kind=dbl), parameter :: pi = 3.141592653589793239d0 ! complex(kind=dbl) :: temp1,temp2,temp3 integer :: mu ! ! If it is not the right graph and the right cut, this default ! value will be returned. ! check = (0.0d0,0.0d0) ! temp1 = 0.0d0 temp2 = 0.0d0 temp3 = 0.0d0 ! IF (graphnumber.EQ.12) THEN ! IF ( (cutindex(1).EQ.5).AND.(cutindex(2).EQ.4) & .AND.(cutindex(3).EQ.1) ) THEN DO mu = 1,3 temp1 = temp1 + kc(5,mu)*kc(5,mu) temp2 = temp2 + kc(4,mu)*kc(4,mu) END DO ELSE RETURN END IF ! ELSE IF (graphnumber.EQ.11) THEN ! IF ( (cutindex(1).EQ.5).AND.(cutindex(2).EQ.4) & .AND.(cutindex(3).EQ.1) ) THEN DO mu = 1,3 temp1 = temp1 + kc(5,mu)*kc(5,mu) temp2 = temp2 + kc(4,mu)*kc(4,mu) END DO ELSE RETURN END IF ! ELSE IF (graphnumber.EQ.10) THEN ! IF ( (cutindex(1).EQ.7).AND.(cutindex(2).EQ.5) & .AND.(cutindex(3).EQ.4).AND.(cutindex(4).EQ.1) ) THEN DO mu = 1,3 temp1 = temp1 + kc(7,mu)*kc(7,mu) temp2 = temp2 + kc(6,mu)*kc(6,mu) temp3 = temp3 + kc(1,mu)*kc(1,mu) END DO ELSE RETURN END IF ! ELSE IF (graphnumber.EQ.9) THEN ! IF ( (cutindex(1).EQ.8).AND.(cutindex(2).EQ.7) & .AND.(cutindex(3).EQ.3).AND.(cutindex(4).EQ.5) ) THEN DO mu = 1,3 temp1 = temp1 + kc(8,mu)*kc(8,mu) temp2 = temp2 + kc(6,mu)*kc(6,mu) temp3 = temp3 + kc(1,mu)*kc(1,mu) END DO ELSE RETURN END IF ! ELSE IF (graphnumber.EQ.8) THEN ! IF ( (cutindex(1).EQ.8).AND.(cutindex(2).EQ.6) & .AND.(cutindex(3).EQ.4).AND.(cutindex(4).EQ.1) ) THEN DO mu = 1,3 temp1 = temp1 + kc(5,mu)*kc(5,mu) temp2 = temp2 + kc(8,mu)*kc(8,mu) temp3 = temp3 + kc(1,mu)*kc(1,mu) END DO ELSE RETURN END IF ! ELSE IF (graphnumber.EQ.7) THEN ! IF ( (cutindex(1).EQ.5).AND.(cutindex(2).EQ.4) & .AND.(cutindex(3).EQ.1).AND.(cutindex(4).EQ.6) ) THEN DO mu = 1,3 temp1 = temp1 + kc(5,mu)*kc(5,mu) temp2 = temp2 + kc(7,mu)*kc(7,mu) temp3 = temp3 + kc(1,mu)*kc(1,mu) END DO ELSE RETURN END IF ! ELSE IF (graphnumber.EQ.6) THEN ! IF ( (cutindex(1).EQ.8).AND.(cutindex(2).EQ.5) & .AND.(cutindex(3).EQ.3).AND.(cutindex(4).EQ.6) ) THEN DO mu = 1,3 temp1 = temp1 + kc(5,mu)*kc(5,mu) temp2 = temp2 + kc(7,mu)*kc(7,mu) temp3 = temp3 + kc(1,mu)*kc(1,mu) END DO ELSE RETURN END IF ! ELSE IF (graphnumber.EQ.5) THEN ! IF ( (cutindex(1).EQ.8).AND.(cutindex(2).EQ.7) & .AND.(cutindex(3).EQ.3).AND.(cutindex(4).EQ.1) ) THEN DO mu = 1,3 temp1 = temp1 + kc(5,mu)*kc(5,mu) temp2 = temp2 + kc(8,mu)*kc(8,mu) temp3 = temp3 + kc(1,mu)*kc(1,mu) END DO ELSE RETURN END IF ! ELSE IF (graphnumber.EQ.4) THEN ! IF ( (cutindex(1).EQ.6).AND.(cutindex(2).EQ.4) & .AND.(cutindex(3).EQ.1).AND.(cutindex(4).EQ.7) ) THEN DO mu = 1,3 temp1 = temp1 + kc(7,mu)*kc(7,mu) temp2 = temp2 + kc(5,mu)*kc(5,mu) temp3 = temp3 + kc(1,mu)*kc(1,mu) END DO ELSE RETURN END IF ! ELSE IF (graphnumber.EQ.3) THEN ! IF ( (cutindex(1).EQ.5).AND.(cutindex(2).EQ.4) & .AND.(cutindex(3).EQ.1).AND.(cutindex(4).EQ.6) ) THEN DO mu = 1,3 temp1 = temp1 + kc(7,mu)*kc(7,mu) temp2 = temp2 + kc(5,mu)*kc(5,mu) temp3 = temp3 + kc(1,mu)*kc(1,mu) END DO ELSE RETURN END IF ! ELSE IF (graphnumber.EQ.2) THEN ! IF ( (cutindex(1).EQ.7).AND.(cutindex(2).EQ.6) & .AND.(cutindex(3).EQ.1).AND.(cutindex(4).EQ.4) ) THEN DO mu = 1,3 temp1 = temp1 + kc(7,mu)*kc(7,mu) temp2 = temp2 + kc(5,mu)*kc(5,mu) temp3 = temp3 + kc(1,mu)*kc(1,mu) END DO ELSE RETURN END IF ! ELSE IF (graphnumber.EQ.1) THEN ! IF ( (cutindex(1).EQ.5).AND.(cutindex(2).EQ.4) & .AND.(cutindex(3).EQ.1).AND.(cutindex(4).EQ.7) ) THEN DO mu = 1,3 temp1 = temp1 + kc(7,mu)*kc(7,mu) temp2 = temp2 + kc(5,mu)*kc(5,mu) temp3 = temp3 + kc(1,mu)*kc(1,mu) END DO ELSE RETURN END IF ! ELSE write(nout,*)'Problem with graph number in checkcalc.' stop END IF ! IF (graphnumber.LE.10) THEN ! ! Here is an infrared sensitive check integral: ! check = temp1 * (temp1 + mm**2) check = check * temp2 * (temp2 + mm**2) check = check * (temp3 + mm**2)**3 check = (mm**5/pi**6) /check ! ! Here is a nice smooth check integral: ! CHECK = (TEMP1 + MM**2)**3 ! CHECK = CHECK * (TEMP2 + MM**2)**3 ! CHECK = CHECK * (TEMP3 + (2.0D0*MM)**2)**3 ! CHECK = (512.0D0 * MM**9 / PI**6) /CHECK ! ELSE IF (graphnumber.LE.12) THEN ! check = (temp1 + mm**2)**3 check = check * (temp2 + mm**2)**3 check = (16.0d0 * mm**6 / pi**4) /check ! ELSE write(nout,*)'We were expecting graphnumbers 1,...,12.' stop END IF ! check = jacdeform * jacnewpoint * check ! ! Weight according to the number of points devoted to the current ! graph. ! check = check * groupsizegraph(graphnumber)/groupsizetotal ! RETURN END subroutine checkcalc ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! function density(graphnumber,k,qs,qsigns,maptypes,nmaps,order) ! use beowulf_parameters implicit none ! In: integer :: graphnumber real(kind=dbl) :: k(0:3*size-1,0:3) integer :: qs(maxmaps,0:size),qsigns(maxmaps,0:size) character(len=6) :: maptypes(maxmaps) integer :: nmaps,order ! Out: real(kind=dbl) :: density ! ! Density of Monte Carlo points as a function of |K(p)|'s. ! ! 29 June 1993 ! 12 July 1993 ! 17 July 1994 ! 4 May 1996 ! 21 November 1996 ! 5 December 1996 ! 5 February 1997 ! 15 December 1998 ! 23 December 1998 ! 9 February 1999 ! 10 March 1999 ! 20 August 1999 ! 21 December 2000 ! 20 March 2001 ! 1 February 2002 ! 7 December 2002 ! integer :: nloops integer :: nloops1,nprops1,nverts1,cutmax1 integer :: nloops2,nprops2,nverts2,cutmax2 common /sizes/ nloops1,nprops1,nverts1,cutmax1, & nloops2,nprops2,nverts2,cutmax2 integer :: groupsize(maxgraphs,maxmaps) integer :: groupsizegraph(maxgraphs) integer :: groupsizetotal common /montecarlo/groupsize,groupsizegraph,groupsizetotal ! integer :: mapnumber,l,mu real(kind=dbl) :: p1(3),p2(3),ell1(3),absp1,absp2,absp3 real(kind=dbl) :: temp1,temp2,temp3,p1sq,p2sq,p3sq character(len=6) :: maptype integer :: qsign(0:size),q(0:size) real(kind=dbl) :: rho3,rho2to3d,rho2to3e,rho2to2t,rho2to2s,rho2to1a,rho2to1b real(kind=dbl) :: rhothree,rholoop ! IF (order.EQ.1) THEN nloops = nloops1 ELSE IF (order.EQ.2) THEN nloops = nloops2 ELSE write(nout,*)'Order must be 1 or 2.',order END IF ! IF (order.EQ.1) THEN ! ! We deal with the case of a Born graph first. ! density = 0.0d0 DO mapnumber = 1,nmaps ! DO l = 0,nloops q(l) = qs(mapnumber,l) qsign(l) = qsigns(mapnumber,l) END DO p1sq = 0.0d0 p2sq = 0.0d0 p3sq = 0.0d0 DO mu = 1,3 temp1 = qsign(1)*k(q(1),mu) temp2 = qsign(2)*k(q(2),mu) temp3 = - temp1 - temp2 p1sq = p1sq + temp1**2 p2sq = p2sq + temp2**2 p3sq = p3sq + temp3**2 END DO absp1 = sqrt(p1sq) absp2 = sqrt(p2sq) absp3 = sqrt(p3sq) rhothree = rho3(absp1,absp2,absp3) density = density & + rhothree*groupsize(graphnumber,mapnumber) ! END DO ! ! Alternative for IF (ORDER.EQ.1) THEN ! ELSE IF (order.EQ.2) THEN ! ! We tackle the case of an order alpha_s^2 graph. ! We construct the density as a sum. ! density = 0.0d0 DO mapnumber = 1,nmaps ! maptype = maptypes(mapnumber) DO l = 0,nloops q(l) = qs(mapnumber,l) qsign(l) = qsigns(mapnumber,l) END DO ! ! First, we need the kinematic variables for this map. ! p1sq = 0.0d0 p2sq = 0.0d0 p3sq = 0.0d0 DO mu = 1,3 ell1(mu) = qsign(1)*k(q(1),mu) temp1 = qsign(2)*k(q(2),mu) temp2 = qsign(3)*k(q(3),mu) temp3 = - temp1 - temp2 p1(mu) = temp1 p1sq = p1sq + temp1**2 p2(mu) = temp2 p2sq = p2sq + temp2**2 p3sq = p3sq + temp3**2 END DO absp1 = sqrt(p1sq) absp2 = sqrt(p2sq) absp3 = sqrt(p3sq) ! ! Now, there are two factors, one for the 'final state momenta' and ! one for the 'loop momentum.' ! rhothree = rho3(absp1,absp2,absp3) ! IF (maptype.EQ.'t2to3d') THEN rholoop = rho2to3d(p1,p2,ell1) ELSE IF (maptype.EQ.'t2to3e') THEN rholoop = rho2to3e(p1,p2,ell1) ELSE IF (maptype.EQ.'t2to2t') THEN rholoop = rho2to2t(p1,p2,ell1) ELSE IF (maptype.EQ.'t2to2s') THEN rholoop = rho2to2s(p1,p2,ell1) ELSE IF (maptype.EQ.'t2to1a') THEN rholoop = rho2to1a(p1,p2,ell1) ELSE IF (maptype.EQ.'t2to1b') THEN rholoop = rho2to1b(p1,p2,ell1) ELSE write(nout,*)'Bad maptype in density.' stop END IF ! density = density & + rhothree*rholoop*groupsize(graphnumber,mapnumber) ! ! Close DO MAPNUMBER = 1,NMAPS ! END DO ! ! Close for IF (ORDER.EQ.1) THEN ... ELSE IF (ORDER.EQ.2) THEN ! ELSE write(nout,*)'Order should have been 1 or 2 in density.' stop END IF ! RETURN END function density ! !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! Subroutines associated with NEWPOINT and DENSITY ! CHOOSEx and RHOx where x = 3, 2to2T, 2to2S, 2to3D, 2to3E, 2to1A,2to1B !2345678901234567890123456789012345678901234567890123456789012345678901234567890 ! subroutine choose3(p1,p2,p3,ok) ! use beowulf_parameters implicit none ! Out: real(kind=dbl) :: p1(3),p2(3),p3(3) logical :: ok ! ! Generates momenta P1(mu),P2(mu),P3(mu) for a three body final ! state with a distribution in momentum fractions x1,x2,x3 ! proportional to ! ! [max(1-x1,1-x2,1-x3)]^B/[(1-x1)*(1-x2)*(1-x3)]^B. ! ! 28 December 2000 ! 16 January 2001 ! real(kind=dbl) :: badnesslimit,cancellimit,thrustcut common /limits/ badnesslimit,cancellimit,thrustcut real(kind=dbl), parameter :: onethird = 0.3333333333333333333d0 real(kind=dbl), parameter :: twothirds = 0.6666666666666666667d0 real(kind=dbl), parameter :: pi = 3.141592653589793239d0 ! ! The parameter E3PAR should match between CHOOSE3 and RHO3. ! real(kind=dbl), parameter :: e3par = 1.5d0 ! ! The parameters A, B, and C need to match between CHOOSE3 and RHO3. ! CHOOSE3 uses A, while RHO3 uses B and C. The relation is ! B = 1 - 1/A and then C is the normalization factor and is