bin/BCVEGPY/totfun.F

0001
0002 c**********************************************************************

0003
0004       double precision function totfun(zup,papawt)
0005         implicit double precision (a-h,o-z)
0006         implicit integer (i-n)
0007
0008 #include "inclcon.h"
0009
0010 c...pythia common block.

0011       common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
0012         common/pypars/mstp(200),parp(200),msti(200),pari(200)
0013         parameter (maxnup=500)
0014       common/hepeup/nup,idprup,xwgtup,scalup,aqedup,aqcdup,idup(maxnup),
0015      &istup(maxnup),mothup(2,maxnup),icolup(2,maxnup),pup(5,maxnup),
0016      &vtimup(maxnup),spinup(maxnup)
0017       save /hepeup/
0018
0019 c...user transformation.

0020         double complex colmat,bundamp
0021         common/upcom/ecm,pmbc,pmb,pmc,fbcc,pmomup(5,8),
0022      &  colmat(10,64),bundamp(4),pmomzero(5,8)
0023 c...transform the bound state information.

0024         common/counter/ibcstate,nev
0025         common/rconst/pi
0026         common/ptpass/ptmin,ptmax,crossmax,etamin,etamax,
0027      &  smin,smax,ymin,ymax,psetamin,psetamax
0028         common/confine/ptcut,etacut
0029         common/colflow/amp2cf(10),smatval
0030
0031 c...parameters transformtion used in totfun()

0032       common/funtrans/nq2,npdfu
0033
0034 c...to get the subprocess cross-section.

0035       common/subopen/subfactor,subenergy,isubonly
0036
0037 c...generate---switch for full events.

0038         logical generate
0039         common/genefull/generate
0040
0041 c...to get the distribution of an extra factor z=(2(k1+k2).p_bc)/shat.

0042       common/extraz/zfactor,zmin,zmax
0043         common/outpdf/ioutpdf,ipdfnum
0044         common/intinip/iinip
0045         common/intinif/iinif
0046
0047 c...for transform the subprocess information, i.e.,  whether using

0048 c...the subprocess q\bar{q}->bc+b+\bar{c} to generate events.

0049       common/qqbar/iqqbar,iqcode
0050
0051       dimension xpp(-25:25),xppbar(-25:25),zup(7),pboo(4),pc(4),pl(4)
0052       data conv/3.8938573d+8/ !pb
0053       common /ppp/ pp(4,40),guv(4)
0054
0055 c------------------------------------------------

0056
0057       totfun=0.0d0
0058       phase = 0.
0059
0060 c------------------------------------------------

0061       if(isubonly.eq.1) then
0062          x1=subfactor
0063            x2=subfactor
0064         else
0065          taumin =((pmbc+pmb+pmc)/ecm)**2
0066            taumax =1.0d0
0067            tau=(taumax-taumin)*zup(6)+taumin
0068            yymin= dlog(dsqrt(tau))
0069            yymax=-dlog(dsqrt(tau))
0070            yy   =(yymax-yymin)*zup(7)+yymin
0071            x1 =dsqrt(tau)*exp(yy)
0072            x2 =dsqrt(tau)*exp(-yy)
0073         end if
0074
0075 c-------------------------------------------------

0076
0077 c... gluon 1, in lab

0078       pup(1,1)= 0.0d0
0079       pup(2,1)= 0.0d0
0080       pup(3,1)= ecm*x1/2.0d0
0081       pup(4,1)= ecm*x1/2.0d0
0082       pup(5,1)= 0.0d0
0083
0084 c... gluon 2, in lab

0085       pup(1,2)= 0.0d0
0086       pup(2,2)= 0.0d0
0087       pup(3,2)=-ecm*x2/2.0d0
0088       pup(4,2)= ecm*x2/2.0d0
0089       pup(5,2)= 0.0d0
0090
0091 c...change momtuma of the final particals into lab coordinate system.

0092 c...the original one getting from phase_gen is result for c.m. system.

0093         do i=1,4
0094           pboo(i)=pup(i,1)+pup(i,2)
0095       end do
0096
0097       do 101, i=1,3
0098             do j=1,4
0099               pc(j)=pmomup(j,i+2)
0100           end do
0101                 call lorentz(pboo,pc,pl)
0102             do j=1,4
0103               pmomup(j,i+2)=pl(j)
0104           end do
0105 101   continue
0106
0107 c...set up kinematics of the out going particles: bc, b and c~.

