bin/BCVEGPY/s_free.F

0001 c**********************************************************************

0002
0003 c...helicity amplitude of the free quark part: g+g->\bar{b}+c+b+\bar{c}.

0004
0005         subroutine freehelicity
0006         implicit double precision (a-h,o-z)
0007         implicit integer (i-n)
0008         double complex bfun,colmat,inpup,sppup,yup,xup,bundamp
0009       common/upcom/ecm,pmbc,pmb,pmc,fbcc,pmomup(5,8),
0010      &  colmat(10,64),bundamp(4),pmomzero(5,8)
0011       common/tryup/bfun(9,4,64),yup(100),xup(2),idp,idq0,idk1,idk2,
0012      &  idb1,idb2,idc1,idc2
0013 c...for transform the subprocess information, i.e.,  whether using

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

0015       common/qqbar/iqqbar,iqcode
0016 c...whether getting the color-octet component contributions. here only

0017 c...for gg->bc+b+~c (bc in color-octet s-wave states.).

0018       common/coloct/ioctet
0019
0020       dimension propup(14,4)
0021
0022 c...the four subroutine are used to get bfun(). first->original from cb and

0023 c...cc; second->inter-change k_1 and k_2; third->inter-change q_b1,q_b2 and

0024 c...q_c1,q_c2; four->inter-change both. here, bfun() stands for the basic

0025 c...feynman diagrams which is free from color factors and proporgators.

0026 c...the results are get from mathmatical programs: helicitybasic.nb

0027 c...and texbasic.nb. here for easy programming, we use the number idk1,idk2,

0028 c...idb1,idb2,idc1,idc2,idq0,idp to identify the particles correspondingly.

0029
0030       idq0=8
0031         idp=3
0032
0033       do 700, ii=1,4
0034          if (ii .eq. 1) then
0035           idk1=1
0036           idk2=2
0037           idb1=4
0038           idb2=6
0039           idc1=7
0040           idc2=5
0041          end if
0042
0043          if (ii .eq. 2) then
0044           idk1=2
0045           idk2=1
0046           idb1=4
0047           idb2=6
0048           idc1=7
0049           idc2=5
0050        end if
0051
0052        if (ii .eq. 3) then
0053           idk1=1
0054           idk2=2
0055           idb1=7
0056           idb2=5
0057           idc1=4
0058           idc2=6
0059        end if
0060
0061        if (ii .eq. 4) then
0062           idk1=2
0063           idk2=1
0064           idb1=7
0065           idb2=5
0066           idc1=4
0067           idc2=6
0068        end if
0069
0070 c...define some basic spinnor products that are needed in helicity amplitude.

