program superCHIC implicit double precision(a-y) implicit complex*16(z) double precision pt1(2),pt2(2),ptdif(2),ptx(2),q(4,20),g(4,4) double precision pboo(4),pcm(4),plb(4) double precision n1(4),n2(4),n1rot(4),n2rot(4),e0(4),e0rot(4) double precision q5(4),q6(4),q7(4),q8(4),q9(4),q6c(4),q7c(4) double precision tvec(4),uvec(4),pgtot(4) complex*16 zq6(4),zq5(4),zq7(4),zq8(4) complex*16 echiplus(4),echiminus(4),echi0(4),echi1(4),cechi1(4) complex*16 epsi0(4),epsiplus(4),epsiminus(4),epsi1(4) &,cepsi1(4),epsi1c(4) complex*16 echi2p2(4,4),echi2p1(4,4),echi20(4,4),echi2m1(4,4) &,echi2m2(4,4),echi2(4,4),cechi2(4,4) integer h,i,j,l,k,chiflag,nhist,ntotal,pflag,eflag &,psiflag,icut,cross,ncut,allflag,xflag,ID(20),ID0,ID1,ID2 common/mom/q common/mompt/pt1,pt2,ptx common/vars/s,rts,mchi,ptcut,etacut,ecut,rcut ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c SuperCHIC Monte Carlo generator for central c c exclusive production: c c c c p(1) p(2) --> p(3) + X(5) + p(4) c c c c xflag=0 : Chi_c meson production c c c c xflag=1 : Chi_b meson production c c c c chiflag=0 chi(0++) c c c c chiflag=1 chi(1++) c c c c chiflag=2 chi(2++) c c c c chiflag=3 : generate all c c c c xflag=2 : digamma CEP :- c c c c X(5) = gamma(12)gamma(13) c c c c Calculates total cross section integrated over a c c given chic rapidity interval (cross = 1) with c c a phenomenological fit for the rapidity dependence, c c or dsigma/dy(chic) at y=0 (cross = 2). c c c c chi_b/chi_c decays: c c c c X(5) --> gamma(6) psi/ups(7) c c c c psi/ups(7) --> mu(8) mu(9) 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 cms energies: c c c c eflag=0 rts= 1.96 TeV c c eflag=1 rts= 14 TeV c c c c note: changing rts from these values manually will c c result in incorrect normalisation due to survival c c factor energy dependence c c c c This is version 1.11: 14 Jan 2010 c c c c Comments etc to: lh367@cam.ac.uk c c c c For more details see, e.g.: c c c c L.A. Harland-Lang, V.A. Khoze, M.G. Ryskin c c and W.J. Stirling arXiv:hep-ph/0909.4748 (chi_c/b) c c c c V.A. Khoze, A.D. Martin, M.G. Ryskin c c and W.J. Stirling Eur. Phys. J. C 38 (2005) 475 c c [arXiv:hepph/0409037] (digamma) c c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c set flags xflag=2 ! chi_c=0, chi_b=1, gamgam=2 chiflag=0 ! chi_(0,1,2) or all = 3 eflag=0 ! 1.96 TeV = 0 or 14 TeV = 1 c total or differential cross section (see preamble) cross=1 c chi cuts: pt(muons) and y(chi), implemented in subroutine icut ptcut=0d0 ! ptmin of final state muons ycut=2d0 ! |rapidity| max of chic etacut=0.6d6 ! |pseudorapidity| max of final state muons c digamma cuts: ecut=5d0 rcut=1d0 ! |pseudorapidity| max of photons c initialise counters etc ntotal=1000000 nhist=1 sum=0d0 sum1=0d0 sum2=0d0 ncut=0 c beam energy etc if(eflag.eq.0)then rts=1.96d3 ! GeV elseif(eflag.eq.1)then rts=14d3 endif ebeam=rts/2d0 s=rts*rts zi=(0d0,1d0) rt2=dsqrt(2d0) pi=dacos(-1d0) c Proton form factor (as in ~ exp(bt)), in GeV^-2 b0=2d0 ! Elastic slope b=2d0*b0 bp=rts/dsqrt(2.0d0) alphap=0.1d0 ! Pomeron slope delta=0.2d0 ! Pomeron intercept c open(30,file="erec.dat") c initialise histograms do i=1,nhist call histo3(i) enddo if(xflag.eq.0.or.xflag.eq.1)then if(chiflag.eq.3)then allflag=0 else allflag=1 endif else allflag=1 endif c set masses and other parameters c chi_c values if(xflag.eq.0)then rtso=60d0 if(eflag.eq.0)then sv0=0.098d0 !eikonal survival factor^2, chi_c0 sv1=0.146d0 ! " " " , chi_c1 elseif(eflag.eq.1)then sv0=0.057d0 sv1=0.086d0 endif c J/psi ID(7)=443 ! PDG number mpsi=3.097d0 brpsimu=5.93d-2 ! Br(J/psi --> mu+mu-) errbrpsi=0.06d-2 ccc chi_c(0++) ID0=10441 mchi0=3.415d0 beff0=6.64d0 ! effective slope neff0=1.59d0 ! normalisation brchipsi0=1.14d-2 ! Br(chi_c --> J/psi) errbrchi0=0.08d-2 c 'bare' central sigma at 60 Gev in nb bare0=2.2d3/1.74d0 ccc chi_c(1++) ID1=20443 mchi1=3.511d0 beff1=4.57d0 neff1=0.36d0 brchipsi1=34.1d-2 errbrchi1=1.5d-2 bare1=2.2d3/3.511d0**2/1.74d0 ccc chi_c(2++) ID2=445 mchi2=3.566d0 c chi2 survival factor set below beff2=5.9d0 neff2=0.714d0 brchipsi2=19.4d-2 errbrchi2=0.8d-2 bare2=2.2d3/1.74d0 c rapidity dependence (~same for 3 chi states) ng=30d0+8.1d-3*(rts-1.96d3) c chi_b values elseif(xflag.eq.1)then rtso=173d0 if(eflag.eq.0)then sv0=0.087d0 sv1=0.142d0 elseif(eflag.eq.1)then sv0=0.049d0 sv1=0.082d0 endif c Upsilon ID(7)=553 mpsi=9.46d0 brpsimu=2.48d-2 errbrpsi=0.05d-2 ccc chi_b(0++) ID0=10551 mchi0=9.589d0 beff0=6.1d0 neff0=2.32d0 brchipsi0=3d-2 !!! Note branching ratio experimentally < 6% only errbrchi0=1d-2 c 'bare' central sigma at 60 Gev in pb bare0=2.68d3/3.05d0 ccc chi_b(1++) ID1=20553 mchi1=9.893d0 beff1=4.5d0 neff1=0.548d0 brchipsi1=35d-2 errbrchi1=8d-2 bare1=2.68d3/(3.05d0*9.983d0**2) ccc chi_b(2++) ID2=555 mchi2=9.912d0 beff2=5.4d0 neff2=0.85d0 brchipsi2=22d-2 errbrchi2=4d-2 bare2=2.68d3/3.05d0 ng=14d0+3.5d-3*(rts-1.96d3) 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) if(allflag.eq.0)then call r2455(ranc) if(ranc.lt.0.3333)then chiflag=0 elseif(ranc.lt.0.6666)then chiflag=1 else chiflag=2 endif endif if(xflag.eq.0.or.xflag.eq.1)then ! chi_c/b if(chiflag.eq.0)then ID(5)=ID0 mchi=mchi0 sv=sv0 beff=beff0 neff=neff0 brchipsi=brchipsi0 bare=bare0 errbrchi=errbrchi0 elseif(chiflag.eq.1)then ID(5)=ID1 mchi=mchi1 sv=sv1 beff=beff1 neff=neff1 brchipsi=brchipsi1 bare=bare1 errbrchi=errbrchi1 elseif(chiflag.eq.2)then ID(5)=ID2 mchi=mchi2 beff=beff2 neff=neff2 brchipsi=brchipsi2 bare=bare2 errbrchi=errbrchi2 endif mx=mchi mmu=0.106d0 elseif(xflag.eq.2)then ! digamma c integrate over digamma system invariant mass call r2455(rm) c mgam integration range mmax=70d0 ! note high mmax values strongly suppressed mmin=2d0*ecut c mmin=10d0 my=(1d0/(mmin)**3-1d0/(mmax)**3)/3d0*rm-1d0/(3d0*mmin**3) mgam=1d0/(-my*3d0)**(1d0/3d0) weightgam=mgam**4*(1d0/(mmin)**3-1d0/(mmax)**3)/3d0 mx=mgam beff=4d0 sv=0.048d0 bare=1d0 ng=29.5d0 endif bb=beff xoo=mx/rtso xo=mx/rts ccccccccccccccccccccccccccccccccccccccccccccccccccc c incoming protons ID(1)=2212 q(1,1)=0d0 q(2,1)=0d0 q(3,1)=ebeam q(4,1)=ebeam if(eflag.eq.0)then ID(2)=-2212 else ID(2)=2212 endif q(1,2)=0d0 q(2,2)=0d0 q(3,2)=-ebeam q(4,2)=ebeam c Transverse momenta c Note: minor change so proton 1 does not always have p_x=0 phi1=2d0*pi*ran0 phi2=2d0*pi*ran1 pt1sq=-dlog(ran2)/(bb) pt2sq=-dlog(ran3)/(bb) pt1(1)=dsqrt(pt1sq)*dsin(phi1) pt1(2)=dsqrt(pt1sq)*dcos(phi1) pt2(1)=dsqrt(pt2sq)*dsin(phi2) pt2(2)=dsqrt(pt2sq)*dcos(phi2) ptx(1)=-pt1(1)-pt2(1) ptx(2)=-pt1(2)-pt2(2) ptxsq=ptx(1)**2+ptx(2)**2 ptdif(1)=pt1(1)-pt2(1) ptdif(2)=pt1(2)-pt2(2) ptw=(ptdif(1))**2+(ptdif(2))**2 weight=(dexp(bb*pt1sq)*dexp(bb*pt2sq))/bb**2 if(xflag.eq.0.or.xflag.eq.1)then if(chiflag.eq.2)then c chic(2++): pt(chi) dependent survival factor call chi2sv(dsqrt(ptxsq),xflag,eflag,sv) c fit to exact p_t dependence weight=weight*(pt1sq*pt2sq+0.105d0*pt1(2)*pt2(2)) &/(pt1sq*pt2sq) endif endif cccccccccccccccccccccccccccccccccccccccccccccccccccc c x rapidity c transverse mass rmx=dsqrt(mx**2+ptxsq) if(cross.eq.1)then c total cross section calculation ym=dlog(rts/rmx) if(ycut.gt.ym)then ymax=ym ! Kinematic limit on chi rapidity else ymax=ycut endif ymin=-ymax yx=ymax*ran4*ran5**(1d0/3d0) if(ran6.gt.0.5d0) yx=-yx weighty=ymax**3*4d0/3d0/(ymax**2-yx**2) weight=weight*weighty elseif(cross.eq.2)then c dsigma/dy|y=yx yx=0d0 endif c longitudinal momenta fractions x1=dexp(yx)*rmx/rts x2=dexp(-yx)*rmx/rts ccccccccccccccccccccccccccccccccccccccccccccccccccccc aa1=bp*(1d0-x1) aa2=bp*(1d0-x2) cc1=0.5d0*pt2sq cc2=0.5d0*pt1sq c impose massless on-shell condition by solving c p1+ + cc1/p2- = aa1 c p2- + cc2/p1+ = aa2 root1sq=(cc1-cc2-aa1*aa2)**2-4d0*cc2*aa1*aa2 root2sq=(cc2-cc1-aa1*aa2)**2-4d0*cc1*aa1*aa2 if(root1sq.le.0d0.or.root2sq.le.0d0)then weight=0d0 else p1p=(cc2-cc1+aa1*aa2+dsqrt(root1sq))/(2d0*aa2) p2m=(cc1-cc2+aa1*aa2+dsqrt(root2sq))/(2d0*aa1) p1m=pt1sq/(2d0*p1p) p2p=pt2sq/(2d0*p2m) ccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Outgoing protons ID(3)=ID(1) q(1,3)=pt1(1) q(2,3)=pt1(2) q(3,3)=(p1p-p1m)/rt2 q(4,3)=(p1p+p1m)/rt2 ID(4)=ID(2) q(1,4)=pt2(1) q(2,4)=pt2(2) q(3,4)=(p2p-p2m)/rt2 q(4,4)=(p2p+p2m)/rt2 c chi momentum do j=1,4 q(j,5)=q(j,1)+q(j,2)-q(j,3)-q(j,4) enddo c print*,q(1,5),q(2,5),q(3,5),q(4,5) c print*,dsqrt((q(4,5)**2-q(3,5)**2-q(2,5)**2-q(1,5)**2)) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Final state particle 4-vector generation ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c digamma CEP if(xflag.eq.2)then x11=x1 x22=x2 d1=-x1*rts*(1d0-dsqrt(1d0+4d0*pt1sq/(x1*rts)**2)) &/2d0 d2=-x2*rts*(1d0-dsqrt(1d0+4d0*pt2sq/(x2*rts)**2)) &/2d0 c Assign incoming gluon momenta... do l=1,4 pboo(l)=q(l,5) q(l,10)=x1*q(l,1) q(l,11)=x2*q(l,2) enddo q(1,10)=-pt1(1) q(2,10)=-pt1(2) q(1,11)=-pt2(1) q(2,11)=-pt2(2) q(4,10)=q(4,10)+d1 q(4,11)=q(4,11)+d2 c ...and make on-shell approximation pboo(4)=pboo(4)+d1+d2 sh=(pboo(4)**2-pboo(3)**2-pboo(2)**2-pboo(1)**2) mxx=dsqrt(sh) c 2-body phase space for weight generation call r2455(rang1) call r2455(rang2) costg=2d0*rang1-1d0 phig=2d0*pi*rang2 sintg=dsqrt(1d0-costg**2) cphig=dcos(phig) sphig=dsin(phig) pgmod=mxx/2d0 pcm(4)=pgmod pcm(1)=pgmod*sintg*sphig pcm(2)=pgmod*sintg*cphig pcm(3)=pgmod*costg call boost(mxx,pboo,pcm,plb) do k=1,4 q(k,12)=plb(k) q(k,13)=pboo(k)-q(k,12) tvec(k)=-q(k,13)+x22*q(k,2) uvec(k)=-q(k,13)+x11*q(k,1) enddo tvec(1)=tvec(1)-pt2(1) tvec(2)=tvec(2)-pt2(2) tvec(4)=tvec(4)+d2 uvec(1)=uvec(1)-pt1(1) uvec(2)=uvec(2)-pt1(2) uvec(4)=uvec(4)+d1 t=rdot(tvec,tvec) u=rdot(uvec,uvec) c generate photon 4-vectors (exact kinematics) do l=1,4 pboo(l)=q(l,5) enddo costg=2d0*rang1-1d0 phig=2d0*pi*rang2 sintg=dsqrt(1d0-costg**2) cphig=dcos(phig) sphig=dsin(phig) pgmod=mgam/2d0 pcm(4)=pgmod pcm(1)=pgmod*sintg*sphig pcm(2)=pgmod*sintg*cphig pcm(3)=pgmod*costg call boost(mx,pboo,pcm,plb) do k=1,4 q(k,12)=plb(k) q(k,13)=pboo(k)-q(k,12) pgtot(k)=q(k,12)+q(k,13) enddo q(4,10)=q(4,10)-d1 q(4,11)=q(4,11)-d2 call boost(mx,pboo,pcm,plb) do k=1,4 q(k,12)=plb(k) q(k,13)=pboo(k)-q(k,12) pgtot(k)=q(k,12)+q(k,13) enddo ID(12)=22 ID(13)=22 elseif(xflag.eq.0.or.xflag.eq.1)then ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c chi_(c/b) polarization + decay products ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c generate chi polarization vectors in pt(chi)=0 frame... c longitudinal vector e0(1)=0d0 e0(2)=0d0 e0(3)=q(4,5)/mchi e0(4)=dsqrt(q(4,5)**2-mchi**2)/mchi c real transverse basis vectors n1(1)=1d0 n1(2)=0d0 n1(3)=0d0 n1(4)=0d0 n2(1)=0d0 n2(2)=1d0 n2(3)=0d0 n2(4)=0d0 c ...and rotate to lab basis do k=1,4 pboo(k)=q(k,5) enddo call rotate(e0,5,e0rot) call rotate(n1,5,n1rot) call rotate(n2,5,n2rot) c generate polarization vectors do j=1,4 echiplus(j)=-(n1rot(j)+zi*n2rot(j))/rt2 echiminus(j)=(n1rot(j)-zi*n2rot(j))/rt2 echi0(j)=e0rot(j) enddo c various checks c qmod=dsqrt(q(1,5)**2+q(2,5)**2+q(3,5)**2) c emod=sqrt(abs(echi0(1))**2+abs(echi0(2))**2+abs(echi0(3))**2) c print*,q(1,5)/qmod,q(2,5)/qmod,q(3,5)/qmod c print*,echi0(1)/emod,echi0(2)/emod,echi0(3)/emod c print*,abs(echi0(4))**2-abs(echi0(3))**2-abs(echi0 c &(2))**2-abs(echi0(1))**2 c print*,cdabs(echiminus(4)*q(4,5)-echiminus(3)*q(3,5)-echiminus(2)* c &q(2,5)-echiminus(1)*q(1,5)) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 ID(6)=22 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc 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 ID(8)=-13 ID(9)=13 endif ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c place cuts call cut(icut) if(icut.eq.0)then weight=0d0 ncut=ncut+1 goto 700 endif ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Weight evaluation ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(xflag.eq.2)then c (g(10)g(11) --> gam(12)gam(13)) subprocess (J_z=0) amplitude: zamppp=1d0 zammpp=1d0 zapmmp=-0.5d0*((sh**2+t**2)/u**2)*(dlog(-t/sh))**2 & -((t-sh)/u)*dlog(-t/sh)-1d0-zi*pi*(((t**2+sh**2)/u**2) & *dlog(-t/sh)+(t-sh)/u) zapppp=-0.5d0*((t**2+u**2)/sh**2)*(dlog(t/u)**2+pi**2) & -((t-u)/sh)*dlog(t/u)-1d0 c average over gluon helicities at amplitude level: zamp1=0.5d0*(zapppp+zammpp) zamp2=0.5d0*(zamppp+zamppp) mu=0.5d0*mgam ! hard scale mu = mgam/2 if(mu.gt.4.75d0)then nf=5d0 lambdacap=0.22d0 qf=11d0/9d0 else nf=4d0 lambdacap=0.3d0 qf=10d0/9d0 endif b=(33d0-2d0*nf)/(12d0*pi) b1=(153d0-19d0*nf)/(2d0*pi*(33d0-2d0*nf)) alphas=1d0/(b*dlog(mu**2/lambdacap**2)) alphas=alphas*(1d0-alphas*b1*dlog(dlog(mu**2/lambdacap**2))) weight=weight*2d0*(cdabs(zamp1)**2+cdabs(zamp2)**2) weight=weight*(8d0*0.5d0*(qf)*(1d0/137d0)*alphas)**2 weight=weight/(2d0*16d0*pi*mgam**2) call lumw(mgam,wl) !'effective luminosity' weight=weight*wl/mgam*2d0 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc elseif(xflag.eq.0.or.xflag.eq.1)then c gg --> chi(1++) vertex if(chiflag.eq.1)then c generate chi polarization at random call r2455(ranchi) do j=1,4 if(ranchi.lt.0.3333)then ! lambda=0 echi1(j)=echi0(j) cechi1(j)=echi0(j) pflag=1 elseif(ranchi.lt.0.6666)then ! lambda=+1 echi1(j)=echiplus(j) cechi1(j)=conjg(echiplus(j)) pflag=2 else ! lambda=-1 echi1(j)=echiminus(j) cechi1(j)=conjg(echiminus(j)) pflag=3 endif enddo weight=weight*(cdabs(echi1(1)*dcmplx(ptdif(2))-echi1(2) &*ptdif(1)))**2 c re-weight as summing over polarizations weight=weight*3d0 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc c gg --> chi(2++) vertex elseif(chiflag.eq.2)then c polarization tensors with Jz=+2,+1,0,-1,-2 do l=1,4 do k=1,4 echi2p2(l,k)=echiplus(l)*echiplus(k) echi2p1(l,k)=(echiplus(l)*echi0(k)+echi0(l)*echiplus(k))/rt2 echi20(l,k)=(echiplus(l)*echiminus(k)+2d0*echi0(l)*echi0(k) & +echiminus(l)*echiplus(k))/dsqrt(6d0) echi2m1(l,k)=(echiminus(l)*echi0(k)+echi0(l)*echiminus(k))/rt2 echi2m2(l,k)=echiminus(l)*echiminus(k) enddo enddo c print*,echi20(4,4)-echi20(3,3)-echi20(2,2)-echi20(1,1) c print*,echi20(2,4)*q(4,5)-echi20(2,3)*q(3,5)-echi20(2,2)*q(2,5) c &-echi20(2,1)*q(1,5) c generate chi(2++) polarization states at random call r2455(ranchi) do j=1,4 do l=1,4 if(ranchi.lt.0.2)then ! lambda=+2 echi2(j,l)=echi2m2(j,l) cechi2(j,l)=conjg(echi2m2(j,l)) pflag=1 elseif(ranchi.lt.0.4)then ! lambda=+1 echi2(j,l)=echi2m1(j,l) cechi2(j,l)=conjg(echi2m1(j,l)) pflag=2 elseif(ranchi.lt.0.6)then ! lambda=0 echi2(j,l)=echi20(j,l) cechi2(j,l)=conjg(echi20(j,l)) pflag=3 elseif(ranchi.lt.0.8)then ! lambda=-1 echi2(j,l)=echi2p1(j,l) cechi2(j,l)=conjg(echi2p1(j,l)) pflag=4 else ! lambda=-2 echi2(j,l)=echi2p2(j,l) cechi2(j,l)=conjg(echi2p2(j,l)) pflag=5 endif enddo enddo weight=weight*2d0*(cdabs(s*(echi2(1,1)*pt1(1)*pt2(1) &+echi2(1,2)*pt1(1)*pt2(2)+echi2(2,1)*pt1(2)*pt2(1) &+echi2(2,2)*pt1(2)*pt2(2))-2d0*(pt1(1)*pt2(1)+pt1(2)*pt2(2)) &*(echi2(3,3)*q(3,1)*q(3,2)-echi2(3,4)*q(3,1)*q(4,2) &-echi2(4,3)*q(4,1)*q(3,2)+echi2(4,4)*q(4,1)*q(4,2))))**2/s**2 c re-weight as summing over chi polarizations weight=weight*5d0 endif endif ccccccccccccccccccccccccccccccccccccccccccccccccccccccc t1=-rts*p1m*rt2 t2=-rts*p2p*rt2 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc c vertex factor fnt=dexp((bb*t1+bb*t2)) c energy dependence for chi_c/b fnx=(x1*x2/xoo**2)**(-2d0*delta) c matrix element rapidity dependence (phenomenological fit) c note: fny=1 by definition for y(chic)=0 fny=((1d0-x1)*(1d0-x2)/(1d0-xo)**2)**NG c Enhanced rescattering effects, S_{enh}^2. NB: value subject to change if(eflag.eq.0)then senh=0.308d0*mx**(0.273d0) elseif(eflag.eq.1)then senh=0.1885d0*mx**0.0881d0 endif ccccccccccccccccccccccccccc if(xflag.eq.0.or.xflag.eq.1)then weight=weight*fnt*bare*fnx*fny*neff*sv*senh*brchipsi*brpsimu elseif(xflag.eq.2)then weight=weight*fnt*weightgam*fny*senh*bb**2*sv*389d9 endif if(allflag.eq.0)then weight=weight*3d0 endif endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c chi_(c/b) decays cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(xflag.eq.0.or.xflag.eq.1)then c decay 5(mchi) --> 6(0) + 7(mpsi) c assign J/psi polarization vectors... e0(1)=0d0 e0(2)=0d0 e0(3)=q(4,7)/mpsi e0(4)=dsqrt(q(4,7)**2-mpsi**2)/mpsi call rotate(e0,7,e0rot) call rotate(n1,7,n1rot) call rotate(n2,7,n2rot) do j=1,4 epsiplus(j)=-(n1rot(j)+zi*n2rot(j))/rt2 epsiminus(j)=(n1rot(j)-zi*n2rot(j))/rt2 epsi0(j)=e0rot(j) enddo c ...and generate at random call r2455(ranchi1) if(chiflag.eq.0)then ! J/psi transversely polarized do j=1,4 if(ranchi1.lt.0.5)then epsi1(j)=epsiplus(j) ! lambda=+1 cepsi1(j)=conjg(epsiplus(j)) psiflag=2 else epsi1(j)=epsiminus(j) ! lambda=-1 cepsi1(j)=conjg(epsiminus(j)) psiflag=3 endif enddo elseif(chiflag.eq.1.or.chiflag.eq.2)then do j=1,4 if(ranchi1.lt.0.3333)then ! lambda=0 epsi1(j)=epsi0(j) cepsi1(j)=epsi0(j) psiflag=1 elseif(ranchi1.lt.0.6666)then ! lambda=+1 epsi1(j)=epsiplus(j) cepsi1(j)=conjg(epsiplus(j)) psiflag=2 else ! lambda=-1 epsi1(j)=epsiminus(j) cepsi1(j)=conjg(epsiminus(j)) psiflag=3 endif enddo endif c various vectors... do j=1,4 q5(j)=q(j,5) q6(j)=q(j,6) q7(j)=q(j,7) zq5(j)=q(j,5) zq6(j)=q(j,6) zq7(j)=q(j,7) enddo ccccccccccccccccccccccccccccccccccccccccccccccccccc c calculates (chi and psi polarization dependent) weight for decay cccccccccccccccccccccccccccccccccccccccccccccccccccc c chi1 cccccccccccccccccccccccccccccccccccccccccccccccccccc if(chiflag.eq.1)then call cdot(zq6,echi1,zq6echi1) call cdot(zq6,cechi1,zq6cechi1) call dot(q6,q7,q6q7) call cdot(zq7,echi1,zq7echi1) pnorm=cdabs(zq6echi1)**2+q6q7*(q6q7+2d0*dreal(zq6cechi1 &*zq7echi1))/mpsi**2 call cdot(echi1,epsi1,zechi1epsi1) call cdot(zq6,cepsi1,zq6cepsi1) weight=weight*(cdabs(zq6echi1)**2+cdabs(zq6cepsi1)**2+ &2d0*dreal(zq6cechi1*zechi1epsi1*zq6cepsi1))/pnorm c re-weight as summing over J/psi polarizations weight=weight*3d0 cccccccccccccccccccccccccccccccccccccccccccccccccccc c chi2 cccccccccccccccccccccccccccccccccccccccccccccccccccc elseif(chiflag.eq.2)then c initialise values weight2=0d0 pnorm=0d0 c metric tensor do j=1,4 do k=1,4 g(j,k)=0d0 enddo enddo g(1,1)=-1d0 g(2,2)=-1d0 g(3,3)=-1d0 g(4,4)=1d0 c generate contravariant 4-vectors q6c(4)=q(4,6) q7c(4)=q(4,7) epsi1c(4)=epsi1(4) do j=1,3 q6c(j)=-q(j,6) q7c(j)=-q(j,7) epsi1c(j)=-epsi1(j) enddo call dot(q6,q7,q6q7) call cdot(zq6,cepsi1,zq6cepsi1) call cdot(zq6,epsi1,zq6epsi1) c normalisation do j=1,4 do l=1,4 do k=1,4 do h=1,4 pnorm=pnorm+g(k,j)*(echi2(h,k)*cechi2(j,l)*(q6q7**2* &g(h,l)+mchi**2*q6c(h)*q6c(l))-2d0*q6q7*q6c(h)*q7c(l) &*dreal(echi2(h,k)*cechi2(j,l))) enddo enddo enddo enddo c polarization dependent weight do j=1,4 do l=1,4 do k=1,4 do h=1,4 weight2=weight2+(2d0*dreal(echi2(h,k)* &cechi2(l,j)*q6c(h)*(q6q7*epsi1c(j)-zq6epsi1*q7c(j) &)*(conjg(epsi1c(k))*q7c(l)-conjg(epsi1c(l))*q7c(k))) &-echi2(h,k)*cechi2(l,j)*(q6c(l)*q6c(h)*(mpsi**2*epsi1c(j) &*conjg(epsi1c(k))-q7c(k)*q7c(j))+g(h,l)*(q6q7*epsi1c(j)- &zq6epsi1*q7c(j))*(q6q7*conjg(epsi1c(k)) &-zq6cepsi1*q7c(k)))) enddo enddo enddo enddo weight=weight*weight2/pnorm weight=weight*3d0 endif cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c decay 7(mpsi) --> 8(mmu) + 9(mmu) do k=1,4 q8(k)=q(k,8) zq8(k)=q(k,8) q9(k)=q(k,9) enddo c calculate (polarization dependent) weights for decay call dot(q8,q9,q8q9) call cdot(epsi1,zq8,zq8epsi1) weight=weight*(q8q9-2d0*cdabs(zq8epsi1)**2+mmu**2) &/(2d0*(q8q9+2d0*mmu**2)) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc weight=weight*3d0 endif ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc wthist=weight/dfloat(ntotal) call binit(wthist) cccccc Calculate average weight and error 700 sum=sum+weight sum1=sum1+weight**2 sum2=sum2+weight*sv c Event record c$$$ do k=3,9 c$$$ c$$$ write(30,300)ID(k),q(1,k),q(2,k),q(3,k),q(4,k) c$$$ c$$$ enddo c$$$ c$$$ write(30,301)weight c$$$ write(30,*) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c End of event loop c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc enddo c$$$ 300 format(i10,1x,E17.10,1x,E17.10,1x,E17.10,1x,E17.10) c$$$ 301 format(E17.10) c close(30) sum=sum/dfloat(ntotal) sum1=sum1/dfloat(ntotal) var=dsqrt((sum1-sum**2)/dfloat(ntotal)) sum2=sum2/dfloat(ntotal) c create histograms do i=1,nhist call histo2(i,0) enddo c print information to screen c print*,sum,var if(allflag.eq.1)then if(xflag.eq.0)then if(chiflag.eq.0)then write(6,*) write(6,*)' chi_c(0++) production' elseif(chiflag.eq.1)then write(6,*) write(6,*)' chi_c(1++) production' elseif(chiflag.eq.2)then write(6,*) write(6,*)' chi_c(2++) production' endif if(cross.eq.1)then write(6,640)rts,mchi,mpsi,brchipsi,errbrchi,brpsimu,errbrpsi, . beff,delta, . sum2/sum,senh,ntotal,ntotal-ncut,2d0*ymax,sum,var elseif(cross.eq.2)then write(6,650)rts,mchi,mpsi,brchipsi,errbrchi,brpsimu,errbrpsi, . beff,delta, . sum2/sum,senh,ntotal,ntotal-ncut,yx,sum,var endif elseif(xflag.eq.1)then if(chiflag.eq.0)then write(6,*) write(6,*)' chi_b(0++) production' elseif(chiflag.eq.1)then write(6,*) write(6,*)' chi_b(1++) production' elseif(chiflag.eq.2)then write(6,*) write(6,*)' chi_b(2++) production' endif if(cross.eq.1)then write(6,641)rts,mchi,mpsi,brchipsi,errbrchi,brpsimu,errbrpsi, . beff,delta, . sum2/sum,senh,ntotal,ntotal-ncut,2d0*ymax,sum,var elseif(cross.eq.2)then write(6,651)rts,mchi,mpsi,brchipsi,errbrchi,brpsimu,errbrpsi, . beff,delta, . sum2/sum,senh,ntotal,ntotal-ncut,yx,sum,var endif endif elseif(allflag.eq.0)then if(xflag.eq.0)then write(6,*) write(6,*)' chi_c(0++),chi_c(1++),chi_c(2++) production' if(cross.eq.1)then write(6,660)rts,mpsi,delta,seff,ntotal,ntotal-ncut,2d0*ymax & ,sum,var elseif(cross.eq.2)then write(6,670)rts,mpsi,delta,seff,ntotal,ntotal-ncut,yx & ,sum,var endif elseif(xflag.eq.1)then write(6,*) write(6,*)' chi_b(0++),chi_b(1++),chi_b(2++) production' if(cross.eq.1)then write(6,661)rts,mpsi,delta,seff,ntotal,ntotal-ncut,2d0*ymax & ,sum,var elseif(cross.eq.2)then write(6,671)rts,mpsi,delta,seff,ntotal,ntotal-ncut,yx & ,sum,var endif endif endif if(xflag.eq.2)then write(6,*) write(6,*)' digamma production' if(cross.eq.1)then write(6,540)rts,beff,sv,ntotal,ntotal-ncut . ,2d0*ymax,sum,var elseif(cross.eq.2)then write(6,550)rts,beff,sum2/sum,ntotal,ntotal-ncut . ,yx,sum,var endif endif 540 format(/, . 3x,'collider energy (Gev) : ',f10.3,/, . 3x,'t slope b (Gev^-2) : ',f10.5,/, . 3x,'S^2_eik : ',f10.5,/, . 3x,'number of events : ',i9,/, . 3x,'events passing cuts : ',i9,/, . 3x,'Delta y : ',f10.5,/,/, . 3x,'sigma (fb) : ',f10.5,' +- ',f10.5, . /) 550 format(/, . 3x,'collider energy (Gev) : ',f10.3,/, . 3x,'t slope (Gev^-2) : ',f10.5,/, . 3x,'S^2_eik : ',f10.5,/, . 3x,'number of events : ',i9,/, . 3x,'events passing cuts : ',i9,/,/, . 3x,'dsigma/dy|y=',f3.1,' (fb) : ',f10.5,' +- ',f10.5, . /) 660 format(/, . 3x,'collider energy (Gev) : ',f10.3,/, . 3x,'j/psi mass (Gev) : ',f10.5,/, . 3x,'delta=alpha_p(0)-1 : ',f10.5,/, . 3x,'S^2_enh : ',f10.5,/, . 3x,'number of events : ',i9,/, . 3x,'events passing cuts : ',i9,/, . 3x,'Delta y : ',f10.5,/,/, . 3x,'sigma (nb) : ',f10.5,' +- ',f10.5, . /) 670 format(/, . 3x,'collider energy (Gev) : ',f10.3,/, . 3x,'j/psi mass (Gev) : ',f10.5,/, . 3x,'delta=alpha_p(0)-1 : ',f10.5,/, . 3x,'S^2_enh : ',f10.5,/, . 3x,'number of events : ',i9,/, . 3x,'events passing cuts : ',i9,/,/, . 3x,'dsigma/dy|y=',f3.1,' (nb) : ',f10.5,' +- ',f10.5, . /) 640 format(/, . 3x,'collider energy (Gev) : ',f10.3,/, . 3x,'chi mass (Gev) : ',f10.5,/, . 3x,'j/psi mass (Gev) : ',f10.5,/, . 3x,'br(chic --> psi gamma) : ',f10.5,' +- ',f10.5,/ . 3x,'br(psi --> mu+mu-) : ',f10.5,' +- ',f10.5,/ . 3x,'t slope b_eff (Gev^-2) : ',f10.5,/, . 3x,'delta=alpha_p(0)-1 : ',f10.5,/, . 3x,'S^2_eik : ',f10.5,/, . 3x,'S^2_enh : ',f10.5,/, . 3x,'number of events : ',i9,/, . 3x,'events passing cuts : ',i9,/, . 3x,'Delta y : ',f10.5,/,/, . 3x,'sigma (nb) : ',f10.5,' +- ',f10.5, . /) 650 format(/, . 3x,'collider energy (Gev) : ',f10.3,/, . 3x,'chi mass (Gev) : ',f10.5,/, . 3x,'j/psi mass (Gev) : ',f10.5,/, . 3x,'br(chic --> psi gamma) : ',f10.5,' +- ',f10.5,/ . 3x,'br(psi --> mu+mu-) : ',f10.5,' +- ',f10.5,/ . 3x,'t slope b_eff (Gev^-2) : ',f10.5,/, . 3x,'delta=alpha_p(0)-1 : ',f10.5,/, . 3x,'S^2_eik : ',f10.5,/, . 3x,'S^2_enh : ',f10.5,/, . 3x,'number of events : ',i9,/, . 3x,'events passing cuts : ',i9,/,/, . 3x,'dsigma/dy|y=',f3.1,' (nb) : ',f10.5,' +- ',f10.5, . /) 661 format(/, . 3x,'collider energy (Gev) : ',f10.3,/, . 3x,'Upsilon mass (Gev) : ',f10.5,/, . 3x,'delta=alpha_p(0)-1 : ',f10.5,/, . 3x,'S^2_enh : ',f10.5,/, . 3x,'number of events : ',i9,/, . 3x,'events passing cuts : ',i9,/, . 3x,'Delta y : ',f10.5,/,/, . 3x,'sigma (pb) : ',f10.5,' +- ',f10.5, . /) 671 format(/, . 3x,'collider energy (Gev) : ',f10.3,/, . 3x,'Upsilon mass (Gev) : ',f10.5,/, . 3x,'delta=alpha_p(0)-1 : ',f10.5,/, . 3x,'S^2_enh : ',f10.5,/, . 3x,'number of events : ',i9,/, . 3x,'events passing cuts : ',i9,/,/, . 3x,'dsigma/dy|y=',f3.1,' (pb) : ',f10.5,' +- ',f10.5, . /) 641 format(/, . 3x,'collider energy (Gev) : ',f10.3,/, . 3x,'chi mass (Gev) : ',f10.5,/, . 3x,'Upsilon mass (Gev) : ',f10.5,/, . 3x,'br(chib --> ups gamma) : ',f10.5,' +- ',f10.5,/ . 3x,'br(eps --> mu+mu-) : ',f10.5,' +- ',f10.5,/ . 3x,'t slope b_eff (Gev^-2) : ',f10.5,/, . 3x,'delta=alpha_p(0)-1 : ',f10.5,/, . 3x,'S^2_eik : ',f10.5,/, . 3x,'S^2_enh : ',f10.5,/, . 3x,'number of events : ',i9,/, . 3x,'events passing cuts : ',i9,/, . 3x,'Delta y : ',f10.5,/,/, . 3x,'sigma (pb) : ',f10.5,' +- ',f10.5, . /) 651 format(/, . 3x,'collider energy (Gev) : ',f10.3,/, . 3x,'chi mass (Gev) : ',f10.5,/, . 3x,'Upsilon mass (Gev) : ',f10.5,/, . 3x,'br(chib --> ups gamma) : ',f10.5,' +- ',f10.5,/ . 3x,'br(ups --> mu+mu-) : ',f10.5,' +- ',f10.5,/ . 3x,'t slope b_eff (Gev^-2) : ',f10.5,/, . 3x,'delta=alpha_p(0)-1 : ',f10.5,/, . 3x,'S^2_eik : ',f10.5,/, . 3x,'S^2_enh : ',f10.5,/, . 3x,'number of events : ',i9,/, . 3x,'events passing cuts : ',i9,/,/, . 3x,'dsigma/dy|y=',f3.1,' (pb) : ',f10.5,' +- ',f10.5, . /) stop end ccccccccccccccccccccccccccccccccccccccccccccc c calculates lorentzian dot product for real 4-vectors double precision function rdot(a,b) double precision a(4),b(4) rdot=a(4)*b(4)-a(3)*b(3)-a(2)*b(2)-a(1)*b(1) return end c calculates lorentzian dot product for real 4-vectors subroutine dot(v1,v2,dt) double precision v1(4),v2(4),dt dt=-(v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3) &-v1(4)*v2(4)) return end c calculates lorentzian dot product for complex 4-vectors subroutine cdot(v1,v2,dt) complex*16 v1(4),v2(4) complex*16 dt dt=-(v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3) &-v1(4)*v2(4)) return end c reads in 'effective luminosity' for digamma CEP subroutine lumw(mgam,wl) double precision wl,mgam,m,v1 double precision warr(55),rm(55) integer o,v data warr/ & 1.88180014,2.03930215,1.54512922,1.20778322,0.946732564, & 0.752739935,0.571403819,0.469902617,0.39174939,0.33073910, & 0.26914951,0.23212474,0.20142495,0.17610557,0.15085951, & 0.13341344,0.11921523,0.10696576,0.09341737,0.08453672, & 0.07673550,0.06986549,0.06377273,0.05831993,0.05358671, & 0.04936697,0.04447949,0.04107049,0.03807317,0.03534359, & 0.03292864,0.03068611,0.02862244,0.02673630,0.02505504, & 0.02346584,0.02203998,0.02073150,0.01932107,0.01817036, & 0.01715060,0.01619606,0.01530502,0.01444515,0.01368000, & 0.01297174,0.01226509,0.01165396,0.01106171,0.00693411, & 0.00447138,0.00313886,0.00214952,0.0015098,0.0010868/ data rm/2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15., & 16.,17.,18.,19.,20.,21.,22.,23.,24.,25.,26.,27.,28.,29., & 30.,31.,32.,33.,34.,35.,36.,37.,38.,39.,40.,41.,42.,43., & 44.,45.,46.,47.,48.,49.,50.,60.,70.,80.,90.,100.,110/ if(mgam.lt.50d0)then v=nint((mgam-2d0)/1d0)+1 v1=((mgam-2d0)/1d0)+1d0 if((v1-dble(v)).lt.0d0)then v=v-1 endif else v=nint((mgam-50d0)/10d0)+49 v1=((mgam-50d0)/10d0)+49d0 if((v1-dble(v)).lt.0d0)then v=v-1 endif endif m=(dlog(warr(v+1))-dlog(warr(v)))/ &(rm(v+1)-rm(v)) wl=m*(mgam-rm(v))+dlog(warr(v)) wl=dexp(wl) return end c calculates survival factor for chi(2++) subroutine chi2sv(ptchi,xf,ef,svout) dimension svarrctev(40),svarrbtev(40),svarrblhc(40),svarrclhc(40) double precision svarr(40) double precision m,svout,ptchi,v1 integer value,xf,j,ef data svarrctev/ & 0.3325825,0.3325477,0.3323151,0.3315647,0.3298286,0.3265279, & 0.3210464,0.3128476,0.3016232,0.2874232,0.2707122,0.2523044, & 0.2331958,0.2143598,0.1965831,0.1803875,0.1660326,0.1535684, & 0.1429038,0.1338680,0.1262562,0.1198595,0.1144814,0.1099459, & 0.1061001,0.1028139,9.9978045E-02,9.7502880E-02,9.5315881E-02, & 9.3360610E-02,9.1595337E-02,8.9991853E-02,8.8534869E-02, & 8.7220058E-02,8.6053327E-02,8.5048161E-02,8.4223129E-02, & 8.3599038E-02,8.3195269E-02,8.3026290E-02/ data svarrbtev/ & 0.3008520,0.3008241,0.3006442,0.3000666,0.2987302,0.2961843, & 0.2919343,0.2855230,0.2766347,0.2652014,0.2514664,0.2359743, & 0.2194756,0.2027799,0.1866139,0.1715269,0.1578583,0.1457568, & 0.1352249,0.1261688,0.1184421,0.1118765,0.1063030,0.1015628, & 9.7514302E-02,9.4034046E-02,9.1017500E-02,8.8378668E-02, & 8.6048037E-02,8.3972082E-02,8.2111478E-02,8.0440089E-02, & 7.8943558E-02,7.7617422E-02,7.6465458E-02,7.5496800E-02, & 7.4724257E-02,7.4160472E-02,7.3815629E-02,7.3694818E-02/ data svarrblhc/ & 0.2026696,0.2026499,0.2025217,0.2021093,0.2011547,0.1993349, & 0.1962965,0.1917116,0.1853539,0.1771736,0.1673443,0.1562550, & 0.1444426,0.1324866,0.1209078,0.1101001,0.1003076,9.1637686E-02, & 8.4093176E-02,7.7607617E-02,7.2076514E-02,6.7379661E-02, & 6.3395485E-02,6.0009770E-02,5.7119429E-02,5.4634426E-02, & 5.2477520E-02,5.0584164E-02,4.8901454E-02,4.7387592E-02, & 4.6011236E-02,4.4750430E-02,4.3592654E-02,4.2533256E-02, & 4.1574374E-02,4.0723760E-02,3.9992884E-02,3.9394580E-02, & 3.8941417E-02,3.8643118E-02/ data svarrclhc/ & 0.2296173,0.2295918,0.2294217,0.2288716,0.2275986,0.2251778, & 0.2211568,0.2151411,0.2069038,0.1964808,0.1842119,0.1706945, & 0.1566596,0.1428219,0.1297599,0.1178574,0.1073063,9.8144226E-02, & 9.0305425E-02,8.3665058E-02,7.8073338E-02,7.3377274E-02, & 6.9432534E-02,6.6109546E-02,6.3295096E-02,6.0892198E-02, & 5.8818642E-02,5.7005975E-02,5.5397697E-02,5.3948846E-02, & 5.2624661E-02,5.1400553E-02,5.0261818E-02,4.9202461E-02, & 4.8224945E-02,4.7338821E-02,4.6558995E-02,4.5903802E-02, & 4.5392185E-02,4.5041133E-02/ if(xf.eq.0)then if(ef.eq.0)then do j=1,40 svarr(j)=svarrctev(j) enddo elseif(ef.eq.1)then do j=1,40 svarr(j)=svarrclhc(j) enddo endif elseif(xf.eq.1)then if(ef.eq.0)then do j=1,40 svarr(j)=svarrbtev(j) enddo elseif(ef.eq.1)then do j=1,40 svarr(j)=svarrblhc(j) enddo endif endif if(ptchi.lt.2.5d-2)then value=1 m=svarr(1)/2.5d-2 svout=m*ptchi elseif(ptchi.lt.1.975d0)then value=nint(ptchi/0.05d0) v1=ptchi/0.05d0 if((v1-dble(value)).lt.0d0)then value=value-1 endif m=(svarr(value+1)-svarr(value))/0.05d0 svout=m*(ptchi-0.05d0*value-0.025d0)+svarr(value+1) else svout=svarr(40) endif 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 binning subroutine subroutine binit(wt) implicit double precision(a-y) double precision q(4,20),pt1(2),pt2(2),ptx(2) common/mom/q common/mompt/pt1,pt2,ptx common/vars/s,rts,mchi,ptcut,etacut,ecut,rcut pt1sqr=((q(1,3)**2+q(2,3)**2)) pt2sqr=(q(1,4)**2+q(2,4)**2) ptpsi=dsqrt(q(1,7)**2+q(2,7)**2) ptchi=dsqrt(q(1,5)**2+q(2,5)**2) ptmu1=dsqrt(q(1,8)**2+q(2,8)**2) ptmu2=dsqrt(q(1,9)**2+q(2,9)**2) ptgam=dsqrt(q(1,6)**2+q(2,6)**2) c call histo1(1,10,0d0,1d0,dabs(costst),wt) call histo1(1,40,0d0,1.6d0,ptchi,wt) c call histo1(1,30,0d0,6d0,ptmu1,wt) c call histo1(1,30,0d0,2.3d0,ptmu2,wt) c call histo1(1,30,0d0,1.5d0,ptpsi,wt) c call histo1(1,30,0d0,0.7d0,ptgam,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,ptcut,etacut,ecut,rcut icut=1 c -- insert cuts here if desired... c for example, a minimum energy cut on the final state leptons c if(q(4,8).lt.1d0) return c if(q(4,9).lt.1d0) return c chi_{c/b) cuts ccc pt(muons) ptmu1=dsqrt(q(1,8)**2+q(2,8)**2) ptmu2=dsqrt(q(1,9)**2+q(2,9)**2) 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) c pseudorapidities of muons c eta1=dabs(0.5d0*dlog((pmod1+q(3,8))/(pmod1-q(3,8)))) c eta2=dabs(0.5d0*dlog((pmod2+q(3,9))/(pmod2-q(3,9)))) c if(ptmu1.gt.ptcut.and.ptmu2.gt.ptcut.and.eta1.lt.etacut.and.eta2 c &.lt.etacut)return ccc gamgam cuts pmod1=dsqrt(q(1,12)**2+q(2,12)**2+q(3,12)**2) pmod2=dsqrt(q(1,13)**2+q(2,13)**2+q(3,13)**2) eta1=dabs(0.5d0*dlog((pmod1+q(3,12))/(pmod1-q(3,12)))) eta2=dabs(0.5d0*dlog((pmod2+q(3,13))/(pmod2-q(3,13)))) et1=dsqrt(q(4,12)**2-q(3,12)**2) et2=dsqrt(q(4,13)**2-q(3,13)**2) et=et1+et2 c et=dsqrt(q(1,5)**2+q(2,5)**2) c ecut=0d0 if(et1.gt.ecut.and.et2.gt.ecut.and.eta1.lt.rcut.and.eta2. <.rcut)return cccccccccccccccccccccccccc icut=0 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 c open(50,file="output.dat") 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 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 c close(50) 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 * * 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