0108       do i=3,5
0109           do j=1,5
0110             pup(j,i)=pmomup(j,i)
0111           end do
0112         end do
0113
0114 c...for s-wave bound state, taking non-relativistic approximation.

0115         do i=1,5
0116           pmomup(i,6)=pmb/(pmb+pmc)*pmomup(i,3)
0117           pmomup(i,7)=pmc/(pmb+pmc)*pmomup(i,3)
0118       end do
0119
0120 c...incoming gluon momenta used in s_amp.for.

0121         do i=1,2
0122           do j=1,5
0123             pmomup(j,i)=pup(j,i)
0124           end do
0125         end do
0126
0127 c...this part is from the inner part of pythia subroutine pyp().

0128       ptbc  =dsqrt(pup(1,3)**2+pup(2,3)**2)
0129       pr    =max(1.0d-16,pup(5,3)**2+pup(1,3)**2+pup(2,3)**2)
0130         prs   =max(1.0d-16,pup(1,3)**2+pup(2,3)**2)
0131       eta   =sign(dlog(min((dsqrt(pr+pup(3,3)**2)+dabs(pup(3,3)))
0132      &          /dsqrt(pr),1.0d+20)),pup(3,3))
0133         pseta =sign(dlog(min((dsqrt(prs+pup(3,3)**2)+dabs(pup(3,3)))
0134      &          /dsqrt(prs),1.0d+20)),pup(3,3))
0135
0136 c...other confinement can also be added here.

0137         if(ptbc.lt.ptcut .or. abs(eta).gt.etacut) then
0138           if (generate) then
0139             do ii=1,10
0140               amp2cf(ii)=0.0d0
0141             end do
0142         end if
0143           smatval=0.0d0
0144           return
0145         end if
0146
0147 c...energy scale.

0148         if(nq2.eq.1) q2 =x1*x2*ecm**2/4.0d0
0149       if(nq2.eq.2) q2 =x1*x2*ecm**2
0150         if(nq2.eq.3.or.nq2.eq.8) q2 =ptbc**2.0d0+pup(5,3)**2
0151         if(nq2.eq.4) then
0152           q=0.0d0
0153           do i=3,5
0154             q=q+dsqrt(pup(1,i)**2+pup(2,i)**2+pup(5,i)**2)
0155           end do
0156           q2=q**2
0157         end if
0158         if(nq2.eq.5) then
0159           q=0.0d0
0160           do i=3,5
0161             q=q+dsqrt(pup(1,i)**2.0d0+pup(2,i)**2.0d0+pup(5,i)**2.0d0)
0162           end do
0163           q2=(q/3.0d0)**2.0d0
0164         end if
0165         if(nq2.eq.6.or.nq2.eq.7) then
0166           q2=pmb**2+pup(1,4)**2+pup(2,4)**2
0167         end if
0168 c...this is the energy scale used in gouz's program

0169         if(nq2.eq.9) then
0170           q2=4.0d0*pmb**2
0171         end if
0172         alps  = 0.00
0173         alps2 = 0.00
0174 c...get the value of alphas. all are in leading order.

0175         if(ioutpdf.eq.1) then
0176            if(ipdfnum.eq.100) alps =alpgrv(q2,1)*4*pi
0177            if(ipdfnum.eq.200) alps =alpmsrt(dsqrt(q2),0.220d0,0)
0178            if(ipdfnum.eq.300) alps =alpcteq(q2,1)*4*pi
0179         else
0180            alps =pyalps(q2)
0181       end if
0182
0183 c...two energy scale for alphas.

0184 c...alphas^4=alphas^2(\mu_b**2)*alphas^2(\mu_c**2).

0185         if(nq2.eq.6.or.nq2.eq.8) then
0186            alps1=alps
0187            if(ioutpdf.eq.1) then
0188               q2=4.0d0*pmc**2.0d0
0189                   if(ipdfnum.eq.100) alps2 =alpgrv(q2,1)*4*pi
0190               if(ipdfnum.eq.200) alps2 =alpmsrt(dsqrt(q2),0.22d0,0)
0191               if(ipdfnum.eq.300) alps2 =alpcteq(q2,1)*4*pi
0192            else
0193               alps2 =pyalps(4.0d0*pmc**2.0d0)
0194          end if
0195            alps =dsqrt(alps1*alps2)
0196         end if
0197
0198 c...store scale choice and alphas.

0199       scalup=dsqrt(q2)
0200       aqcdup=alps
0201
0202         if(isubonly.eq.0) then
0203          if(ioutpdf.eq.0) then
0204 c...tevatron