0071          yup(1)=sppup(idk1,idc2,idq0)
0072          yup(2)=sppup(idk1,idb2,idq0)
0073          yup(3)=sppup(idk1,idq0,idb2)
0074          yup(4)=sppup(idk1,idq0,idc2)
0075          yup(5)=sppup(idk1,idq0,idc1)
0076          yup(6)=sppup(idk1,idq0,idb1)
0077          yup(7)=sppup(idk1,idk2,idq0)
0078          yup(8)=sppup(idk1,idq0,idk2)
0079          yup(9)=sppup(idk1,idc1,idq0)
0080
0081          yup(10)=sppup(idk2,idb1,idq0)
0082          yup(11)=sppup(idk2,idc2,idq0)
0083          yup(12)=sppup(idk2,idk1,idq0)
0084          yup(13)=sppup(idk2,idq0,idb2)
0085          yup(14)=sppup(idk2,idq0,idc2)
0086          yup(15)=sppup(idk2,idq0,idc1)
0087          yup(16)=sppup(idk2,idb2,idq0)
0088          yup(17)=sppup(idk2,idq0,idb1)
0089          yup(18)=sppup(idk2,idc1,idq0)
0090
0091          yup(19)=sppup(idq0,idk2,idq0)
0092          yup(20)=sppup(idq0,idb1,idc1)
0093          yup(21)=sppup(idq0,idk1,idb1)
0094          yup(22)=sppup(idq0,idc2,idb1)
0095          yup(23)=sppup(idq0,idb1,idc2)
0096          yup(24)=sppup(idq0,idb2,idc2)
0097          yup(25)=sppup(idq0,idb1,idk1)
0098          yup(26)=sppup(idq0,idk2,idb1)
0099          yup(27)=sppup(idq0,idc1,idq0)
0100          yup(28)=sppup(idq0,idc2,idq0)
0101          yup(29)=sppup(idq0,idb1,idq0)
0102          yup(30)=sppup(idq0,idb2,idq0)
0103          yup(31)=sppup(idq0,idk1,idq0)
0104          yup(32)=sppup(idq0,idc1,idb1)
0105          yup(33)=sppup(idq0,idk2,idc2)
0106
0107          yup(34)=sppup(idc2,idq0,idc1)
0108          yup(35)=sppup(idb2,idq0,idb1)
0109          yup(36)=sppup(idb2,idc2,idq0)
0110          yup(37)=sppup(idb1,idq0,idc2)
0111          yup(38)=sppup(idq0,idk1,idb2)
0112          yup(39)=sppup(idq0,idb2,idc1)
0113          yup(40)=sppup(idq0,idk2,idc1)
0114          yup(41)=sppup(idq0,idk2,idb2)
0115          yup(42)=sppup(idq0,idc1,idb2)
0116          xup(1)=inpup(idq0,idk1)*inpup(idq0,idk2)
0117          xup(2)=dconjg(xup(1))
0118
0119          if(iqqbar.eq.1) then
0120         yup(43)=sppup(idb1,idq0,idc1)
0121           yup(44)=sppup(idc2,idb1,idk2)
0122         yup(45)=sppup(idc2,idk1,idk2)
0123         yup(46)=sppup(idk1,idb1,idk2)
0124           yup(47)=sppup(idk1,idk2,idb1)
0125         yup(48)=sppup(idb2,idk1,idk2)
0126           yup(49)=sppup(idk1,idb2,idc1)
0127           yup(50)=sppup(idk1,idb2,idk2)
0128
0129           yup(51)=sppup(idb2,idk1,idb1)
0130         yup(52)=sppup(idc2,idk2,idb1)
0131           yup(53)=sppup(idk2,idb1,idk1)
0132           yup(54)=sppup(idk2,idk1,idc1)
0133           yup(55)=sppup(idc2,idb1,idk1)
0134           yup(56)=sppup(idc2,idk2,idk1)
0135           yup(57)=sppup(idk2,idk1,idb1)
0136           yup(58)=sppup(idk2,idb2,idc1)
0137           yup(59)=sppup(idb2,idk2,idk1)
0138           yup(60)=sppup(idk2,idb2,idk1)
0139
0140           yup(61)=sppup(idb2,idb1,idq0)
0141           yup(62)=sppup(idb2,idb1,idc1)
0142           yup(63)=sppup(idb2,idk2,idc1)
0143           yup(64)=sppup(idb2,idb1,idc2)
0144           yup(65)=sppup(idc1,idq0,idb1)
0145           yup(66)=sppup(idk1,idk2,idb2)
0146           yup(67)=sppup(idc1,idb2,idk2)
0147           yup(68)=sppup(idc2,idb2,idk2)
0148           yup(69)=sppup(idb2,idq0,idc1)
0149           yup(70)=sppup(idc1,idb1,idb2)
0150
0151           yup(71)=sppup(idb2,idq0,idc2)
0152           yup(72)=sppup(idc1,idk2,idk1)
0153           yup(73)=sppup(idc2,idb2,idk1)
0154           yup(74)=sppup(idc1,idb1,idk2)
0155           yup(75)=sppup(idq0,idb2,idb1)
0156           yup(76)=sppup(idc1,idk2,idb1)
0157           yup(77)=sppup(idb2,idk2,idc2)
0158           yup(78)=sppup(idc2,idk2,idb2)
0159           yup(79)=sppup(idq0,idk1,idc1)
0160           yup(80)=sppup(idq0,idk1,idc2)
0161
0162           yup(81)=sppup(idc1,idb1,idk1)
0163           yup(82)=sppup(idk1,idb2,idb1)
0164           yup(83)=sppup(idc2,idb2,idb1)
0165           yup(84)=sppup(idc1,idb2,idb1)
0166           yup(85)=sppup(idc2,idb1,idb2)
0167           yup(86)=sppup(idk1,idb2,idc2)
0168          end if
0169
0170        if(iqqbar.eq.0) then
0171           if (ii.eq.1) call bfirst
0172           if (ii.eq.2) call bsecond
0173           if (ii.eq.3) call bthird
0174           if (ii.eq.4) call bfourth
0175          end if
0176
0177          if(iqqbar.eq.1) then
0178           if (ii.eq.1) then
0179             do kk=1,64
0180                   bfun(1,ii,kk)=0.0d0
0181             bfun(2,ii,kk)=0.0d0
0182               bfun(3,ii,kk)=0.0d0
0183               bfun(4,ii,kk)=0.0d0
0184             end do
0185                 call bfirst
0186           end if
0187           if (ii.eq.3) then
0188             do kk=1,64
0189                   bfun(1,ii,kk)=0.0d0
0190             bfun(2,ii,kk)=0.0d0
0191               bfun(3,ii,kk)=0.0d0
0192               bfun(4,ii,kk)=0.0d0
0193             end do
0194                 call bthird
0195           end if
0196          end if
0197
0198 c...the needed basic propagators. propup(n,1) from cb and cc sets. other

0199 c...propagators can be obtained from them by interchange gluons and quarks,

0200 c...which can be achieved by interchanging the identity number of the

0201 c...particles (i.e. idk1,idk2,...). here the correspondence of the diagram

0202 c...labels, mcb1,mcb2,..., to the original feynman diagrams are shown in paper.

0203        if(iqqbar.eq.0) then
0204 c...propup(1,1)<-mcb1,propup(2,1)<-mcb2,propup(3,1)<-mcb3,propup(4,1)<-mcb4

0205        if (ii.eq.1 .or. ii.eq.2) then
0206          propup(1,ii)=1/(-2*dotup(idb1,idk2))/(2*dotup(idc1,idc1)
0207      &   +2*dotup(idc2,idc1)-2*dotup(idk1,idc2)-2*dotup(idk1,idc1))
0208      &   /(-2*dotup(idk1,idc2))
0209            propup(2,ii)=1/(-2*dotup(idb1,idk2))/(2*dotup(idc1,idc1)
0210      &   +2*dotup(idc2,idc1)-2*dotup(idk1,idc2)-2*dotup(idk1,idc1))
0211      &   /(-2*dotup(idc1,idk1))
0212            propup(3,ii)=1/(-2*dotup(idb2,idk2))/(2*dotup(idc1,idc1)
0213      &   +2*dotup(idc2,idc1)-2*dotup(idk1,idc2)-
0214      &     2*dotup(idk1,idc1))/(-2*dotup(idc2,idk1))
0215            propup(4,ii)=1/(-2*dotup(idb2,idk2))/(2*dotup(idc1,idc1)
0216      &   +2*dotup(idc2,idc1)-2*dotup(idk1,idc2)-
0217      &   2*dotup(idk1,idc1))/(-2*dotup(idc1,idk1))
0218          end if
0219
0220 c...propup(5,1)->mcc1,propup(6,1)->mcc3,propup(7,1)->mcc5

