*********************************************************************** * * * FOURJPHACT * * A. Ballestrero * * Vers 0.2 (March 2000) * * * The parameters ppmax, irad, icom for the call to Pythia subroutine * py4jet are set at 0.d0, 0, 0 in a data statement in fxn function. * They can of course be changed if needed * *********************************************************************** implicit real*8 (a-b,d-h,o-z) implicit double complex (c) character*60 string DIMENSION region(18),y1(50,10),iv(32) COMMON/abqrea/rmz,gamz,s2w,rmz2,alfas,alfainv,pi,e_cm,ycut,xm(6) COMMON/abqpar/beta,rlim,s_col,smin,emcoupl,estrinf1,estrinf2 COMMON/abqopz/isr,icos,icut,ialfasrun,if1,if2,iff1,iid,iycut COMMON /abqran/ idum COMMON/abqsalv/iv,iy,idum2 COMMON/abqcou/fer,fel,zer,zel,f3l,f4l,f5l,f6l,z3l,z4l,z5l,z6l, & f3r,f4r,f5r,f6r,z3r,z4r,z5r,z6r COMMON/abqcom/czipr,ccz,czero,cuno,cim c with vegas COMMON/abresl/resl(10),standdevl(10) c for cuts COMMON/abqcuts/e_min(3:6),e_max(3:6),thbeam_min(3:6), & thbeam_max(3:6),thsep_min(6),thsep_max(6),rm_min(6),rm_max(6), & beamcut_min(3:6),beamcut_max(3:6),sepcut_min(6),sepcut_max(6), & rm_min2(6),rm_max2(6),pt_min(3:6),pt_max(3:6) c for distributions PARAMETER (ndismax=50,nintmax=50,nbinmax=500,nitmax=10) DIMENSION test(ndismax),devbin(ndismax,nbinmax) DIMENSION rnorm(ndismax,nbinmax,nitmax) DIMENSION nbintot(ndismax),nbintest(ndismax) COMMON/abdist/distr_estrinf(ndismax,nintmax+1), & bin_width(ndismax,nintmax), & distr_local(ndismax,nbinmax,nitmax), & dev_local(ndismax,nbinmax,nitmax), & tail_local(ndismax,nitmax) COMMON/abidis/ncallbin(ndismax,nbinmax,nitmax), & nbin_number(ndismax, & nintmax),nbin_sum(ndismax,nintmax),nsubint(ndismax),ndistr, & idistr,it1,init COMMON/abcdis/string(ndismax) c for event generation COMMON/abflat/rmaxfxn,rmaxfxn_1it,rmaxfxn_2it,scalemax COMMON/abifla/itmx,novermax,iflat,iseed,istorvegas,istormom, & iterm,ipyth,interf COMMON/abfla2/irepeat,nevent,nflevts c for oneshot common/abonei/ioneshot,icompprob common/abones/prob(20) common/abmass/rmd,rmu,rms,rmc,rmb,rmt EXTERNAL fxn,gammln DATA rmz/91.1867d0/,gamz/2.4974d0/, & rmd/0.33D0/,rmu/0.33D0/,rms/0.50D0/,rmc/1.50D0/, & rmb/4.8D0/,rmt/175D0/, & alfainv/128.07d0/,alfas/0.1230d0/, & s2w/0.231030912451068d0/,gf/1.16637d-05/, & pi/3.141592653589793238462643d0/, & rme/0.51099906d-3/ data czero/(0.d0,0.d0)/,cuno/(1.d0,0.d0)/,cim/(0.d0,1.d0)/ read*, e_cm read *, ioneshot if (ioneshot.eq.1) then read*, icompprob else read*, if1 ! flavour index 3 4 (1=d,2=u,3=s,4=c,5=b,6=t) read*, if2 !flavour index 5 6 (1=d,2=u,3=s,4=c,5=b,6=t,21=gg) endif read*, ialfasrun ! READ*,isr ! yes/no ISR READ*,iycut ! iycut=0 no ycut; =1 durham; =2 jade; =3 cambridge IF (iycut.gt.0) then READ*, ycut ENDIF READ*,icut ! yes/no cuts IF(icut.EQ.1)THEN READ*,e_min ! 4 energy lower cuts (GeV) READ*,e_max ! 4 energy upper cuts (GeV) READ*,rm_min ! 6 invariant mass lower limits (GeV) ! (34, 35, 36, 45, 46, 56) READ*,rm_max ! 6 invariant mass upper limits (GeV) READ*,pt_min ! 4 transverse momenta lower cuts (GeV) READ*,pt_max ! 4 transverse momenta upper cuts (GeV) READ*,icos ! angular cuts in deg (0) or cos (1) READ*,thbeam_min! 4 particle-beam angle lower (in degrees) cuts READ*,thbeam_max! 4 particle-beam angle upper (in degrees) cuts READ*,thsep_min ! 6 particle-particle angle lower (in degs) cuts READ*,thsep_max ! 6 particle-particle angle upper (in degs) cuts ENDIF ************************************************************* c this part has to be repeated in input file if ioneshot=1 c and icompprob =1 has been chosen ************************************************************* 100 READ*,idistr ! yes/no distributions IF(idistr.EQ.1)THEN READ*,ndistr ! number of distributions DO i=1,ndistr READ*,nsubint(i) ! number of sub-intervals with different ! binning READ*,(distr_estrinf(i,j),j=1,nsubint(i)+1) !lower limits of ! each subint (which coincide with the upper limit !of the previous one) + upper limit of the last subint READ*,(nbin_number(i,j),j=1,nsubint(i)) !number of bins in !each subint ENDDO !i ENDIF READ*,iflat ! yes/no flat event generation IF(iflat.EQ.1)THEN READ*,scalemax ! scale factor for the maximum READ*,istorvegas ! yes/no VEGAS data stored READ*,irepeat ! 0=normal;1=repeat 2nd only;2=rep. with ! nflevts fixed ! 1=repeat only second iteration using vegas ! data stored in previous run, ! 2=repeat, generating an exact number ! nflevts of flat events using vegas data ! stored in previous run IF(irepeat.eq.2)THEN READ*, nflevts ! number of events to be generated as ! described above END IF READ*,istormom ! yes/no momenta of flat events written ! in .dat files READ*,ipyth ! yes/no call to Pythia ENDIF READ*,acc ! integration accuracy READ*,iterm ! yes/no thermalization READ*,ncall_term ! thermalization calls per iteration READ*,itmx_term ! thermalization iterations READ*,ncall ! integration calls per iteration READ*,itmx ! integration iterations ! (2 for flat event generation) *controls if (ioneshot.eq.1.and.icompprob.eq.0) then open (unit=34,status='old',file='abprob.dat') read (34,*) prob endif if (ioneshot.ne.1)then if(if1.ne.1.and.if1.ne.2.and.if1.ne.3.and.if1.ne.4.and. & if1.ne.5) then print*, ' IPROC NOT ALLOWED ! wrong if1' stop endif if(if2.ne.1.and.if2.ne.2.and.if2.ne.3.and.if2.ne.4.and. & if2.ne.5.and.if2.ne.21) then print*, ' IPROC NOT ALLOWED ! wrong if2' stop endif endif IF(idistr.EQ.1)THEN IF(ndistr.GT.ndismax)THEN PRINT*,'ERROR:' PRINT*,'Number of distributions cannot exceed' PRINT*,'NDISMAX given as parameter' STOP ENDIF DO i=1,ndistr IF(nsubint(i).GT.nintmax)THEN PRINT*,'ERROR:' PRINT*,'Number of sub-intervals cannot exceed' PRINT*,'NINTMAX given as parameter' STOP ENDIF DO j=1,nsubint(i) nbintest(i)=nbintest(i)+nbin_number(i,j) ENDDO !j IF(nbintest(i).GT.nbinmax)THEN PRINT*,'ERROR:' PRINT*,'total number of bins cannot exceed' PRINT*,'NBINMAX given as parameter' STOP ENDIF ENDDO !i ENDIF IF(itmx_term.GT.nitmax)THEN PRINT*,'ERROR:' PRINT*,'Number of thermalization iterations cannot exceed' PRINT*,'NITMAX given as parameter' STOP ENDIF IF(itmx.GT.nitmax)THEN PRINT*,'ERROR:' PRINT*,'Number of integration iterations cannot exceed' PRINT*,'NITMAX given as parameter' STOP ENDIF IF(iflat.EQ.1.AND.itmx.NE.2)THEN PRINT*,'ERROR' PRINT*,'Flat events generation needs ITMX =2' stop ENDIF * controls end cosh if (ioneshot.eq.1.and.icompprob.eq.1) then ichann=20 cosh eventuali altre flag a 0, come iflat etc., che andranno salvate nel frattempo else ichann=1 endif do jjjjj=1,ichann c default value for if2 (=21 for gluon) and c iid=0 (no identical particles) iid=0 if (ioneshot.eq.1.and.icompprob.eq.1) then if (jjjjj.le.5) then if1=jjjjj if2=21 elseif (jjjjj.gt.5) then if (jjjjj.le.10) then if1=1 if2=jjjjj-5 elseif (jjjjj.le.14) then if1=2 if2=jjjjj-9 elseif (jjjjj.le.17) then if1=3 if2=jjjjj-12 elseif (jjjjj.le.19) then if1=4 if2=jjjjj-14 else if1=5 if2=5 endif endif endif coshend s_col=e_cm**2 emcoupl=(1.d0/pi/alfainv)**4 c initialize parameters if (ioneshot.eq.0.or.(ioneshot.eq.1.and.icompprob.eq.1)) then CALL initialize endif rmeisr=rme IF (icut.EQ.1) THEN DO i=3,6 IF(icos.EQ.1)THEN beamcut_min(i)=thbeam_min(i) beamcut_max(i)=thbeam_max(i) ELSE IF(icos.EQ.0)THEN beamcut_min(i)=cos(thbeam_min(i)*pi/180.d0) beamcut_max(i)=cos(thbeam_max(i)*pi/180.d0) ENDIF ENDDO !i DO i=1,6 IF(icos.EQ.1)THEN sepcut_min(i)=thsep_min(i) sepcut_max(i)=thsep_max(i) ELSE IF(icos.EQ.0)THEN sepcut_min(i)=cos(thsep_min(i)*pi/180.d0) sepcut_max(i)=cos(thsep_max(i)*pi/180.d0) ENDIF rm_min2(i)=rm_min(i)**2 rm_max2(i)=rm_max(i)**2 ENDDO !i ENDIF * Estremi di integrazione sulle masse m34 ed m56: IF(icut.EQ.0.OR.rm_min(1).EQ.0.d0)THEN estrinf_34=xm(3)+xm(4) ELSE IF(icut.EQ.1.AND.rm_min(1).NE.0.d0)THEN estrinf_34=rm_min(1) ENDIF IF(icut.EQ.0.OR.rm_min(6).EQ.0.d0)THEN estrinf_56=xm(5)+xm(6) ELSE IF(icut.EQ.1.AND.rm_min(6).NE.0.d0)THEN estrinf_56=rm_min(6) ENDIF IF(isr.EQ.1)THEN rl=log(s_col/rmeisr/rmeisr) alfa_me=1.d0/137.0359895d0 beta=2.d0*alfa_me*(rl-1.d0)/pi gamma=exp(gammln(1.d0+beta/2.d0)) ge=0.5772156649d0 rlim=exp(beta*(0.75d0-ge)/2.d0)*(beta/2.d0)/gamma ENDIF CALL printer * Routine Vegas parameter: IF(isr.EQ.0)THEN ndim=7 ELSE ndim=9 ENDIF * Integration limits: DO i=1,ndim region(i)=0.d0 region(ndim+i)=1.d0 ENDDO !i * estremi di integrazione sulle due masse invarianti estrinf1=estrinf_34 estrinf2=estrinf_56 smin=(xm(3)+xm(4)+xm(5)+xm(6))**2 * Thermalization: IF(iterm.EQ.1.AND.irepeat.EQ.0)THEN init=0 PRINT*,' ' PRINT*,'Thermalization' CALL vegas(region,ndim,fxn,init,ncall_term,itmx_term,nprn, & avgi1,sd1,rchi2a,acc,y1,it1,ndo1,si1,swgt1,schi1) ENDIF !thermalization * initialization of the variables used if IFLAT=1 IF(iflat.EQ.1)THEN * n. random generato con RAN2(ISEED) iseed=19753179 * numero di punti che eccedono il massimo fissato novermax=0 * n. di momenti generati nevent=0 * maximum function value at first and second vegas iteration rmaxfxn=0.d0 rmaxfxn_1it=0.d0 rmaxfxn_2it=0.d0 IF(istormom.EQ.1)THEN OPEN(unit=23,file='ABMOM.DAT',status='new') ENDIF ENDIF * inizializzazione di variabili usate se IDISTR=1 IF(idistr.EQ.1)THEN DO m=1,ndistr test(m)=0.d0 DO i=1,nbinmax DO j=1,nitmax distr_local(m,i,j)=0.d0 dev_local(m,i,j)=0.d0 ENDDO !j ENDDO !i ENDDO !m DO m=1,ndistr DO i=1,nitmax tail_local(m,i)=0.d0 ENDDO !i ENDDO !m DO m=1,ndistr DO i=1,nintmax nbin_sum(m,i)=0 ENDDO !i ENDDO !m DO m=1,ndistr DO i=1,nsubint(m) bin_width(m,i)=(distr_estrinf(m,i+1)-distr_estrinf(m,i))/ & nbin_number(m,i) ENDDO !i ENDDO !m ENDIF * integration loop IF(iterm.EQ.1)THEN init=1 ELSE init=0 ENDIF IF(iflat.EQ.0)THEN CALL vegas(region,ndim,fxn,init,ncall,itmx,nprn,avgi1, & sd1,rchi2a,acc,y1,it1,ndo1,si1,swgt1,schi1) ELSE IF(iflat.EQ.1)THEN acc=0.d0 DO nit=1,2 IF(nit.EQ.2)init=2 IF((nit.EQ.1.AND.irepeat.EQ.0).OR.(nit.EQ.2))THEN IF(irepeat.GE.1)THEN it1=2 OPEN(unit=24,file='abvegas.dat',status='old') READ(24,*)idum,idum2,iy DO i=1,32 READ(24,*)iv(i) ENDDO !i READ(24,*)avgi1,sd1,rchi2a DO i=1,50 DO j=1,10 READ(24,*)y1(i,j) ENDDO !j ENDDO !i READ(24,*)ndo1,si1,swgt1,schi1 READ(24,*)rmaxfxn READ(24,*)avgi_tot,sd_tot CLOSE(24) c eventually avoid vegas printing if(irepeat.GT.1) nprn=-1 c ENDIF CALL vegas(region,ndim,fxn,init,ncall,nit,nprn,avgi1, & sd1,rchi2a,acc,y1,it1,ndo1,si1,swgt1,schi1) IF(nit.EQ.1)THEN rmaxfxn_1it=rmaxfxn ENDIF IF(istorvegas.EQ.1.AND.nit.EQ.1)THEN OPEN(unit=24,file='abvegas.dat',status='new') WRITE(24,*)idum,idum2,iy DO i=1,32 WRITE(24,*)iv(i) ENDDO !i WRITE(24,*)avgi1,sd1,rchi2a DO i=1,50 DO j=1,10 WRITE(24,*)y1(i,j) ENDDO !j ENDDO !i WRITE(24,*)ndo1,si1,swgt1,schi1 ENDIF IF(istorvegas.EQ.1.AND.nit.EQ.2.and.irepeat.eq.0)THEN rmaxfxn=rmaxfxn_2it WRITE(24,*)rmaxfxn WRITE(24,*)avgi1,sd1 CLOSE(24) ENDIF ENDIF ENDDO !nit ENDIF *aggend if (irepeat.ne.2) then avgi_tot=avgi1 sd_tot=sd1 endif it=it1-1 PRINT*,' ' PRINT 151 151 FORMAT('------------------------------------------------------') PRINT*,' ' PRINT 201,avgi_tot,sd_tot 201 FORMAT(' Sigma = ',d13.7,' +/-',d10.3,' (pb)') cosh if (icompprob.eq.1) then prob(JJJJJ)=avgi_tot endif coshend *agg IF(iflat.EQ.1)THEN PRINT*,'Informations about flat events generation:' PRINT*,' ---------------------- ' IF(irepeat.EQ.0)THEN PRINT 203,rmaxfxn_1it 203 FORMAT(' Maximum after first VEGAS iteration = ',d9.3) ENDIF if (irepeat.ne.2) then PRINT 204,rmaxfxn_2it 204 FORMAT(' Maximum after second VEGAS iteration = ',d9.3) else PRINT*,'Maximum = ',rmaxfxn_2it endif PRINT 205,NEVENT 205 FORMAT(' Flat events number = ',i9) PRINT 206,novermax 206 FORMAT(' number of function values over maximum = ',i9) ENDIF cosh enddo !jjjjj if (ichann.eq.20) then c print*, ' ' c print*, 'cross sections of all channels' c print*, prob c print*, ' ' probtot=0.d0 do i=1,20 c print*, 'prob(i)', i, prob(i) probtot=probtot+prob(i) enddo do i=1,20 prob(i)=prob(i)/probtot enddo c print*, 'prob' c print*, prob c print*, ' ' do i=2,20 prob(i)=prob(i)+prob(i-1) enddo open (unit=35,status='new',file='abprob.dat') write (35,*) prob c print*, 'summed prob' c print*, prob c print*, ' ' c print*, ' ' c print*, 'total cross section of all channels' c print*,probtot icompprob=0 goto 100 endif coshend *aggend * Distributions: IF(idistr.EQ.1)THEN OPEN(unit=22,file='ABDIS.DAT',status='new') WRITE(22,*)'N_iteration=',it1-1 WRITE(22,*)'cross-section:',avgi_tot,'+/-',sd_tot,' (pb) ' WRITE(22,*) it1=it1-1 DO m=1,ndistr WRITE(22,*)string(m) DO n=1,nsubint(m) nbintot(m)=nbintot(m)+nbin_number(m,n) ENDDO !n DO i=1,nbintot(m) DO j=1,it1 IF(ncallbin(m,i,j).ge.2)THEN dev_local(m,i,j)=(ncallbin(m,i,j)*dev_local(m,i,j)- & distr_local(m,i,j)**2)/(ncallbin(m,i,j)-1) ENDIF ENDDO !j ENDDO !i DO i=1,nbinmax devbin(m,i)=0.d0 ENDDO !i DO i=1,nbintot(m) DO j=1,it1 IF(dev_local(m,i,j).ne.0.d0)THEN devbin(m,i)=devbin(m,i)+1.d0/dev_local(m,i,j) ENDIF ENDDO !j ENDDO !i DO i=1,nbintot(m) DO j=1,it1 rnorm(m,i,j)=dev_local(m,i,j)*devbin(m,i) ENDDO !j ENDDO !i DO i=1,nbintot(m) IF(devbin(m,i).gt.0.d0)THEN devbin(m,i)=1.d0/sqrt(devbin(m,i)) ENDIF ENDDO !i DO i=1,nbintot(m) IF(rnorm(m,i,it1).NE.0.d0)THEN test(m)=test(m)+distr_local(m,i,it1) distr_local(m,i,it1)=distr_local(m,i,it1)/rnorm(m,i,it1) ENDIF ENDDO !i DO i=1,nbintot(m) DO j=1,it1-1 IF(rnorm(m,i,j).NE.0.d0)THEN distr_local(m,i,it1)=distr_local(m,i,it1)+ & distr_local(m,i,j)/rnorm(m,i,j) ENDIF ENDDO !j ENDDO !i k=0 DO i=1,nsubint(m) delrm=bin_width(m,i) rmed=distr_estrinf(m,i)-delrm/2.d0 IF(i.eq.1)THEN k=1 ELSE k=k+nbin_number(m,i-1) ENDIF DO j=k,nbin_number(m,i)+k-1 rmed=rmed+delrm write(22,*)rmed,distr_local(m,j,it1)/bin_width(m,i), & devbin(m,j)/bin_width(m,i) ENDDO !j ENDDO !i WRITE(22,*)'test(m)=',test(m) WRITE(22,*)'tail_local(m,it1)=',tail_local(m,it1) WRITE(22,*)'test=',test(m)+tail_local(m,it1)-resl(it1) ENDDO !m CLOSE(22) ENDIF STOP END subroutine initialize implicit real*8 (a-b,d-h,o-z) implicit double complex (c) COMMON/abqrea/rmz,gamz,s2w,rmz2,alfas,alfainv,pi,e_cm,ycut,xm(6) COMMON/abqopz/isr,icos,icut,ialfasrun,if1,if2,iff1,iid,iycut COMMON/abqcou/fer,fel,zer,zel,f3l,f4l,f5l,f6l,z3l,z4l,z5l,z6l, & f3r,f4r,f5r,f6r,z3r,z4r,z5r,z6r COMMON/abqcom/czipr,ccz,czero,cuno,cim c for oneshot common/abmass/rmd,rmu,rms,rmc,rmb,rmt c initialize particles masses and identities if (if1.eq.1) then xm(3)=rmd elseif (if1.eq.2) then xm(3)=rmu elseif (if1.eq.3) then xm(3)=rms elseif (if1.eq.4) then xm(3)=rmc elseif (if1.eq.5) then xm(3)=rmb elseif (if1.eq.6) then xm(3)=rmt else print*, ' IF1 NOT ALLOWED ! ' stop endif if (if2.ne.21) then if (if2.eq.1) then xm(5)=rmd elseif (if2.eq.2) then xm(5)=rmu elseif (if2.eq.3) then xm(5)=rms elseif (if2.eq.4) then xm(5)=rmc elseif (if2.eq.5) then xm(5)=rmb elseif (if2.eq.6) then xm(5)=rmt else print*, ' IF1 NOT ALLOWED ! ' stop endif else xm(5)=0.d0 if (if1.eq.1.or.if1.eq.3.or.if1.eq.5) then iff1=2 else iff1=1 endif endif xm(4)=xm(3) xm(6)=xm(5) rmz2=rmz**2 czipr=(1.d0,0.d0) ccz=cim*gamz*rmz ** valori di accoppiamenti cr e cl sw=sqrt(s2w) rcw=sqrt(1.-s2w) rcotw=rcw/sw ** fotone-elettrone left e right fel=-1.d0 fer=-1.d0 ** Zeta-elettrone left e right zer=(+1.d0*s2w)/rcw/sw zel=(-.5d0+1.d0*s2w)/rcw/sw zvr=0.d0 zvl=.5d0/rcw/sw ** Zeta-quarkup left e right fqul=2.d0/3.d0 fqur=fqul zqur=(-2.d0/3.d0*s2w)/rcw/sw zqul=(.5d0-2.d0/3.d0*s2w)/rcw/sw ** Zeta-quarkdown left e right fqdl=-1.d0/3.d0 fqdr=fqdl zqdr=(1.d0/3.d0*s2w)/rcw/sw zqdl=(-.5d0+1.d0/3.d0*s2w)/rcw/sw ** W left (right=0) wcl=1.d0/sw/sqrt(2.d0) if(mod(if1,2).eq.0) then f3l=fqul f3r=fqur f4l=fqul f4r=fqur z3l=zqul z3r=zqur z4l=zqul z4r=zqur else f3l=fqdl f3r=fqdr f4l=fqdl f4r=fqdr z3l=zqdl z3r=zqdr z4l=zqdl z4r=zqdr endif if (if2.ne.21)then if(mod(if2,2).eq.0) then f5l=fqul f5r=fqur f6l=fqul f6r=fqur z5l=zqul z5r=zqur z6l=zqul z6r=zqur else f5l=fqdl f5r=fqdr f6l=fqdl f6r=fqdr z5l=zqdl z5r=zqdr z6l=zqdl z6r=zqdr endif endif if (if1.eq.if2) then iid=1 else iid=0 endif return end REAL*8 FUNCTION fxn(x,wgt) IMPLICIT REAL*8 (a-h,o-z) CHARACTER*60 string REAL*4 singlep(0:3) DIMENSION p1(0:3),p2(0:3),p3(0:3),p4(0:3),p5(0:3),p6(0:3) DIMENSION p(4,6),x(9) DIMENSION cc(3:6,3:6),y(3:6,3:6),rmod(3:6),cosbeam(3:6),pt(3:6) c for gg dimension qgg(4,6) c for cambridge dimension v(3:6,3:6),costheta(3:6,3:6),ee(3:6) COMMON/abqrea/rmz,gamz,s2w,rmz2,alfas,alfainv,pi,e_cm,ycut,xm(6) COMMON/abqpar/beta,rlim,s_col,smin,emcoupl,estrinf1,estrinf2 COMMON/abqopz/isr,icos,icut,ialfasrun,if1,if2,iff1,iid,iycut COMMON /abqran/ idum COMMON/abstat/ncall_eff c for cuts COMMON/abqcuts/e_min(3:6),e_max(3:6),thbeam_min(3:6), & thbeam_max(3:6),thsep_min(6),thsep_max(6),rm_min(6),rm_max(6), & beamcut_min(3:6),beamcut_max(3:6),sepcut_min(6),sepcut_max(6), & rm_min2(6),rm_max2(6),pt_min(3:6),pt_max(3:6) c for distributions PARAMETER (ndismax=50,nintmax=50,nbinmax=500,nitmax=10) DIMENSION distr_var(ndismax) COMMON/abdist/distr_estrinf(ndismax,nintmax+1), & bin_width(ndismax,nintmax), & distr_local(ndismax,nbinmax,nitmax), & dev_local(ndismax,nbinmax,nitmax), & tail_local(ndismax,nitmax) COMMON/abidis/ncallbin(ndismax,nbinmax,nitmax), & nbin_number(ndismax, & nintmax),nbin_sum(ndismax,nintmax),nsubint(ndismax),ndistr, & idistr,it1,init COMMON/abcdis/string(ndismax) c for event generation COMMON/abflat/rmaxfxn,rmaxfxn_1it,rmaxfxn_2it,scalemax COMMON/abifla/itmx,novermax,iflat,iseed,istorvegas,istormom, & iterm,ipyth,interf COMMON/abfla2/irepeat,nevent,nflevts c with Pythia C...HEPEVT commonblock. PARAMETER (NMXHEP=4000) COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP), &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP) c for oneshot common/abonei/ioneshot,icompprob common/abones/prob(20) DATA ncall_eff/0/ c no mass to the electrons in this phase space DATA rme/0.d0/ c parameters for Pythia py4jet call DATA ppmax/0.d0/,irad/0/,icom/0/ cosh if (ioneshot.eq.1.and.icompprob.eq.0) then xprob=ran2(idum) if (xprob.lt.prob(1)) then if1=1 if2=21 weigpr=prob(1) elseif (xprob.lt.prob(2)) then if1=2 if2=21 weigpr=prob(2)-prob(1) elseif (xprob.lt.prob(3)) then if1=3 if2=21 weigpr=prob(3)-prob(2) elseif (xprob.lt.prob(4)) then if1=4 if2=21 weigpr=prob(4)-prob(3) elseif (xprob.lt.prob(5)) then if1=5 if2=21 weigpr=prob(5)-prob(4) elseif (xprob.lt.prob(6)) then if1=1 if2=1 weigpr=prob(6)-prob(5) elseif (xprob.lt.prob(7)) then if1=1 if2=2 weigpr=prob(7)-prob(6) elseif (xprob.lt.prob(8)) then if1=1 if2=3 weigpr=prob(8)-prob(7) elseif (xprob.lt.prob(9)) then if1=1 if2=4 weigpr=prob(9)-prob(8) elseif (xprob.lt.prob(10)) then if1=1 if2=5 weigpr=prob(10)-prob(9) elseif (xprob.lt.prob(11)) then if1=2 if2=2 weigpr=prob(11)-prob(10) elseif (xprob.lt.prob(12)) then if1=2 if2=3 weigpr=prob(12)-prob(11) elseif (xprob.lt.prob(13)) then if1=2 if2=4 weigpr=prob(13)-prob(12) elseif (xprob.lt.prob(14)) then if1=2 if2=5 weigpr=prob(14)-prob(13) elseif (xprob.lt.prob(15)) then if1=3 if2=3 weigpr=prob(15)-prob(14) elseif (xprob.lt.prob(16)) then if1=3 if2=4 weigpr=prob(16)-prob(15) elseif (xprob.lt.prob(17)) then if1=3 if2=5 weigpr=prob(17)-prob(16) elseif (xprob.lt.prob(18)) then if1=4 if2=4 weigpr=prob(18)-prob(17) elseif (xprob.lt.prob(19)) then if1=4 if2=5 weigpr=prob(19)-prob(18) elseif (xprob.lt.1.d0) then if1=5 if2=5 weigpr=prob(20)-prob(19) else print*, 'UNEXPECTED PROBABILITY' stop endif call initialize * Estremi di integrazione sulle masse m34 ed m56: IF(icut.EQ.0.OR.rm_min(1).EQ.0.d0)THEN estrinf1=xm(3)+xm(4) ELSE IF(icut.EQ.1.AND.rm_min(1).NE.0.d0)THEN estrinf1=rm_min(1) ENDIF IF(icut.EQ.0.OR.rm_min(6).EQ.0.d0)THEN estrinf2=xm(5)+xm(6) ELSE IF(icut.EQ.1.AND.rm_min(6).NE.0.d0)THEN estrinf2=rm_min(6) ENDIF smin=(xm(3)+xm(4)+xm(5)+xm(6))**2 endif coshend el=e_cm/2.d0 IF(isr.EQ.1)THEN x1=1.d0-x(3)**(2.d0/beta) x2=1.d0-x(4)**(2.d0/beta) IF(x1.LE.0.d0.OR.x2.LE.0.d0.OR.x1.GT.1.d0.OR.x2.GT.1.d0)THEN fxn=0.d0 RETURN ENDIF IF(x1.EQ.1.d0)THEN dx1=rlim ELSE dx1=rlim+(1.d0-x1)**(1.d0-beta/2.d0)*(-beta*(1.d0+x1)/4.d0 & +(beta**2/32.d0)*(-4.d0*(1.d0+x1)*log(1.d0-x1)+3.d0* & (1.d0+x1)*log(x1)-4.d0*log(x1)/(1.d0-x1)-5.d0-x1)) ENDIF IF(x2.EQ.1.d0)THEN dx2=rlim ELSE dx2=rlim+(1.d0-x2)**(1.d0-beta/2.d0)*(-beta*(1.d0+x2)/4.d0 & +(beta**2/32.d0)*(-4.d0*(1.d0+x2)*log(1.d0-x2)+3.d0* & (1.d0+x2)*log(x2)-4.d0*log(x2)/(1.d0-x2)-5.d0-x2)) ENDIF str_fun=dx1*dx2*4.d0/beta/beta IF(str_fun.LE.0.d0)THEN fxn=0.d0 RETURN ENDIF IF(((x1*el)**2-rme**2).LT.0.d0.OR.((x2*el)**2-rme**2). & LT.0.d0)THEN fxn=0.d0 RETURN ENDIF s=2.d0*rme**2+2.d0*x1*x2*el**2+2.d0*sqrt((x1*el)**2-rme**2) & *sqrt((x2*el)**2-rme**2) IF(s.LT.smin)THEN fxn=0.d0 RETURN ENDIF bcm=(sqrt((x1*el)**2-rme**2)-sqrt((x2*el)**2-rme**2))/ & ((x1+x2)*el) IF((1.d0-bcm**2).LT.0.d0)THEN fxn=0.d0 RETURN ENDIF gcm=1.d0/sqrt(1.d0-bcm**2) * per ottenere l'elettrone sul mass-shell fisso l'energia e ricavo * l'impulso: p1(0)=x1*el p1(1)=0.d0 p1(2)=0.d0 p1(3)=max(sqrt((x1*el)**2-rme**2),0.d0) p2(0)=x2*el p2(1)=0.d0 p2(2)=0.d0 p2(3)=-max(sqrt((x2*el)**2-rme**2),0.d0) ELSE s=s_col bcm=0.d0 gcm=1.d0 p1(0)=el p1(1)=0.d0 p1(2)=0.d0 p1(3)=max(sqrt(el**2-rme**2),0.d0) p2(0)=el p2(1)=0.d0 p2(2)=0.d0 p2(3)=-max(sqrt(el**2-rme**2),0.d0) ENDIF IF(x(1).EQ.1.d0.OR.x(1).EQ.0.d0.OR.x(2).EQ. & 1.d0.OR.x(2).EQ.0.d0)THEN fxn=0.d0 RETURN ENDIF IF(icut.EQ.0.OR.rm_max(1).GE.e_cm)THEN estrsup1=sqrt(s)-estrinf2 IF(estrsup1.LE.estrinf1)THEN fxn=0.d0 RETURN ENDIF ELSE IF(icut.EQ.1.AND.rm_max(1).LT.e_cm)THEN IF(rm_max(1).GT.(sqrt(s)-estrinf2))THEN estrsup1=sqrt(s)-estrinf2 IF(estrsup1.LE.estrinf1)THEN fxn=0.d0 RETURN ENDIF ELSE estrsup1=rm_max(1) ENDIF ENDIF * (34) and (56) moments in the lab frame: rj34=estrsup1-estrinf1 xm(1)=x(1)*rj34+estrinf1 IF(icut.EQ.0.OR.rm_max(6).GE.e_cm)THEN estrsup2=sqrt(s)-xm(1) ELSE IF(icut.EQ.1.AND.rm_max(6).LT.e_cm)THEN IF(rm_max(6).GT.(sqrt(s)-xm(1)))THEN estrsup2=sqrt(s)-xm(1) ELSE estrsup2=rm_max(6) ENDIF ENDIF rj56=estrsup2-estrinf2 xm(2)=rj56*x(2)+estrinf2 app=((s-xm(1)**2-xm(2)**2)**2-4.d0*xm(1)**2*xm(2)**2)/(4.d0*s) IF(app.LE.0.d0)THEN fxn=0.d0 RETURN ENDIF p34_cm=sqrt(app) IF(isr.EQ.0)THEN c34=x(3)*2.d0-1.d0 c3=x(4)*2.d0-1.d0 ph3=x(5)*2.d0*pi c5=x(6)*2.d0-1.d0 ph5=x(7)*2.d0*pi ELSE c34=x(5)*2.d0-1.d0 c3=x(6)*2.d0-1.d0 ph3=x(7)*2.d0*pi c5=x(8)*2.d0-1.d0 ph5=x(9)*2.d0*pi ENDIF app=1.d0-c34**2 IF(app.LT.0.d0)THEN fxn=0.d0 RETURN ENDIF s34=sqrt(app) p(4,1)=sqrt(p34_cm**2+xm(1)**2) p(1,1)=p34_cm*s34 p(3,1)=p34_cm*c34 p(4,2)=sqrt(p34_cm**2+xm(2)**2) p(1,2)=-p(1,1) p(3,2)=-p(3,1) * (3) and (4) moments in the (34) frame: app=(xm(1)**2-xm(3)**2-xm(4)**2)**2-4.d0*(xm(3)*xm(4))**2 IF(app.LE.0.d0)THEN fxn=0.d0 RETURN ENDIF p3_34=sqrt(app)/(2.d0*xm(1)) app=1.d0-c3**2 IF(app.LT.0.d0)THEN fxn=0.d0 RETURN ENDIF s3=sqrt(app) p(4,3)=sqrt(p3_34**2+xm(3)**2) p(1,3)=p3_34*s3*cos(ph3) p(2,3)=p3_34*s3*sin(ph3) p(3,3)=p3_34*c3 p(4,4)=sqrt(p3_34**2+xm(4)**2) p(1,4)=-p(1,3) p(2,4)=-p(2,3) p(3,4)=-p(3,3) p(2,1)=0.d0 p(2,2)=0.d0 * boost in the lab frame: p3(0)=(p(4,1)*p(4,3)+p(1,1)*p(1,3)+p(2,1)*p(2,3) & +p(3,1)*p(3,3))/xm(1) trasf=(p(4,3)+p3(0))/(xm(1)+p(4,1)) p3(1)=p(1,3)+trasf*p(1,1) p3(2)=p(2,3)+trasf*p(2,1) p3(3)=p(3,3)+trasf*p(3,1) p4(0)=(p(4,1)*p(4,4)+p(1,1)*p(1,4)+p(2,1)*p(2,4) & +p(3,1)*p(3,4))/xm(1) trasf=(p(4,4)+p4(0))/(xm(1)+p(4,1)) p4(1)=p(1,4)+trasf*p(1,1) p4(2)=p(2,4)+trasf*p(2,1) p4(3)=p(3,4)+trasf*p(3,1) * (5) and (6) moments in the (56) frame: app=(xm(2)**2-xm(5)**2-xm(6)**2)**2-4.d0*(xm(5)*xm(6))**2 IF(app.LE.0.d0)THEN fxn=0.d0 RETURN ENDIF p5_56=sqrt(app)/(2.d0*xm(2)) app=1.d0-c5**2 IF(app.LT.0.d0)THEN fxn=0.d0 RETURN ENDIF s5=sqrt(app) p(4,5)=sqrt(p5_56**2+xm(5)**2) p(1,5)=p5_56*s5*cos(ph5) p(2,5)=p5_56*s5*sin(ph5) p(3,5)=p5_56*c5 p(4,6)=sqrt(p5_56**2+xm(6)**2) p(1,6)=-p(1,5) p(2,6)=-p(2,5) p(3,6)=-p(3,5) * (5) and (6) moments in the lab frame: p5(0)=(p(4,2)*p(4,5)+p(1,2)*p(1,5)+p(2,2)*p(2,5) & +p(3,2)*p(3,5))/xm(2) trasf=(p(4,5)+p5(0))/(xm(2)+p(4,2)) p5(1)=p(1,5)+trasf*p(1,2) p5(2)=p(2,5)+trasf*p(2,2) p5(3)=p(3,5)+trasf*p(3,2) p6(0)=(p(4,2)*p(4,6)+p(1,2)*p(1,6)+p(2,2)*p(2,6) & +p(3,2)*p(3,6))/xm(2) trasf=(p(4,6)+p6(0))/(xm(2)+p(4,2)) p6(1)=p(1,6)+trasf*p(1,2) p6(2)=p(2,6)+trasf*p(2,2) p6(3)=p(3,6)+trasf*p(3,2) * boost nel sistema del collider IF(isr.EQ.1)THEN p3(0)=gcm*(p3(0)+bcm*p3(3)) p3(1)=p3(1) p3(2)=p3(2) p3(3)=p3(3)/gcm+bcm*p3(0) p4(0)=gcm*(p4(0)+bcm*p4(3)) p4(1)=p4(1) p4(2)=p4(2) p4(3)=p4(3)/gcm+bcm*p4(0) p5(0)=gcm*(p5(0)+bcm*p5(3)) p5(1)=p5(1) p5(2)=p5(2) p5(3)=p5(3)/gcm+bcm*p5(0) p6(0)=gcm*(p6(0)+bcm*p6(3)) p6(1)=p6(1) p6(2)=p6(2) p6(3)=p6(3)/gcm+bcm*p6(0) ENDIF * rotazioni per dare un angolo fi (solo per distribuzioni,generazione eventi..) fi=ran2(iseed)*2.d0*pi aux=p3(1)*cos(fi)-p3(2)*sin(fi) p3(2)=p3(1)*sin(fi)+p3(2)*cos(fi) p3(1)=aux aux=p4(1)*cos(fi)-p4(2)*sin(fi) p4(2)=p4(1)*sin(fi)+p4(2)*cos(fi) p4(1)=aux aux=p5(1)*cos(fi)-p5(2)*sin(fi) p5(2)=p5(1)*sin(fi)+p5(2)*cos(fi) p5(1)=aux aux=p6(1)*cos(fi)-p6(2)*sin(fi) p6(2)=p6(1)*sin(fi)+p6(2)*cos(fi) p6(1)=aux * Cuts: IF (icut.EQ.1) THEN rmod(3)=sqrt(p3(1)**2+p3(2)**2+p3(3)**2) rmod(4)=sqrt(p4(1)**2+p4(2)**2+p4(3)**2) rmod(5)=sqrt(p5(1)**2+p5(2)**2+p5(3)**2) rmod(6)=sqrt(p6(1)**2+p6(2)**2+p6(3)**2) * Energy cuts: IF(p3(0).LE.e_min(3).OR.p4(0).LE.e_min(4).OR.p5(0).LE. & e_min(5).OR.p6(0).LE.e_min(6).OR.p3(0).GE.e_max(3).OR. & p4(0).GE.e_max(4).OR.p5(0).GE.e_max(5).OR.p6(0).GE. & e_max(6))THEN fxn=0.d0 RETURN ENDIF * beam-angle cuts: IF(beamcut_min(3).NE.1.d0.OR.beamcut_min(4).NE.1.d0.OR. & beamcut_min(5).NE.1.d0.OR.beamcut_min(6).NE.1.d0.OR. & beamcut_max(3).NE.-1.d0.OR.beamcut_max(4).NE.-1.d0.OR. & beamcut_max(5).NE.-1.d0.OR.beamcut_max(6).NE.-1.d0)THEN cosbeam(3)=p3(3)/rmod(3) cosbeam(4)=p4(3)/rmod(4) cosbeam(5)=p5(3)/rmod(5) cosbeam(6)=p6(3)/rmod(6) IF(cosbeam(3).GE.beamcut_min(3).OR.cosbeam(4).GE. & beamcut_min(4).OR.cosbeam(5).GE.beamcut_min(5).OR. & cosbeam(6).GE.beamcut_min(6).OR.cosbeam(3).LE. & beamcut_max(3).OR.cosbeam(4).LE.beamcut_max(4).OR. & cosbeam(5).LE.beamcut_max(5).OR.cosbeam(6).LE. & beamcut_max(6))THEN fxn=0.d0 RETURN ENDIF ENDIF * angular separation cuts: IF(sepcut_min(1).NE.1.d0.OR.sepcut_min(2).NE.1.d0.OR. & sepcut_min(3).NE.1.d0.OR.sepcut_min(4).NE.1.d0.OR. & sepcut_min(5).NE.1.d0.OR.sepcut_min(6).NE.1.d0.OR. & sepcut_max(1).NE.-1.d0.OR.sepcut_max(2).NE.-1.d0.OR. & sepcut_max(3).NE.-1.d0.OR.sepcut_max(4).NE.-1.d0.OR. & sepcut_max(5).NE.-1.d0.OR.sepcut_max(6).NE.-1.d0)THEN cc(3,4)=(p3(1)*p4(1)+p3(2)*p4(2)+p3(3)*p4(3))/rmod(3)/rmod(4) cc(3,5)=(p3(1)*p5(1)+p3(2)*p5(2)+p3(3)*p5(3))/rmod(3)/rmod(5) cc(3,6)=(p3(1)*p6(1)+p3(2)*p6(2)+p3(3)*p6(3))/rmod(3)/rmod(6) cc(4,5)=(p4(1)*p5(1)+p4(2)*p5(2)+p4(3)*p5(3))/rmod(4)/rmod(5) cc(4,6)=(p4(1)*p6(1)+p4(2)*p6(2)+p4(3)*p6(3))/rmod(4)/rmod(6) cc(5,6)=(p5(1)*p6(1)+p5(2)*p6(2)+p5(3)*p6(3))/rmod(5)/rmod(6) IF(cc(3,4).GE.sepcut_min(1).OR.cc(3,5).GE.sepcut_min(2).OR. & cc(3,6).GE.sepcut_min(3).OR.cc(4,5).GE.sepcut_min(4).OR. & cc(4,6).GE.sepcut_min(5).OR.cc(5,6).GE.sepcut_min(6).OR. & cc(3,4).LE.sepcut_max(1).OR.cc(3,5).LE.sepcut_max(2).OR. & cc(3,6).LE.sepcut_max(3).OR.cc(4,5).LE.sepcut_max(4).OR. & cc(4,6).LE.sepcut_max(5).OR.cc(5,6).LE.sepcut_max(6))THEN fxn=0.d0 RETURN ENDIF ENDIF * invariant mass cuts: IF(rm_min(1).GT.0.d0.OR.rm_min(2).GT.0.d0.OR.rm_min(3).GT. & 0.d0.OR.rm_min(4).GT.0.d0.OR.rm_min(5).GT.0.d0.OR. & rm_min(6).GT.0.d0.OR.rm_max(1).LT.e_cm.OR.rm_max(2).LT. & e_cm.OR.rm_max(3).LT.e_cm.OR.rm_max(4).LT.e_cm.OR. & rm_max(5).LT.e_cm.OR.rm_max(6).LT.e_cm)THEN y(3,4)=xm(3)**2+xm(4)**2+2.d0*(p3(0)*p4(0)- & p3(1)*p4(1)-p3(2)*p4(2)-p3(3)*p4(3)) y(3,5)=xm(3)**2+xm(5)**2+2.d0*(p3(0)*p5(0)- & p3(1)*p5(1)-p3(2)*p5(2)-p3(3)*p5(3)) y(3,6)=xm(3)**2+xm(6)**2+2.d0*(p3(0)*p6(0)- & p3(1)*p6(1)-p3(2)*p6(2)-p3(3)*p6(3)) y(4,5)=xm(4)**2+xm(5)**2+2.d0*(p4(0)*p5(0)- & p4(1)*p5(1)-p4(2)*p5(2)-p4(3)*p5(3)) y(4,6)=xm(4)**2+xm(6)**2+2.d0*(p4(0)*p6(0)- & p4(1)*p6(1)-p4(2)*p6(2)-p4(3)*p6(3)) y(5,6)=xm(5)**2+xm(6)**2+2.d0*(p5(0)*p6(0)- & p5(1)*p6(1)-p5(2)*p6(2)-p5(3)*p6(3)) IF(y(3,4).LE.rm_min2(1).OR.y(3,5).LE.rm_min2(2).OR.y(3,6). & LE.rm_min2(3).OR.y(4,5).LE.rm_min2(4).OR.y(4,6).LE. & rm_min2(5).OR.y(5,6).LE.rm_min2(6).OR.y(3,4).GE. & rm_max2(1).OR.y(3,5).GE.rm_max2(2).OR.y(3,6).GE. & rm_max2(3).OR.y(4,5).GE.rm_max2(4).OR.y(4,6).GE. & rm_max2(5).OR.y(5,6).GE.rm_max2(6))THEN fxn=0.d0 RETURN ENDIF ENDIF * Pt cuts: IF(pt_min(3).GT.0.d0.OR.pt_min(4).GT.0.d0.OR.pt_min(5).GT. & 0.d0.OR.pt_min(6).GT.0.d0.OR.pt_max(3).LT.e_cm.OR.pt_max(4). & LT.e_cm.OR.pt_max(5).LT.e_cm.OR.pt_max(6).LT.e_cm)THEN pt(3)=sqrt(p3(1)**2+p3(2)**2) pt(4)=sqrt(p4(1)**2+p4(2)**2) pt(5)=sqrt(p5(1)**2+p5(2)**2) pt(6)=sqrt(p6(1)**2+p6(2)**2) IF(pt(3).LE.pt_min(3).OR.pt(4).LE.pt_min(4).OR.pt(5).LE. & pt_min(5).OR.pt(6).LE.pt_min(6).OR.pt(3).GE.pt_max(3).OR. & pt(4).GE.pt_max(4).OR.pt(5).GE.pt_max(5).OR.pt(6).GE. & pt_max(6))THEN fxn=0.d0 RETURN ENDIF ENDIF ENDIF !icut * here include additional cuts if (iycut.GT.0) then p3modulo=sqrt(p3(1)**2+p3(2)**2+p3(3)**2) p4modulo=sqrt(p4(1)**2+p4(2)**2+p4(3)**2) p5modulo=sqrt(p5(1)**2+p5(2)**2+p5(3)**2) p6modulo=sqrt(p6(1)**2+p6(2)**2+p6(3)**2) rprodspt34=p3(1)*p4(1)+p3(2)*p4(2)+p3(3)*p4(3) rprodspt35=p3(1)*p5(1)+p3(2)*p5(2)+p3(3)*p5(3) rprodspt36=p3(1)*p6(1)+p3(2)*p6(2)+p3(3)*p6(3) rprodspt45=p5(1)*p4(1)+p5(2)*p4(2)+p5(3)*p4(3) rprodspt46=p6(1)*p4(1)+p6(2)*p4(2)+p6(3)*p4(3) rprodspt56=p5(1)*p6(1)+p5(2)*p6(2)+p5(3)*p6(3) rcoseno34=rprodspt34/p3modulo/p4modulo rcoseno35=rprodspt35/p3modulo/p5modulo rcoseno36=rprodspt36/p3modulo/p6modulo rcoseno45=rprodspt45/p4modulo/p5modulo rcoseno46=rprodspt46/p4modulo/p6modulo rcoseno56=rprodspt56/p5modulo/p6modulo if (iycut.eq.1) then ydur34=2.d0*min(p3(0)**2,p4(0)**2)*(1-rcoseno34)/s ydur35=2.d0*min(p3(0)**2,p5(0)**2)*(1-rcoseno35)/s ydur36=2.d0*min(p3(0)**2,p6(0)**2)*(1-rcoseno36)/s ydur45=2.d0*min(p4(0)**2,p5(0)**2)*(1-rcoseno45)/s ydur46=2.d0*min(p4(0)**2,p6(0)**2)*(1-rcoseno46)/s ydur56=2.d0*min(p5(0)**2,p6(0)**2)*(1-rcoseno56)/s if ( (ydur34.lt.ycut).or.(ydur35.lt.ycut).or. & (ydur36.lt.ycut).or.(ydur45.lt.ycut).or. & (ydur46.lt.ycut).or.(ydur56.lt.ycut)) then fxn=0.d0 return endif endif if (iycut.eq.2) then yjade34=2.d0*p3(0)*p4(0)*(1-rcoseno34)/s yjade35=2.d0*p3(0)*p5(0)*(1-rcoseno35)/s yjade36=2.d0*p3(0)*p6(0)*(1-rcoseno36)/s yjade45=2.d0*p4(0)*p5(0)*(1-rcoseno45)/s yjade46=2.d0*p4(0)*p6(0)*(1-rcoseno46)/s yjade56=2.d0*p5(0)*p6(0)*(1-rcoseno56)/s if ( (yjade34.lt.ycut).or.(yjade35.lt.ycut).or. & (yjade36.lt.ycut).or.(yjade45.lt.ycut).or. & (yjade46.lt.ycut).or.(yjade56.lt.ycut)) then fxn=0.d0 return endif endif if (iycut.eq.3) then ee(3)=p3(0) ee(4)=p4(0) ee(5)=p5(0) ee(6)=p6(0) costheta(3,4)=rcoseno34 costheta(3,5)=rcoseno35 costheta(3,6)=rcoseno36 costheta(4,5)=rcoseno45 costheta(4,6)=rcoseno46 costheta(5,6)=rcoseno56 do i=3,5 do j=i+1,6 v(i,j)=2.d0*(1.d0-costheta(i,j)) enddo enddo rmax=max(v(3,4),v(3,5),v(3,6),v(4,5),v(4,6),v(5,6)) do k=1,3 vmin=min(v(3,4),v(3,5),v(3,6),v(4,5),v(4,6),v(5,6)) iimin=0 jjmin=0 do i=3,5 do j=i+1,6 if (vmin.eq.v(i,j)) then iimin=i jjmin=j endif enddo enddo if (iimin.eq.0.or.jjmin.eq.0) then print*, 'iimin and/or jjmin not found' stop else if (ee(iimin).lt.ee(jjmin)) then iimark=iimin else iimark=jjmin endif ytest=ee(iimark)**2*vmin if (ytest.lt.ycut) then fxn=0.d0 return else do i=3,iimark-1 v(i,iimark)=rmax+1.d0 enddo do i=iimark+1,6 v(iimark,i)=rmax+1.d0 enddo endif endif enddo endif endif IF(idistr.EQ.1)THEN * here include parton distributions * include 'abdis.dis' * here commented is an example of distributions c c string(1)=' Distribution: 34 invariant mass' c distr_var(1)=sqrt((p3(0)+p4(0))**2-(p3(1)+p4(1))**2- c & (p3(2)+p4(2))**2-(p3(3)+p4(3))**2) c string(2)='Distribution: particle 5 energy ' c distr_var(2)=p5(0) c * end abdis ENDIF fxn=4.d0*p3_34*p5_56*p34_cm*(pi)**3*rj34*rj56/sqrt(s) cosh if (ioneshot.eq.1.and.icompprob.eq.0) then fxn=fxn/weigpr endif coshend * ISR introduction: IF(isr.EQ.1)THEN fxn=fxn*str_fun ENDIF if (if2.ne.21) then fxn=fxn*ee_4q(p1,p2,p3,p4,p5,p6) fxn=fxn*emcoupl/s/2.d0 c alfas in matrix element and previous result for alfa_em**4 fxn=fxn*alfainv**2 else qgg(3,2)=p1(1) qgg(1,2)=p1(2) qgg(2,2)=p1(3) qgg(4,2)=p1(0) qgg(3,1)=p2(1) qgg(1,1)=p2(2) qgg(2,1)=p2(3) qgg(4,1)=p2(0) qgg(3,3)=p5(1) qgg(1,3)=p5(2) qgg(2,3)=p5(3) qgg(4,3)=p5(0) qgg(3,4)=p6(1) qgg(1,4)=p6(2) qgg(2,4)=p6(3) qgg(4,4)=p6(0) qgg(3,5)=p3(1) qgg(1,5)=p3(2) qgg(2,5)=p3(3) qgg(4,5)=p3(0) qgg(3,6)=p4(1) qgg(1,6)=p4(2) qgg(2,6)=p4(3) qgg(4,6)=p4(0) call s_2e2q2g(qgg,xm(3),iff1,res) fxn=fxn*res if (ialfasrun.eq.1) then scale=(p5(0)+p6(0))**2-(p5(1)+p6(1))**2-(p5(2)+p6(2)) & **2-(p5(3)+p6(3))**2 if (scale.le.0.d0) then fxn=0.d0 else scale=sqrt(scale) endif alfasgg=rals(.233d0,scale,5) else alfasgg=alfas endif fxn=fxn*emcoupl/s/2.d0 c previous result for alfa_em**4 fxn=fxn*(alfainv* alfasgg)**2 endif fxn=fxn*0.38937966d+9 !sigma(pb) *agg IF(iflat.EQ.1)THEN IF((iterm.EQ.0).OR.(iterm.EQ.1.AND.init.GE.1))THEN iter=it1 ENDIF IF(iter.EQ.1.AND. & ((iterm.EQ.0).OR.(iterm.EQ.1.AND.init.GE.1)))THEN rmaxfxn=max(rmaxfxn,fxn*wgt) ENDIF IF(iter.EQ.2.AND. & ((iterm.EQ.0).OR.(iterm.EQ.1.AND.init.GE.1)))THEN rmaxfxn_2it=max(rmaxfxn_2it,fxn*wgt) ENDIF IF(iter.EQ.2.AND. & ((iterm.EQ.0).OR.(iterm.EQ.1.AND.init.GE.1)))THEN IF(fxn*wgt.GT.rmaxfxn*scalemax)THEN novermax=novermax+1 ENDIF ry=scalemax*rmaxfxn*ran2(iseed) IF(ry.LE.fxn*wgt)THEN * numero di eventi generati NEVENT=NEVENT+1 NEVHEP=NEVENT IF(ipyth.EQ.1)THEN NHEP=4 DO i=1,4 ISTHEP(i)=1 ENDDO !i if (if2.eq.21) then DO i=0,3 IF(i.EQ.0)j=4 IF(i.NE.0)j=i phep(j,1)=p3(i) phep(j,2)=p4(i) phep(j,3)=p5(i) phep(j,4)=p6(i) ENDDO !i phep(5,1)=xm(3) phep(5,2)=xm(4) phep(5,3)=xm(5) phep(5,4)=xm(6) idhep(1)=if1 idhep(2)=-if1 idhep(3)=if2 idhep(4)=if2 else rmm34=(p3(0)+p4(0))**2-(p3(1)+p4(1))**2- 1 (p3(2)+p4(2))**2-(p3(3)+p4(3))**2 rmm56=(p5(0)+p6(0))**2-(p5(1)+p6(1))**2- 1 (p5(2)+p6(2))**2-(p5(3)+p6(3))**2 if (iid.eq.0)then if (rmm34.gt.rmm56) then DO i=0,3 IF(i.EQ.0)j=4 IF(i.NE.0)j=i phep(j,1)=p3(i) phep(j,2)=p4(i) phep(j,3)=p5(i) phep(j,4)=p6(i) ENDDO !i phep(5,1)=xm(3) phep(5,2)=xm(4) phep(5,3)=xm(5) phep(5,4)=xm(6) idhep(1)=if1 idhep(2)=-if1 idhep(3)=if2 idhep(4)=-if2 else DO i=0,3 IF(i.EQ.0)j=4 IF(i.NE.0)j=i phep(j,1)=p5(i) phep(j,2)=p6(i) phep(j,3)=p3(i) phep(j,4)=p4(i) ENDDO !i phep(5,1)=xm(5) phep(5,2)=xm(6) phep(5,3)=xm(3) phep(5,4)=xm(4) idhep(1)=if2 idhep(2)=-if2 idhep(3)=if1 idhep(4)=-if1 endif else idhep(1)=if1 idhep(2)=-if1 idhep(3)=if1 idhep(4)=-if1 rmm36=(p3(0)+p6(0))**2-(p3(1)+p6(1))**2- 1 (p3(2)+p6(2))**2-(p3(3)+p6(3))**2 rmm54=(p5(0)+p4(0))**2-(p5(1)+p4(1))**2- 1 (p5(2)+p4(2))**2-(p5(3)+p4(3))**2 rmmmin=min(rmm34,rmm56,rmm36,rmm54) if (rmmmin.eq.rmm34) then DO i=0,3 IF(i.EQ.0)j=4 IF(i.NE.0)j=i phep(j,1)=p5(i) phep(j,2)=p6(i) phep(j,3)=p3(i) phep(j,4)=p4(i) ENDDO phep(5,1)=xm(5) phep(5,2)=xm(6) phep(5,3)=xm(3) phep(5,4)=xm(4) elseif (rmmmin.eq.rmm56) then DO i=0,3 IF(i.EQ.0)j=4 IF(i.NE.0)j=i phep(j,1)=p3(i) phep(j,2)=p4(i) phep(j,3)=p5(i) phep(j,4)=p6(i) ENDDO !i phep(5,1)=xm(3) phep(5,2)=xm(4) phep(5,3)=xm(5) phep(5,4)=xm(6) elseif (rmmmin.eq.rmm36) then DO i=0,3 IF(i.EQ.0)j=4 IF(i.NE.0)j=i phep(j,1)=p5(i) phep(j,2)=p4(i) phep(j,3)=p3(i) phep(j,4)=p6(i) ENDDO !i phep(5,1)=xm(5) phep(5,2)=xm(4) phep(5,3)=xm(3) phep(5,4)=xm(6) elseif (rmmmin.eq.rmm54) then DO i=0,3 IF(i.EQ.0)j=4 IF(i.NE.0)j=i phep(j,1)=p3(i) phep(j,2)=p6(i) phep(j,3)=p5(i) phep(j,4)=p4(i) ENDDO !i phep(5,1)=xm(3) phep(5,2)=xm(6) phep(5,3)=xm(5) phep(5,4)=xm(4) else print*, 'ERROR IN FINDING MINIMAL INVARIAN MASS' stop endif endif endif * call to Pythia call py4jet(ppmax,irad,icom) ctest if (nevent.lt.10) then call pylist(1) endif ctestend * reinitialize HEPEVT do i=1,4000 do j=1,5 phep(j,i)=0.d0 end do idhep(i)=0 end do ENDIF IF(istormom.EQ.1)THEN nout=23 DO m=0,3 singlep(m)=p3(m) ENDDO !m WRITE(nout,*)singlep DO m=0,3 singlep(m)=p4(m) ENDDO !m WRITE(nout,*)singlep DO m=0,3 singlep(m)=p5(m) ENDDO !m WRITE(nout,*)singlep DO m=0,3 singlep(m)=p6(m) ENDDO !m WRITE(nout,*)singlep ENDIF ENDIF ENDIF ENDIF *aggend * numero di chiamate effettive dopo avere superato i tagli ncall_eff=ncall_eff+1 * distribuzioni pesate * distribuzioni pesate IF(idistr.EQ.1.AND.((iterm.EQ.1.AND.init.NE.0).OR.iterm.EQ.0))THEN delhist=fxn*wgt devhist=delhist**2 DO i=1,ndistr DO j=1,nsubint(i)-1 IF(j.GE.2)THEN nbin_sum(i,j)=nbin_sum(i,j-1)+nbin_number(i,j-1) ENDIF IF(nbin_number(i,j).NE.0.AND.distr_var(i).LT.distr_estrinf & (i,j+1).AND.distr_var(i).GE.distr_estrinf(i,j))THEN IF(j.EQ.1)THEN nbin=int((distr_var(i)-distr_estrinf(i,j))/ & bin_width(i,j))+1 ELSE nbin=int((distr_var(i)-distr_estrinf(i,j))/ & bin_width(i,j))+(1+nbin_sum(i,j)) ENDIF distr_local(i,nbin,it1)=distr_local(i,nbin,it1)+delhist dev_local(i,nbin,it1)=dev_local(i,nbin,it1)+devhist ncallbin(i,nbin,it1)=ncallbin(i,nbin,it1)+1 ENDIF ENDDO !j j=nsubint(i) nbin_sum(i,j)=nbin_sum(i,j-1)+nbin_number(i,j-1) IF(nbin_number(i,j).NE.0.AND.distr_var(i).LE.distr_estrinf & (i,j+1).AND.distr_var(i).GE.distr_estrinf(i,j))THEN nbin=int((distr_var(i)-distr_estrinf(i,j))/bin_width(i,j))+ & (1+nbin_sum(i,j)) IF(nbin.EQ.(nbin_sum(i,j)+nbin_number(i,j)+1))THEN nbin=nbin-1 ENDIF distr_local(i,nbin,it1)=distr_local(i,nbin,it1)+delhist dev_local(i,nbin,it1)=dev_local(i,nbin,it1)+devhist ncallbin(i,nbin,it1)=ncallbin(i,nbin,it1)+1 ENDIF IF(distr_var(i).GT.distr_estrinf(i,nsubint(i)+1).OR. & distr_var(i).LT.distr_estrinf(i,1))THEN tail_local(i,it1)=tail_local(i,it1)+delhist ENDIF ENDDO !i ENDIF RETURN END REAL*8 FUNCTION gammln(xx) IMPLICIT NONE REAL*8 xx INTEGER j REAL*8 ser,stp,tmp,xa,y,cof(6) SAVE cof,stp DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, *-.5395239384953d-5,2.5066282746310005d0/ xa=xx y=xa tmp=xa+5.5d0 tmp=(xa+0.5d0)*log(tmp)-tmp ser=1.000000000190015d0 DO 11 j=1,6 y=y+1.d0 ser=ser+cof(j)/y 11 CONTINUE gammln=tmp+log(stp*ser/xa) RETURN END C (C) Copr. 1986-92 Numerical Recipes Software #>,1')5c). SUBROUTINE printer IMPLICIT REAL*8 (a-h,o-z) character*3 cha3_1, cha3_2 COMMON/abqrea/rmz,gamz,s2w,rmz2,alfas,alfainv,pi,e_cm,ycut,xm(6) COMMON/abqopz/isr,icos,icut,ialfasrun,if1,if2,iff1,iid,iycut c for cuts COMMON/abqcuts/e_min(3:6),e_max(3:6),thbeam_min(3:6), & thbeam_max(3:6),thsep_min(6),thsep_max(6),rm_min(6),rm_max(6), & beamcut_min(3:6),beamcut_max(3:6),sepcut_min(6),sepcut_max(6), & rm_min2(6),rm_max2(6),pt_min(3:6),pt_max(3:6) c for event generation COMMON/abflat/rmaxfxn,rmaxfxn_1it,rmaxfxn_2it,scalemax COMMON/abifla/itmx,novermax,iflat,iseed,istorvegas,istormom, & iterm,ipyth,interf COMMON/abfla2/irepeat,nevent,nflevts c for oneshot common/abonei/ioneshot,icompprob print*, ' ' print*, ' ' if (ioneshot.eq.1.and.icompprob.eq.0) then print*, ' e+e- -> all 4 partons' else if (if1.eq.1) cha3_1='dd~' if (if1.eq.2) cha3_1='uu~' if (if1.eq.3) cha3_1='ss~' if (if1.eq.4) cha3_1='cc~' if (if1.eq.5) cha3_1='bb~' if (if1.eq.6) cha3_1='tt~' if (if2.eq.21) then print*, ' e+e- -> '//cha3_1//'gg' else if (if2.eq.1) cha3_2='dd~' if (if2.eq.2) cha3_2='uu~' if (if2.eq.3) cha3_2='ss~' if (if2.eq.4) cha3_2='cc~' if (if2.eq.5) cha3_2='bb~' if (if2.eq.6) cha3_2='tt~' print*, ' e+e- -> '//cha3_1//cha3_2 endif endif print*, ' ' print*, ' **************' print*, ' ' PRINT*,'INPUT' PRINT 101,e_cm 101 FORMAT(' cm energy = ',d13.7,' GeV') PRINT*,' ' PRINT*,'DATA' PRINT 471,s2w 471 FORMAT(' s2w = ',d13.7) PRINT 481,alfainv 481 FORMAT(' 1/alfa_em = ',d13.7) if (ialfasrun.ne.1) then PRINT 485,alfas 485 FORMAT(' alfas = ',d13.7) endif PRINT*,' ' PRINT*,'OPTIONS' IF (isr.EQ.1) THEN PRINT*,'Born + QED' ELSE PRINT*,'Born only' ENDIF if (ialfasrun.eq.1) then print*,'alfas running' endif PRINT*,' ' if (iycut.eq.0) then print *, 'NO YCUT ' elseif (iycut.eq.1) then print 497, ycut 497 Format ('DURHAM ycut = ', d15.5) elseif (iycut.eq.2) then print 498, ycut 498 Format ('JADE ycut = ', d15.5) elseif (iycut.eq.3) then print 499, ycut 499 Format ('CAMBRIDGE ycut = ', d15.5) else print*, 'error: IYCUT NOT VALID' stop endif PRINT*,' ' IF (icut.EQ.1) THEN PRINT*,'Cuts :' PRINT*,'-----------------' PRINT 501,e_min(3),e_min(4),e_min(5),e_min(6) 501 FORMAT(' ENERGY_MIN(3,4,5,6) =( ',f6.2, & ',',f6.2,',',f6.2,',',f6.2,' ) GeV') PRINT 511,e_max(3),e_max(4),e_max(5),e_max(6) 511 FORMAT(' ENERGY_MAX(3,4,5,6) =( ',f7.2, & ',',f7.2,',',f7.2,',',f7.2,' ) GeV') PRINT 801,rm_min(1),rm_min(2),rm_min(3),rm_min(4),rm_min(5), & rm_min(6) 801 FORMAT(' MASS_MIN(34,35,36,45,46,56) =( ',f6.2, & ',',f6.2,',',f6.2,',',f6.2,',',f6.2,',',f6.2,' ) GeV') PRINT 851,rm_max(1),rm_max(2),rm_max(3),rm_max(4),rm_max(5), & rm_max(6) 851 FORMAT(' MASS_MAX(34,35,36,45,46,56) =( ',f7.2, & ',',f7.2,',',f7.2,',',f7.2,',',f7.2,',',f7.2,' ) GeV') PRINT 901,pt_min(3),pt_min(4),pt_min(5),pt_min(6) 901 FORMAT(' PT_MIN(3,4,5,6) =( ',f6.2, & ',',f6.2,',',f6.2,',',f6.2,' ) GeV') PRINT 911,pt_max(3),pt_max(4),pt_max(5),pt_max(6) 911 FORMAT(' PT_MAX(3,4,5,6) =( ',f7.2, & ',',f7.2,',',f7.2,',',f7.2,' ) GeV') IF (icos.EQ.0)THEN PRINT 701,thbeam_min(3),thbeam_min(4),thbeam_min(5), & thbeam_min(6) 701 FORMAT(' THBEAM_MIN(3,4,5,6) =( ',f6.2, & ',',f6.2,',',f6.2,',',f6.2,' ) deg') PRINT 702,thbeam_max(3),thbeam_max(4),thbeam_max(5), & thbeam_max(6) 702 FORMAT(' THBEAM_MAX(3,4,5,6) =( ',f6.2, & ',',f6.2,',',f6.2,',',f6.2,' ) deg') ELSE IF(icos.EQ.1)THEN PRINT 711,thbeam_min(3),thbeam_min(4),thbeam_min(5), & thbeam_min(6) 711 FORMAT(' COSBEAM_MAX(3,4,5,6) =( ',f6.2, & ',',f6.2,',',f6.2,',',f6.2,' )') PRINT 712,thbeam_max(3),thbeam_max(4),thbeam_max(5), & thbeam_max(6) 712 FORMAT(' COSBEAM_MIN(3,4,5,6) =( ',f6.2, & ',',f6.2,',',f6.2,',',f6.2,' )') ENDIF IF(icos.EQ.0)THEN PRINT 601,thsep_min(1),thsep_min(2),thsep_min(3),thsep_min(4), & thsep_min(5),thsep_min(6) 601 FORMAT(' THSEP_MIN(34,35,36,45,46,56) =( ',f6.2, & ',',f6.2,',',f6.2,',',f6.2,',',f6.2,',',f6.2,' ) deg') PRINT 602,thsep_max(1),thsep_max(2),thsep_max(3),thsep_max(4), & thsep_max(5),thsep_max(6) 602 FORMAT(' THSEP_MAX(34,35,36,45,46,56) =( ',f6.2, & ',',f6.2,',',f6.2,',',f6.2,',',f6.2,',',f6.2,' ) deg') ELSE IF(icos.EQ.1)THEN PRINT 611,thsep_min(1),thsep_min(2),thsep_min(3),thsep_min(4), & thsep_min(5),thsep_min(6) 611 FORMAT(' COSSEP_MAX(34,35,36,45,46,56) =( ',f6.2, & ',',f6.2,',',f6.2,',',f6.2,',',f6.2,',',f6.2,' )') PRINT 612,thsep_max(1),thsep_max(2),thsep_max(3),thsep_max(4), & thsep_max(5),thsep_max(6) 612 FORMAT(' COSSEP_MIN(34,35,36,45,46,56) =( ',f6.2, & ',',f6.2,',',f6.2,',',f6.2,',',f6.2,',',f6.2,' )') ENDIF ENDIF PRINT 951 951 FORMAT('------------------------------------------------------') PRINT*,' ' *agg IF(iflat.EQ.1)THEN PRINT*,'Flat events generation' IF(istorvegas.EQ.1)THEN cdaf ricordarsi di usare abvegas.dat per dati PRINT*,'VEGAS data stored in ABVEGAS.DAT' ENDIF IF(irepeat.EQ.1)THEN PRINT*,'second VEGAS iteration repeated' ENDIF PRINT 952,scalemax 952 FORMAT(' Maximum scale factor = ',d9.3) IF(istormom.EQ.1)THEN cdaf ricordarsi di usare abmom.dat per dati PRINT*,'Flat events stored in ABMOM.DAT' ENDIF IF(ipyth.EQ.1)THEN PRINT*,'with calls to Pythia' ENDIF ENDIF *aggend RETURN END SUBROUTINE vegas(region,ndim,fxn,init,ncall,itmx,nprn,tgral,sd, * chi2a,acc,xi,it,ndo,si,swgt,schi) Cc ho aggiunto acc INTEGER init,itmx,ncall,ndim,nprn,ndmx,mxdim,ncall_eff REAL*8 tgral,chi2a,sd,acc,region(2*ndim),fxn,alph,tiny PARAMETER (alph=1.5,ndmx=50,mxdim=10,tiny=1.e-30) C PARAMETER (ALPH=1.5,NDMX=100,MXDIM=10,TINY=1.e-30) EXTERNAL fxn CU USES fxn,ran2,rebin INTEGER i,idum,it,j,k,mds,nd,ndo,ng,npg,ia(mxdim),kg(mxdim) REAL*8 calls,dv2g,dxg,f,f2,f2b,fb,rc,ti,tsi,wgt,xjac,xn,xnd,xo, *d(ndmx,mxdim),di(ndmx,mxdim),dt(mxdim),dx(mxdim),r(ndmx),x(mxdim), *xi(ndmx,mxdim),xin(ndmx),ran2 REAL*8 schi,si,swgt,resl,standdevl COMMON/abresl/resl(10),standdevl(10) COMMON /abqran/ idum COMMON/abstat/ncall_eff COMMON/abfla2/irepeat,nevent,nflevts DATA mds/1/ SAVE IF(init.LE.0)THEN mds=1 ndo=1 C** it=1 C** DO 11 j=1,ndim xi(1,j)=1. 11 CONTINUE ENDIF IF (init.LE.1)THEN si=0. swgt=0. schi=0. C** it=1 C** ENDIF IF (init.LE.2)THEN nd=ndmx ng=1 IF(mds.NE.0)THEN ng=(ncall/2.+0.25)**(1./ndim) mds=1 IF((2*ng-ndmx).ge.0)then mds=-1 npg=ng/ndmx+1 nd=ng/npg ng=npg*nd ENDIF ENDIF k=ng**ndim npg=max(ncall/k,2) calls=npg*k dxg=1./ng dv2g=(calls*dxg**ndim)**2/npg/npg/(npg-1.) xnd=nd dxg=dxg*xnd xjac=1./calls DO 12 j=1,ndim dx(j)=region(j+ndim)-region(j) xjac=xjac*dx(j) 12 CONTINUE IF(nd.NE.ndo)THEN DO 13 i=1,nd r(i)=1. 13 CONTINUE DO 14 j=1,ndim CALL rebin(ndo/xnd,nd,r,xin,xi(1,j)) 14 CONTINUE ndo=nd ENDIF IF(nprn.GE.0) WRITE(*,200) ndim,calls,it,itmx ENDIF DO 28 it=it,itmx C** IF(it.GE.2.AND.acc*abs(tgral).ge.sd) RETURN C** ti=0. tsi=0. DO 16 j=1,ndim kg(j)=1 DO 15 i=1,nd d(i,j)=0. di(i,j)=0. 15 CONTINUE 16 CONTINUE 10 CONTINUE c aggiunta per far scegliere a caso cella in cui valutare funzione c per generazione piatta con numero determinato di eventi nflevts IF (irepeat.EQ.2) THEN DO WHILE (nevent.LT.nflevts) wgt=xjac DO j=1,ndim kg(j)=ran2(idum)*ng+1 xn=(kg(j)-ran2(idum))*dxg+1. ia(j)=max(min(int(xn),ndmx),1) IF(ia(j).GT.1)THEN xo=xi(ia(j),j)-xi(ia(j)-1,j) rc=xi(ia(j)-1,j)+(xn-ia(j))*xo ELSE xo=xi(ia(j),j) rc=(xn-ia(j))*xo ENDIF x(j)=region(j)+rc*dx(j) wgt=wgt*xo*xnd END DO f=wgt*fxn(x,wgt) END DO RETURN ENDIF fb=0. f2b=0. DO 19 k=1,npg wgt=xjac DO 17 j=1,ndim xn=(kg(j)-ran2(idum))*dxg+1. ia(j)=max(min(int(xn),ndmx),1) IF(ia(j).GT.1)THEN xo=xi(ia(j),j)-xi(ia(j)-1,j) rc=xi(ia(j)-1,j)+(xn-ia(j))*xo ELSE xo=xi(ia(j),j) rc=(xn-ia(j))*xo ENDIF x(j)=region(j)+rc*dx(j) wgt=wgt*xo*xnd 17 CONTINUE f=wgt*fxn(x,wgt) f2=f*f fb=fb+f f2b=f2b+f2 DO 18 j=1,ndim di(ia(j),j)=di(ia(j),j)+f IF(mds.GE.0) d(ia(j),j)=d(ia(j),j)+f2 18 CONTINUE 19 CONTINUE f2b=sqrt(f2b*npg) f2b=(f2b-fb)*(f2b+fb) IF (f2b.LE.0.) f2b=tiny ti=ti+fb tsi=tsi+f2b IF(mds.LT.0)THEN DO 21 j=1,ndim d(ia(j),j)=d(ia(j),j)+f2b 21 CONTINUE ENDIF DO 22 k=ndim,1,-1 kg(k)=mod(kg(k),ng)+1 IF(kg(k).NE.1) GOTO 10 22 CONTINUE tsi=tsi*dv2g wgt=1./tsi si=si+dble(wgt)*dble(ti) schi=schi+dble(wgt)*dble(ti)**2 swgt=swgt+dble(wgt) tgral=si/swgt chi2a=max((schi-si*tgral)/(it-.99d0),0.d0) sd=sqrt(1./swgt) tsi=sqrt(tsi) IF(nprn.GE.0)THEN C** aggiunta di ncall_eff e sua inizializzazione dopo ogni iterazione C** ho modificato anche il FORMAT 201 WRITE(*,201) it,ncall_eff,it,ti,tsi,tgral,sd,chi2a ncall_eff=0 resl(it)=ti standdevl(it)=tsi IF(nprn.NE.0)THEN DO 23 j=1,ndim WRITE(*,202) j,(xi(i,j),di(i,j),i=1+nprn/2,nd,nprn) 23 CONTINUE ENDIF ENDIF c** era qui l'aggiunta DO 25 j=1,ndim xo=d(1,j) xn=d(2,j) d(1,j)=(xo+xn)/2. dt(j)=d(1,j) DO 24 i=2,nd-1 rc=xo+xn xo=xn xn=d(i+1,j) d(i,j)=(rc+xn)/3. dt(j)=dt(j)+d(i,j) 24 CONTINUE d(nd,j)=(xo+xn)/2. dt(j)=dt(j)+d(nd,j) 25 CONTINUE DO 27 j=1,ndim rc=0. DO 26 i=1,nd IF(d(i,j).lt.tiny) d(i,j)=tiny r(i)=((1.-d(i,j)/dt(j))/(log(dt(j))-log(d(i,j))))**alph rc=rc+r(i) 26 CONTINUE CALL rebin(rc/xnd,nd,r,xin,xi(1,j)) 27 CONTINUE 28 CONTINUE RETURN 200 FORMAT(/' input parameters for vegas: ndim=',i3,' ncall=', *f12.0/28x,' it=',i5,' itmx=',i5) 201 FORMAT(/' iteration no.',i3,':',12x,'effective ncall=',i11/ *' iteration no.',i3,': ','integral =',g14.7,'+/- ',g9.2/ *' all iterations: integral =',g14.7,'+/- ',g9.3,' chi**2/it' *'n =',g9.2) 202 FORMAT(/' data for axis ',i2/' X delta i ', *' x delta i ',' x delta i ',/(1x, *f7.5,1x,g11.4,5x,f7.5,1x,g11.4,5x,f7.5,1x,g11.4)) END C (C) Copr. 1986-92 Numerical Recipes Software #>,1')5c). FUNCTION ran2(idum) INTEGER idum,im1,im2,imm1,ia1,ia2,iq1,iq2,ir1,ir2,ntab,ndiv REAL*8 ran2,am,eps,rnmx PARAMETER (im1=2147483563,im2=2147483399,am=1./im1,imm1=im1-1, *ia1=40014,ia2=40692,iq1=53668,iq2=52774,ir1=12211,ir2=3791, *ntab=32,ndiv=1+imm1/ntab,eps=1.2e-7,rnmx=1.-eps) INTEGER idum2,j,k,iv(ntab),iy COMMON/abqsalv/iv,iy,idum2 DATA idum2/123456789/, iv/ntab*0/, iy/0/ IF (idum.LE.0) THEN idum=max(-idum,1) idum2=idum DO 11 j=ntab+8,1,-1 k=idum/iq1 idum=ia1*(idum-k*iq1)-k*ir1 IF (idum.LT.0) idum=idum+im1 IF (j.LE.ntab) iv(j)=idum 11 CONTINUE iy=iv(1) ENDIF k=idum/iq1 idum=ia1*(idum-k*iq1)-k*ir1 IF (idum.LT.0) idum=idum+im1 k=idum2/iq2 idum2=ia2*(idum2-k*iq2)-k*ir2 IF (idum2.LT.0) idum2=idum2+im2 j=1+iy/ndiv iy=iv(j)-idum2 iv(j)=idum IF(iy.LT.1)iy=iy+imm1 ran2=min(am*iy,rnmx) RETURN END C (C) Copr. 1986-92 Numerical Recipes Software #>,1')5c). SUBROUTINE rebin(rc,nd,r,xin,xi) INTEGER nd REAL*8 rc,r(*),xi(*),xin(*) INTEGER i,k REAL*8 dr,xn,xo k=0 xn=0. dr=0. DO 11 i=1,nd-1 1 IF(rc.GT.dr)THEN k=k+1 dr=dr+r(k) xo=xn xn=xi(k) GOTO 1 ENDIF dr=dr-rc xin(i)=xn-(xn-xo)*dr/r(k) 11 CONTINUE DO 12 i=1,nd-1 xi(i)=xin(i) 12 CONTINUE xi(nd)=1. RETURN END C (C) Copr. 1986-92 Numerical Recipes Software #>,1')5c). real*8 function ee_4q(p1,p2,p3,p4,p5,p6) * impulsi p1(0:3),p2(0:3),p3(0:3),p4(0:3),p5(0:3),p6(0:3) * 1 e+, 2 e-, * implicit real*8 (a-b,d-h,o-z) implicit double complex (c) dimension p1(0:3),p2(0:3),p3(0:3),p4(0:3),p5(0:3),p6(0:3) dimension pau(0:3),p12(0:3),p14(0:3),p16(0:3),p23(0:3),p25(0:3), & p34(0:3),p36(0:3),p45(0:3),p56(0:3), & p123(0:3),p125(0:3),p134(0:3),p136(0:3),p145(0:3),p156(0:3), & p235(0:3),p345(0:3),p356(0:3) dimension cn_qcd(2,2,2,2,2,2),cn_qcd_sc(2,2,2,2,2,2) structure/res/ double complex id(0:1) end structure record/res/cn7(2,2,2,2,2,2),cn8(2,2,2,2,2,2), & cn9(2,2,2,2,2,2),cn10(2,2,2,2,2,2), structure/polcom/ double complex e(0:3),ek0,v end structure record/polcom/c12f(2,2),c12z(2,2), & c34f(2,2),c56f(2,2) structure/tu/ double complex a(2,2),b(2,2),c(2,2),d(2,2) end structure record/tu/l3_12fz(2,2),l3_56f(2,2), & l5_12fz(2,2),l5_34f(2,2), & v12f(0:3),v12z(0:3),v56f(0:3),v34f(0:3),tapptt(2,2,2,2) structure/td/ double complex a(2,2),b(2,2),c(2,2),d(2,2) end structure record/td/ r4_12fz(2,2),r4_56f(2,2),r6_12fz(2,2),r6_34f(2,2) COMMON/abqrea/rmz,gamz,s2w,rmz2,alfas,alfainv,pi,e_cm,ycut,xm(6) COMMON/abqopz/isr,icos,icut,ialfasrun,if1,if2,iff1,iid,iycut COMMON/abqcou/fer,fel,zer,zel,f3l,f4l,f5l,f6l,z3l,z4l,z5l,z6l, & f3r,f4r,f5r,f6r,z3r,z4r,z5r,z6r COMMON/abqcom/czipr,ccz,czero,cuno,cim DATA rmme/0.d0/ * pk0 -- p=p1 p1k0=p1(0)-p1(1) * pk0 -- p=p2 p2k0=p2(0)-p2(1) * pk0 -- p=p3 p3k0=p3(0)-p3(1) * pk0 -- p=p4 p4k0=p4(0)-p4(1) * pk0 -- p=p5 p5k0=p5(0)-p5(1) * pk0 -- p=p6 p6k0=p6(0)-p6(1) * Impulsi dei propagatori do m=0,3 p12(m)=-p1(m)-p2(m) p14(m)=-p1(m)+p4(m) p16(m)=-p1(m)+p6(m) p23(m)=-p2(m)+p3(m) p34(m)=p3(m)+p4(m) p25(m)=-p2(m)+p5(m) p36(m)=p3(m)+p6(m) p45(m)=p5(m)+p4(m) p56(m)=p5(m)+p6(m) p123(m)=p12(m)+p3(m) p125(m)=-p1(m)+p25(m) p134(m)=-p1(m)+p34(m) p136(m)=p16(m)+p3(m) p145(m)=-p1(m)+p45(m) p156(m)=p5(m)+p16(m) p235(m)=p25(m)+p3(m) p345(m)=p34(m)+p5(m) p356(m)=p36(m)+p5(m) end do if (ialfasrun.eq.1) then rm34=p34(0)*p34(0)-p34(1)*p34(1)-p34(2)*p34(2)-p34(3)*p34(3) rm34=sqrt(max(rm34,0.d0)) rm56=p56(0)*p56(0)-p56(1)*p56(1)-p56(2)*p56(2)-p56(3)*p56(3) rm56=sqrt(max(rm56,0.d0)) if (iid.eq.1) then rm45=p45(0)*p45(0)-p45(1)*p45(1)-p45(2)*p45(2)-p45(3)*p45(3) rm45=sqrt(max(rm45,0.d0)) rm36=p36(0)*p36(0)-p36(1)*p36(1)-p36(2)*p36(2)-p36(3)*p36(3) rm36=sqrt(max(rm36,0.d0)) endif endif * pk0 -- p=p123 p123k0=p123(0)-p123(1) * p.q -- p.q=p123q,p=p123,q=p123 p123q=p123(0)*p123(0)-p123(1)*p123(1)-p123(2)*p123(2)-p123 & (3)*p123(3) * pk0 -- p=p125 p125k0=p125(0)-p125(1) * p.q -- p.q=p125q,p=p125,q=p125 p125q=p125(0)*p125(0)-p125(1)*p125(1)-p125(2)*p125(2)-p125 & (3)*p125(3) * pk0 -- p=p134 p134k0=p134(0)-p134(1) * p.q -- p.q=p134q,p=p134,q=p134 p134q=p134(0)*p134(0)-p134(1)*p134(1)-p134(2)*p134(2)-p134 & (3)*p134(3) * pk0 -- p=p145 p145k0=p145(0)-p145(1) * p.q -- p.q=p145q,p=p145,q=p145 p145q=p145(0)*p145(0)-p145(1)*p145(1)-p145(2)*p145(2)-p145 & (3)*p145(3) * pk0 -- p=p156 p156k0=p156(0)-p156(1) * p.q -- p.q=p156q,p=p156,q=p156 p156q=p156(0)*p156(0)-p156(1)*p156(1)-p156(2)*p156(2)-p156 & (3)*p156(3) * pk0 -- p=p235 p235k0=p235(0)-p235(1) * p.q -- p.q=p235q,p=p235,q=p235 p235q=p235(0)*p235(0)-p235(1)*p235(1)-p235(2)*p235(2)-p235 & (3)*p235(3) * pk0 -- p=p136 p136k0=p136(0)-p136(1) * p.q -- p.q=p136q,p=p136,q=p136 p136q=p136(0)*p136(0)-p136(1)*p136(1)-p136(2)*p136(2)-p136 & (3)*p136(3) * pk0 -- p=p345 p345k0=p345(0)-p345(1) * p.q -- p.q=p345q,p=p345,q=p345 p345q=p345(0)*p345(0)-p345(1)*p345(1)-p345(2)*p345(2)-p345 & (3)*p345(3) * pk0 -- p=p356 p356k0=p356(0)-p356(1) * p.q -- p.q=p356q,p=p356,q=p356 p356q=p356(0)*p356(0)-p356(1)*p356(1)-p356(2)*p356(2)-p356 & (3)*p356(3) * p.q -- p.q=p1p2,p=p1,q=p2 p1p2=p1(0)*p2(0)-p1(1)*p2(1)-p1(2)*p2(2)-p1(3)*p2(3) s=2.d0*p1p2+2.d0*rmme**2 ss=sqrt(s) E=ss/2.d0 * vettori di polarizzazione e+e- quqd=p1p2 * il - qui e in cdw e' perche' i propagatori bosonici portano un segno - * rispetto a quelli fermionici cdz=-1.d0/(s*czipr-rmz2+ccz) df=-1.d0/s fac1=(df*fer) fac2=(df*fel) cfac1z=(cdz*zer) cfac2z=(cdz*zel) * T -- qu=p1,qd=p2,v=0,a=v12f(0).a,b=v12f(0).b,c=v12f(0).c,d=v12f(0).d,cr=f * ac1,cl=fac2,nsum=0 eps_0=-p1(2)*p2(3)+p2(2)*p1(3) ceps_0=eps_0*cim ceps_1=p1(3)*cim ceps_2=p2(3)*cim auxa=-quqd+p1k0*p2(0)+p2k0*p1(0) v12f(0).a(1,1)=fac1*(auxa+ceps_0) v12f(0).a(2,2)=fac2*(auxa-ceps_0) v12f(0).b(1,2)=-fac2*(p2(2)+ceps_2) v12f(0).b(2,1)=fac1*(p2(2)-ceps_2) v12f(0).c(1,2)=fac1*(p1(2)+ceps_1) v12f(0).c(2,1)=fac2*(-p1(2)+ceps_1) v12f(0).d(1,1)=fac2 v12f(0).d(2,2)=fac1 * T -- qu=p1,qd=p2,v=1,a=v12f(1).a,b=v12f(1).b,c=v12f(1).c,d=v12f(1).d,cr=f * ac1,cl=fac2,nsum=0 auxa=-quqd+p1k0*p2(1)+p2k0*p1(1) v12f(1).a(1,1)=fac1*(auxa+ceps_0) v12f(1).a(2,2)=fac2*(auxa-ceps_0) v12f(1).b(1,2)=-fac2*(p2(2)+ceps_2) v12f(1).b(2,1)=fac1*(p2(2)-ceps_2) v12f(1).c(1,2)=fac1*(p1(2)+ceps_1) v12f(1).c(2,1)=fac2*(-p1(2)+ceps_1) v12f(1).d(1,1)=fac2 v12f(1).d(2,2)=fac1 * T -- qu=p1,qd=p2,v=2,a=v12f(2).a,b=v12f(2).b,c=v12f(2).c,d=v12f(2).d,cr=f * ac1,cl=fac2,nsum=0 eps_0=-p1k0*p2(3)+p2k0*p1(3) ceps_0=eps_0*cim auxa=p1k0*p2(2)+p2k0*p1(2) v12f(2).a(1,1)=fac1*(auxa+ceps_0) v12f(2).a(2,2)=fac2*(auxa-ceps_0) v12f(2).b(1,2)=-fac2*p2k0 v12f(2).b(2,1)=fac1*p2k0 v12f(2).c(1,2)=fac1*p1k0 v12f(2).c(2,1)=-fac2*p1k0 * T -- qu=p1,qd=p2,v=3,a=v12f(3).a,b=v12f(3).b,c=v12f(3).c,d=v12f(3).d,cr=f * ac1,cl=fac2,nsum=0 eps_0=p1k0*p2(2)-p2k0*p1(2) ceps_0=eps_0*cim ceps_1=p1k0*cim ceps_2=p2k0*cim auxa=+p1k0*p2(3)+p2k0*p1(3) v12f(3).a(1,1)=fac1*(auxa+ceps_0) v12f(3).a(2,2)=fac2*(auxa-ceps_0) v12f(3).b(1,2)=-fac2*ceps_2 v12f(3).b(2,1)=-fac1*ceps_2 v12f(3).c(1,2)=fac1*ceps_1 v12f(3).c(2,1)=fac2*ceps_1 do m=0,3 * mline -- res=c12f(&).e(m),abcd=v12f(m).,m1=(-rmme),m2=rmme,den=0 do iut=1,2 do jut=1,2 c12f(iut,jut).e(m)=v12f(m).a(iut,jut)+(-rmme)*v12f(m).b(iut & ,jut)+rmme*v12f(m).c(iut,jut)+(-rmme)*rmme*v12f(m).d(iut,jut & ) enddo enddo end do * T -- qu=p1,qd=p2,v=0,a=v12z(0).a,b=v12z(0).b,c=v12z(0).c,d=v12z(0).d,cr=c * fac1z,cl=cfac2z,nsum=0 eps_0=-p1(2)*p2(3)+p2(2)*p1(3) ceps_0=eps_0*cim ceps_1=p1(3)*cim ceps_2=p2(3)*cim auxa=-quqd+p1k0*p2(0)+p2k0*p1(0) v12z(0).a(1,1)=cfac1z*(auxa+ceps_0) v12z(0).a(2,2)=cfac2z*(auxa-ceps_0) v12z(0).b(1,2)=-cfac2z*(p2(2)+ceps_2) v12z(0).b(2,1)=cfac1z*(p2(2)-ceps_2) v12z(0).c(1,2)=cfac1z*(p1(2)+ceps_1) v12z(0).c(2,1)=cfac2z*(-p1(2)+ceps_1) v12z(0).d(1,1)=cfac2z v12z(0).d(2,2)=cfac1z * T -- qu=p1,qd=p2,v=1,a=v12z(1).a,b=v12z(1).b,c=v12z(1).c,d=v12z(1).d,cr=c * fac1z,cl=cfac2z,nsum=0 auxa=-quqd+p1k0*p2(1)+p2k0*p1(1) v12z(1).a(1,1)=cfac1z*(auxa+ceps_0) v12z(1).a(2,2)=cfac2z*(auxa-ceps_0) v12z(1).b(1,2)=-cfac2z*(p2(2)+ceps_2) v12z(1).b(2,1)=cfac1z*(p2(2)-ceps_2) v12z(1).c(1,2)=cfac1z*(p1(2)+ceps_1) v12z(1).c(2,1)=cfac2z*(-p1(2)+ceps_1) v12z(1).d(1,1)=cfac2z v12z(1).d(2,2)=cfac1z * T -- qu=p1,qd=p2,v=2,a=v12z(2).a,b=v12z(2).b,c=v12z(2).c,d=v12z(2).d,cr=c * fac1z,cl=cfac2z,nsum=0 eps_0=-p1k0*p2(3)+p2k0*p1(3) ceps_0=eps_0*cim auxa=p1k0*p2(2)+p2k0*p1(2) v12z(2).a(1,1)=cfac1z*(auxa+ceps_0) v12z(2).a(2,2)=cfac2z*(auxa-ceps_0) v12z(2).b(1,2)=-cfac2z*p2k0 v12z(2).b(2,1)=cfac1z*p2k0 v12z(2).c(1,2)=cfac1z*p1k0 v12z(2).c(2,1)=-cfac2z*p1k0 * T -- qu=p1,qd=p2,v=3,a=v12z(3).a,b=v12z(3).b,c=v12z(3).c,d=v12z(3).d,cr=c * fac1z,cl=cfac2z,nsum=0 eps_0=p1k0*p2(2)-p2k0*p1(2) ceps_0=eps_0*cim ceps_1=p1k0*cim ceps_2=p2k0*cim auxa=+p1k0*p2(3)+p2k0*p1(3) v12z(3).a(1,1)=cfac1z*(auxa+ceps_0) v12z(3).a(2,2)=cfac2z*(auxa-ceps_0) v12z(3).b(1,2)=-cfac2z*ceps_2 v12z(3).b(2,1)=-cfac1z*ceps_2 v12z(3).c(1,2)=cfac1z*ceps_1 v12z(3).c(2,1)=cfac2z*ceps_1 do m=0,3 * mline -- res=c12z(&).e(m),abcd=v12z(m).,m1=(-rmme),m2=rmme,den=0 do iut=1,2 do jut=1,2 c12z(iut,jut).e(m)=v12z(m).a(iut,jut)+(-rmme)*v12z(m).b(iut & ,jut)+rmme*v12z(m).c(iut,jut)+(-rmme)*rmme*v12z(m).d(iut,jut & ) enddo enddo end do **** moltiplico per i propagatori bosonici (G_mn - p_m p_n/ Mq)/dp do i1=1,2 do i2=1,2 c12z(i1,i2).v=c12z(i1,i2).e(0)*p12(0) do i=1,3 c12z(i1,i2).v=c12z(i1,i2).v- & c12z(i1,i2).e(i)*p12(i) enddo enddo enddo do i1=1,2 do i2=1,2 do mu=0,3 c12z(i1,i2).e(mu)=c12z(i1,i2).e(mu)-p12(mu) & *c12z(i1,i2).v/rmz2 enddo enddo enddo **** do i1=1,2 do i2=1,2 * pk0 -- p=c12f(i1,i2).e c12f(i1,i2).ek0=c12f(i1,i2).e(0)-c12f(i1,i2).e(1) end do end do do i1=1,2 do i2=1,2 * pk0 -- p=c12z(i1,i2).e c12z(i1,i2).ek0=c12z(i1,i2).e(0)-c12z(i1,i2).e(1) end do end do do ide=0,iid if (ide.eq.1) then * p3<->p5 per rifare tutti i diagrammi con questo scambio per particelle id. do mu=0,3 pau(mu)=p3(mu) p3(mu)=p5(mu) p5(mu)=pau(mu) pau(mu)=p25(mu) p25(mu)=p23(mu) p23(mu)=pau(mu) pau(mu)=p34(mu) p34(mu)=p45(mu) p45(mu)=pau(mu) pau(mu)=p36(mu) p36(mu)=p56(mu) p56(mu)=pau(mu) pau(mu)=p123(mu) p123(mu)=p125(mu) p125(mu)=pau(mu) pau(mu)=p134(mu) p134(mu)=p145(mu) p145(mu)=pau(mu) pau(mu)=p156(mu) p156(mu)=p136(mu) p136(mu)=pau(mu) end do pkau=p3k0 p3k0=p5k0 p5k0=pkau pkau=p3q p3q=p5q p5q=pkau pkau=p123k0 p123k0=p125k0 p125k0=pkau pkau=p123q p123q=p125q p125q=pkau pkau=p134k0 p134k0=p145k0 p145k0=pkau pkau=p134q p134q=p145q p145q=pkau pkau=p156k0 p156k0=p136k0 p136k0=pkau pkau=p156q p156q=p136q p136q=pkau xmau=xm(3) xm(3)=xm(5) xm(5)=xmau endif * quqd -- p=p5,q=p6 quqd=p5(0)*p6(0)-p5(1)*p6(1)-p5(2)*p6(2)-p5(3)*p6(3) df=-1.d0/(xm(5)**2+xm(6)**2+2.d0*quqd) fac1=(df*f5r) fac2=(df*f5l) * T -- qu=p5,qd=p6,v=0,a=v56f(0).a,b=v56f(0).b,c=v56f(0).c,d=v56f(0).d,cr=f * ac1,cl=fac2,nsum=0 eps_0=-p5(2)*p6(3)+p6(2)*p5(3) ceps_0=eps_0*cim ceps_1=p5(3)*cim ceps_2=p6(3)*cim auxa=-quqd+p5k0*p6(0)+p6k0*p5(0) v56f(0).a(1,1)=fac1*(auxa+ceps_0) v56f(0).a(2,2)=fac2*(auxa-ceps_0) v56f(0).b(1,2)=-fac2*(p6(2)+ceps_2) v56f(0).b(2,1)=fac1*(p6(2)-ceps_2) v56f(0).c(1,2)=fac1*(p5(2)+ceps_1) v56f(0).c(2,1)=fac2*(-p5(2)+ceps_1) v56f(0).d(1,1)=fac2 v56f(0).d(2,2)=fac1 * T -- qu=p5,qd=p6,v=1,a=v56f(1).a,b=v56f(1).b,c=v56f(1).c,d=v56f(1).d,cr=f * ac1,cl=fac2,nsum=0 auxa=-quqd+p5k0*p6(1)+p6k0*p5(1) v56f(1).a(1,1)=fac1*(auxa+ceps_0) v56f(1).a(2,2)=fac2*(auxa-ceps_0) v56f(1).b(1,2)=-fac2*(p6(2)+ceps_2) v56f(1).b(2,1)=fac1*(p6(2)-ceps_2) v56f(1).c(1,2)=fac1*(p5(2)+ceps_1) v56f(1).c(2,1)=fac2*(-p5(2)+ceps_1) v56f(1).d(1,1)=fac2 v56f(1).d(2,2)=fac1 * T -- qu=p5,qd=p6,v=2,a=v56f(2).a,b=v56f(2).b,c=v56f(2).c,d=v56f(2).d,cr=f * ac1,cl=fac2,nsum=0 eps_0=-p5k0*p6(3)+p6k0*p5(3) ceps_0=eps_0*cim auxa=p5k0*p6(2)+p6k0*p5(2) v56f(2).a(1,1)=fac1*(auxa+ceps_0) v56f(2).a(2,2)=fac2*(auxa-ceps_0) v56f(2).b(1,2)=-fac2*p6k0 v56f(2).b(2,1)=fac1*p6k0 v56f(2).c(1,2)=fac1*p5k0 v56f(2).c(2,1)=-fac2*p5k0 * T -- qu=p5,qd=p6,v=3,a=v56f(3).a,b=v56f(3).b,c=v56f(3).c,d=v56f(3).d,cr=f * ac1,cl=fac2,nsum=0 eps_0=p5k0*p6(2)-p6k0*p5(2) ceps_0=eps_0*cim ceps_1=p5k0*cim ceps_2=p6k0*cim auxa=+p5k0*p6(3)+p6k0*p5(3) v56f(3).a(1,1)=fac1*(auxa+ceps_0) v56f(3).a(2,2)=fac2*(auxa-ceps_0) v56f(3).b(1,2)=-fac2*ceps_2 v56f(3).b(2,1)=-fac1*ceps_2 v56f(3).c(1,2)=fac1*ceps_1 v56f(3).c(2,1)=fac2*ceps_1 do m=0,3 * mline -- res=c56f(&).e(m),abcd=v56f(m).,m1=xm(5),m2=(-xm(6)),den=0 do iut=1,2 do jut=1,2 c56f(iut,jut).e(m)=v56f(m).a(iut,jut)+xm(5)*v56f(m).b(iut, & jut)+(-xm(6))*v56f(m).c(iut,jut)+xm(5)*(-xm(6))*v56f(m).d & (iut,jut) enddo enddo end do do i5=1,2 do i6=1,2 * pk0 -- p=c56f(i5,i6).e c56f(i5,i6).ek0=c56f(i5,i6).e(0)-c56f(i5,i6).e(1) end do end do * quqd -- p=p3,q=p4 quqd=p3(0)*p4(0)-p3(1)*p4(1)-p3(2)*p4(2)-p3(3)*p4(3) df=-1.d0/(xm(3)**2+xm(4)**2+2.d0*quqd) fac1=(df*f3r) fac2=(df*f3l) * T -- qu=p3,qd=p4,v=0,a=v34f(0).a,b=v34f(0).b,c=v34f(0).c,d=v34f(0).d,cr=f * ac1,cl=fac2,nsum=0 eps_0=-p3(2)*p4(3)+p4(2)*p3(3) ceps_0=eps_0*cim ceps_1=p3(3)*cim ceps_2=p4(3)*cim auxa=-quqd+p3k0*p4(0)+p4k0*p3(0) v34f(0).a(1,1)=fac1*(auxa+ceps_0) v34f(0).a(2,2)=fac2*(auxa-ceps_0) v34f(0).b(1,2)=-fac2*(p4(2)+ceps_2) v34f(0).b(2,1)=fac1*(p4(2)-ceps_2) v34f(0).c(1,2)=fac1*(p3(2)+ceps_1) v34f(0).c(2,1)=fac2*(-p3(2)+ceps_1) v34f(0).d(1,1)=fac2 v34f(0).d(2,2)=fac1 * T -- qu=p3,qd=p4,v=1,a=v34f(1).a,b=v34f(1).b,c=v34f(1).c,d=v34f(1).d,cr=f * ac1,cl=fac2,nsum=0 auxa=-quqd+p3k0*p4(1)+p4k0*p3(1) v34f(1).a(1,1)=fac1*(auxa+ceps_0) v34f(1).a(2,2)=fac2*(auxa-ceps_0) v34f(1).b(1,2)=-fac2*(p4(2)+ceps_2) v34f(1).b(2,1)=fac1*(p4(2)-ceps_2) v34f(1).c(1,2)=fac1*(p3(2)+ceps_1) v34f(1).c(2,1)=fac2*(-p3(2)+ceps_1) v34f(1).d(1,1)=fac2 v34f(1).d(2,2)=fac1 * T -- qu=p3,qd=p4,v=2,a=v34f(2).a,b=v34f(2).b,c=v34f(2).c,d=v34f(2).d,cr=f * ac1,cl=fac2,nsum=0 eps_0=-p3k0*p4(3)+p4k0*p3(3) ceps_0=eps_0*cim auxa=p3k0*p4(2)+p4k0*p3(2) v34f(2).a(1,1)=fac1*(auxa+ceps_0) v34f(2).a(2,2)=fac2*(auxa-ceps_0) v34f(2).b(1,2)=-fac2*p4k0 v34f(2).b(2,1)=fac1*p4k0 v34f(2).c(1,2)=fac1*p3k0 v34f(2).c(2,1)=-fac2*p3k0 * T -- qu=p3,qd=p4,v=3,a=v34f(3).a,b=v34f(3).b,c=v34f(3).c,d=v34f(3).d,cr=f * ac1,cl=fac2,nsum=0 eps_0=p3k0*p4(2)-p4k0*p3(2) ceps_0=eps_0*cim ceps_1=p3k0*cim ceps_2=p4k0*cim auxa=+p3k0*p4(3)+p4k0*p3(3) v34f(3).a(1,1)=fac1*(auxa+ceps_0) v34f(3).a(2,2)=fac2*(auxa-ceps_0) v34f(3).b(1,2)=-fac2*ceps_2 v34f(3).b(2,1)=-fac1*ceps_2 v34f(3).c(1,2)=fac1*ceps_1 v34f(3).c(2,1)=fac2*ceps_1 do m=0,3 * mline -- res=c34f(&).e(m),abcd=v34f(m).,m1=xm(3),m2=(-xm(4)),den=0 do iut=1,2 do jut=1,2 c34f(iut,jut).e(m)=v34f(m).a(iut,jut)+xm(3)*v34f(m).b(iut, & jut)+(-xm(4))*v34f(m).c(iut,jut)+xm(3)*(-xm(4))*v34f(m).d & (iut,jut) enddo enddo end do do i3=1,2 do i4=1,2 * pk0 -- p=c34f(i3,i4).e c34f(i3,i4).ek0=c34f(i3,i4).e(0)-c34f(i3,i4).e(1) end do end do *****attaccamento di c12f(2) a 3 * quqd -- p=p3,q=p123 quqd=p3(0)*p123(0)-p3(1)*p123(1)-p3(2)*p123(2)-p3(3)*p123( & 3) do i1=1,2 do i2=1,2 * T -- qu=p3,qd=p123,v=c12f(i1,i2).e,a=l3_12fz(i1,i2).a,b=l3_12fz(i1,i2).b, * c=l3_12fz(i1,i2).c,d=l3_12fz(i1,i2).d,cr=f3r,cl=f3l,nsum=0 ceps_0=-c12f(i1,i2).ek0*(p3(2)*p123(3)-p123(2)*p3(3))+p3k0 & *(c12f(i1,i2).e(2)*p123(3)-p123(2)*c12f(i1,i2).e(3))-p123 & k0*(c12f(i1,i2).e(2)*p3(3)-p3(2)*c12f(i1,i2).e(3)) ceps_0=ceps_0*cim ceps_1=-c12f(i1,i2).e(3)*p3k0+p3(3)*c12f(i1,i2).ek0 ceps_1=ceps_1*cim ceps_2=-c12f(i1,i2).e(3)*p123k0+p123(3)*c12f(i1,i2).ek0 ceps_2=ceps_2*cim cvqu=c12f(i1,i2).e(0)*p3(0)-c12f(i1,i2).e(1)*p3(1)-c12f(i1 & ,i2).e(2)*p3(2)-c12f(i1,i2).e(3)*p3(3) cvqd=c12f(i1,i2).e(0)*p123(0)-c12f(i1,i2).e(1)*p123(1)-c12 & f(i1,i2).e(2)*p123(2)-c12f(i1,i2).e(3)*p123(3) cauxa=-c12f(i1,i2).ek0*quqd+p3k0*cvqd+p123k0*cvqu cauxb=-c12f(i1,i2).ek0*p123(2)+p123k0*c12f(i1,i2).e(2) cauxc=+c12f(i1,i2).ek0*p3(2)-p3k0*c12f(i1,i2).e(2) l3_12fz(i1,i2).a(1,1)=f3r*(cauxa+ceps_0) l3_12fz(i1,i2).a(2,2)=f3l*(cauxa-ceps_0) l3_12fz(i1,i2).b(1,2)=f3l*(cauxb-ceps_2) l3_12fz(i1,i2).b(2,1)=f3r*(-cauxb-ceps_2) l3_12fz(i1,i2).c(1,2)=f3r*(cauxc+ceps_1) l3_12fz(i1,i2).c(2,1)=f3l*(-cauxc+ceps_1) l3_12fz(i1,i2).d(1,1)=f3l*c12f(i1,i2).ek0 l3_12fz(i1,i2).d(2,2)=f3r*c12f(i1,i2).ek0 end do end do *****attaccamento di c12z(2) a 3 do i1=1,2 do i2=1,2 * T -- qu=p3,qd=p123,v=c12z(i1,i2).e,a=l3_12fz(i1,i2).a,b=l3_12fz(i1,i2).b, * c=l3_12fz(i1,i2).c,d=l3_12fz(i1,i2).d,cr=z3r,cl=z3l,nsum=1 ceps_0=-c12z(i1,i2).ek0*(p3(2)*p123(3)-p123(2)*p3(3))+p3k0 & *(c12z(i1,i2).e(2)*p123(3)-p123(2)*c12z(i1,i2).e(3))-p123 & k0*(c12z(i1,i2).e(2)*p3(3)-p3(2)*c12z(i1,i2).e(3)) ceps_0=ceps_0*cim ceps_1=-c12z(i1,i2).e(3)*p3k0+p3(3)*c12z(i1,i2).ek0 ceps_1=ceps_1*cim ceps_2=-c12z(i1,i2).e(3)*p123k0+p123(3)*c12z(i1,i2).ek0 ceps_2=ceps_2*cim cvqu=c12z(i1,i2).e(0)*p3(0)-c12z(i1,i2).e(1)*p3(1)-c12z(i1 & ,i2).e(2)*p3(2)-c12z(i1,i2).e(3)*p3(3) cvqd=c12z(i1,i2).e(0)*p123(0)-c12z(i1,i2).e(1)*p123(1)-c12 & z(i1,i2).e(2)*p123(2)-c12z(i1,i2).e(3)*p123(3) cauxa=-c12z(i1,i2).ek0*quqd+p3k0*cvqd+p123k0*cvqu cauxb=-c12z(i1,i2).ek0*p123(2)+p123k0*c12z(i1,i2).e(2) cauxc=+c12z(i1,i2).ek0*p3(2)-p3k0*c12z(i1,i2).e(2) l3_12fz(i1,i2).a(1,1)=l3_12fz(i1,i2).a(1,1)+z3r*(cauxa+cep & s_0) l3_12fz(i1,i2).a(2,2)=l3_12fz(i1,i2).a(2,2)+z3l*(cauxa-cep & s_0) l3_12fz(i1,i2).b(1,2)=l3_12fz(i1,i2).b(1,2)+z3l*(cauxb-cep & s_2) l3_12fz(i1,i2).b(2,1)=l3_12fz(i1,i2).b(2,1)+z3r*(-cauxb-ce & ps_2) l3_12fz(i1,i2).c(1,2)=l3_12fz(i1,i2).c(1,2)+z3r*(cauxc+cep & s_1) l3_12fz(i1,i2).c(2,1)=l3_12fz(i1,i2).c(2,1)+z3l*(-cauxc+ce & ps_1) l3_12fz(i1,i2).d(1,1)=l3_12fz(i1,i2).d(1,1)+z3l*c12z(i1,i2 & ).ek0 l3_12fz(i1,i2).d(2,2)=l3_12fz(i1,i2).d(2,2)+z3r*c12z(i1,i2 & ).ek0 end do end do *****attaccamento di c12f(2) a 4 * quqd -- p=p356,q=p4 quqd=p356(0)*p4(0)-p356(1)*p4(1)-p356(2)*p4(2)-p356(3)*p4( & 3) do i1=1,2 do i2=1,2 * T -- qu=p356,qd=p4,v=c12f(i1,i2).e,a=r4_12fz(i1,i2).a,b=r4_12fz(i1,i2).b, * c=r4_12fz(i1,i2).c,d=r4_12fz(i1,i2).d,cr=f4r,cl=f4l,nsum=0 ceps_0=-c12f(i1,i2).ek0*(p356(2)*p4(3)-p4(2)*p356(3))+p356 & k0*(c12f(i1,i2).e(2)*p4(3)-p4(2)*c12f(i1,i2).e(3))-p4k0*( & c12f(i1,i2).e(2)*p356(3)-p356(2)*c12f(i1,i2).e(3)) ceps_0=ceps_0*cim ceps_1=-c12f(i1,i2).e(3)*p356k0+p356(3)*c12f(i1,i2).ek0 ceps_1=ceps_1*cim ceps_2=-c12f(i1,i2).e(3)*p4k0+p4(3)*c12f(i1,i2).ek0 ceps_2=ceps_2*cim cvqu=c12f(i1,i2).e(0)*p356(0)-c12f(i1,i2).e(1)*p356(1)-c12 & f(i1,i2).e(2)*p356(2)-c12f(i1,i2).e(3)*p356(3) cvqd=c12f(i1,i2).e(0)*p4(0)-c12f(i1,i2).e(1)*p4(1)-c12f(i1 & ,i2).e(2)*p4(2)-c12f(i1,i2).e(3)*p4(3) cauxa=-c12f(i1,i2).ek0*quqd+p356k0*cvqd+p4k0*cvqu cauxb=-c12f(i1,i2).ek0*p4(2)+p4k0*c12f(i1,i2).e(2) cauxc=+c12f(i1,i2).ek0*p356(2)-p356k0*c12f(i1,i2).e(2) r4_12fz(i1,i2).a(1,1)=f4r*(cauxa+ceps_0) r4_12fz(i1,i2).a(2,2)=f4l*(cauxa-ceps_0) r4_12fz(i1,i2).b(1,2)=f4l*(cauxb-ceps_2) r4_12fz(i1,i2).b(2,1)=f4r*(-cauxb-ceps_2) r4_12fz(i1,i2).c(1,2)=f4r*(cauxc+ceps_1) r4_12fz(i1,i2).c(2,1)=f4l*(-cauxc+ceps_1) r4_12fz(i1,i2).d(1,1)=f4l*c12f(i1,i2).ek0 r4_12fz(i1,i2).d(2,2)=f4r*c12f(i1,i2).ek0 end do end do *****attaccamento di c12z(2) a 4 do i1=1,2 do i2=1,2 * T -- qu=p356,qd=p4,v=c12z(i1,i2).e,a=r4_12fz(i1,i2).a,b=r4_12fz(i1,i2).b, * c=r4_12fz(i1,i2).c,d=r4_12fz(i1,i2).d,cr=z4r,cl=z4l,nsum=1 ceps_0=-c12z(i1,i2).ek0*(p356(2)*p4(3)-p4(2)*p356(3))+p356 & k0*(c12z(i1,i2).e(2)*p4(3)-p4(2)*c12z(i1,i2).e(3))-p4k0*( & c12z(i1,i2).e(2)*p356(3)-p356(2)*c12z(i1,i2).e(3)) ceps_0=ceps_0*cim ceps_1=-c12z(i1,i2).e(3)*p356k0+p356(3)*c12z(i1,i2).ek0 ceps_1=ceps_1*cim ceps_2=-c12z(i1,i2).e(3)*p4k0+p4(3)*c12z(i1,i2).ek0 ceps_2=ceps_2*cim cvqu=c12z(i1,i2).e(0)*p356(0)-c12z(i1,i2).e(1)*p356(1)-c12 & z(i1,i2).e(2)*p356(2)-c12z(i1,i2).e(3)*p356(3) cvqd=c12z(i1,i2).e(0)*p4(0)-c12z(i1,i2).e(1)*p4(1)-c12z(i1 & ,i2).e(2)*p4(2)-c12z(i1,i2).e(3)*p4(3) cauxa=-c12z(i1,i2).ek0*quqd+p356k0*cvqd+p4k0*cvqu cauxb=-c12z(i1,i2).ek0*p4(2)+p4k0*c12z(i1,i2).e(2) cauxc=+c12z(i1,i2).ek0*p356(2)-p356k0*c12z(i1,i2).e(2) r4_12fz(i1,i2).a(1,1)=r4_12fz(i1,i2).a(1,1)+z4r*(cauxa+cep & s_0) r4_12fz(i1,i2).a(2,2)=r4_12fz(i1,i2).a(2,2)+z4l*(cauxa-cep & s_0) r4_12fz(i1,i2).b(1,2)=r4_12fz(i1,i2).b(1,2)+z4l*(cauxb-cep & s_2) r4_12fz(i1,i2).b(2,1)=r4_12fz(i1,i2).b(2,1)+z4r*(-cauxb-ce & ps_2) r4_12fz(i1,i2).c(1,2)=r4_12fz(i1,i2).c(1,2)+z4r*(cauxc+cep & s_1) r4_12fz(i1,i2).c(2,1)=r4_12fz(i1,i2).c(2,1)+z4l*(-cauxc+ce & ps_1) r4_12fz(i1,i2).d(1,1)=r4_12fz(i1,i2).d(1,1)+z4l*c12z(i1,i2 & ).ek0 r4_12fz(i1,i2).d(2,2)=r4_12fz(i1,i2).d(2,2)+z4r*c12z(i1,i2 & ).ek0 end do end do *****attaccamento di c56f a 3 * quqd -- p=p3,q=p356 quqd=p3(0)*p356(0)-p3(1)*p356(1)-p3(2)*p356(2)-p3(3)*p356( & 3) do i5=1,2 do i6=1,2 * T -- qu=p3,qd=p356,v=c56f(i5,i6).e,a=l3_56f(i5,i6).a,b=l3_56f(i5,i6).b,c= * l3_56f(i5,i6).c,d=l3_56f(i5,i6).d,cr=f3r,cl=f3l,nsum=0 ceps_0=-c56f(i5,i6).ek0*(p3(2)*p356(3)-p356(2)*p3(3))+p3k0 & *(c56f(i5,i6).e(2)*p356(3)-p356(2)*c56f(i5,i6).e(3))-p356 & k0*(c56f(i5,i6).e(2)*p3(3)-p3(2)*c56f(i5,i6).e(3)) ceps_0=ceps_0*cim ceps_1=-c56f(i5,i6).e(3)*p3k0+p3(3)*c56f(i5,i6).ek0 ceps_1=ceps_1*cim ceps_2=-c56f(i5,i6).e(3)*p356k0+p356(3)*c56f(i5,i6).ek0 ceps_2=ceps_2*cim cvqu=c56f(i5,i6).e(0)*p3(0)-c56f(i5,i6).e(1)*p3(1)-c56f(i5 & ,i6).e(2)*p3(2)-c56f(i5,i6).e(3)*p3(3) cvqd=c56f(i5,i6).e(0)*p356(0)-c56f(i5,i6).e(1)*p356(1)-c56 & f(i5,i6).e(2)*p356(2)-c56f(i5,i6).e(3)*p356(3) cauxa=-c56f(i5,i6).ek0*quqd+p3k0*cvqd+p356k0*cvqu cauxb=-c56f(i5,i6).ek0*p356(2)+p356k0*c56f(i5,i6).e(2) cauxc=+c56f(i5,i6).ek0*p3(2)-p3k0*c56f(i5,i6).e(2) l3_56f(i5,i6).a(1,1)=f3r*(cauxa+ceps_0) l3_56f(i5,i6).a(2,2)=f3l*(cauxa-ceps_0) l3_56f(i5,i6).b(1,2)=f3l*(cauxb-ceps_2) l3_56f(i5,i6).b(2,1)=f3r*(-cauxb-ceps_2) l3_56f(i5,i6).c(1,2)=f3r*(cauxc+ceps_1) l3_56f(i5,i6).c(2,1)=f3l*(-cauxc+ceps_1) l3_56f(i5,i6).d(1,1)=f3l*c56f(i5,i6).ek0 l3_56f(i5,i6).d(2,2)=f3r*c56f(i5,i6).ek0 end do end do *****attaccamento di c56f a 4 * quqd -- p=p123,q=p4 quqd=p123(0)*p4(0)-p123(1)*p4(1)-p123(2)*p4(2)-p123(3)*p4( & 3) do i5=1,2 do i6=1,2 * T -- qu=p123,qd=p4,v=c56f(i5,i6).e,a=r4_56f(i5,i6).a,b=r4_56f(i5,i6).b,c= * r4_56f(i5,i6).c,d=r4_56f(i5,i6).d,cr=f4r,cl=f4l,nsum=0 ceps_0=-c56f(i5,i6).ek0*(p123(2)*p4(3)-p4(2)*p123(3))+p123 & k0*(c56f(i5,i6).e(2)*p4(3)-p4(2)*c56f(i5,i6).e(3))-p4k0*( & c56f(i5,i6).e(2)*p123(3)-p123(2)*c56f(i5,i6).e(3)) ceps_0=ceps_0*cim ceps_1=-c56f(i5,i6).e(3)*p123k0+p123(3)*c56f(i5,i6).ek0 ceps_1=ceps_1*cim ceps_2=-c56f(i5,i6).e(3)*p4k0+p4(3)*c56f(i5,i6).ek0 ceps_2=ceps_2*cim cvqu=c56f(i5,i6).e(0)*p123(0)-c56f(i5,i6).e(1)*p123(1)-c56 & f(i5,i6).e(2)*p123(2)-c56f(i5,i6).e(3)*p123(3) cvqd=c56f(i5,i6).e(0)*p4(0)-c56f(i5,i6).e(1)*p4(1)-c56f(i5 & ,i6).e(2)*p4(2)-c56f(i5,i6).e(3)*p4(3) cauxa=-c56f(i5,i6).ek0*quqd+p123k0*cvqd+p4k0*cvqu cauxb=-c56f(i5,i6).ek0*p4(2)+p4k0*c56f(i5,i6).e(2) cauxc=+c56f(i5,i6).ek0*p123(2)-p123k0*c56f(i5,i6).e(2) r4_56f(i5,i6).a(1,1)=f4r*(cauxa+ceps_0) r4_56f(i5,i6).a(2,2)=f4l*(cauxa-ceps_0) r4_56f(i5,i6).b(1,2)=f4l*(cauxb-ceps_2) r4_56f(i5,i6).b(2,1)=f4r*(-cauxb-ceps_2) r4_56f(i5,i6).c(1,2)=f4r*(cauxc+ceps_1) r4_56f(i5,i6).c(2,1)=f4l*(-cauxc+ceps_1) r4_56f(i5,i6).d(1,1)=f4l*c56f(i5,i6).ek0 r4_56f(i5,i6).d(2,2)=f4r*c56f(i5,i6).ek0 end do end do *****attaccamento di c12f(2) a 5 * quqd -- p=p5,q=p125 quqd=p5(0)*p125(0)-p5(1)*p125(1)-p5(2)*p125(2)-p5(3)*p125( & 3) do i1=1,2 do i2=1,2 * T -- qu=p5,qd=p125,v=c12f(i1,i2).e,a=l5_12fz(i1,i2).a,b=l5_12fz(i1,i2).b, * c=l5_12fz(i1,i2).c,d=l5_12fz(i1,i2).d,cr=f5r,cl=f5l,nsum=0 ceps_0=-c12f(i1,i2).ek0*(p5(2)*p125(3)-p125(2)*p5(3))+p5k0 & *(c12f(i1,i2).e(2)*p125(3)-p125(2)*c12f(i1,i2).e(3))-p125 & k0*(c12f(i1,i2).e(2)*p5(3)-p5(2)*c12f(i1,i2).e(3)) ceps_0=ceps_0*cim ceps_1=-c12f(i1,i2).e(3)*p5k0+p5(3)*c12f(i1,i2).ek0 ceps_1=ceps_1*cim ceps_2=-c12f(i1,i2).e(3)*p125k0+p125(3)*c12f(i1,i2).ek0 ceps_2=ceps_2*cim cvqu=c12f(i1,i2).e(0)*p5(0)-c12f(i1,i2).e(1)*p5(1)-c12f(i1 & ,i2).e(2)*p5(2)-c12f(i1,i2).e(3)*p5(3) cvqd=c12f(i1,i2).e(0)*p125(0)-c12f(i1,i2).e(1)*p125(1)-c12 & f(i1,i2).e(2)*p125(2)-c12f(i1,i2).e(3)*p125(3) cauxa=-c12f(i1,i2).ek0*quqd+p5k0*cvqd+p125k0*cvqu cauxb=-c12f(i1,i2).ek0*p125(2)+p125k0*c12f(i1,i2).e(2) cauxc=+c12f(i1,i2).ek0*p5(2)-p5k0*c12f(i1,i2).e(2) l5_12fz(i1,i2).a(1,1)=f5r*(cauxa+ceps_0) l5_12fz(i1,i2).a(2,2)=f5l*(cauxa-ceps_0) l5_12fz(i1,i2).b(1,2)=f5l*(cauxb-ceps_2) l5_12fz(i1,i2).b(2,1)=f5r*(-cauxb-ceps_2) l5_12fz(i1,i2).c(1,2)=f5r*(cauxc+ceps_1) l5_12fz(i1,i2).c(2,1)=f5l*(-cauxc+ceps_1) l5_12fz(i1,i2).d(1,1)=f5l*c12f(i1,i2).ek0 l5_12fz(i1,i2).d(2,2)=f5r*c12f(i1,i2).ek0 end do end do *****attaccamento di c12z(2) a 5 do i1=1,2 do i2=1,2 * T -- qu=p5,qd=p125,v=c12z(i1,i2).e,a=l5_12fz(i1,i2).a,b=l5_12fz(i1,i2).b, * c=l5_12fz(i1,i2).c,d=l5_12fz(i1,i2).d,cr=z5r,cl=z5l,nsum=1 ceps_0=-c12z(i1,i2).ek0*(p5(2)*p125(3)-p125(2)*p5(3))+p5k0 & *(c12z(i1,i2).e(2)*p125(3)-p125(2)*c12z(i1,i2).e(3))-p125 & k0*(c12z(i1,i2).e(2)*p5(3)-p5(2)*c12z(i1,i2).e(3)) ceps_0=ceps_0*cim ceps_1=-c12z(i1,i2).e(3)*p5k0+p5(3)*c12z(i1,i2).ek0 ceps_1=ceps_1*cim ceps_2=-c12z(i1,i2).e(3)*p125k0+p125(3)*c12z(i1,i2).ek0 ceps_2=ceps_2*cim cvqu=c12z(i1,i2).e(0)*p5(0)-c12z(i1,i2).e(1)*p5(1)-c12z(i1 & ,i2).e(2)*p5(2)-c12z(i1,i2).e(3)*p5(3) cvqd=c12z(i1,i2).e(0)*p125(0)-c12z(i1,i2).e(1)*p125(1)-c12 & z(i1,i2).e(2)*p125(2)-c12z(i1,i2).e(3)*p125(3) cauxa=-c12z(i1,i2).ek0*quqd+p5k0*cvqd+p125k0*cvqu cauxb=-c12z(i1,i2).ek0*p125(2)+p125k0*c12z(i1,i2).e(2) cauxc=+c12z(i1,i2).ek0*p5(2)-p5k0*c12z(i1,i2).e(2) l5_12fz(i1,i2).a(1,1)=l5_12fz(i1,i2).a(1,1)+z5r*(cauxa+cep & s_0) l5_12fz(i1,i2).a(2,2)=l5_12fz(i1,i2).a(2,2)+z5l*(cauxa-cep & s_0) l5_12fz(i1,i2).b(1,2)=l5_12fz(i1,i2).b(1,2)+z5l*(cauxb-cep & s_2) l5_12fz(i1,i2).b(2,1)=l5_12fz(i1,i2).b(2,1)+z5r*(-cauxb-ce & ps_2) l5_12fz(i1,i2).c(1,2)=l5_12fz(i1,i2).c(1,2)+z5r*(cauxc+cep & s_1) l5_12fz(i1,i2).c(2,1)=l5_12fz(i1,i2).c(2,1)+z5l*(-cauxc+ce & ps_1) l5_12fz(i1,i2).d(1,1)=l5_12fz(i1,i2).d(1,1)+z5l*c12z(i1,i2 & ).ek0 l5_12fz(i1,i2).d(2,2)=l5_12fz(i1,i2).d(2,2)+z5r*c12z(i1,i2 & ).ek0 end do end do *****attaccamento di c12f(2) a 6 * quqd -- p=p345,q=p6 quqd=p345(0)*p6(0)-p345(1)*p6(1)-p345(2)*p6(2)-p345(3)*p6( & 3) do i1=1,2 do i2=1,2 * T -- qu=p345,qd=p6,v=c12f(i1,i2).e,a=r6_12fz(i1,i2).a,b=r6_12fz(i1,i2).b, * c=r6_12fz(i1,i2).c,d=r6_12fz(i1,i2).d,cr=f6r,cl=f6l,nsum=0 ceps_0=-c12f(i1,i2).ek0*(p345(2)*p6(3)-p6(2)*p345(3))+p345 & k0*(c12f(i1,i2).e(2)*p6(3)-p6(2)*c12f(i1,i2).e(3))-p6k0*( & c12f(i1,i2).e(2)*p345(3)-p345(2)*c12f(i1,i2).e(3)) ceps_0=ceps_0*cim ceps_1=-c12f(i1,i2).e(3)*p345k0+p345(3)*c12f(i1,i2).ek0 ceps_1=ceps_1*cim ceps_2=-c12f(i1,i2).e(3)*p6k0+p6(3)*c12f(i1,i2).ek0 ceps_2=ceps_2*cim cvqu=c12f(i1,i2).e(0)*p345(0)-c12f(i1,i2).e(1)*p345(1)-c12 & f(i1,i2).e(2)*p345(2)-c12f(i1,i2).e(3)*p345(3) cvqd=c12f(i1,i2).e(0)*p6(0)-c12f(i1,i2).e(1)*p6(1)-c12f(i1 & ,i2).e(2)*p6(2)-c12f(i1,i2).e(3)*p6(3) cauxa=-c12f(i1,i2).ek0*quqd+p345k0*cvqd+p6k0*cvqu cauxb=-c12f(i1,i2).ek0*p6(2)+p6k0*c12f(i1,i2).e(2) cauxc=+c12f(i1,i2).ek0*p345(2)-p345k0*c12f(i1,i2).e(2) r6_12fz(i1,i2).a(1,1)=f6r*(cauxa+ceps_0) r6_12fz(i1,i2).a(2,2)=f6l*(cauxa-ceps_0) r6_12fz(i1,i2).b(1,2)=f6l*(cauxb-ceps_2) r6_12fz(i1,i2).b(2,1)=f6r*(-cauxb-ceps_2) r6_12fz(i1,i2).c(1,2)=f6r*(cauxc+ceps_1) r6_12fz(i1,i2).c(2,1)=f6l*(-cauxc+ceps_1) r6_12fz(i1,i2).d(1,1)=f6l*c12f(i1,i2).ek0 r6_12fz(i1,i2).d(2,2)=f6r*c12f(i1,i2).ek0 end do end do *****attaccamento di c12z(2) a 6 do i1=1,2 do i2=1,2 * T -- qu=p345,qd=p6,v=c12z(i1,i2).e,a=r6_12fz(i1,i2).a,b=r6_12fz(i1,i2).b, * c=r6_12fz(i1,i2).c,d=r6_12fz(i1,i2).d,cr=z6r,cl=z6l,nsum=1 ceps_0=-c12z(i1,i2).ek0*(p345(2)*p6(3)-p6(2)*p345(3))+p345 & k0*(c12z(i1,i2).e(2)*p6(3)-p6(2)*c12z(i1,i2).e(3))-p6k0*( & c12z(i1,i2).e(2)*p345(3)-p345(2)*c12z(i1,i2).e(3)) ceps_0=ceps_0*cim ceps_1=-c12z(i1,i2).e(3)*p345k0+p345(3)*c12z(i1,i2).ek0 ceps_1=ceps_1*cim ceps_2=-c12z(i1,i2).e(3)*p6k0+p6(3)*c12z(i1,i2).ek0 ceps_2=ceps_2*cim cvqu=c12z(i1,i2).e(0)*p345(0)-c12z(i1,i2).e(1)*p345(1)-c12 & z(i1,i2).e(2)*p345(2)-c12z(i1,i2).e(3)*p345(3) cvqd=c12z(i1,i2).e(0)*p6(0)-c12z(i1,i2).e(1)*p6(1)-c12z(i1 & ,i2).e(2)*p6(2)-c12z(i1,i2).e(3)*p6(3) cauxa=-c12z(i1,i2).ek0*quqd+p345k0*cvqd+p6k0*cvqu cauxb=-c12z(i1,i2).ek0*p6(2)+p6k0*c12z(i1,i2).e(2) cauxc=+c12z(i1,i2).ek0*p345(2)-p345k0*c12z(i1,i2).e(2) r6_12fz(i1,i2).a(1,1)=r6_12fz(i1,i2).a(1,1)+z6r*(cauxa+cep & s_0) r6_12fz(i1,i2).a(2,2)=r6_12fz(i1,i2).a(2,2)+z6l*(cauxa-cep & s_0) r6_12fz(i1,i2).b(1,2)=r6_12fz(i1,i2).b(1,2)+z6l*(cauxb-cep & s_2) r6_12fz(i1,i2).b(2,1)=r6_12fz(i1,i2).b(2,1)+z6r*(-cauxb-ce & ps_2) r6_12fz(i1,i2).c(1,2)=r6_12fz(i1,i2).c(1,2)+z6r*(cauxc+cep & s_1) r6_12fz(i1,i2).c(2,1)=r6_12fz(i1,i2).c(2,1)+z6l*(-cauxc+ce & ps_1) r6_12fz(i1,i2).d(1,1)=r6_12fz(i1,i2).d(1,1)+z6l*c12z(i1,i2 & ).ek0 r6_12fz(i1,i2).d(2,2)=r6_12fz(i1,i2).d(2,2)+z6r*c12z(i1,i2 & ).ek0 end do end do *****attaccamento di c34f a 5 * quqd -- p=p5,q=p345 quqd=p5(0)*p345(0)-p5(1)*p345(1)-p5(2)*p345(2)-p5(3)*p345( & 3) do i3=1,2 do i4=1,2 * T -- qu=p5,qd=p345,v=c34f(i3,i4).e,a=l5_34f(i3,i4).a,b=l5_34f(i3,i4).b,c= * l5_34f(i3,i4).c,d=l5_34f(i3,i4).d,cr=f5r,cl=f5l,nsum=0 ceps_0=-c34f(i3,i4).ek0*(p5(2)*p345(3)-p345(2)*p5(3))+p5k0 & *(c34f(i3,i4).e(2)*p345(3)-p345(2)*c34f(i3,i4).e(3))-p345 & k0*(c34f(i3,i4).e(2)*p5(3)-p5(2)*c34f(i3,i4).e(3)) ceps_0=ceps_0*cim ceps_1=-c34f(i3,i4).e(3)*p5k0+p5(3)*c34f(i3,i4).ek0 ceps_1=ceps_1*cim ceps_2=-c34f(i3,i4).e(3)*p345k0+p345(3)*c34f(i3,i4).ek0 ceps_2=ceps_2*cim cvqu=c34f(i3,i4).e(0)*p5(0)-c34f(i3,i4).e(1)*p5(1)-c34f(i3 & ,i4).e(2)*p5(2)-c34f(i3,i4).e(3)*p5(3) cvqd=c34f(i3,i4).e(0)*p345(0)-c34f(i3,i4).e(1)*p345(1)-c34 & f(i3,i4).e(2)*p345(2)-c34f(i3,i4).e(3)*p345(3) cauxa=-c34f(i3,i4).ek0*quqd+p5k0*cvqd+p345k0*cvqu cauxb=-c34f(i3,i4).ek0*p345(2)+p345k0*c34f(i3,i4).e(2) cauxc=+c34f(i3,i4).ek0*p5(2)-p5k0*c34f(i3,i4).e(2) l5_34f(i3,i4).a(1,1)=f5r*(cauxa+ceps_0) l5_34f(i3,i4).a(2,2)=f5l*(cauxa-ceps_0) l5_34f(i3,i4).b(1,2)=f5l*(cauxb-ceps_2) l5_34f(i3,i4).b(2,1)=f5r*(-cauxb-ceps_2) l5_34f(i3,i4).c(1,2)=f5r*(cauxc+ceps_1) l5_34f(i3,i4).c(2,1)=f5l*(-cauxc+ceps_1) l5_34f(i3,i4).d(1,1)=f5l*c34f(i3,i4).ek0 l5_34f(i3,i4).d(2,2)=f5r*c34f(i3,i4).ek0 end do end do *****attaccamento di c34f a 6 * quqd -- p=p125,q=p6 quqd=p125(0)*p6(0)-p125(1)*p6(1)-p125(2)*p6(2)-p125(3)*p6( & 3) do i3=1,2 do i4=1,2 * T -- qu=p125,qd=p6,v=c34f(i3,i4).e,a=r6_34f(i3,i4).a,b=r6_34f(i3,i4).b,c= * r6_34f(i3,i4).c,d=r6_34f(i3,i4).d,cr=f6r,cl=f6l,nsum=0 ceps_0=-c34f(i3,i4).ek0*(p125(2)*p6(3)-p6(2)*p125(3))+p125 & k0*(c34f(i3,i4).e(2)*p6(3)-p6(2)*c34f(i3,i4).e(3))-p6k0*( & c34f(i3,i4).e(2)*p125(3)-p125(2)*c34f(i3,i4).e(3)) ceps_0=ceps_0*cim ceps_1=-c34f(i3,i4).e(3)*p125k0+p125(3)*c34f(i3,i4).ek0 ceps_1=ceps_1*cim ceps_2=-c34f(i3,i4).e(3)*p6k0+p6(3)*c34f(i3,i4).ek0 ceps_2=ceps_2*cim cvqu=c34f(i3,i4).e(0)*p125(0)-c34f(i3,i4).e(1)*p125(1)-c34 & f(i3,i4).e(2)*p125(2)-c34f(i3,i4).e(3)*p125(3) cvqd=c34f(i3,i4).e(0)*p6(0)-c34f(i3,i4).e(1)*p6(1)-c34f(i3 & ,i4).e(2)*p6(2)-c34f(i3,i4).e(3)*p6(3) cauxa=-c34f(i3,i4).ek0*quqd+p125k0*cvqd+p6k0*cvqu cauxb=-c34f(i3,i4).ek0*p6(2)+p6k0*c34f(i3,i4).e(2) cauxc=+c34f(i3,i4).ek0*p125(2)-p125k0*c34f(i3,i4).e(2) r6_34f(i3,i4).a(1,1)=f6r*(cauxa+ceps_0) r6_34f(i3,i4).a(2,2)=f6l*(cauxa-ceps_0) r6_34f(i3,i4).b(1,2)=f6l*(cauxb-ceps_2) r6_34f(i3,i4).b(2,1)=f6r*(-cauxb-ceps_2) r6_34f(i3,i4).c(1,2)=f6r*(cauxc+ceps_1) r6_34f(i3,i4).c(2,1)=f6l*(-cauxc+ceps_1) r6_34f(i3,i4).d(1,1)=f6l*c34f(i3,i4).ek0 r6_34f(i3,i4).d(2,2)=f6r*c34f(i3,i4).ek0 end do end do **** 8 diagrammi di qcd (4) * accoppiamenti right e left con fotone sono gli stessi r3r5=f3r*f5r **** Diagramma CN7 **** do i1=1,2 do i2=1,2 do i5=1,2 do i6=1,2 * TT -- aa=tapptt(i1,i2,i5,i6).a,bb=tapptt(i1,i2,i5,i6).b,cc=tapptt(i1,i2,i * 5,i6).c,dd=tapptt(i1,i2,i5,i6).d,a1=l3_56f(i5,i6).a,b1=l3_56f(i5,i6).b,c1= * l3_56f(i5,i6).c,d1=l3_56f(i5,i6).d,a2=r4_12fz(i1,i2).a,b2=r4_12fz(i1,i2).b * ,c2=r4_12fz(i1,i2).c,d2=r4_12fz(i1,i2).d,prq=p356q,m=xm(3) tapptt(i1,i2,i5,i6).a(1,1)=l3_56f(i5,i6).a(1,1)*r4_12fz(i1 & ,i2).a(1,1)+l3_56f(i5,i6).c(1,2)*p356q*r4_12fz(i1,i2).b(2 & ,1) tapptt(i1,i2,i5,i6).b(1,1)=xm(3)*(l3_56f(i5,i6).d(1,1)*r4_ & 12fz(i1,i2).a(1,1)+l3_56f(i5,i6).b(1,2)*r4_12fz(i1,i2).b( & 2,1)) tapptt(i1,i2,i5,i6).c(1,1)=xm(3)*(l3_56f(i5,i6).a(1,1)*r4_ & 12fz(i1,i2).d(1,1)+l3_56f(i5,i6).c(1,2)*r4_12fz(i1,i2).c( & 2,1)) tapptt(i1,i2,i5,i6).d(1,1)=l3_56f(i5,i6).d(1,1)*p356q*r4_1 & 2fz(i1,i2).d(1,1)+l3_56f(i5,i6).b(1,2)*r4_12fz(i1,i2).c(2 & ,1) tapptt(i1,i2,i5,i6).a(1,2)=xm(3)*(l3_56f(i5,i6).a(1,1)*r4_ & 12fz(i1,i2).b(1,2)+l3_56f(i5,i6).c(1,2)*r4_12fz(i1,i2).a( & 2,2)) tapptt(i1,i2,i5,i6).b(1,2)=l3_56f(i5,i6).d(1,1)*p356q*r4_1 & 2fz(i1,i2).b(1,2)+l3_56f(i5,i6).b(1,2)*r4_12fz(i1,i2).a(2 & ,2) tapptt(i1,i2,i5,i6).c(1,2)=l3_56f(i5,i6).a(1,1)*r4_12fz(i1 & ,i2).c(1,2)+l3_56f(i5,i6).c(1,2)*p356q*r4_12fz(i1,i2).d(2 & ,2) tapptt(i1,i2,i5,i6).d(1,2)=xm(3)*(l3_56f(i5,i6).d(1,1)*r4_ & 12fz(i1,i2).c(1,2)+l3_56f(i5,i6).b(1,2)*r4_12fz(i1,i2).d( & 2,2)) tapptt(i1,i2,i5,i6).a(2,1)=xm(3)*(l3_56f(i5,i6).c(2,1)*r4_ & 12fz(i1,i2).a(1,1)+l3_56f(i5,i6).a(2,2)*r4_12fz(i1,i2).b( & 2,1)) tapptt(i1,i2,i5,i6).b(2,1)=l3_56f(i5,i6).b(2,1)*r4_12fz(i1 & ,i2).a(1,1)+l3_56f(i5,i6).d(2,2)*p356q*r4_12fz(i1,i2).b(2 & ,1) tapptt(i1,i2,i5,i6).c(2,1)=l3_56f(i5,i6).c(2,1)*p356q*r4_1 & 2fz(i1,i2).d(1,1)+l3_56f(i5,i6).a(2,2)*r4_12fz(i1,i2).c(2 & ,1) tapptt(i1,i2,i5,i6).d(2,1)=xm(3)*(l3_56f(i5,i6).b(2,1)*r4_ & 12fz(i1,i2).d(1,1)+l3_56f(i5,i6).d(2,2)*r4_12fz(i1,i2).c( & 2,1)) tapptt(i1,i2,i5,i6).a(2,2)=l3_56f(i5,i6).c(2,1)*p356q*r4_1 & 2fz(i1,i2).b(1,2)+l3_56f(i5,i6).a(2,2)*r4_12fz(i1,i2).a(2 & ,2) tapptt(i1,i2,i5,i6).b(2,2)=xm(3)*(l3_56f(i5,i6).b(2,1)*r4_ & 12fz(i1,i2).b(1,2)+l3_56f(i5,i6).d(2,2)*r4_12fz(i1,i2).a( & 2,2)) tapptt(i1,i2,i5,i6).c(2,2)=xm(3)*(l3_56f(i5,i6).c(2,1)*r4_ & 12fz(i1,i2).c(1,2)+l3_56f(i5,i6).a(2,2)*r4_12fz(i1,i2).d( & 2,2)) tapptt(i1,i2,i5,i6).d(2,2)=l3_56f(i5,i6).b(2,1)*r4_12fz(i1 & ,i2).c(1,2)+l3_56f(i5,i6).d(2,2)*p356q*r4_12fz(i1,i2).d(2 & ,2) end do end do end do end do do i1=1,2 do i2=1,2 do i5=1,2 do i6=1,2 * mline -- res=cn7(i1,i2,&,i5,i6).id(ide),abcd=tapptt(i1,i2,i5,i6).,m1=xm(3 * ),m2=(-xm(4)),den=((p356q-xm(3)**2)*p356k0*r3r5) do iut=1,2 do jut=1,2 cn7(i1,i2,iut,jut,i5,i6).id(ide)=(tapptt(i1,i2,i5,i6).a(iu & t,jut)+xm(3)*tapptt(i1,i2,i5,i6).b(iut,jut)+(-xm(4))*tapp & tt(i1,i2,i5,i6).c(iut,jut)+xm(3)*(-xm(4))*tapptt(i1,i2,i5 & ,i6).d(iut,jut))/((p356q-xm(3)**2)*p356k0*r3r5) enddo enddo end do end do end do end do **** Diagramma CN8 **** do i1=1,2 do i2=1,2 do i5=1,2 do i6=1,2 * TT -- aa=tapptt(i1,i2,i5,i6).a,bb=tapptt(i1,i2,i5,i6).b,cc=tapptt(i1,i2,i * 5,i6).c,dd=tapptt(i1,i2,i5,i6).d,a1=l3_12fz(i1,i2).a,b1=l3_12fz(i1,i2).b,c * 1=l3_12fz(i1,i2).c,d1=l3_12fz(i1,i2).d,a2=r4_56f(i5,i6).a,b2=r4_56f(i5,i6) * .b,c2=r4_56f(i5,i6).c,d2=r4_56f(i5,i6).d,prq=p123q,m=xm(3) tapptt(i1,i2,i5,i6).a(1,1)=l3_12fz(i1,i2).a(1,1)*r4_56f(i5 & ,i6).a(1,1)+l3_12fz(i1,i2).c(1,2)*p123q*r4_56f(i5,i6).b(2 & ,1) tapptt(i1,i2,i5,i6).b(1,1)=xm(3)*(l3_12fz(i1,i2).d(1,1)*r4 & _56f(i5,i6).a(1,1)+l3_12fz(i1,i2).b(1,2)*r4_56f(i5,i6).b( & 2,1)) tapptt(i1,i2,i5,i6).c(1,1)=xm(3)*(l3_12fz(i1,i2).a(1,1)*r4 & _56f(i5,i6).d(1,1)+l3_12fz(i1,i2).c(1,2)*r4_56f(i5,i6).c( & 2,1)) tapptt(i1,i2,i5,i6).d(1,1)=l3_12fz(i1,i2).d(1,1)*p123q*r4_ & 56f(i5,i6).d(1,1)+l3_12fz(i1,i2).b(1,2)*r4_56f(i5,i6).c(2 & ,1) tapptt(i1,i2,i5,i6).a(1,2)=xm(3)*(l3_12fz(i1,i2).a(1,1)*r4 & _56f(i5,i6).b(1,2)+l3_12fz(i1,i2).c(1,2)*r4_56f(i5,i6).a( & 2,2)) tapptt(i1,i2,i5,i6).b(1,2)=l3_12fz(i1,i2).d(1,1)*p123q*r4_ & 56f(i5,i6).b(1,2)+l3_12fz(i1,i2).b(1,2)*r4_56f(i5,i6).a(2 & ,2) tapptt(i1,i2,i5,i6).c(1,2)=l3_12fz(i1,i2).a(1,1)*r4_56f(i5 & ,i6).c(1,2)+l3_12fz(i1,i2).c(1,2)*p123q*r4_56f(i5,i6).d(2 & ,2) tapptt(i1,i2,i5,i6).d(1,2)=xm(3)*(l3_12fz(i1,i2).d(1,1)*r4 & _56f(i5,i6).c(1,2)+l3_12fz(i1,i2).b(1,2)*r4_56f(i5,i6).d( & 2,2)) tapptt(i1,i2,i5,i6).a(2,1)=xm(3)*(l3_12fz(i1,i2).c(2,1)*r4 & _56f(i5,i6).a(1,1)+l3_12fz(i1,i2).a(2,2)*r4_56f(i5,i6).b( & 2,1)) tapptt(i1,i2,i5,i6).b(2,1)=l3_12fz(i1,i2).b(2,1)*r4_56f(i5 & ,i6).a(1,1)+l3_12fz(i1,i2).d(2,2)*p123q*r4_56f(i5,i6).b(2 & ,1) tapptt(i1,i2,i5,i6).c(2,1)=l3_12fz(i1,i2).c(2,1)*p123q*r4_ & 56f(i5,i6).d(1,1)+l3_12fz(i1,i2).a(2,2)*r4_56f(i5,i6).c(2 & ,1) tapptt(i1,i2,i5,i6).d(2,1)=xm(3)*(l3_12fz(i1,i2).b(2,1)*r4 & _56f(i5,i6).d(1,1)+l3_12fz(i1,i2).d(2,2)*r4_56f(i5,i6).c( & 2,1)) tapptt(i1,i2,i5,i6).a(2,2)=l3_12fz(i1,i2).c(2,1)*p123q*r4_ & 56f(i5,i6).b(1,2)+l3_12fz(i1,i2).a(2,2)*r4_56f(i5,i6).a(2 & ,2) tapptt(i1,i2,i5,i6).b(2,2)=xm(3)*(l3_12fz(i1,i2).b(2,1)*r4 & _56f(i5,i6).b(1,2)+l3_12fz(i1,i2).d(2,2)*r4_56f(i5,i6).a( & 2,2)) tapptt(i1,i2,i5,i6).c(2,2)=xm(3)*(l3_12fz(i1,i2).c(2,1)*r4 & _56f(i5,i6).c(1,2)+l3_12fz(i1,i2).a(2,2)*r4_56f(i5,i6).d( & 2,2)) tapptt(i1,i2,i5,i6).d(2,2)=l3_12fz(i1,i2).b(2,1)*r4_56f(i5 & ,i6).c(1,2)+l3_12fz(i1,i2).d(2,2)*p123q*r4_56f(i5,i6).d(2 & ,2) end do end do end do end do do i1=1,2 do i2=1,2 do i5=1,2 do i6=1,2 * mline -- res=cn8(i1,i2,&,i5,i6).id(ide),abcd=tapptt(i1,i2,i5,i6).,m1=xm(3 * ),m2=(-xm(4)),den=((p123q-xm(3)**2)*p123k0*r3r5) do iut=1,2 do jut=1,2 cn8(i1,i2,iut,jut,i5,i6).id(ide)=(tapptt(i1,i2,i5,i6).a(iu & t,jut)+xm(3)*tapptt(i1,i2,i5,i6).b(iut,jut)+(-xm(4))*tapp & tt(i1,i2,i5,i6).c(iut,jut)+xm(3)*(-xm(4))*tapptt(i1,i2,i5 & ,i6).d(iut,jut))/((p123q-xm(3)**2)*p123k0*r3r5) enddo enddo end do end do end do end do **** Diagramma CN9 **** do i1=1,2 do i2=1,2 do i3=1,2 do i4=1,2 * TT -- aa=tapptt(i1,i2,i3,i4).a,bb=tapptt(i1,i2,i3,i4).b,cc=tapptt(i1,i2,i * 3,i4).c,dd=tapptt(i1,i2,i3,i4).d,a1=l5_34f(i3,i4).a,b1=l5_34f(i3,i4).b,c1= * l5_34f(i3,i4).c,d1=l5_34f(i3,i4).d,a2=r6_12fz(i1,i2).a,b2=r6_12fz(i1,i2).b * ,c2=r6_12fz(i1,i2).c,d2=r6_12fz(i1,i2).d,prq=p345q,m=xm(5) tapptt(i1,i2,i3,i4).a(1,1)=l5_34f(i3,i4).a(1,1)*r6_12fz(i1 & ,i2).a(1,1)+l5_34f(i3,i4).c(1,2)*p345q*r6_12fz(i1,i2).b(2 & ,1) tapptt(i1,i2,i3,i4).b(1,1)=xm(5)*(l5_34f(i3,i4).d(1,1)*r6_ & 12fz(i1,i2).a(1,1)+l5_34f(i3,i4).b(1,2)*r6_12fz(i1,i2).b( & 2,1)) tapptt(i1,i2,i3,i4).c(1,1)=xm(5)*(l5_34f(i3,i4).a(1,1)*r6_ & 12fz(i1,i2).d(1,1)+l5_34f(i3,i4).c(1,2)*r6_12fz(i1,i2).c( & 2,1)) tapptt(i1,i2,i3,i4).d(1,1)=l5_34f(i3,i4).d(1,1)*p345q*r6_1 & 2fz(i1,i2).d(1,1)+l5_34f(i3,i4).b(1,2)*r6_12fz(i1,i2).c(2 & ,1) tapptt(i1,i2,i3,i4).a(1,2)=xm(5)*(l5_34f(i3,i4).a(1,1)*r6_ & 12fz(i1,i2).b(1,2)+l5_34f(i3,i4).c(1,2)*r6_12fz(i1,i2).a( & 2,2)) tapptt(i1,i2,i3,i4).b(1,2)=l5_34f(i3,i4).d(1,1)*p345q*r6_1 & 2fz(i1,i2).b(1,2)+l5_34f(i3,i4).b(1,2)*r6_12fz(i1,i2).a(2 & ,2) tapptt(i1,i2,i3,i4).c(1,2)=l5_34f(i3,i4).a(1,1)*r6_12fz(i1 & ,i2).c(1,2)+l5_34f(i3,i4).c(1,2)*p345q*r6_12fz(i1,i2).d(2 & ,2) tapptt(i1,i2,i3,i4).d(1,2)=xm(5)*(l5_34f(i3,i4).d(1,1)*r6_ & 12fz(i1,i2).c(1,2)+l5_34f(i3,i4).b(1,2)*r6_12fz(i1,i2).d( & 2,2)) tapptt(i1,i2,i3,i4).a(2,1)=xm(5)*(l5_34f(i3,i4).c(2,1)*r6_ & 12fz(i1,i2).a(1,1)+l5_34f(i3,i4).a(2,2)*r6_12fz(i1,i2).b( & 2,1)) tapptt(i1,i2,i3,i4).b(2,1)=l5_34f(i3,i4).b(2,1)*r6_12fz(i1 & ,i2).a(1,1)+l5_34f(i3,i4).d(2,2)*p345q*r6_12fz(i1,i2).b(2 & ,1) tapptt(i1,i2,i3,i4).c(2,1)=l5_34f(i3,i4).c(2,1)*p345q*r6_1 & 2fz(i1,i2).d(1,1)+l5_34f(i3,i4).a(2,2)*r6_12fz(i1,i2).c(2 & ,1) tapptt(i1,i2,i3,i4).d(2,1)=xm(5)*(l5_34f(i3,i4).b(2,1)*r6_ & 12fz(i1,i2).d(1,1)+l5_34f(i3,i4).d(2,2)*r6_12fz(i1,i2).c( & 2,1)) tapptt(i1,i2,i3,i4).a(2,2)=l5_34f(i3,i4).c(2,1)*p345q*r6_1 & 2fz(i1,i2).b(1,2)+l5_34f(i3,i4).a(2,2)*r6_12fz(i1,i2).a(2 & ,2) tapptt(i1,i2,i3,i4).b(2,2)=xm(5)*(l5_34f(i3,i4).b(2,1)*r6_ & 12fz(i1,i2).b(1,2)+l5_34f(i3,i4).d(2,2)*r6_12fz(i1,i2).a( & 2,2)) tapptt(i1,i2,i3,i4).c(2,2)=xm(5)*(l5_34f(i3,i4).c(2,1)*r6_ & 12fz(i1,i2).c(1,2)+l5_34f(i3,i4).a(2,2)*r6_12fz(i1,i2).d( & 2,2)) tapptt(i1,i2,i3,i4).d(2,2)=l5_34f(i3,i4).b(2,1)*r6_12fz(i1 & ,i2).c(1,2)+l5_34f(i3,i4).d(2,2)*p345q*r6_12fz(i1,i2).d(2 & ,2) end do end do end do end do do i1=1,2 do i2=1,2 do i3=1,2 do i4=1,2 * mline -- res=cn9(i1,i2,i3,i4,&).id(ide),abcd=tapptt(i1,i2,i3,i4).,m1=xm(5 * ),m2=(-xm(6)),den=((p345q-xm(5)**2)*p345k0*r3r5) do iut=1,2 do jut=1,2 cn9(i1,i2,i3,i4,iut,jut).id(ide)=(tapptt(i1,i2,i3,i4).a(iu & t,jut)+xm(5)*tapptt(i1,i2,i3,i4).b(iut,jut)+(-xm(6))*tapp & tt(i1,i2,i3,i4).c(iut,jut)+xm(5)*(-xm(6))*tapptt(i1,i2,i3 & ,i4).d(iut,jut))/((p345q-xm(5)**2)*p345k0*r3r5) enddo enddo end do end do end do end do **** Diagramma CN10 **** do i1=1,2 do i2=1,2 do i3=1,2 do i4=1,2 * TT -- aa=tapptt(i1,i2,i3,i4).a,bb=tapptt(i1,i2,i3,i4).b,cc=tapptt(i1,i2,i * 3,i4).c,dd=tapptt(i1,i2,i3,i4).d,a1=l5_12fz(i1,i2).a,b1=l5_12fz(i1,i2).b,c * 1=l5_12fz(i1,i2).c,d1=l5_12fz(i1,i2).d,a2=r6_34f(i3,i4).a,b2=r6_34f(i3,i4) * .b,c2=r6_34f(i3,i4).c,d2=r6_34f(i3,i4).d,prq=p125q,m=xm(5) tapptt(i1,i2,i3,i4).a(1,1)=l5_12fz(i1,i2).a(1,1)*r6_34f(i3 & ,i4).a(1,1)+l5_12fz(i1,i2).c(1,2)*p125q*r6_34f(i3,i4).b(2 & ,1) tapptt(i1,i2,i3,i4).b(1,1)=xm(5)*(l5_12fz(i1,i2).d(1,1)*r6 & _34f(i3,i4).a(1,1)+l5_12fz(i1,i2).b(1,2)*r6_34f(i3,i4).b( & 2,1)) tapptt(i1,i2,i3,i4).c(1,1)=xm(5)*(l5_12fz(i1,i2).a(1,1)*r6 & _34f(i3,i4).d(1,1)+l5_12fz(i1,i2).c(1,2)*r6_34f(i3,i4).c( & 2,1)) tapptt(i1,i2,i3,i4).d(1,1)=l5_12fz(i1,i2).d(1,1)*p125q*r6_ & 34f(i3,i4).d(1,1)+l5_12fz(i1,i2).b(1,2)*r6_34f(i3,i4).c(2 & ,1) tapptt(i1,i2,i3,i4).a(1,2)=xm(5)*(l5_12fz(i1,i2).a(1,1)*r6 & _34f(i3,i4).b(1,2)+l5_12fz(i1,i2).c(1,2)*r6_34f(i3,i4).a( & 2,2)) tapptt(i1,i2,i3,i4).b(1,2)=l5_12fz(i1,i2).d(1,1)*p125q*r6_ & 34f(i3,i4).b(1,2)+l5_12fz(i1,i2).b(1,2)*r6_34f(i3,i4).a(2 & ,2) tapptt(i1,i2,i3,i4).c(1,2)=l5_12fz(i1,i2).a(1,1)*r6_34f(i3 & ,i4).c(1,2)+l5_12fz(i1,i2).c(1,2)*p125q*r6_34f(i3,i4).d(2 & ,2) tapptt(i1,i2,i3,i4).d(1,2)=xm(5)*(l5_12fz(i1,i2).d(1,1)*r6 & _34f(i3,i4).c(1,2)+l5_12fz(i1,i2).b(1,2)*r6_34f(i3,i4).d( & 2,2)) tapptt(i1,i2,i3,i4).a(2,1)=xm(5)*(l5_12fz(i1,i2).c(2,1)*r6 & _34f(i3,i4).a(1,1)+l5_12fz(i1,i2).a(2,2)*r6_34f(i3,i4).b( & 2,1)) tapptt(i1,i2,i3,i4).b(2,1)=l5_12fz(i1,i2).b(2,1)*r6_34f(i3 & ,i4).a(1,1)+l5_12fz(i1,i2).d(2,2)*p125q*r6_34f(i3,i4).b(2 & ,1) tapptt(i1,i2,i3,i4).c(2,1)=l5_12fz(i1,i2).c(2,1)*p125q*r6_ & 34f(i3,i4).d(1,1)+l5_12fz(i1,i2).a(2,2)*r6_34f(i3,i4).c(2 & ,1) tapptt(i1,i2,i3,i4).d(2,1)=xm(5)*(l5_12fz(i1,i2).b(2,1)*r6 & _34f(i3,i4).d(1,1)+l5_12fz(i1,i2).d(2,2)*r6_34f(i3,i4).c( & 2,1)) tapptt(i1,i2,i3,i4).a(2,2)=l5_12fz(i1,i2).c(2,1)*p125q*r6_ & 34f(i3,i4).b(1,2)+l5_12fz(i1,i2).a(2,2)*r6_34f(i3,i4).a(2 & ,2) tapptt(i1,i2,i3,i4).b(2,2)=xm(5)*(l5_12fz(i1,i2).b(2,1)*r6 & _34f(i3,i4).b(1,2)+l5_12fz(i1,i2).d(2,2)*r6_34f(i3,i4).a( & 2,2)) tapptt(i1,i2,i3,i4).c(2,2)=xm(5)*(l5_12fz(i1,i2).c(2,1)*r6 & _34f(i3,i4).c(1,2)+l5_12fz(i1,i2).a(2,2)*r6_34f(i3,i4).d( & 2,2)) tapptt(i1,i2,i3,i4).d(2,2)=l5_12fz(i1,i2).b(2,1)*r6_34f(i3 & ,i4).c(1,2)+l5_12fz(i1,i2).d(2,2)*p125q*r6_34f(i3,i4).d(2 & ,2) end do end do end do end do do i1=1,2 do i2=1,2 do i3=1,2 do i4=1,2 * mline -- res=cn10(i1,i2,i3,i4,&).id(ide),abcd=tapptt(i1,i2,i3,i4).,m1=xm( * 5),m2=(-xm(6)),den=((p125q-xm(5)**2)*p125k0*r3r5) do iut=1,2 do jut=1,2 cn10(i1,i2,i3,i4,iut,jut).id(ide)=(tapptt(i1,i2,i3,i4).a(i & ut,jut)+xm(5)*tapptt(i1,i2,i3,i4).b(iut,jut)+(-xm(6))*tap & ptt(i1,i2,i3,i4).c(iut,jut)+xm(5)*(-xm(6))*tapptt(i1,i2,i & 3,i4).d(iut,jut))/((p125q-xm(5)**2)*p125k0*r3r5) enddo enddo end do end do end do end do end do !fine del do per particelle identiche ******************qui mettere le somme dei diagrammi******************** if (ialfasrun.eq.1) then als34=rals(.233d0,rm34,5) als56=rals(.233d0,rm56,5) if (iid.eq.1) then als54=rals(.233d0,rm45,5) als36=rals(.233d0,rm36,5) endif else als34=alfas als56=alfas if (iid.eq.1) then als54=alfas als36=alfas endif endif res=0.d0 rcolfac=2.d0 rcolfacint_sc=-2.d0/3.d0 DO i1=1,2 do i2=1,2 do i3=1,2 do i4=1,2 do i5=1,2 do i6=1,2 cn_qcd(i1,i2,i3,i4,i5,i6)=als56* & (cn7(i1,i2,i3,i4,i5,i6).id(0)+cn8(i1,i2,i3,i4,i5,i6).id(0))+ & als34* & (cn9(i1,i2,i3,i4,i5,i6).id(0)+cn10(i1,i2,i3,i4,i5,i6).id(0)) res=res+rcolfac*(dreal(cn_qcd(i1,i2,i3,i4,i5,i6))**2 & +dimag(cn_qcd(i1,i2,i3,i4,i5,i6))**2) if (iid.eq.1) then cn_qcd_sc(i1,i2,i5,i4,i3,i6)=als36* & (cn7(i1,i2,i3,i4,i5,i6).id(1)+cn8(i1,i2,i3,i4,i5,i6).id(1))+ & als54* & (cn9(i1,i2,i3,i4,i5,i6).id(1)+cn10(i1,i2,i3,i4,i5,i6).id(1)) endif enddo enddo enddo enddo enddo ENDDO !i1 if (iid.eq.1) then DO i1=1,2 do i2=1,2 do i3=1,2 do i4=1,2 do i5=1,2 do i6=1,2 res=res+rcolfac*(dreal(cn_qcd_sc(i1,i2,i3,i4,i5,i6))**2 & +dimag(cn_qcd_sc(i1,i2,i3,i4,i5,i6))**2) res=res-2.d0*rcolfacint_sc*dreal(cn_qcd(i1,i2,i3,i4,i5,i6) & *conjg(cn_qcd_sc(i1,i2,i3,i4,i5,i6))) enddo enddo enddo enddo enddo ENDDO !i1 endif IF (iid.EQ.1) THEN rden=16.d0 ELSE rden=4.d0 ENDIF ee_4q=res/p1k0/p2k0/p3k0/p4k0/p5k0/p6k0/rden RETURN END cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Computes the matrix element for ee-->qq+gg averaged over incoming c polarizations and colors. Couplings (4*pi*a_em*4*pi*a_s)**2 not c included. c c c Factor 1/2 for two identical particles in the final state c is included, therefore THIS AMPLITUDE SQUARED HAS TO BE INTEGRATED c OVER THE FULL FOUR-PARTICLE PHASE SPACE. c c c Authors: A. Ballestrero, E. Maina and S. Moretti. U. of Torino c e_mail:ballestrero,maina,moretti@to.infn.it c Phone +39 11 6527222 (A. B.) c 7203 (E. M.) c c Date : Dec. 2, 1992. c c it=1 --> q is up-like c it=2 --> q is down-like c rm=quark mass c c qin(4,6) are the momenta: c qin(*,1),qin(*,2) e-,e+ c qin(*,3) gluon c qin(*,4) gluon c qin(*,5) quark c qin(*,6) antiquark c c qin(4,*)= energy 1,2,3=x,y,z c c Momentum conservation: c qin(i,1)+qin(i,2)=qin(i,3)+qin(i,4)+qin(i,5)+qin(i,6) i=1,4 c c The Z0 mass and width and the weak mixing angle c are in a parameter statement. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine s_2e2q2g(qin,rm,it,res) c implicit real*8 (a-d,f-h,o,q-t,v-z) implicit complex*16 (e,p,u) dimension uq1(2,2,2),uq2(2,2,2), 1 uq3(2,2,2),uq4(2,2,2),uq5(2,2,2),uq6(2,2,2) dimension p_z(4,2,2),p_g(4,2,2) dimension e_z(2,2,2,2,2),e_g(2,2,2,2,2) dimension p5(2,2,2,0:4),p6(2,2,2,0:4),pe(2,2,2),e(2,2,2,0:4,3) dimension u2(2,2,2,0:4,3,0:4,3), 1 u3(2,2,0:4,3,2),u1(2,2,2,0:4,3) dimension q(4,0:6),am(6),s(0:6,0:6),nd(0:4),nmax(0:4),d1(0:4) dimension nsc(4),ne(0:4),qin(4,6) dimension u(2,2),uu(2,2) parameter (rmz=91.1867d0,gamz=2.4974d0,s2w=0.231030912451068d0, 1 p_null=(0.,0.)) data nd/-1,0,1,-1,-1/,nmax/1,0,2,2,2/ c do i=1,6 do j=1,4 q(j,i)=qin(j,i) end do end do sw=sqrt(s2w) cw=sqrt(1.-s2w) am(1)=0.d0 am(2)=0.d0 am(3)=0.d0 am(4)=0.d0 am(5)=rm am(6)=rm c g_l=-1.d0 g_r=g_l z_l=(-.5d0+1.d0*s2w)/cw/sw z_r=(+1.d0*s2w)/cw/sw if(it.eq.1) then g5_l=2.d0/3.d0 g5_r=g5_l z5_r=(-2.d0/3.d0*s2w)/cw/sw z5_l=(.5d0-2.d0/3.d0*s2w)/cw/sw elseif(it.eq.2) then g5_l=-1.d0/3.d0 g5_r=g5_l z5_r=(1.d0/3.d0*s2w)/cw/sw z5_l=(-.5d0+1.d0/3.d0*s2w)/cw/sw else write(6,*)'it not allowed' end if c qf1=sqrt(q(4,1)+q(3,1)) qf2=sqrt(q(4,2)+q(3,2)) qf3=sqrt(q(4,3)+q(3,3)) qf4=sqrt(q(4,4)+q(3,4)) c uq1(1,2,2)=qf1*(1,0) uq1(2,2,2)=(q(1,1)*(1,0)+q(2,1)*(0,1))/qf1 uq1(2,1,1)=uq1(1,2,2) uq1(1,1,1)=-conjg(uq1(2,2,2)) c uq2(2,1,2)=-qf2*(1,0) uq2(1,1,2)=-(-q(1,2)*(1,0)-q(2,2)*(0,1))/qf2 uq2(1,2,1)=uq2(2,1,2) uq2(2,2,1)=-conjg(uq2(1,1,2)) c uq3(1,2,2)=qf3*(1,0) uq3(2,2,2)=(q(1,3)*(1,0)+q(2,3)*(0,1))/qf3 uq3(2,1,1)=uq3(1,2,2) uq3(1,1,1)=-conjg(uq3(2,2,2)) c uq4(1,2,2)=qf4*(1,0) uq4(2,2,2)=(q(1,4)*(1,0)+q(2,4)*(0,1))/qf4 uq4(2,1,1)=uq4(1,2,2) uq4(1,1,1)=-conjg(uq4(2,2,2)) c p_z(4,2,1)=z_l*(uq2(1,1,2)*uq1(1,1,1)+uq2(2,1,2)*uq1(2,1,1)) p_z(3,2,1)=-z_l*(uq2(1,1,2)*uq1(1,1,1)-uq2(2,1,2)*uq1(2,1,1)) p_z(2,2,1)=-z_l*(-(0,1)*(uq2(1,1,2)*uq1(2,1,1) 1 -uq2(2,1,2)*uq1(1,1,1))) p_z(1,2,1)=-z_l*(uq2(1,1,2)*uq1(2,1,1)+uq2(2,1,2)*uq1(1,1,1)) c p_z(4,1,2)=z_r*(uq2(1,2,1)*uq1(1,2,2)+uq2(2,2,1)*uq1(2,2,2)) p_z(3,1,2)=z_r*(uq2(1,2,1)*uq1(1,2,2)-uq2(2,2,1)*uq1(2,2,2)) p_z(2,1,2)=z_r*(-(0,1)*(uq2(1,2,1)*uq1(2,2,2) 1 -uq2(2,2,1)*uq1(1,2,2))) p_z(1,1,2)=z_r*(uq2(1,2,1)*uq1(2,2,2)+uq2(2,2,1)*uq1(1,2,2)) c p_g(4,2,1)=g_l*(uq2(1,1,2)*uq1(1,1,1)+uq2(2,1,2)*uq1(2,1,1)) p_g(3,2,1)=-g_l*(uq2(1,1,2)*uq1(1,1,1)-uq2(2,1,2)*uq1(2,1,1)) p_g(2,2,1)=-g_l*(-(0,1)*(uq2(1,1,2)*uq1(2,1,1) 1 -uq2(2,1,2)*uq1(1,1,1))) p_g(1,2,1)=-g_l*(uq2(1,1,2)*uq1(2,1,1)+uq2(2,1,2)*uq1(1,1,1)) c p_g(4,1,2)=g_r*(uq2(1,2,1)*uq1(1,2,2)+uq2(2,2,1)*uq1(2,2,2)) p_g(3,1,2)=g_r*(uq2(1,2,1)*uq1(1,2,2)-uq2(2,2,1)*uq1(2,2,2)) p_g(2,1,2)=g_r*(-(0,1)*(uq2(1,2,1)*uq1(2,2,2) 1 -uq2(2,2,1)*uq1(1,2,2))) p_g(1,1,2)=g_r*(uq2(1,2,1)*uq1(2,2,2)+uq2(2,2,1)*uq1(1,2,2)) c sn=q(4,3)*q(4,4) do k=1,3 sn=sn-q(k,3)*q(k,4) end do snorm=1.d0/sqrt(sn) c s_hat=q(4,1)*q(4,2) do k=1,3 s_hat=s_hat-q(k,1)*q(k,2) end do s_hat=2.d0*s_hat up1z=(s_hat-rmz*rmz+rmz*gamz*(0,1)) rp1g=s_hat z_rat=z5_r/z5_l do i=1,2 j=3-i e_z(1,1,1,i,j)=(p_z(4,i,j)+p_z(3,i,j))/up1z*z5_l e_z(2,1,1,i,j)=(p_z(1,i,j)+p_z(2,i,j)*(0,1))/up1z*z5_l e_z(1,2,1,i,j)=(p_z(1,i,j)-p_z(2,i,j)*(0,1))/up1z*z5_l e_z(2,2,1,i,j)=(p_z(4,i,j)-p_z(3,i,j))/up1z*z5_l c e_z(1,1,2,i,j)=e_z(2,2,1,i,j)*z_rat e_z(2,1,2,i,j)=-e_z(2,1,1,i,j)*z_rat e_z(1,2,2,i,j)=-e_z(1,2,1,i,j)*z_rat e_z(2,2,2,i,j)=e_z(1,1,1,i,j)*z_rat c e_g(1,1,1,i,j)=(p_g(4,i,j)+p_g(3,i,j))/rp1g*g5_l e_g(2,1,1,i,j)=(p_g(1,i,j)+p_g(2,i,j)*(0,1))/rp1g*g5_l e_g(1,2,1,i,j)=(p_g(1,i,j)-p_g(2,i,j)*(0,1))/rp1g*g5_l e_g(2,2,1,i,j)=(p_g(4,i,j)-p_g(3,i,j))/rp1g*g5_l c e_g(1,1,2,i,j)=e_g(2,2,1,i,j) e_g(2,1,2,i,j)=-e_g(2,1,1,i,j) e_g(1,2,2,i,j)=-e_g(1,2,1,i,j) e_g(2,2,2,i,j)=e_g(1,1,1,i,j) end do c do i=1,2 j=3-i do l=1,2 do m=1,2 do n=1,2 e(n,m,l,2,j)=e_g(n,m,l,i,j)+e_z(n,m,l,i,j) end do end do end do end do c qm5=sqrt(q(4,5)*q(4,5)-am(5)*am(5)) q5=sqrt(qm5+q(3,5)) f5p=sqrt((q(4,5)/qm5+1.d0)/2.d0) f5m=sqrt((q(4,5)/qm5-1.d0)/2.d0) c uq5(1,1,2)=q5*f5m*(1,0) uq5(2,1,2)=(q(1,5)*(1,0)-q(2,5)*(0,1))/q5*f5m uq5(1,2,2)=q5*f5p*(1,0) uq5(2,2,2)=(q(1,5)*(1,0)-q(2,5)*(0,1))/q5*f5p c uq5(1,1,1)=-conjg(uq5(2,2,2)) uq5(2,1,1)=uq5(1,2,2) uq5(1,2,1)=-conjg(uq5(2,1,2)) uq5(2,2,1)=uq5(1,1,2) c qm6=sqrt(q(4,6)*q(4,6)-am(6)*am(6)) q6=sqrt(qm6+q(3,6)) f6p=sqrt((q(4,6)/qm6+1.d0)/2.d0) f6m=sqrt((q(4,6)/qm6-1.d0)/2.d0) c uq6(1,1,2)=-(-q(1,6)*(1,0)+q(2,6)*(0,1))/q6*f6p uq6(2,1,2)=-q6*f6p*(1,0) uq6(1,2,2)=(-q(1,6)*(1,0)+q(2,6)*(0,1))/q6*f6m uq6(2,2,2)=q6*f6m*(1,0) c uq6(1,1,1)=uq6(2,2,2) uq6(2,1,1)=-conjg(uq6(1,2,2)) uq6(1,2,1)=uq6(2,1,2) uq6(2,2,1)=-conjg(uq6(1,1,2)) c do i=1,4 q(i,2)=q(i,2)+q(i,1) end do c do i=1,4 q(i,0)=q(i,3)+q(i,4) end do c do i=1,6 s(i,i)=am(i)*am(i) end do c sq=q(4,0)*q(4,0) do k=1,3 sq=sq-q(k,0)*q(k,0) end do s(0,0)=sq c sq=q(4,2)*q(4,2) do k=1,3 sq=sq-q(k,2)*q(k,2) end do s(2,2)=sq c do i=0,6 do j=i+1,6 sq=q(4,i)*q(4,j) do k=1,3 sq=sq-q(k,i)*q(k,j) end do s(i,j)=sq s(j,i)=sq end do end do c do i=2,4 p5(1,1,1,i)=(q(4,5)-nd(i)*q(4,i)+(q(3,5)-nd(i)*q(3,i)))*(1,0) p5(2,1,1,i)=(q(1,5)-nd(i)*q(1,i))*(1,0) 1 +(q(2,5)-nd(i)*q(2,i))*(0,1) p5(1,2,1,i)=conjg(p5(2,1,1,i)) p5(2,2,1,i)=(q(4,5)-nd(i)*q(4,i)-(q(3,5)-nd(i)*q(3,i)))*(1,0) c p5(1,1,2,i)=p5(2,2,1,i) p5(2,1,2,i)=-p5(2,1,1,i) p5(1,2,2,i)=-p5(1,2,1,i) p5(2,2,2,i)=p5(1,1,1,i) end do do i=2,4 p6(1,1,1,i)=(-q(4,6)+nd(i)*q(4,i) 1 +(-q(3,6)+nd(i)*q(3,i)))*(1,0) p6(2,1,1,i)=(-q(1,6)+nd(i)*q(1,i))*(1,0) 1 +(-q(2,6)+nd(i)*q(2,i))*(0,1) p6(1,2,1,i)=conjg(p6(2,1,1,i)) p6(2,2,1,i)=(-q(4,6)+nd(i)*q(4,i) 1 -(-q(3,6)+nd(i)*q(3,i)))*(1,0) c p6(1,1,2,i)=p6(2,2,1,i) p6(2,1,2,i)=-p6(2,1,1,i) p6(1,2,2,i)=-p6(1,2,1,i) p6(2,2,2,i)=p6(1,1,1,i) end do do k=1,2 do i=1,2 do j=1,2 p5(j,i,k,0)=p6(j,i,k,2) end do end do end do c pe(1,1,1)=-(q(4,3)-q(4,4)+(q(3,3)-q(3,4)))*(1,0) pe(2,1,1)=-((q(1,3)-q(1,4))*(1,0)+(q(2,3)-q(2,4))*(0,1)) pe(1,2,1)=conjg(pe(2,1,1)) pe(2,2,1)=-(q(4,3)-q(4,4)-(q(3,3)-q(3,4)))*(1,0) c pe(1,1,2)=pe(2,2,1) pe(1,2,2)=-pe(1,2,1) pe(2,1,2)=-pe(2,1,1) pe(2,2,2)=pe(1,1,1) c do j1=1,2 do j2=1,2 e(j1,j2,1,3,2)=uq4(j1,2,2)*conjg(uq3(j2,2,2))*snorm end do end do do j1=1,2 do j2=1,2 e(j1,j2,2,3,2)=uq3(j1,1,1)*conjg(uq4(j2,1,1))*snorm end do end do do j1=1,2 do j2=1,2 e(j1,j2,2,3,1)=uq4(j1,1,1)*conjg(uq3(j2,1,1))*snorm end do end do do j1=1,2 do j2=1,2 e(j1,j2,1,3,1)=uq3(j1,2,2)*conjg(uq4(j2,2,2))*snorm end do end do c do j1=1,2 do j2=1,2 e(j1,j2,1,4,2)=uq3(j1,2,2)*conjg(uq4(j2,2,2))*snorm end do end do do j1=1,2 do j2=1,2 e(j1,j2,2,4,2)=uq4(j1,1,1)*conjg(uq3(j2,1,1))*snorm end do end do do j1=1,2 do j2=1,2 e(j1,j2,2,4,1)=uq3(j1,1,1)*conjg(uq4(j2,1,1))*snorm end do end do do j1=1,2 do j2=1,2 e(j1,j2,1,4,1)=uq4(j1,2,2)*conjg(uq3(j2,2,2))*snorm end do end do c do k=1,2 do i=1,2 do j=1,2 e(j,i,k,0,1)=-pe(j,i,k)/(2.d0*s(3,4)) end do end do end do c nsc(1)=0 nsc(2)=2 nsc(3)=3 nsc(4)=4 c do isc=1,4 i=nsc(isc) d1(i)=-2.d0*nd(i)*s(5,i)+s(i,i) amm=am(5) do iu=1,2 do j=1,nmax(i) do k=1,2 nc=3-k u(1,k)=uq5(1,nc,iu)*e(1,1,nc,i,j) 1 +uq5(2,nc,iu)*e(2,1,nc,i,j) u(2,k)=uq5(1,nc,iu)*e(1,2,nc,i,j) 1 +uq5(2,nc,iu)*e(2,2,nc,i,j) end do do k=1,2 nc=3-k u1(1,k,iu,i,j)=(u(1,nc)*p5(1,1,nc,i)+u(2,nc)*p5(2,1,nc,i) 1 +u(1,k)*amm)/d1(i) u1(2,k,iu,i,j)=(u(1,nc)*p5(1,2,nc,i)+u(2,nc)*p5(2,2,nc,i) 1 +u(2,k)*amm)/d1(i) end do end do end do end do c do i=2,4 do i1=2,4 if (i1.eq.i) go to 100 d2=d1(i)-2.d0*nd(i1)*s(5,i1)+s(i1,i1) 1 +2.d0*nd(i)*nd(i1)*s(i,i1) amm=am(5) do j=1,nmax(i) do j1=1,nmax(i1) do iu=1,2 do k=1,2 nc=3-k u(1,k)=u1(1,nc,iu,i,j)*e(1,1,nc,i1,j1) 1 +u1(2,nc,iu,i,j)*e(2,1,nc,i1,j1) u(2,k)=u1(1,nc,iu,i,j)*e(1,2,nc,i1,j1) 1 +u1(2,nc,iu,i,j)*e(2,2,nc,i1,j1) end do nres=9-i-i1 do k=1,2 nc=3-k u2(1,k,iu,i,j,i1,j1)=(u(1,nc)*p6(1,1,nc,nres) 1 +u(2,nc)*p6(2,1,nc,nres)+u(1,k)*amm)/d2 u2(2,k,iu,i,j,i1,j1)=(u(1,nc)*p6(1,2,nc,nres) 1 +u(2,nc)*p6(2,2,nc,nres)+u(2,k)*amm)/d2 end do end do end do end do 100 continue end do end do c do isc=1,4 i=nsc(isc) do j=1,nmax(i) do j1=1,2 do k=1,2 nc=3-k u3(1,k,i,j,j1)=e(1,1,nc,i,j)*uq6(1,nc,j1) 1 +e(1,2,nc,i,j)*uq6(2,nc,j1) u3(2,k,i,j,j1)=e(2,1,nc,i,j)*uq6(1,nc,j1) 1 +e(2,2,nc,i,j)*uq6(2,nc,j1) end do end do end do end do rsum=0.d0 rdif=0.d0 ne(0)=1 do i2=1,nmax(2) ne(2)=i2 do i5=1,2 do i6=1,2 call gg3(ne,i5,i6,u1,u3,u_3g) do i3=1,nmax(3) ne(3)=i3 do i4=1,nmax(4) ne(4)=i4 call gg(ne,i5,i6,u2,u3,u_sum,u_dif) dsum=abs(u_sum) if(i3.eq.i4) u_dif=u_dif+u_3g ddif=abs(u_dif) rs=28.d0/3.d0*dsum*dsum rd=12.d0*ddif*ddif rsum=rsum+rs rdif=rdif+rd end do end do end do end do end do res=rsum+rdif c res=res/8.d0 return end subroutine gg(ne,ichi1,ichi2,u2,u3,u_sum,u_dif) implicit real*8 (a-d,f-h,o,q-t,v-z) implicit complex*16 (e,p,u) dimension np(3,6),ne(0:4),usum(2) dimension u2(2,2,2,0:4,3,0:4,3),u3(2,2,0:4,3,2) data np/2,3,4, 2 3,4,2, 3 3,2,4, 4 2,4,3, 5 4,3,2, 6 4,2,3/ c ind=0 do k1=1,2 usum(k1)=(0.d0,0.d0) do i=1,3 ind=ind+1 n1=np(1,ind) n2=np(2,ind) n3=np(3,ind) us=u2(1,1,ichi1,n1,ne(n1),n2,ne(n2)) 1 *u3(1,2,n3,ne(n3),ichi2) 1 +u2(2,1,ichi1,n1,ne(n1),n2,ne(n2)) 1 *u3(2,2,n3,ne(n3),ichi2) 1 +u2(1,2,ichi1,n1,ne(n1),n2,ne(n2)) 1 *u3(1,1,n3,ne(n3),ichi2) 1 +u2(2,2,ichi1,n1,ne(n1),n2,ne(n2)) 1 *u3(2,1,n3,ne(n3),ichi2) usum(k1)=usum(k1)+us end do end do u_sum=(usum(1)+usum(2))/2.d0 u_dif=(usum(1)-usum(2))/2.d0 return end subroutine gg3(ne,ichi1,ichi2,u1,u3,u_3g) implicit real*8 (a-d,f-h,o,q-t,v-z) implicit complex*16 (e,p,u) dimension np3(2,2),ne(0:4) dimension u3(2,2,0:4,3,2),u1(2,2,2,0:4,3) data np3/0,2, 2 2,0/ c u_3g=(0.d0,0.d0) do i=1,2 n1=np3(1,i) n2=np3(2,i) us=u1(1,1,ichi1,n1,ne(n1))*u3(1,2,n2,ne(n2),ichi2) 1 +u1(2,1,ichi1,n1,ne(n1))*u3(2,2,n2,ne(n2),ichi2) 1 +u1(1,2,ichi1,n1,ne(n1))*u3(1,1,n2,ne(n2),ichi2) 1 +u1(2,2,ichi1,n1,ne(n1))*u3(2,1,n2,ne(n2),ichi2) u_3g=u_3g+us end do return end real*8 function rals(qcdl,rs,nf) implicit real*8(a-h,o-z) * data pi/3.141592653589793238462643d0/,bqm/4.7/ * * questa riga si puo' modificare o togliere. Comunque non ha senso andare * a scale molto inferiori a 2 GeV con 4 nfe, ed il programma si pianta per * rs dell'ordine di qcdl. if (rs.lt.2.d0) rs=2.d0 * if(rs.lt.bqm) then ex1= 2.d0/25.d0 ex2= 963.d0/14375.d0 rat= bqm/qcdl qcdlp= qcdl*rat**ex1*(2.d0*log(rat))**ex2 nfe= nf-1 else qcdlp= qcdl nfe= nf endif * qcdb0= 11.d0-2.d0/3.d0*nfe qcdb1= 102.d0-38.d0/3.d0*nfe qcdb2= 0.5d0*(2857.d0+nfe*(-5033.d0/9.d0+325.d0/27.d0*nfe)) qcda= 2.d0*log(rs/qcdlp) qcdal= log(qcda) rals= 4.d0*pi/qcdb0/qcda*(1.d0-qcdb1/qcdb0**2/qcda*qcdal+ # (qcdb1/qcdb0**2/qcda)**2*( (qcdal-0.5d0)**2+ # qcdb2*qcdb0/qcdb1**2-5.d0/4.d0)) * return end