0205 c...evaluate parton distribution for (g1<--p,g2<---p~)

0206            if(npdfu.eq.1) then
0207              call pypdfu(2212,x1,q2,xpp)
0208            call pypdfu(-2212,x2,q2,xppbar)
0209            end if
0210 c...lhc

0211 c...evaluate parton distribution for (g1<--p,g2<---p)

0212            if(npdfu.eq.2) then
0213              call pypdfu(2212,x1,q2,xpp)
0214            call pypdfu(2212,x2,q2,xppbar)
0215            end if
0216          else
0217            if(ipdfnum.eq.100) then
0218              call grv98pa(1, x1, q2, uv, dv, us, ds, ss, gl1)
0219            call grv98pa(1, x2, q2, uv, dv, us, ds, ss, gl2)
0220              if(iqqbar.eq.0) then
0221                    xpp(21)=gl1
0222                    xppbar(21)=gl2
0223              else
0224                if(iqcode.eq.1) then
0225                  if(npdfu.eq.1) then     !tevatron
0226                             xpp(iqcode)=uv+us
0227                     xpp(-iqcode)=us
0228                             xppbar(iqcode)=uv+us
0229                             xppbar(-iqcode)=us
0230                      end if
0231                  if(npdfu.eq.2) then     !lhc
0232                             xpp(iqcode)=uv+us
0233                     xpp(-iqcode)=us
0234                             xppbar(iqcode)=us
0235                             xppbar(-iqcode)=uv+us
0236                      end if
0237                    end if
0238                    if(iqcode.eq.2) then
0239                  if(npdfu.eq.1) then     !tevatron (u-p,~u-~p and u-~p,~u-p)
0240                             xpp(iqcode)=dv+ds
0241                     xpp(-iqcode)=ds
0242                             xppbar(iqcode)=dv+ds
0243                             xppbar(-iqcode)=ds
0244                      end if
0245                  if(npdfu.eq.2) then     !lhc (u-p,~u-p and u-p,~u-p)
0246                             xpp(iqcode)=dv+ds
0247                     xpp(-iqcode)=ds
0248                             xppbar(iqcode)=ds
0249                             xppbar(-iqcode)=dv+ds
0250                      end if
0251                    end if
0252                    if(iqcode.eq.3) then      ! the same for tevatron or lhc
0253                      xpp(3)=2.0d0*ss
0254                      xppbar(3)=2.0d0*ss
0255                    end if
0256            end if
0257            end if
0258            if(ipdfnum.eq.200) then
0259              qq=dsqrt(q2)
0260                  call mrstlo(x1,qq,1,upv,dnv,usea,dsea,str,chm,bot,glu1)
0261              call mrstlo(x2,qq,1,upv,dnv,usea,dsea,str,chm,bot,glu2)
0262              if(iqqbar.eq.0) then
0263                    xpp(21)=glu1
0264                    xppbar(21)=glu2
0265              else
0266                if(iqcode.eq.1) then
0267                  if(npdfu.eq.1) then     !tevatron (u-p,~u-~p and u-~p,~u-p)
0268                             xpp(iqcode)=upv+usea
0269                     xpp(-iqcode)=usea
0270                             xppbar(iqcode)=upv+usea
0271                             xppbar(-iqcode)=usea
0272                      end if
0273                  if(npdfu.eq.2) then     !lhc (u-p,~u-p and u-p,~u-p)
0274                             xpp(iqcode)=upv+usea
0275                     xpp(-iqcode)=usea
0276                             xppbar(iqcode)=usea
0277                             xppbar(-iqcode)=upv+usea
0278                      end if
0279                    end if
0280                    if(iqcode.eq.2) then
0281                  if(npdfu.eq.1) then     !tevatron (u-p,~u-~p and u-~p,~u-p)
0282                             xpp(iqcode)=dnv+dsea
0283                     xpp(-iqcode)=dsea
0284                             xppbar(iqcode)=dnv+dsea
0285                             xppbar(-iqcode)=dsea
0286                      end if
0287                  if(npdfu.eq.2) then     !lhc (u-p,~u-p and u-p,~u-p)
0288                             xpp(iqcode)=dnv+dsea
0289                     xpp(-iqcode)=dsea
0290                             xppbar(iqcode)=dsea
0291                             xppbar(-iqcode)=dnv+dsea
0292                      end if
0293                    end if
0294                    if(iqcode.eq.3) then
0295                      xpp(3)=str*2.0d0
0296                      xppbar(3)=str*2.0d0
0297              end if
0298                  end if
0299            end if
0300            if(ipdfnum.eq.300) then
0301                  qq=dsqrt(q2)
0302 c...cteq6l.