0221          propup(5,ii)=1/(2*dotup(idb1,idb1)+2*dotup(idb1,idb2))
0222      &   /(-2*dotup(idc2,idk2))/(2*dotup(idk1,idk2)-
0223      &   2*dotup(idk1,idc2)-2*dotup(idc2,idk2))
0224          propup(6,ii)=1/(2*dotup(idb1,idb1)+2*dotup(idb1,idb2))
0225      &   /(-2*dotup(idc2,idk2))/(-2*dotup(idc1,idk1))
0226          propup(7,ii)=1/(2*dotup(idb1,idb1)+2*dotup(idb1,idb2))
0227      &   /(-2*dotup(idc1,idk1))/(2*dotup(idk1,idk2)
0228      &     -2*dotup(idc1,idk1)-2*dotup(idc1,idk2))
0229
0230 c...propup(8,1)->mcc7, propup(9,1)->mcc8

0231          if (ii.eq.1 .or. ii.eq.3) then
0232            propup(8,ii)=1/(2*dotup(idb1,idb1)+2*dotup(idb1,idb2))
0233      &   /(2*dotup(idk1,idk2))/(2*dotup(idk1,idk2)
0234      &     -2*dotup(idk1,idc2)-2*dotup(idc2,idk2))
0235            propup(9,ii)=1/(2*dotup(idb1,idb1)+2*dotup(idb1,idb2))
0236      &   /(2*dotup(idk1,idk2))/(2*dotup(idk1,idk2)
0237      &     -2*dotup(idc1,idk1)-2*dotup(idc1,idk2))
0238          end if
0239
0240 c...propup(10,1)->mco1, propup(11,1)->mco2

0241          propup(10,ii)=1/(2*dotup(idb1,idb1)+2*dotup(idb1,idb2))
0242      &   /(-2*dotup(idc2,idk1))/(2*dotup(idb1,idb1)+2*dotup(idb1,idb2)
0243      &     -2*dotup(idb1,idk2)-2*dotup(idb2,idk2))
0244          propup(11,ii)=1/(2*dotup(idb1,idb1)+2*dotup(idb1,idb2))
0245      &   /(-2*dotup(idc1,idk1))/(2*dotup(idb1,idb1)+2*dotup(idb1,idb2)
0246      &   -2*dotup(idb1,idk2)-2*dotup(idb2,idk2))
0247
0248 c...propup(12,1)->moo1

0249          if (ii.eq.1 .or. ii.eq.2) then
0250            propup(12,ii)=1/(2*dotup(idb1,idb1)+2*dotup(idb1,idb2))/
0251      &   (2*dotup(idc1,idc1)+2*dotup(idc2,idc1))/(2*dotup(idb1,idb1)
0252      &     +2*dotup(idb1,idb2)-2*dotup(idb1,idk1)-2*dotup(idb2,idk1))
0253          end if
0254
0255 c...propup(13,1)->moo3, propup(14,1)->moo4

0256          if (ii.eq.1) then
0257            propup(13,ii)=1/(2*dotup(idb1,idb1)+2*dotup(idb1,idb2))/
0258      &   (2*dotup(idk1,idk2))/(2*dotup(idc1,idc1)+2*dotup(idc2,idc1))
0259            propup(14,ii)=1/(2*dotup(idb1,idb1)+2*dotup(idb1,idb2))/
0260      &   (2*dotup(idc1,idc1)+2*dotup(idc2,idc1))
0261          end if
0262          end if
0263
0264        if(iqqbar.eq.1) then
0265            if (ii.eq.1 .or. ii.eq.3) then
0266              propup(1,ii)=1/(2*dotup(idk1,idk2))/(2*dotup(idk1,idk2)-
0267      &           2*dotup(idk1,idb1)-2*dotup(idk2,idb1))/
0268      &         (2*dotup(idc1,idc1)+2*dotup(idc1,idc2))
0269                  propup(2,ii)=1/(2.*dotup(idk1,idk2))/(2*dotup(idk1,idk2)-
0270      &           2*dotup(idk1,idb2)-2*dotup(idk2,idb2))/
0271      &         (2*dotup(idc1,idc1)+2*dotup(idc1,idc2))
0272                  propup(4,ii)=1/(-2*dotup(idk1,idb1)-2*dotup(idk1,idb2)+
0273      &           2*dotup(idb1,idb1)+2*dotup(idb1,idb2))/
0274      &         (2*dotup(idc1,idc1)+2*dotup(idc1,idc2))/
0275      &         (2*dotup(idb1,idb1)+2*dotup(idb1,idb2))
0276              if(ii.eq.1) then
0277                    propup(3,ii)=1/(2*dotup(idk1,idk2))/
0278      &         (2*dotup(idc1,idc1)+2*dotup(idc1,idc2))/
0279      &         (2*dotup(idb1,idb1)+2*dotup(idb1,idb2))
0280              end if
0281            end if
0282          end if
0283 700   continue
0284
0285 c...five collectted matrix elements according to the color factor c_1ij,

0286 c...c_2ij, c_3ij, c_4ij and c_5ij. note here the dsqrt part in the denominator

0287 c...is needed, for all the quarks are massive. the following result is get

0288 c...from the mathmatical program: matrixtotal.nb. here kk=1-64 stands for the

0289 c...helicity of the amplitude.

0290
0291       if(iqqbar.eq.0) then
0292 c0--------

