program CHIGEN implicit double precision(a-y) implicit complex*16(z) double precision q(4,20) double precision pboo(4),pcm(4),plb(4) double precision rho1chi(3),rho2chi(5) integer r,l,h,i,j,k,m,iset,nhist,ntotal,eflag,ichi,icut,ncut, &id(20),id0,id1,id2,pdf,num,match,mode,iord,iap,d,p,o,pol character prefix*50,fsp*10,beam*10 common/mom/q common/vars/s,rts,mchi,yx common/cuts/etaelmax,etaelmin,ptelmin,ptphmin,ptchimin,ptpsimin &,ptpsimax common/flags/pdf,eflag common/alph/alfas common/match/mm,nn ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c CHIGEN Monte Carlo generator for inclusive colour c c singlet chi_(cJ) --> J/\psi + gamma production c c c c p(1) p(2) --> chi(5) --> gamma(6) + psi(7) c c c c --> gamma(6) + el+(8) + el-(9) + X c c c c Momenta for each event in array q(i,j), where j is c c the particle label and i is the 4-momentum c c component, with: c c c c 1,2 = transverse components c c 3 = beam component c c 4 = energy c c c c PDG number for ith particle in arrary ID(i) c c c c This version: 18 May 2011 c c c c Comments etc to: lh367@cam.ac.uk c c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c physics parameters rts=7d3 ! centre of mass energy mode=200 ! parton choice, see below kfac=2d0 ! k-factor mpsi=3.096916d0 ! pdg 2010 value mchi0=3.41475d0 ! pdg 2010 value mchi1=3.51066d0 ! pdg 2010 value mchi2=3.55620d0 ! pdg 2010 value mmu=0.105658d0 ! pdg 2010 value - el = muon as default ID0=10441 ! pdg numbers ID1=20443 ID2=445 ID(6)=22 ID(7)=443 ID(8)=-13 ID(9)=13 brat0=0.0116d0 ! chi0 --> psi gamma pdg 2010 value brat1=0.3440d0 ! chi1 --> psi gamma pdg 2010 value brat2=0.1950d0 ! chi2 --> psi gamma pdg 2010 value bratj=0.0593d0 ! psi --> el+ el- pdg 2010 value c other parameters ebeam=rts/2d0 s=rts*rts zi=(0d0,1d0) rt2=dsqrt(2d0) pi=dacos(-1d0) bp=rts/dsqrt(2.0d0) c initialise counters etc ntotal=1000000 nhist=6 sum=0d0 sum1=0d0 sum2=0d0 ncut=0 num=1 c initialise histograms do i=1,nhist call histo3(i) enddo c select collider iap=1 if(iap.eq.0)then beam='proton-proton' elseif(iap.eq.1)then beam='proton-antiproton' endif c select partons iset=0 ! no error sets if (mode.eq.200) then ! MSTW 2008 LO prefix = "mstw2008lo" else if (mode.eq.201) then ! MSTW 2008 NLO prefix = "mstw2008nlo" else if (mode.eq.202) then ! MSTW 2008 NNLO prefix = "mstw2008nnlo" else write(6,*) "Error: mode = ",mode stop end if c -- initialise alphas iord=0 FR2 = 1.D0 MUR = 1.D0 ASMUR = 0.5D0 MC = 1.4D0 MB = 4.75D0 MT = 1.D10 CALL INITALPHAS(IORD, FR2, MUR, ASMUR, MC, MB, MT) c chi/psi polarisations on (pol=1) or off (pol=2) pol=2 c (independent) polarization factors xsi1=1d0 ! chi1(long.)/chi1(tot) xsi2a=0d0 ! chi2(2+/-)/chi2(tot) xsi2b=1d0 ! chi2(1+/-)/chi2(tot) ccc note: fractions must add up to 1! if(xsi1.gt.1d0.or.(xsi2a+xsi2b).gt.1d0)then print*,'Error: Polarisation fractions greater than unity!' stop endif c cuts: applied in subroutine 'icut' etaelmax=5d0 ! max eta of leptons, photon etaelmin=2d0 ! min eta of leptons, photon ptelmin=0d0 ! min pt of leptons ptphmin=0d0 ! min pt of photons ptchimin=0d0 ! min pt of chi ptpsimin=0d0 ! min pt of psi ptpsimax=5d5 ! max pt of psi ptmax=20d0**2 ! max ptsq of chi ptmin=0d0**2 ! min ptsq of chi c 2-->1 to 2-->2 matching parameters match=2 if(match.eq.1)then mm=3.9d0**2 nn=0.2d0 klo=1d0 bb=0.18d0 elseif(match.eq.2)then mm=3.2d0**2 nn=0.5d0 klo=1d0 bb=0.3d0 elseif(match.eq.3)then mm=2.5d0**2 nn=0.8d0 klo=1d0 bb=0.4d0 elseif(match.eq.4)then mm=2d0**2 nn=1.5d0 klo=1d0 bb=0.5d0 endif c chi_c0 total width (normalisation) gam0=10.5d-3 !PDG 2009 value c set the desired chi_c (=0,1,2) ichi=1 if(ichi.eq.0) then mchi=mchi0 brat=brat0 fsp='chi0' elseif(ichi.eq.1) then mchi=mchi1 brat=brat1 fsp='chi1' elseif(ichi.eq.2) then mchi=mchi2 brat=brat2 fsp='chi2' endif ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Start of event loop c c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do i=1,ntotal weight=0d0 call r2455(ran0) call r2455(ran1) call r2455(ran2) call r2455(ran3) call r2455(ran4) call r2455(ran5) call r2455(ran6) ccccccccccccccccccccccccccccccccccccccccccccccccccc c incoming protons ID(1)=2212 q(1,1)=0d0 q(2,1)=0d0 q(3,1)=ebeam q(4,1)=ebeam ID(2)=2212 q(1,2)=0d0 q(2,2)=0d0 q(3,2)=-ebeam q(4,2)=ebeam cccccccccccccccccccccccccccccccccccccccccccccccccccc c final state phase space call r2455(ranx1) call r2455(ranx2) call r2455(ranx3) call r2455(ranx4) phix=2d0*pi*ranx1 ptxsq=ptmin+(ptmax-ptmin)*ranx2 ptxx=dsqrt(ptxsq)*dcos(phix) ptxy=dsqrt(ptxsq)*dsin(phix) ccccccccccccccccc ccc chi1 smearing ccccccccccccccccccc call r2455(rphiq) call r2455(rptq) phiq=2d0*pi*rphiq qtsq=mm*rptq qtx=dsqrt(qtsq)*dcos(phiq) qty=dsqrt(qtsq)*dsin(phiq) ptx=ptxx+qtx pty=ptxy+qty ptxsq=ptx**2+pty**2 if(dsqrt(ptxsq).lt.1d-30)then weight=0d0 goto 700 endif rmx=dsqrt(ptxsq+mchi**2) ! m_\perp ymax=dlog(rts/rmx) ! chi rapidity ymin=-ymax yx=ymin+(ymax-ymin)*ranx3 ygmax=(dlog((rts-rmx*dexp(yx))/dsqrt(ptxsq))) ygmin=-(dlog((rts-rmx*dexp(-yx))/dsqrt(ptxsq))) if(ygmax.lt.ygmin)then weight=0d0 goto 700 endif yg=ygmin+(ygmax-ygmin)*ranx4 c nlo (2->2) momentum fractions x1c1=(dsqrt(ptxsq)*dexp(yg)+rmx*dexp(yx))/rts x2c1=(dsqrt(ptxsq)*dexp(-yg)+rmx*dexp(-yx))/rts weightloc1=2d0*pi*(ptmax-ptmin) weightloc1=weightloc1*mm*dexp(-bb*qtsq)*bb weightloc1=weightloc1*(ymax-ymin)*(ygmax-ygmin) c generate fusing gluon momenta q(1,10)=0d0 q(2,10)=0d0 q(3,10)=x1c1*rts/2d0 q(4,10)=x1c1*rts/2d0 q(1,11)=0d0 q(2,11)=0d0 q(3,11)=-x2c1*rts/2d0 q(4,11)=x2c1*rts/2d0 c chi momentum q(1,5)=ptx q(2,5)=pty q(3,5)=rmx*(dexp(yx)-dexp(-yx))/2d0 q(4,5)=rmx*(dexp(yx)+dexp(-yx))/2d0 do k=1,4 q(k,12)=q(k,10)+q(k,11)-q(k,5) enddo c Mandelstam variables thc1=-2d0*(q(4,12)*q(4,10)-q(3,12)*q(3,10)-q(2,12)*q(2,10)- & q(1,12)*q(1,10)) uhc1=-2d0*(q(4,12)*q(4,11)-q(3,12)*q(3,11)-q(2,12)*q(2,11)- & q(1,12)*q(1,11)) shc1=2d0*(q(4,11)*q(4,10)-q(3,11)*q(3,10)-q(2,11)*q(2,10)- & q(1,11)*q(1,10)) ccccccccccccccccccccccccc ptxsq=ptxx**2+ptxy**2 if(dsqrt(ptxsq).lt.1d-30)then weight=0d0 goto 700 endif rmx=dsqrt(ptxsq+mchi**2) ! m_\perp ymax=dlog(rts/rmx) ! chi rapidity ymin=-ymax yx=ymin+(ymax-ymin)*ranx3 ygmax=(dlog((rts-rmx*dexp(yx))/dsqrt(ptxsq))) ygmin=-(dlog((rts-rmx*dexp(-yx))/dsqrt(ptxsq))) if(ygmax.lt.ygmin)then weight=0d0 goto 700 endif yg=ygmin+(ygmax-ygmin)*ranx4 c nlo (2->2) momentum fractions x1=(dsqrt(ptxsq)*dexp(yg)+rmx*dexp(yx))/rts x2=(dsqrt(ptxsq)*dexp(-yg)+rmx*dexp(-yx))/rts c lo (2->1) momentum fractions x11=dexp(yx)*rmx/rts x22=dexp(-yx)*rmx/rts weightnlo=2d0*pi*(ptmax-ptmin) weightnlo=weightnlo*(ymax-ymin)*(ygmax-ygmin) weightlo=dexp(-bb*ptxsq)*bb weightlo=weightlo*(ptmax-ptmin)*(ymax-ymin) c generate fusing gluon momenta q(1,10)=0d0 q(2,10)=0d0 q(3,10)=x1*rts/2d0 q(4,10)=x1*rts/2d0 q(1,11)=0d0 q(2,11)=0d0 q(3,11)=-x2*rts/2d0 q(4,11)=x2*rts/2d0 c chi momentum q(1,5)=ptxx q(2,5)=ptxy q(3,5)=rmx*(dexp(yx)-dexp(-yx))/2d0 q(4,5)=rmx*(dexp(yx)+dexp(-yx))/2d0 do k=1,4 q(k,12)=q(k,10)+q(k,11)-q(k,5) enddo c Mandelstam variables th=-2d0*(q(4,12)*q(4,10)-q(3,12)*q(3,10)-q(2,12)*q(2,10)- & q(1,12)*q(1,10)) uh=-2d0*(q(4,12)*q(4,11)-q(3,12)*q(3,11)-q(2,12)*q(2,11)- & q(1,12)*q(1,11)) sh=2d0*(q(4,11)*q(4,10)-q(3,11)*q(3,10)-q(2,11)*q(2,10)- & q(1,11)*q(1,10)) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c chi decays cccccccccccccccccccccccccccccccccccccccccccccccccccccccc c -- decay 5(mchi) --> 6(0) + 7(mpsi) call r2455(ran7) call r2455(ran8) phi6=2d0*pi*ran7 cost6=2d0*ran8-1d0 sint6=dsqrt(1d0-cost6**2) cphi6=dcos(phi6) sphi6=dsin(phi6) pmod=(mchi**2-mpsi**2)/(2d0*mchi) pcm(4)=pmod pcm(1)=pmod*sint6*sphi6 pcm(2)=pmod*sint6*cphi6 pcm(3)=pmod*cost6 do k=1,4 pboo(k)=q(k,5) enddo call boost(mchi,pboo,pcm,plb) do k=1,4 q(k,6)=plb(k) q(k,7)=q(k,5)-q(k,6) enddo do k=1,4 pcm(k)=q(k,6) enddo do k=1,3 pboo(k)=-q(k,7) enddo pboo(4)=q(4,7) call boost(mpsi,pboo,pcm,plb) c photon in psi r.f. do k=1,4 q(k,16)=plb(k) enddo cccccccccccccccccccccccccccccccccccccccccccccccccccc c -- decay 7(mpsi) --> 8+(mmu) + 9-(mmu) call r2455(ran9) call r2455(ran10) cost8=2d0*ran9-1d0 phi8=2d0*pi*ran10 sint8=dsqrt(1d0-cost8**2) cphi8=dcos(phi8) sphi8=dsin(phi8) pcm(4)=mpsi/2d0 pcmod=dsqrt(pcm(4)**2-mmu**2) pcm(1)=pcmod*sint8*sphi8 pcm(2)=pcmod*sint8*cphi8 pcm(3)=pcmod*cost8 do k=1,4 pboo(k)=q(k,7) enddo call boost(mpsi,pboo,pcm,plb) do k=1,4 q(k,8)=plb(k) q(k,9)=q(k,7)-q(k,8) enddo do k=1,3 pboo(k)=-q(k,7) enddo pboo(4)=q(4,7) do k=1,4 pcm(k)=q(k,8) enddo call boost(mpsi,pboo,pcm,plb) c lepton in psi r.f. do k=1,4 q(k,18)=plb(k) enddo ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c place cuts call cut(icut) if(icut.eq.0)then weight=0d0 ncut=ncut+1 goto 700 endif ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c polarization effects cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do k=1,3 pboo(k)=-q(k,5) enddo pboo(4)=q(4,5) do k=1,4 pcm(k)=q(k,7) enddo call boost(mchi,pboo,pcm,plb) c psi in chi r.f. (dir. of helicity axis) do k=1,4 q(k,17)=plb(k) enddo do k=1,4 pcm(k)=q(k,6) enddo call boost(mchi,pboo,pcm,plb) c photon in chi r.f. (dir. of helicity axis) do k=1,4 q(k,15)=plb(k) enddo cost=(q(1,18)*q(1,16)+q(2,18)*q(2,16)+q(3,18)*q(3,16))/ &dsqrt((q(1,18)**2+q(2,18)**2+q(3,18)**2)* &(q(1,16)**2+q(2,16)**2+q(3,16)**2)) costa=(q(1,17)*q(1,5)+q(2,17)*q(2,5)+q(3,17)*q(3,5))/ &dsqrt((q(1,17)**2+q(2,17)**2+q(3,17)**2)* &(q(1,5)**2+q(2,5)**2+q(3,5)**2)) if(pol.eq.1)then cccccccccccccccccccccccccccccccccccccccccccccc c chi_0 decay cccccccccccccccccccccccccccccccccccccccccccccc if(ichi.eq.0)then weightpol=(1d0+cost**2)*3d0/4d0 endif cccccccccccccccccccccccccccccccccccccccccccccc c chi_1 decay ccccccccccccccccccccccccccccccccccccccccccccc if(ichi.eq.1)then rho1chi(1)=xsi1 rho1chi(2)=(1d0-rho1chi(1))/2d0 rho1chi(3)=rho1chi(2) weightm0=0.5d0-0.5d0*cost**2+cost**2*costa**2 weightmp=0.5d0-0.5d0*cost**2*costa**2 weightpol=rho1chi(1)*weightm0+2d0*rho1chi(2)*weightmp weightpol=weightpol*9d0/4d0 endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccc c chi2 decay ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(ichi.eq.2)then rho2chi(1)=xsi2a/2d0 rho2chi(2)=xsi2b/2d0 rho2chi(3)=1d0-xsi2a-xsi2b rho2chi(4)=rho2chi(2) rho2chi(5)=rho2chi(1) weightm0=0.25d0+0.3d0*costa**2-0.45d0*costa**4+0.25d0*cost**2 &-1.5d0*cost**2*costa**2+1.35d0*cost**2*costa**4 weightm1=0.3d0-0.3d0*costa**2+0.3d0*costa**4+ &0.6d0*cost**2*costa**2-0.9d0*cost**2*costa**4 weightm2=0.225d0+0.15d0*costa**2-0.075d0*costa**4-0.075d0*cost**2 &+0.15d0*cost**2*costa**2+0.225d0*cost**2*costa**4 weightpol=2d0*rho2chi(1)*weightm2+2d0*rho2chi(2)*weightm1+ &rho2chi(3)*weightm0 weightpol=weightpol*15d0/4d0 endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc else weightpol=1d0 endif ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Weight evaluation ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc mu=rmx/2d0 ! hard scale alfas=alphas(mu) c PDFs... glu1lo= GetOnePDF(prefix,iset,x11,mu,0) glu2lo= GetOnePDF(prefix,iset,x22,mu,0) call GetAllPDFs(prefix,iset,x1,mu,upv1,dnv1,upb1,dnb1,str1,sbar1, & chm1,cbar1,bot1,bbar1,glu1nlo,phot) up1=upv1+upb1 dn1=dnv1+dnb1 if(iap.eq.0)then call GetAllPDFs(prefix,iset,x2,mu,upv2,dnv2,upb2,dnb2,str2,sbar2, & chm2,cbar2,bot2,bbar2,glu2nlo,phot) up2=upv2+upb2 dn2=dnv2+dnb2 elseif(iap.eq.1)then call GetAllPDFs(prefix,iset,x2,mu,upv2,dnv2,up2,dn2,sbar2,str2, & cbar2,chm2,bbar2,bot2,glu2nlo,phot) upb2=upv2+up2 dnb2=dnv2+dn2 endif q1=dn1+up1+str1+chm1+bot1 q1bar=dnb1+upb1+sbar1+cbar1+bbar1 q2=dn2+up2+str2+chm2+bot2 q2bar=dnb2+upb2+sbar2+cbar2+bbar2 ccccccccccccccccccc call GetAllPDFs(prefix,iset,x1c1,mu,upv1,dnv1,upb1,dnb1,str1 &,sbar1,chm1,cbar1,bot1,bbar1,glu1nloc1,phot) up1=upv1+upb1 dn1=dnv1+dnb1 if(iap.eq.0)then call GetAllPDFs(prefix,iset,x2c1,mu,upv2,dnv2,upb2,dnb2,str2 &,sbar2,chm2,cbar2,bot2,bbar2,glu2nloc1,phot) up2=upv2+upb2 dn2=dnv2+dnb2 elseif(iap.eq.1)then call GetAllPDFs(prefix,iset,x2c1,mu,upv2,dnv2,up2,dn2,sbar2,str2, & cbar2,chm2,bbar2,bot2,glu2nloc1,phot) upb2=upv2+up2 dnb2=dnv2+dn2 endif q1c1=dn1+up1+str1+chm1+bot1 q1barc1=dnb1+upb1+sbar1+cbar1+bbar1 q2c1=dn2+up2+str2+chm2+bot2 q2barc1=dnb2+upb2+sbar2+cbar2+bbar2 cccccccccccccccccccccccccccccccccccccccccccccccccc c subroutines for (2->2) subprocess cross sections if(ichi.eq.0)then call chi0nlo(sh,th,uh,out) call chi0nloq(sh,uh,th,outqg) call chi0nloq(sh,th,uh,outgq) call chi0nloq(uh,sh,th,outqbq) call chi0nloq(th,sh,uh,outqqb) elseif(ichi.eq.1)then call chi1nlo(sh,th,uh,out) call chi1nloq(sh,uh,th,outqg) call chi1nloq(sh,th,uh,outgq) call chi1nloq(uh,sh,th,outqbq) call chi1nloq(th,sh,uh,outqqb) call chi1nlo(shc1,thc1,uhc1,outc1) call chi1nloq(shc1,uhc1,thc1,outqgc1) call chi1nloq(shc1,thc1,uhc1,outgqc1) call chi1nloq(uhc1,shc1,thc1,outqbqc1) call chi1nloq(thc1,shc1,uhc1,outqqbc1) elseif(ichi.eq.2)then call chi2nlo(sh,th,uh,out) call chi2nloq(sh,uh,th,outqg) call chi2nloq(sh,th,uh,outgq) call chi2nloq(uh,sh,th,outqbq) call chi2nloq(th,sh,uh,outqqb) endif weightgg=out*glu1nlo*glu2nlo/(x1*x2) weightqg=outqg*glu1nlo*(q2+q2bar) & /(x1*x2) weightgq=outgq*glu2nlo*(q1+q1bar) & /(x1*x2) weightqbqp=outqbq*q2bar*q1 weightqbqp=weightqbqp*(-8d0*uh**2/(3d0*sh**2)) weightqqb=outqqb*q1bar*q2 weightqqb=weightqqbp*(-8d0*th**2/(3d0*sh**2)) weightnlo=weightnlo*(weightgg+weightqg+weightgq+ & weightqbq+weightqqb) weightnlo=weightnlo*x1*x2/pi cccccccccccc weightgg=outc1*glu1nloc1*glu2nloc1/(x1c1*x2c1) weightqg=outqgc1*glu1nloc1*(q2c1+q2barc1) & /(x1c1*x2c1) weightgq=outgqc1*glu2nloc1*(q1c1+q1barc1) & /(x1c1*x2c1) weightqbqp=outqbqc1*q2barc1*q1c1 weightqbqp=weightqbqp*(-8d0*uhc1**2/(3d0*shc1**2)) weightqqb=outqqbc1*q1barc1*q2c1 weightqqb=weightqqbp*(-8d0*thc1**2/(3d0*shc1**2)) weightloc1=weightloc1*(weightgg+weightqg+weightgq+ & weightqbq+weightqqb) weightloc1=weightloc1*x1c1*x2c1/pi cccccccccccccccccccccccccccccccccccccccccccccccccc c (2->1) subprocess cross sections rp=0.006d0*mchi**2 gam0=96d0*alfas**2*rp/mchi0**4 if(ichi.eq.0)then weightlo=weightlo*pi**2*gam0/(8d0*mchi*s) elseif(ichi.eq.2)then weightlo=weightlo*pi**2*5d0/(8d0*mchi*s) weightlo=weightlo*128d0*alfas**2*rp/(5d0*mchi**4) endif weightlo=weightlo*glu1lo*glu2lo/(x11*x22)*klo ccccccccccccccccccccccccccccccccccccccccccccccccc if(ichi.eq.1)then weighttot=fw(ptxsq)*weightnlo+(1d0-fw(ptxsq))*weightloc1 else ! matching between lo and nlo weighttot=fw(ptxsq)*weightnlo+(1d0-fw(ptxsq))*weightlo endif weight=brat*bratj*389.379d3 weight=weight*weighttot weight=weight*kfac weight=weight*weightpol 666 wthist=weight/dfloat(ntotal) call binit(wthist) 700 sum=sum+weight sum1=sum1+weight**2 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c End of event loop c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc enddo sum=sum/dfloat(ntotal) sum1=sum1/dfloat(ntotal) var=dsqrt((sum1-sum**2)/dfloat(ntotal)) EFF=1d0*(ntotal-ncut)/ntotal c create histograms DO 6 I=1,NHIST IF(I.EQ.1) WRITE(6,101) IF(I.EQ.2) WRITE(6,102) IF(I.EQ.3) WRITE(6,103) IF(I.EQ.4) WRITE(6,104) IF(I.EQ.5) WRITE(6,105) IF(I.EQ.6) WRITE(6,106) IF(I.EQ.7) WRITE(6,107) CALL HISTO2(I,0) 6 CONTINUE 101 FORMAT(//,' CHI TRANSVERSE MOMENTUM ',/) 102 FORMAT(//,' CHI RAPIDITY ',/) 103 FORMAT(//,' LEPTON TRANSVERSE MOMENTUM ',/) 104 FORMAT(//,' LEPTON RAPIDITY ',/) 105 FORMAT(//,' PHOTON TRANSVERSE MOMENTUM ',/) 106 FORMAT(//,' PHOTON RAPIDITY ',/) 107 FORMAT(//,' ',/) ccccccc print information to screen c print*,sum,var,sum2 WRITE(6,10)sum,var write(6,15)etaelmax,etaelmin,ptelmin WRITE(6,12)ntotal,ncut,eff 10 FORMAT(//,' SIGMA = ',G16.7,' +/-',G16.7,' NB ') 12 FORMAT(/,' [IN,OUT,EFF] : ',2I8,F7.3) 15 FORMAT(/,' [ETAELMAX,ETAELMIN,PTELMIN] : ',3F7.2) WRITE(6,11)fsp,rts/1D3,MCHI,BRAT,prefix 11 FORMAT(/,' Input parameters :',/ . ,' final state : ',A10,/ . ,' sqrt(S) : ',F8.4,' TeV',/ . ,' chi mass : ',F8.4,' GeV',/ . ,' chi br rat : ',F8.4,/ . ,' pdfs : ',A50,/) stop end ccccccccccccccccccccccccccccccccccccccccccccc include 'mstwpdf.f' include 'alphaS.f' c matching function function fw(x) implicit double precision(a-y) common/match/mm,nn fw=1d0-(1d0+dexp(mm*nn))/(dexp(nn*x)+dexp(mm*nn)) return end c binning subroutine subroutine binit(wt) implicit double precision(a-y) double precision q(4,20) common/mom/q common/vars/s,rts,mchi,yx common/vars1/ptgam,etagam,ptel2,ptel1,etael1,etael2 common/vars2/ptchi,etachi,ptpsi,etapsi call histo1(1,60,0d0,10d0,ptchi,wt) call histo1(2,24,0d0,8d0,etachi,wt) call histo1(3,20,0d0,3d0,ptel1,wt) call histo1(3,20,0d0,3d0,ptel2,wt) call histo1(4,20,0d0,6d0,etael1,wt) call histo1(4,20,0d0,6d0,etael2,wt) call histo1(5,20,0d0,1.3d0,ptgam,wt) call histo1(6,20,0d0,6d0,etagam,wt) c$$$ cost=(q(1,18)*q(1,16)+q(2,18)*q(2,16)+q(3,18)*q(3,16))/ c$$$ &dsqrt((q(1,18)**2+q(2,18)**2+q(3,18)**2)* c$$$ &(q(1,16)**2+q(2,16)**2+q(3,16)**2)) c$$$ c$$$ costa=(q(1,17)*q(1,5)+q(2,17)*q(2,5)+q(3,17)*q(3,5))/ c$$$ &dsqrt((q(1,17)**2+q(2,17)**2+q(3,17)**2)* c$$$ &(q(1,5)**2+q(2,5)**2+q(3,5)**2)) c$$$ c$$$ call histo1(1,20,-1d0,1d0,cost,wt) c$$$ call histo1(2,20,-1d0,1d0,costa,wt) return end subroutine cut(icut) implicit double precision(a-y) double precision q(4,20) integer icut common/mom/q common/vars/s,rts,mchi,yx common/vars1/ptgam,etagam,ptel2,ptel1,etael1,etael2 common/vars2/ptchi,etachi,ptpsi,etapsi common/cuts/etaelmax,etaelmin,ptelmin,ptphmin,ptchimin,ptpsimin &,ptpsimax icut=0 C -- photon cuts ptgam=dsqrt(q(1,6)**2+q(2,6)**2) if(ptgam.lt.ptphmin) return etagam=(0.5d0*dlog((q(4,6)+q(3,6))/(q(4,6)-q(3,6)))) if(etagam.gt.etaelmax) return if(etagam.lt.etaelmin) return c if(q(4,6).lt.0.4d0)return c -- lepton cuts ptel1=dsqrt(q(1,8)**2+q(2,8)**2) if(ptel1.lt.ptelmin) return ptel2=dsqrt(q(1,9)**2+q(2,9)**2) if(ptel2.lt.ptelmin) return pmod1=dsqrt(q(1,8)**2+q(2,8)**2+q(3,8)**2) pmod2=dsqrt(q(1,9)**2+q(2,9)**2+q(3,9)**2) etael1=(0.5d0*dlog((pmod1+q(3,8))/(pmod1-q(3,8)))) etael2=(0.5d0*dlog((pmod2+q(3,9))/(pmod2-q(3,9)))) if(etael1.gt.etaelmax) return if(etael1.lt.etaelmin) return if(etael2.gt.etaelmax) return if(etael2.lt.etaelmin) return c -- j/psi variables ptpsi=dsqrt(q(1,7)**2+q(2,7)**2) pmodpsi=dsqrt(q(1,7)**2+q(2,7)**2+q(3,7)**2) etapsi=.5d0*dlog((pmodpsi+q(3,7)) . /(pmodpsi-q(3,7))) c if(etapsi.lt.etaelmin) return c if(etapsi.gt.etaelmax) return if(ptpsi.lt.ptpsimin) return if(ptpsi.gt.ptpsimax) return c -- chi variables ptchi=dsqrt(q(1,5)**2+q(2,5)**2) pmodchi=dsqrt(q(1,5)**2+q(2,5)**2+q(3,5)**2) etachi=.5d0*dlog((pmodchi+q(3,5)) . /(pmodchi-q(3,5))) if(ptchi.lt.ptchimin)return icut=1 return end c prints histograms subroutine histo1(ih,ib,x0,x1,x,w) implicit real*8(a-h,o-z) character*1 regel(30),blank,star dimension h(20,100),hx(20),io(20),iu(20),ii(20) dimension y0(20),y1(20),ic(20) data regel / 30*' ' /,blank /' ' /,star /'*'/ save if(ih.eq.1)then open(10,file="out1.dat") elseif(ih.eq.2)then open(20,file="out2.dat") elseif(ih.eq.3)then open(30,file="out3.dat") elseif(ih.eq.4)then open(40,file="out4.dat") elseif(ih.eq.5)then open(50,file="out5.dat") elseif(ih.eq.6)then open(60,file="out6.dat") endif c print*,hx(ih) y0(ih)=x0 y1(ih)=x1 ic(ih)=ib if(x.lt.x0) goto 11 if(x.gt.x1) goto 12 ix=idint((x-x0)/(x1-x0)*dble(ib))+1 h(ih,ix)=h(ih,ix)+w if(h(ih,ix).gt.hx(ih)) hx(ih)=h(ih,ix) ii(ih)=ii(ih)+1 return 11 iu(ih)=iu(ih)+1 return 12 io(ih)=io(ih)+1 return entry histo2(ih,il) ib1=ic(ih) x01=y0(ih) x11=y1(ih) bsize=(x11-x01)/dble(ib1) hx(ih)=hx(ih)*(1.d0+1.d-06) if(il.eq.0) write(6,21) ih,ii(ih),iu(ih),io(ih) if(il.eq.1) write(6,22) ih,ii(ih),iu(ih),io(ih) 21 format(' no.',i3,' lin : inside,under,over ',3i6) 22 format(' no.',i3,' log : inside,under,over ',3i6) if(ii(ih).eq.0) goto 28 write(6,23) 23 format(35(1h ),3(10h----+----i)) do 27 iv=1,ib1 z=(dble(iv)-0.5d0)/dble(ib1)*(x11-x01)+x01 if(il.eq.1) goto 24 iz=idint(h(ih,iv)/hx(ih)*30.)+1 goto 25 24 iz=-1 if(h(ih,iv).gt.0.d0) .iz=idint(dlog(h(ih,iv))/dlog(hx(ih))*30.)+1 25 if(iz.gt.0.and.iz.le.30) regel(iz)=star write(6,26) z,h(ih,iv)/bsize,(regel(i),i=1,30) c write(50,*)z,h(ih,iv)/bsize ! Print histogram to file if(ih.eq.1)then write(10,*)z,h(ih,iv)/bsize ! Print histogram to file elseif(ih.eq.2)then write(20,*)z,h(ih,iv)/bsize elseif(ih.eq.3)then write(30,*)z,h(ih,iv)/bsize elseif(ih.eq.4)then write(40,*)z,h(ih,iv)/bsize elseif(ih.eq.5)then write(50,*)z,h(ih,iv)/bsize elseif(ih.eq.6)then write(60,*)z,h(ih,iv)/bsize endif 26 format(1h ,2g15.6,4h i,30a1,1hi) 36 format(1h ,2g15.6) if(iz.gt.0.and.iz.le.30) regel(iz)=blank 27 continue write(6,23) return 28 write(6,29) 29 format(' no entries inside histogram') return entry histo3(ih) do 31 i=1,100 31 h(ih,i)=0. hx(ih)=0. io(ih)=0 iu(ih)=0 ii(ih)=0 close(10) close(20) close(30) close(40) close(50) close(60) return end c rotates from pt(q(l))=0 frame to general lab frame subroutine rotate(vin,l,vout) double precision vin(4),vout(4) double precision rmatrix(4,4),q(4,20) double precision sy,cy,sx,cx common/mom/q integer i,j,k,l do k=1,4 vout(k)=0d0 enddo sx=q(1,l)/dsqrt(q(1,l)**2+q(2,l)**2+q(3,l)**2) cx=dsqrt(1d0-sx**2) sy=q(2,l)/dsqrt(q(2,l)**2+q(3,l)**2) cy=q(3,l)/dsqrt(q(2,l)**2+q(3,l)**2) rmatrix(1,1)=cx rmatrix(1,2)=0d0 rmatrix(1,3)=sx rmatrix(1,4)=0d0 rmatrix(2,1)=-sx*sy rmatrix(2,2)=cy rmatrix(2,3)=sy*cx rmatrix(2,4)=0d0 rmatrix(3,1)=-sx*cy rmatrix(3,2)=-sy rmatrix(3,3)=cy*cx rmatrix(3,4)=0d0 rmatrix(4,1)=0d0 rmatrix(4,2)=0d0 rmatrix(4,3)=0d0 rmatrix(4,4)=1d0 do i=1,4 do j=1,4 vout(i)=vout(i)+rmatrix(i,j)*vin(j) enddo enddo return end c boosting subroutine subroutine boost(p,pboo,pcm,plb) real*8 pboo(4),pcm(4),plb(4),p,fact plb(4)=(pboo(4)*pcm(4)+pboo(3)*pcm(3) & +pboo(2)*pcm(2)+pboo(1)*pcm(1))/p fact=(plb(4)+pcm(4))/(p+pboo(4)) do 10 j=1,3 10 plb(j)=pcm(j)+fact*pboo(j) return end c subprocess cross sections subroutine chi0nlo(sh,th,uh,out) implicit double precision (a-y) common/vars/s,rts,mchi,yx common/alph/alfas q=sh*th*uh p=sh*th+th*uh+uh*sh rp=0.006d0*mchi**2 pi=dacos(-1d0) out=(9d0*mchi**4*p**4*(mchi**8-2d0*mchi**4*p+p**2)- &6d0*mchi**2*p**3*q*(2d0*mchi**8-5d0*mchi**4*p+p**2)- &p**2*q**2*(mchi**8+2d0*mchi**4*p-p**2)+2d0*mchi**2*p*q**3* &(mchi**4-p)+6d0*mchi**4*q**4) out=out/(q*(q-mchi**2*p)**4) out=out*4d0*pi*alfas**3*rp/(mchi**3*sh**2) return end subroutine chi0nloq(sh,th,uh,out) implicit double precision (a-y) common/vars/s,rts,mchi,yx common/alph/alfas rp=0.006d0*mchi**2 pi=dacos(-1d0) out=-(th-3d0*mchi**2)**2*(sh**2+uh**2)/(th*(th-mchi**2)**4) out=out*8d0*pi*alfas**3*rp/(9d0*mchi**3*sh**2) return end subroutine chi1nlo(sh,th,uh,out) implicit double precision (a-y) common/vars/s,rts,mchi,yx common/alph/alfas q=sh*th*uh p=sh*th+th*uh+uh*sh rp=0.006d0*mchi**2 pi=dacos(-1d0) out=mchi**2*p**2*(mchi**4-4d0*p)+ &2d0*q*(-mchi**8+5d0*mchi**4*p+p**2)-15d0*mchi**2*q**2 out=out*p**2/(q-mchi**2*p)**4 out=out*12d0*pi*alfas**3*rp/(mchi**3*sh**2) return end subroutine chi1nloq(sh,th,uh,out) implicit double precision (a-y) common/vars/s,rts,mchi,yx common/alph/alfas rp=0.006d0*mchi**2 pi=dacos(-1d0) out=(-th*(sh**2+uh**2)-4d0*mchi**2*sh*uh)/(th-mchi**2)**4 out=out*16d0*pi*alfas**3*rp/(3d0*mchi**3*sh**2) return end subroutine chi2nlo(sh,th,uh,out) implicit double precision (a-y) common/vars/s,rts,mchi,yx common/alph/alfas q=sh*th*uh p=sh*th+th*uh+uh*sh rp=0.006d0*mchi**2 pi=dacos(-1d0) out=12d0*mchi**4*p**4*(mchi**8-2d0*mchi**4*p+p**2)- &3d0*mchi**2*p**3*q*(8d0*mchi**8-mchi**4*p+4d0*p**2)+ &2d0*p**2*q**2*(-7d0*mchi**8+43d0*mchi**4*p+p**2)+ &mchi**2*p*q**3*(16d0*mchi**4-61d0*p)+12d0*mchi**4*q**4 out=out/(q*(q-mchi**2*p)**4) out=out*4d0*pi*alfas**3*rp/(mchi**3*sh**2) return end subroutine chi2nloq(sh,th,uh,out) implicit double precision (a-y) common/vars/s,rts,mchi,yx common/alph/alfas rp=0.006d0*mchi**2 pi=dacos(-1d0) out=(th-mchi**2)**2*(th**2+6d0*mchi**4)- &2d0*sh*uh*(th**2-6d0*mchi**2*(th-mchi**2)) out=-out/(th*(th-mchi**2)**4) out=out*16d0*pi*alfas**3*rp/(9d0*mchi**3*sh**2) return end * * subtractive mitchell-moore generator * ronald kleiss - october 2, 1987 * * the algorithm is N(i)=[ N(i-24) - N(i-55) ]mod M, * implemented in a cirucular array with identifcation * of NR(i+55) and nr(i), such that effectively: * N(1) <--- N(32) - N(1) * N(2) <--- N(33) - N(2) .... * .... N(24) <--- N(55) - N(24) * N(25) <--- N(1) - N(25) .... * .... N(54) <--- N(30) - N(54) * N(55) <--- N(31) - N(55) * * in this version M =2**30 and RM=1/M=2.D0**(-30.D0) * * the array NR has been initialized by putting NR(i)=i * and subsequently running the algorithm 100,000 times. * subroutine R2455(Ran) IMPLICIT REAL*8(A-H,O-Z) DIMENSION N(55) DATA N/ . 980629335, 889272121, 422278310,1042669295, 531256381, . 335028099, 47160432, 788808135, 660624592, 793263632, . 998900570, 470796980, 327436767, 287473989, 119515078, . 575143087, 922274831, 21914605, 923291707, 753782759, . 254480986, 816423843, 931542684, 993691006, 343157264, . 272972469, 733687879, 468941742, 444207473, 896089285, . 629371118, 892845902, 163581912, 861580190, 85601059, . 899226806, 438711780, 921057966, 794646776, 417139730, . 343610085, 737162282,1024718389, 65196680, 954338580, . 642649958, 240238978, 722544540, 281483031,1024570269, . 602730138, 915220349, 651571385, 405259519, 145115737/ DATA M/1073741824/ DATA RM/0.9313225746154785D-09/ DATA K/55/,L/31/ IF(K.EQ.55) THEN K=1 ELSE K=K+1 ENDIF IF(L.EQ.55) THEN L=1 ELSE L=L+1 ENDIF J=N(L)-N(K) IF(J.LT.0) J=J+M N(K)=J RAN=J*RM END