0303              if(iqqbar.eq.0) then
0304                xpp(21)   =ctq6pdf(0,x1,qq)
0305                    xppbar(21)=ctq6pdf(0,x2,qq)
0306              else
0307                if(npdfu.eq.1) then         !tevatron
0308                      xpp(iqcode)=ctq6pdf(iqcode,x1,qq)
0309                      xpp(-iqcode)=ctq6pdf(-iqcode,x1,qq)
0310                      xppbar(iqcode)=ctq6pdf(iqcode,x2,qq)
0311                      xppbar(-iqcode)=ctq6pdf(-iqcode,x2,qq)
0312                end if
0313                if(npdfu.eq.2) then         !lhc
0314                      xpp(iqcode)=ctq6pdf(iqcode,x1,qq)
0315                      xpp(-iqcode)=ctq6pdf(-iqcode,x1,qq)
0316                      xppbar(iqcode)=ctq6pdf(-iqcode,x2,qq)
0317                      xppbar(-iqcode)=ctq6pdf(iqcode,x2,qq)
0318                end if
0319                  end if
0320            end if
0321          end if
0322         end if
0323
0324 c...this ensure the rightness of the extrapolation of the pdf.

0325 c...(by using the pdfs, sometimes it will get negative value)

0326         if(xpp(21).lt.1.0d-16) xpp(21)=0.0d0
0327         if(xppbar(21).lt.1.0d-16) xppbar(21)=0.0d0
0328
0329         if(iqqbar.eq.1) then
0330          if(xpp(iqcode).lt.1.0d-16) xpp(iqcode)=0.0d0
0331          if(xppbar(iqcode).lt.1.0d-16) xppbar(iqcode)=0.0d0
0332          if(xpp(-iqcode).lt.1.0d-16) xpp(-iqcode)=0.0d0
0333          if(xppbar(-iqcode).lt.1.0d-16) xppbar(-iqcode)=0.0d0
0334         end if
0335
0336 c...for the sub-process, taking the constant alphas value.

0337         if(isubonly.eq.1) alps=0.20d0
0338
0339 c...if not generate s-wave, go to the part for p-wave.

0340         if(ibcstate.eq.3) goto 1005
0341         if(ibcstate.eq.4) goto 1005
0342         if(ibcstate.eq.5) goto 1005
0343         if(ibcstate.eq.6) goto 1005
0344
0345 c...common factor for s-wave states.

0346         phase =papawt*alps**4.0d0/(2.0d0**11*pi*3.0d0*dotup(1,2))
0347
0348 c...first to get the square of the amplitude.

0349 c...note 1) the momenta inputed into this subroutine is pmomup(j,i):

0350 c...bc+: i=3, b: i=4, ~c: i=5;

0351 c...j=1: p_x; j=2: p_y; j=3: p_z; j=4: e; j=5: mass;    2) all the momenta

0352 c...now are in the lab system, you may directly get the momenta in c.m.s

0353 c...of gluon-gluon subsystem before running lorentz transformation used

0354 c...above. the cross-section will not change under the cordinate

0355 c...transformation. 3) sigscl and sigvct, after calling this subroutine,

0356 c...is not the final cross-section, we need to add some coefficients.

0357 c...there we don't need an extra q in xsection_bcy, because 4) the momenta

0358 c...are now also stored in pup(j,i), which is transformed according to

0359 c...a pythia common block.

0360         call xsection(sigscl,sigvct)
0361
0362
0363 c...get the right phase for the subprocess: q+~q->bc+b+~c from

0364 c...the same precedure of subprocess g+g->bc+b+~c.

0365         if(iqqbar.eq.1) then
0366           phase=phase*(2.0d0**6/3.0d0**2)
0367         end if
0368
0369 c...the correct cross-section for the subprocess to a particular particle

0370 c...momenta.

0371         sigscl=conv*phase*sigscl
0372         sigvct=conv*phase*sigvct
0373         sigcross = 0.0
0374         if(ibcstate.eq.1) sigcross=sigscl
0375         if(ibcstate.eq.2) sigcross=sigvct
0376
0377 c...get the cross-section for s-wave.