0293       if(ioctet.eq.0) then
0294          do 900, kk=1,64
0295          colmat(1,kk)=(-((bfun(6,1,kk) + (bfun(9,1,kk) - bfun(9,2,kk))*
0296      -   (dotup(1,2) - dotup(5,1) - dotup(5,2)))*propup(8,1)) -
0297      -   bfun(8,3,kk)*propup(9,3) - (4*bfun(5,3,kk) - 4*bfun(5,4,kk) +
0298      -   bfun(9,3,kk)-bfun(9,4,kk))*(dotup(1,2)-dotup(4,1)-dotup(4,2))*
0299      -   propup(9,3) + bfun(8,4,kk)*(propup(7,4) + propup(9,3)) +
0300      -   2*bfun(5,1,kk)*dotup(5,1)*propup(10,1) +(-bfun(2,4,kk)
0301      -   + bfun(4,4,kk))*propup(11,4) + 2*bfun(5,4,kk)*dotup(4,2)
0302      -   *propup(11,4) +(-bfun(2,1,kk) + bfun(4,1,kk))*propup(12,2) +
0303      -   2*bfun(5,1,kk)*dotup(4,4)*propup(12,2)+2*bfun(5,1,kk)*
0304      -   dotup(4,6) *propup(12,2) + 2*bfun(5,1,kk)*dotup(5,1)*
0305      -   propup(12,2) +2*bfun(5,1,kk)*dotup(7,1)*propup(12,2) -
0306      -   bfun(3,1,kk)*(propup(10,1) + propup(12,2)) +bfun(1,1,kk)*
0307      -   (propup(1,1) + propup(10,1) + propup(12,2)) - (bfun(6,1,kk)
0308      -   - bfun(8,1,kk) + bfun(8,2,kk))*propup(13,1) +4*bfun(5,1,kk)
0309      -   *dotup(1,2)*propup(13,1)-4*bfun(5,2,kk)*dotup(1,2)*propup(13,1)
0310      -   +bfun(9,1,kk)*dotup(5,1)*propup(13,1) -bfun(9,2,kk)*dotup(5,1)
0311      -   *propup(13,1) + bfun(9,1,kk)*dotup(5,2)*propup(13,1) -
0312      -   bfun(9,2,kk)*dotup(5,2)*propup(13,1) - 4*bfun(5,1,kk)*
0313      -   dotup(7,1)*propup(13,1)+4*bfun(5,2,kk)*dotup(7,1)*propup(13,1)-
0314      -   bfun(9,1,kk)*dotup(7,1)*propup(13,1) + bfun(9,2,kk)*
0315      -   dotup(7,1)*propup(13,1) - 4*bfun(5,1,kk)*dotup(7,2)*
0316      -   propup(13,1) + 4*bfun(5,2,kk)*dotup(7,2)*propup(13,1) -
0317      -   bfun(9,1,kk)*dotup(7,2)*propup(13,1) + bfun(9,2,kk)
0318      -   *dotup(7,2)*propup(13,1) + bfun(6,2,kk)*(propup(5,2)
0319      -   + propup(8,1) + propup(13,1)) - ((2*bfun(5,1,kk) -
0320      -   4*bfun(5,2,kk)+bfun(9,1,kk)+bfun(9,2,kk))*propup(14,1))/2.)
0321      -   /dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*dotup(8,5))
0322
0323          colmat(2,kk)=((bfun(9,1,kk) - bfun(9,2,kk))*
0324      -   (dotup(1,2) - dotup(5,1) - dotup(5,2))*propup(8,1) -
0325      -   bfun(8,4,kk)*propup(9,3) + 4*bfun(5,3,kk)*dotup(1,2)*
0326      -   propup(9,3) -4*bfun(5,4,kk)*dotup(1,2)*propup(9,3) +
0327      -   bfun(9,3,kk)*dotup(1,2)*propup(9,3)-bfun(9,4,kk)*dotup(1,2)
0328      -   *propup(9,3) -4*bfun(5,3,kk)*dotup(4,1)*propup(9,3) +
0329      -   4*bfun(5,4,kk)*dotup(4,1)*propup(9,3) - bfun(9,3,kk)
0330      -   *dotup(4,1)*propup(9,3) + bfun(9,4,kk)*dotup(4,1)
0331      -   *propup(9,3) -4*bfun(5,3,kk)*dotup(4,2)*propup(9,3) +
0332      -   4*bfun(5,4,kk)*dotup(4,2)*propup(9,3)-bfun(9,3,kk)*dotup(4,2)
0333      -   *propup(9,3) +bfun(9,4,kk)*dotup(4,2)*propup(9,3) +
0334      -   bfun(8,3,kk)*(propup(7,3) + propup(9,3)) +2*bfun(5,2,kk)
0335      -   *dotup(5,2)*propup(10,2) +(-bfun(2,3,kk) + bfun(4,3,kk)
0336      -   )*propup(11,3) +2*bfun(5,3,kk)*dotup(4,1)*propup(11,3) +
0337      -   (-bfun(2,2,kk) + bfun(4,2,kk))*propup(12,1) +
0338      -   2*bfun(5,2,kk)*dotup(4,4)*propup(12,1) +2*bfun(5,2,kk)
0339      -   *dotup(4,6)*propup(12,1) + 2*bfun(5,2,kk)*dotup(5,2)
0340      -   *propup(12,1) +2*bfun(5,2,kk)*dotup(7,2)*propup(12,1) -
0341      -   bfun(3,2,kk)*(propup(10,2) + propup(12,1)) +
0342      -   bfun(1,2,kk)*(propup(1,2) + propup(10,2) + propup(12,1)) +
0343      -   (-bfun(8,1,kk) + bfun(8,2,kk))*propup(13,1) -
0344      -   4*bfun(5,1,kk)*dotup(1,2)*propup(13,1) +4*bfun(5,2,kk)*
0345      -   dotup(1,2)*propup(13,1)-bfun(9,1,kk)*dotup(5,1)*propup(13,1) +
0346      -   bfun(9,2,kk)*dotup(5,1)*propup(13,1) - bfun(9,1,kk)
0347      -   *dotup(5,2)*propup(13,1)+bfun(9,2,kk)*dotup(5,2)*propup(13,1)+
0348      -   4*bfun(5,1,kk)*dotup(7,1)*propup(13,1) - 4*bfun(5,2,kk)
0349      -   *dotup(7,1)*propup(13,1) + bfun(9,1,kk)*dotup(7,1)
0350      -   *propup(13,1) - bfun(9,2,kk)*dotup(7,1)*propup(13,1) +
0351      -   4*bfun(5,1,kk)*dotup(7,2)*propup(13,1) - 4*bfun(5,2,kk)
0352      -   *dotup(7,2)*propup(13,1) + bfun(9,1,kk)*dotup(7,2)
0353      -   *propup(13,1) - bfun(9,2,kk)*dotup(7,2)*propup(13,1) -
0354      -   bfun(6,2,kk)*(propup(8,1) + propup(13,1)) +
0355      -   bfun(6,1,kk)*(propup(5,1) + propup(8,1) + propup(13,1)) -
0356      -   ((-4*bfun(5,1,kk) + 2*bfun(5,2,kk) + bfun(9,1,kk) +
0357      -   bfun(9,2,kk))*propup(14,1))/2.)
0358      -   /dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*dotup(8,5))
0359
0360          colmat(3,kk)=(bfun(4,1,kk)*propup(4,1)+bfun(7,2,kk)*propup(6,2)+
0361      -   bfun(7,4,kk)*propup(6,4) - bfun(6,3,kk)*propup(8,3) -
0362      -   (bfun(9,3,kk)-bfun(9,4,kk))*(dotup(1,2)-dotup(6,1)-dotup(6,2))
0363      -   *propup(8,3) +bfun(6,4,kk)*(propup(5,4) + propup(8,3)) -
0364      -   (4*bfun(5,1,kk) - 4*bfun(5,2,kk) + bfun(9,1,kk) -
0365      -   bfun(9,2,kk))*(dotup(1,2) - dotup(7,1) - dotup(7,2))*
0366      -   propup(9,1) - 2*bfun(5,1,kk)*dotup(5,1)*propup(10,1) +
0367      -   (bfun(2,4,kk) - bfun(4,4,kk) - 2*bfun(5,4,kk)*dotup(4,2))*
0368      -   propup(11,4) - bfun(4,1,kk)*propup(12,2) -
0369      -   2*bfun(5,1,kk)*dotup(4,4)*propup(12,2) -
0370      -   2*bfun(5,1,kk)*dotup(4,6)*propup(12,2) -
0371      -   2*bfun(5,1,kk)*dotup(5,1)*propup(12,2) -
0372      -   2*bfun(5,1,kk)*dotup(7,1)*propup(12,2) +
0373      -   bfun(2,1,kk)*(propup(2,1) + propup(12,2)) -
0374      -   bfun(1,1,kk)*(propup(10,1) + propup(12,2)) +
0375      -   bfun(3,1,kk)*(propup(3,1)+propup(10,1)+propup(12,2))+
0376      -   (bfun(6,1,kk) - bfun(6,2,kk))*propup(13,1) -
0377      -   4*bfun(5,1,kk)*dotup(1,2)*propup(13,1) +
0378      -   4*bfun(5,2,kk)*dotup(1,2)*propup(13,1) -
0379      -   bfun(9,1,kk)*dotup(5,1)*propup(13,1) +
0380      -   bfun(9,2,kk)*dotup(5,1)*propup(13,1) -
0381      -   bfun(9,1,kk)*dotup(5,2)*propup(13,1) +
0382      -   bfun(9,2,kk)*dotup(5,2)*propup(13,1) +
0383      -   4*bfun(5,1,kk)*dotup(7,1)*propup(13,1) -
0384      -   4*bfun(5,2,kk)*dotup(7,1)*propup(13,1) +
0385      -   bfun(9,1,kk)*dotup(7,1)*propup(13,1) -
0386      -   bfun(9,2,kk)*dotup(7,1)*propup(13,1) +
0387      -   4*bfun(5,1,kk)*dotup(7,2)*propup(13,1) -
0388      -   4*bfun(5,2,kk)*dotup(7,2)*propup(13,1) +
0389      -   bfun(9,1,kk)*dotup(7,2)*propup(13,1) -
0390      -   bfun(9,2,kk)*dotup(7,2)*propup(13,1) -
0391      -   bfun(8,1,kk)*(propup(9,1) + propup(13,1)) +
0392      -   bfun(8,2,kk)*(propup(7,2) + propup(9,1) + propup(13,1)) +
0393      -   ((2*bfun(5,1,kk) - 4*bfun(5,2,kk) + bfun(9,1,kk) +
0394      -   bfun(9,2,kk))*propup(14,1))/2.)
0395      -   /dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*dotup(8,5))
0396
0397          colmat(4,kk)=(bfun(4,2,kk)*propup(4,2)+bfun(7,1,kk)*propup(6,1)+
0398      -   bfun(7,3,kk)*propup(6,3) - (bfun(6,4,kk) - (bfun(9,3,kk)
0399      -   -bfun(9,4,kk))*(dotup(1,2)-dotup(6,1)-dotup(6,2)))*propup(8,3)+
0400      -   bfun(6,3,kk)*(propup(5,3) + propup(8,3)) - (bfun(8,2,kk)
0401      -   -(4*bfun(5,1,kk)-4*bfun(5,2,kk)+bfun(9,1,kk) -bfun(9,2,kk))*
0402      -   (dotup(1,2) - dotup(7,1) - dotup(7,2)))*propup(9,1) -
0403      -   2*bfun(5,2,kk)*dotup(5,2)*propup(10,2) +bfun(3,2,kk)
0404      -   *(propup(3,2) + propup(10,2)) + (bfun(2,3,kk) -
0405      -   bfun(4,3,kk) - 2*bfun(5,3,kk)*dotup(4,1))*propup(11,3) +
0406      -   bfun(3,2,kk)*propup(12,1) -bfun(4,2,kk)*propup(12,1) -
0407      -   2*bfun(5,2,kk)*dotup(4,4)*propup(12,1) -
0408      -   2*bfun(5,2,kk)*dotup(4,6)*propup(12,1) -
0409      -   2*bfun(5,2,kk)*dotup(5,2)*propup(12,1) -
0410      -   2*bfun(5,2,kk)*dotup(7,2)*propup(12,1) +
0411      -   bfun(2,2,kk)*(propup(2,2) + propup(12,1)) -
0412      -   bfun(1,2,kk)*(propup(10,2) + propup(12,1)) +
0413      -   (-bfun(6,1,kk) + bfun(6,2,kk))*propup(13,1) -
0414      -   4*bfun(5,2,kk)*dotup(1,2)*propup(13,1) -
0415      -   (bfun(8,2,kk)-4*bfun(5,1,kk)*dotup(1,2))*propup(13,1)+
0416      -   bfun(9,1,kk)*dotup(5,1)*propup(13,1) -
0417      -   bfun(9,2,kk)*dotup(5,1)*propup(13,1) +
0418      -   bfun(9,1,kk)*dotup(5,2)*propup(13,1) -
0419      -   bfun(9,2,kk)*dotup(5,2)*propup(13,1) -
0420      -   4*bfun(5,1,kk)*dotup(7,1)*propup(13,1) +
0421      -   4*bfun(5,2,kk)*dotup(7,1)*propup(13,1) -
0422      -   bfun(9,1,kk)*dotup(7,1)*propup(13,1) +
0423      -   bfun(9,2,kk)*dotup(7,1)*propup(13,1) -
0424      -   4*bfun(5,1,kk)*dotup(7,2)*propup(13,1) +
0425      -   4*bfun(5,2,kk)*dotup(7,2)*propup(13,1) -
0426      -   bfun(9,1,kk)*dotup(7,2)*propup(13,1) +
0427      -   bfun(9,2,kk)*dotup(7,2)*propup(13,1) +
0428      -   bfun(8,1,kk)*(propup(7,1) + propup(9,1) + propup(13,1)) +
0429      -   ((-4*bfun(5,1,kk) + 2*bfun(5,2,kk) + bfun(9,1,kk) +
0430      -   bfun(9,2,kk))*propup(14,1))/2.)
0431      &   /dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*dotup(8,5))
0432
0433          colmat(5,kk)=(-(bfun(6,3,kk)*propup(5,3))-bfun(6,4,kk)*
0434      -   propup(5,4)-bfun(8,1,kk)*propup(7,1)-bfun(8,2,kk)*propup(7,2)+
0435      -   (-bfun(1,3,kk) + bfun(3,3,kk) - 2*bfun(5,3,kk)*dotup(6,1))*
0436      -   propup(10,3) + (-bfun(1,4,kk) + bfun(3,4,kk) -
0437      -   2*bfun(5,4,kk)*dotup(6,2))*propup(10,4) + (bfun(2,2,kk) -
0438      -   2*bfun(5,2,kk)*dotup(7,2))*propup(11,2) -bfun(1,2,kk)*
0439      -   propup(12,1) +(bfun(2,2,kk) + bfun(3,2,kk) - 2*bfun(5,2,kk)
0440      -   *(dotup(4,4) + dotup(4,6) + dotup(5,2) +dotup(7,2)))
0441      -   *propup(12,1) - bfun(4,2,kk)*(propup(4,2) +propup(11,2) +
0442      -   propup(12,1)) +(-bfun(1,1,kk) + bfun(3,1,kk))*propup(12,2)+
0443      -   bfun(2,1,kk)*(propup(11,1) + propup(12,2)) -bfun(4,1,kk)*
0444      -   (propup(4,1) + propup(11,1) + propup(12,2)) +(-bfun(5,2,kk)
0445      -   + bfun(9,1,kk) + bfun(9,2,kk))*propup(14,1) - bfun(5,1,kk)*
0446      -   (2*dotup(7,1)*propup(11,1) +2*(dotup(4,4) + dotup(4,6) +
0447      -   dotup(5,1) + dotup(7,1))*propup(12,2) + propup(14,1)))
0448      &   /dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*dotup(8,5))
0449 900    continue
0450        end if
0451 c0---------