0378         if(isubonly.eq.0) then
0379           if(iqqbar.eq.0) then
0380            totfun =sigcross*xpp(21)*xppbar(21)/x1/x2
0381            if(ioutpdf.eq.1 .and. ipdfnum.eq.300) then
0382             totfun=sigcross*xpp(21)*xppbar(21)
0383            end if
0384           else
0385                 totfun =sigcross*(xpp(iqcode)*xppbar(iqcode)+
0386      &          xpp(-iqcode)*xppbar(-iqcode))/x1/x2
0387           if(ioutpdf.eq.1 .and. ipdfnum.eq.300) then
0388               totfun =sigcross*xpp(iqcode)*xppbar(iqcode)
0389             end if
0390           end if
0391         else
0392           totfun =sigcross
0393         end if
0394
0395 c...getting an extra distribution about z=(2(k1+k2).p_bc)/shat.

0396       zfactor=2.0d0*(dotup(1,3)+dotup(2,3))/(x1*x2*ecm**2.0d0)
0397
0398 c...the following is only to eliminate the numerical uncerntainty,

0399 c...which in principle does not needed. however we added here

0400 c...to avoid some very particular cases.

0401         if(totfun.lt.1.0d-16) totfun=1.0d-16
0402
0403         return
0404
0405 c************************************************

0406 c...the following is for the p-wave states.

0407 c************************************************

0408
0409 1005  continue
0410
0411       if(ibcstate.ne.1.and.ibcstate.ne.2) then
0412 c...the p-wave part is adopted from fdc program and here

0413 c...is to make all the variable consistent to fdc.

0414 c...in the lab system.

0415 c...pp(1,i)=e=pup(4,i);

0416 c...pp((2,3,4),i)=(p_x,p_y,p_z)=pup((1,2,3),i).

0417 c...pp(4)=pup(5)--b quark; pp(5)=pup(4)--\bar{c} quark.

0418          do i=1,3
0419            pp(2,i)=pup(1,i)
0420            pp(3,i)=pup(2,i)
0421            pp(4,i)=pup(3,i)
0422            pp(1,i)=pup(4,i)
0423          end do
0424
0425          pp(2,4)=pup(1,5)
0426          pp(3,4)=pup(2,5)
0427          pp(4,4)=pup(3,5)
0428          pp(1,4)=pup(4,5)
0429
0430          pp(2,5)=pup(1,4)
0431          pp(3,5)=pup(2,4)
0432          pp(4,5)=pup(3,4)
0433          pp(1,5)=pup(4,4)
0434
0435 c...this is an overall color factor, where the average over

0436 c...the initial gluons have been included, i.e. the factor

0437 c...(1/2**6) has been included into cfact, i.e.

0438 c...   cfact=(33.0d0/2.0d0)/(2.0d0**6)

0439          cfact=0.2578125d0
0440
0441 c...this is to get constants for the gluon coupling constant g.

0442          g3=dsqrt(4.0d0*pi*alps)
0443         ampofpw = 0.0
0444          if(ibcstate.eq.3) then
0445             ampofpw=amps2_1p1()
0446             phase=papawt*cfact*g3**8/(2.0d0**9*pi**5*dotup(1,2))
0447          end if
0448          if(ibcstate.eq.4) then
0449             ampofpw=amps2_3p0()
0450             phase=papawt*cfact*g3**8/(2.0d0**9*pi**5*dotup(1,2))
0451          end if
0452          if(ibcstate.eq.5) then
0453             ampofpw=amps2_3p1()
0454             phase=papawt*cfact*g3**8/(2.0d0**9*pi**5*dotup(1,2))
0455          end if
0456          if(ibcstate.eq.6) then
0457             ampofpw=amps2_3p2()
0458             phase=papawt*cfact*g3**8/(2.0d0**9*pi**5*dotup(1,2))
0459          end if
0460 c...get the cross-sectin for p-wave.

0461          if(isubonly.eq.1) then
0462            totfun=conv*phase*ampofpw
0463          else
0464            totfun =conv*phase*ampofpw*xpp(21)*xppbar(21)/x1/x2
0465            if(ioutpdf.eq.1 .and. ipdfnum.eq.300) then
0466              totfun=conv*phase*ampofpw*xpp(21)*xppbar(21)
0467            end if
0468          end if
0469         end if
0470
0471 c...getting an extra distribution about z=(2(k1+k2).p_bc)/shat.

0472       zfactor=2.0d0*(dotup(1,3)+dotup(2,3))/(x1*x2*ecm**2.0d0)
0473
0474 c...the following is only to eliminate the numerical uncerntainty,

0475 c...which in principle does not needed. however we added here

0476 c...to avoid some very particular cases.

0477         if(totfun.lt.1.0d-16) totfun=1.0d-16
0478
0479         return
0480         end