0452 c8---------

0453          if(ioctet.eq.1) then
0454           do kk=1,64
0455             colmat(1,kk)=(-(bfun(7,1,kk)*propup(6,1)) +
0456      -      (bfun(6,2,kk) - (bfun(9,1,kk) - bfun(9,2,kk))*
0457      -      (dotup(1,2) - dotup(5,1) - dotup(5,2)))*propup(8,1) -
0458      -      bfun(6,1,kk)*(propup(5,1)+propup(8,1))+(bfun(8,2,kk)-
0459      -      (4*bfun(5,1,kk)-4*bfun(5,2,kk)+bfun(9,1,kk)-bfun(9,2,kk))*
0460      -      (dotup(1,2) - dotup(7,1) - dotup(7,2)))*propup(9,1) -
0461      -      bfun(8,1,kk)*(propup(7,1) + propup(9,1)))
0462      &      /dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*dotup(8,5))
0463                 colmat(2,kk)=(-(bfun(7,2,kk)*propup(6,2)) +
0464      -      (bfun(6,1,kk) + (bfun(9,1,kk) - bfun(9,2,kk))*
0465      -      (dotup(1,2) - dotup(5,1) - dotup(5,2)))*propup(8,1) -
0466      -      bfun(6,2,kk)*(propup(5,2)+propup(8,1))+(bfun(8,1,kk)+
0467      -      (4*bfun(5,1,kk)-4*bfun(5,2,kk)+bfun(9,1,kk)-bfun(9,2,kk))*
0468      -      (dotup(1,2) - dotup(7,1) - dotup(7,2)))*propup(9,1)-
0469      -      bfun(8,2,kk)*(propup(7,2) + propup(9,1)))
0470      &      /dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*dotup(8,5))
0471         colmat(3,kk)=(bfun(7,1,kk)*propup(6,1)+bfun(7,4,kk)*propup(6,4)-
0472      -      (bfun(1,2,kk) + 2*bfun(5,2,kk)*dotup(5,2))*propup(10,2) +
0473      -      (-bfun(1,3,kk) + bfun(3,3,kk) - 2*bfun(5,3,kk)*dotup(6,1))*
0474      -      propup(10,3) + (bfun(2,4,kk) - bfun(4,4,kk) -
0475      -      2*bfun(5,4,kk)*dotup(4,2))*propup(11,4) -
0476      -      (bfun(1,2,kk) - bfun(2,2,kk) + bfun(4,2,kk) +
0477      -      2*bfun(5,2,kk)*(dotup(4,4) + dotup(4,6) + dotup(5,2) +
0478      -      dotup(7,2)))*propup(12,1) +
0479      -      bfun(3,2,kk)*(propup(3,2) + propup(10,2) + propup(12,1)) -
0480      -      bfun(1,1,kk)*propup(12,2) + bfun(3,1,kk)*propup(12,2) -
0481      -      bfun(4,1,kk)*(propup(11,1) + propup(12,2)) +
0482      -      bfun(2,1,kk)*(propup(2,1) + propup(11,1) + propup(12,2)) +
0483      -      2*bfun(5,1,kk)*(-(dotup(7,1)*propup(11,1)) -
0484      -      (dotup(4,4) + dotup(4,6) + dotup(5,1) + dotup(7,1))*
0485      -      propup(12,2)) - (bfun(5,1,kk) + bfun(5,2,kk) -
0486      -      bfun(9,1,kk) - bfun(9,2,kk))*propup(14,1))
0487      &      /dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*dotup(8,5))
0488          colmat(4,kk)=(bfun(7,2,kk)*propup(6,2)+bfun(7,3,kk)*propup(6,3)
0489      -      -(bfun(1,1,kk) + 2*bfun(5,1,kk)*dotup(5,1))*propup(10,1) +
0490      -      (-bfun(1,4,kk) + bfun(3,4,kk) - 2*bfun(5,4,kk)*dotup(6,2))*
0491      -      propup(10,4) + (bfun(2,3,kk) - bfun(4,3,kk) -
0492      -      2*bfun(5,3,kk)*dotup(4,1))*propup(11,3) -
0493      -      bfun(1,2,kk)*propup(12,1) + bfun(3,2,kk)*propup(12,1) -
0494      -      bfun(4,2,kk)*(propup(11,2) + propup(12,1)) +
0495      -      bfun(2,2,kk)*(propup(2,2) + propup(11,2) + propup(12,1)) -
0496      -      2*bfun(5,2,kk)*(dotup(7,2)*propup(11,2) +
0497      -      (dotup(4,4) + dotup(4,6) + dotup(5,2) + dotup(7,2))*
0498      -      propup(12,1))- (bfun(1,1,kk) - bfun(2,1,kk) + bfun(4,1,kk))
0499      -      *propup(12,2)-2*bfun(5,1,kk)*(dotup(4,4)+dotup(4,6)+
0500      -      dotup(5,1)+dotup(7,1))*propup(12,2)+bfun(3,1,kk)*
0501      -      (propup(3,1)+propup(10,1)+propup(12,2)) - (bfun(5,1,kk) +
0502      -    bfun(5,2,kk) - bfun(9,1,kk) - bfun(9,2,kk))*propup(14,1))
0503      +    /dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*dotup(8,5))
0504           colmat(5,kk)=(-((bfun(6,4,kk) - (bfun(9,3,kk) - bfun(9,4,kk))*
0505      -      (dotup(1,2) - dotup(6,1) - dotup(6,2)))*propup(8,3)) +
0506      -      bfun(6,3,kk)*(propup(5,3) + propup(8,3)) - (bfun(8,2,kk) -
0507      -      (4*bfun(5,1,kk)-4*bfun(5,2,kk)+bfun(9,1,kk)-bfun(9,2,kk))*
0508      -      (dotup(1,2) - dotup(7,1) - dotup(7,2)))*propup(9,1) +
0509      -      (bfun(1,4,kk) - bfun(3,4,kk))*propup(10,4) +
0510      -      2*bfun(5,4,kk)*dotup(6,2)*propup(10,4) + (bfun(1,1,kk) -
0511      -      bfun(3,1,kk))*propup(12,2) - bfun(2,1,kk)*(propup(11,1) +
0512      -      propup(12,2)) + bfun(4,1,kk)*(propup(4,1) + propup(11,1) +
0513      -      propup(12,2)) +(-bfun(6,1,kk) + bfun(6,2,kk))*propup(13,1) -
0514      -      bfun(8,2,kk)*propup(13,1) -(4*bfun(5,2,kk)*(dotup(1,2) -
0515      -      dotup(7,1) - dotup(7,2)) - (bfun(9,1,kk) - bfun(9,2,kk))*
0516      -      (dotup(5,1) + dotup(5,2) - dotup(7,1) - dotup(7,2)))*
0517      -      propup(13,1) + bfun(8,1,kk)*(propup(7,1) + propup(9,1) +
0518      -      propup(13,1)) + bfun(5,1,kk)*(2*(dotup(4,4) + dotup(4,6) +
0519      -      dotup(5,1))*propup(12,2) + 2*dotup(7,1)*(propup(11,1) +
0520      -      propup(12,2) - 2*propup(13,1)) + 4*(dotup(1,2) -
0521      -      dotup(7,2))*propup(13,1) - propup(14,1)) +((4*bfun(5,2,kk)
0522      -      - bfun(9,1,kk) - bfun(9,2,kk))*propup(14,1))/2.)
0523      &      /dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*dotup(8,5))
0524                 colmat(6,kk)=(-((bfun(6,3,kk)+(bfun(9,3,kk)-bfun(9,4,kk))*
0525      -      (dotup(1,2) - dotup(6,1) - dotup(6,2)))*propup(8,3)) +
0526      -      bfun(6,4,kk)*(propup(5,4) + propup(8,3)) -
0527      -      (4*bfun(5,1,kk) - 4*bfun(5,2,kk) + bfun(9,1,kk) -
0528      -      bfun(9,2,kk))*(dotup(1,2) - dotup(7,1) - dotup(7,2))*
0529      -      propup(9,1) + (bfun(1,3,kk) - bfun(3,3,kk))*propup(10,3)+
0530      -      2*bfun(5,3,kk)*dotup(6,1)*propup(10,3) +
0531      -      (bfun(1,2,kk) - bfun(3,2,kk))*propup(12,1) -
0532      -      bfun(2,2,kk)*(propup(11,2) + propup(12,1)) +
0533      -      bfun(4,2,kk)*(propup(4,2) + propup(11,2) + propup(12,1))+
0534      -      (bfun(6,1,kk) - bfun(6,2,kk))*propup(13,1) -
0535      -      4*bfun(5,1,kk)*dotup(1,2)*propup(13,1) -
0536      -      ((bfun(9,1,kk) - bfun(9,2,kk))*(dotup(5,1) +dotup(5,2) -
0537      -      dotup(7,1) - dotup(7,2)) - 4*bfun(5,1,kk)*(dotup(7,1)
0538      -      + dotup(7,2)))*propup(13,1) - bfun(8,1,kk)*(propup(9,1) +
0539      -      propup(13,1)) + bfun(8,2,kk)*(propup(7,2) + propup(9,1) +
0540      -      propup(13,1)) + bfun(5,2,kk)*(2*(dotup(4,4) + dotup(4,6)+
0541      -      dotup(5,2))*propup(12,1) + 2*dotup(7,2)*(propup(11,2) +
0542      -      propup(12,1) - 2*propup(13,1)) + 4*(dotup(1,2) -
0543      -      dotup(7,1))*propup(13,1)-propup(14,1))-((-4*bfun(5,1,kk)+
0544      -      bfun(9,1,kk) + bfun(9,2,kk))*propup(14,1))/2.)
0545      &      /dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*dotup(8,5))
0546                 colmat(7,kk)=(-(bfun(7,3,kk)*propup(6,3)) +
0547      -      (bfun(6,4,kk) - (bfun(9,3,kk) - bfun(9,4,kk))*
0548      -      (dotup(1,2) - dotup(6,1) - dotup(6,2)))*propup(8,3) -
0549      -      bfun(6,3,kk)*(propup(5,3) + propup(8,3)) +
0550      -      (bfun(8,4,kk) - (4*bfun(5,3,kk) - 4*bfun(5,4,kk) +
0551      -      bfun(9,3,kk) - bfun(9,4,kk))*
0552      -      (dotup(1,2) - dotup(4,1) - dotup(4,2)))*propup(9,3) -
0553      -      bfun(8,3,kk)*(propup(7,3) + propup(9,3)))
0554      &      /dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*dotup(8,5))
0555                 colmat(8,kk)=(-(bfun(7,4,kk)*propup(6,4)) +
0556      -      (bfun(6,3,kk) + (bfun(9,3,kk) - bfun(9,4,kk))*
0557      -      (dotup(1,2) - dotup(6,1) - dotup(6,2)))*propup(8,3) -
0558      -      bfun(6,4,kk)*(propup(5,4) + propup(8,3)) +
0559      -      (bfun(8,3,kk) + (4*bfun(5,3,kk) - 4*bfun(5,4,kk) +
0560      -      bfun(9,3,kk) - bfun(9,4,kk))*
0561      -      (dotup(1,2) - dotup(4,1) - dotup(4,2)))*propup(9,3) -
0562      -      bfun(8,4,kk)*(propup(7,4) + propup(9,3)))
0563      &      /dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*dotup(8,5))
0564                 colmat(9,kk)=(-(bfun(1,1,kk)*propup(1,1))-bfun(2,1,kk)*
0565      -      propup(2,1)-bfun(3,1,kk)*propup(3,1)-bfun(4,1,kk)*
0566      &      propup(4,1))/dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*
0567      &      dotup(8,5))
0568                 colmat(10,kk)=(-(bfun(1,2,kk)*propup(1,2))-bfun(2,2,kk)*
0569      -    propup(2,2) - bfun(3,2,kk)*propup(3,2) - bfun(4,2,kk)*
0570      -      propup(4,2))/dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*
0571      &      dotup(8,5))
0572           end do
0573          end if
0574 c8---------

0575       end if
0576
0577         if(iqqbar.eq.1) then
0578          do kk=1,64
0579         colmat(1,kk)=(propup(2,1)*bfun(2,1,kk)+propup(1,3)*bfun(1,3,kk)
0580      &       -propup(3,1)*bfun(3,1,kk)-propup(4,3)*bfun(4,3,kk))
0581      &       /dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*dotup(8,5))
0582         colmat(2,kk)=(propup(1,1)*bfun(1,1,kk)+propup(2,3)*bfun(2,3,kk)
0583      &       +propup(3,1)*bfun(3,1,kk)+propup(4,3)*bfun(4,3,kk))
0584      &       /dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*dotup(8,5))
0585         colmat(3,kk)=(propup(4,1)*bfun(4,1,kk)+propup(4,3)*bfun(4,3,kk))
0586      &     /dsqrt(dotup(8,6)*dotup(8,7)*dotup(8,4)*dotup(8,5))
0587          end do
0588         end if
0589
0590         end