* *------------------------------------------------------------------------------ * subroutine wto(*rs*,*wmi*,*zmi*,*zgi*,*tqmi*,*hmi*,*alszi*,sig*,esig*) * * computes the cross-section for given sqrt(s),M_W,M_Z,G_Z,M_t,M_H,alpha_s *------------------------------------------------------------------------------ * subroutine wto(rs,wmi,zmi,zgi,tqmi,hmi,alszi,sig,esig) implicit real*8 (a-h,o-z) character*1,oxcm,ocoul,oud,oqcd,oglu,opeak,ofl,ofsr,oral,obin, # oprt,opeakn,om,osm,oww,ozz,opeaka,ostop,rio, # oanom,otrans,omh,ord,otop,ostore,omssm,ockm, # omdist,opglu,osmass,oint,oseed,osinw,oprint character*2,ofs,oev character*4,otype * parameter (nin=98,nout=99,ndimmx=9,ninv=10,liw=ndimmx+2, # lw=ndimmx*(ndimmx-1)/2+12*ndimmx) parameter (itmxmx=100000,im=100,npos=512,ifmax=10000) parameter(nbf=2,nt=1,lwa=(nbf*(3*nbf+13))/2,lwat=(nt*(3*nt+13))/2) * common/wtmp/zrm common/wtai/rio common/wtfb/oev common/wtcw/oww common/wtcz/ozz common/wtmod/om common/wtfs/ofs common/wtud/oud common/wtoi/oint common/wtii/iint common/wtsfh/ip0 common/wtkr/ireg common/wtcb/obin common/wtfls/ofl common/wtqcd/als common/wtgg/oglu common/wtsw/osinw common/wtrms/rsm2 common/wtim/ostop common/wtprt/oprt common/wtfsr/ofsr common/wtopa/delc common/wtps/opeak common/wthiggs/hm common/wtdis/dist common/wtsmod/osm common/wtafsz/g2z common/wtkount/ik common/wtckm/ockm common/wtio/oprint common/wtickm/ickm common/wthqcd/alsh common/wtistrf/isf common/wtbme/bfact common/wthiggo/omh common/wtoral/oral common/wtlmsb/qcdl common/wtqcdz/alsz common/wtaqcd/oqcd common/wtap/opeaka common/wtcqcd/iqcd common/wtcac/oanom common/wtgnum/itmx common/wtsp/psg(4) common/wtfsm/fsm(4) common/wtpqcd/opglu common/wopst/ostore common/wtsf/ix0,it0 common/wtpsn/opeakn common/wtomd/omdist common/wtseed/oseed common/wtcoul/ocoul common/wtmult/inm,jnm common/wttopop/otop common/wtchi/hch(36) common/wtspl/prx0(9) common/wopstm/osmass common/wtshel/otrans common/wtipt/ifz(51) common/wtfsw/tfactsw common/wtmssmo/omssm common/wtcurrent/oxcm common/wtnf/ifl(npos) common/wticuts/iac(4) common/wtrm/rbm2,rcm2 common/wtnpr/ipr,ipr0 common/wtisa/isaa,isab common/wttopt/ios,iosf common/wtspln/prx0n(9) common/wtmd/arrinv(10) common/wtvckm/vckm(3,3) common/wtmatx/colf(8,8) common/wtparc/xap(ninv) common/wtcx/xscmx(npos) common/wthx/xshmx(npos) common/wtochannel/otype common/wtpoints/ipm,irm common/wtcutsw/tcs,s0cut common/wtparh/xaph(ninv) common/wtpmx/xmx(npos,9) common/wtdis2/distm,distp common/wtbac/dg1z,dkg,rlg common/wtcx0/xscmx0(npos) common/wthx0/xshmx0(npos) common/wtrmss/chcm2,chsm2 common/wttc/itc,itcc,itcn common/wtpmxh/xmxh(npos,9) common/wtlb/abp,bbp,abm,bbm common/wtpls/wmu,wlg,zmu,zlg common/wtstor/stry(npos,ifmax) common/wtgd/idtl,idt4,idt7,idt9 common/wtncc/chf2,chfp2,conc(10) common/wtcclr/vupl,vupr,vdpl,vdpr common/wtiaz/alzr,alzi,alzir,alzii common/wtalqed/arezm,aizm,aresz,aisz common/wtcmplx/sw(2),sz(2),rsw(2),rsz(2) common/wtmssmi/am,tbeta,rmu,scalm,bat,bab common/wtnclr/vel,ver,velr,vfl,vfr,vfpl,vfpr common/wtcchannel/chu,chup,chd,chdp,fcuc,fcdc common/wtnchannel/chf,chfp,tif,tifp,fcun,fcdn common/wthapar/rhm,rhm2,rhg,rhmg,shg,shgs,opshgs common/wtsasw/separa,asccsw,sthqsw,eusw,edsw,alsw common/wtee/qch,qch2,vqr,vql,hbe(24),hbo(24),hmp(24) common/wtfmass/em,rmm,tm,rnm,uqm,dqm,cqm,sqm,bqm,tqm,dmy common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi common/wtacc/acg1g,aclg,ackg,acg4g,acktg,acltg,acg5g,acg1z,aclz, # ackz,acg4z,acktz,acltz,acg5z common/wtfmass2/em2,amm2,tm2,rnm2,uqm2,dqm2,cqm2,sqm2,bqm2, # tqm2,dmy2 common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,g8,tfact,pih,alw, # eta,feta,beta,g2,tfacth common/wtacchannel/omchu,opchu,omchup,opchup,omchdp,opchdp, # omchd,opchd,hchup,hchu,hchdp,hchd common/wtapar/ars,s,rwm,rwm2,rwg,rwmg,swg,swgs,opswgs,sth2,cth2, # hsth2,tsth2,scth2,asth2,tth2,rzm,rzm2,rzg,rzmg,szg, # szgs,opszgs,sth4,cth4,ve,vf,vfp,rbqm2,rszw,rszw2, # s0w,s0z common/wtsubreg/dsm,usm,dsp,usp,rl(6),rr(6),srl(6),sdsm,sdsp,vvl1, # vvl2,vvl3,ul,omul,suml common/wtcuts/aim(6),bim(6),ae(4),asa(4),bsa(4),afsa(6),bfsa(6), # ombsa(4),opbsa(4),teq,rae(4),omasa(4),opasa(4), # sg12,cg12,sg13,cg13,sg14,cg14,sg23,cg23,sg24, # cg24,sg34,cg34,sct120,sct130,sct140,sct230, # sct240,sct340,sgam(4),cgam(4) common/wtmssm/ams,shm,shms,bhm,bhms,sbeta,cbeta,salpha,calpha, # sbma,cbma,rshm,rshm2,rshg,rshmg,sshg,sshgs, # opsshgs,rbhm,rbhm2,rbhg,rbhmg,sbhg,sbhgs,opsbhgs, # ram,ram2,rag,ramg,sag,sags,opsags common/wtmssmc/chms,chm,rchn,rchm2,rchg,rchmg,schg,schgs,opschgs * dimension xmx0(9) dimension sclm(4) dimension xmxh0(9) dimension idhep(6) dimension phep(5,6) dimension aphep(5,4) dimension vk(ndimmx) dimension prx(ndimmx) dimension zs(2),zl(2) dimension aprx(ndimmx) dimension prxa(ninv,itmxmx) dimension axstr(100),ir(100) dimension vol(npos,4),volt(4) dimension exc(im),extst(itmxmx) dimension excl(im),qtrf(itmxmx) dimension volh(npos,8),volht(8) dimension sigfl(4,6),esigfl(4,6) dimension excwp(im),cwpst(itmxmx) dimension volhq(npos,12),volhtq(12) dimension volhb(npos,20),volhtb(20) dimension st(nt),fvect(nt),wat(lwat) dimension sigpart(3,6),esigpart(3,6) dimension gbl(npos,9,4),gbu(npos,9,4) dimension gblh(npos,9,8),gbuh(npos,9,8) dimension gblhq(npos,9,12),gbuhq(npos,9,12) dimension gblhb(npos,9,20),gbuhb(npos,9,20) dimension exce(im),exetst(itmxmx),exc4j(itmxmx) dimension wmx(lw),bl(ndimmx),bu(ndimmx),iw(liw) dimension fvecw(nbf),waw(lwa),waz(lwa),fvecz(nbf) dimension pggfz(2),pgglqz(2),p3qz(2),fzz(2),fwz(2) dimension amps(itmxmx),amps1(itmxmx),amps2(itmxmx) dimension siga(4),sigan(4),sigah(8),sigahq(12),sigahb(20) dimension xsmx(npos,4),xsmxn(npos,4),xsmxh(npos,8), # xsmxhq(npos,12),xsmxhb(npos,20) dimension oxsmx(npos,4),oxsmxn(npos,4),oxsmxh(npos,8), # oxsmxhq(npos,12),oxsmxhb(npos,20) dimension exc1(im),extst1(itmxmx),exc2(im),extst2(itmxmx), # exc3(im),exc4(im),exc5(im),exc6(im),exc7(im), # extst3(itmxmx),extst4(itmxmx),extst5(itmxmx), # extst6(itmxmx),excm(im) * external d01gcf,wtoregion,wtoxsc,wtoxsn,wtoxsm43,wtoxsnh19, # wtoxsnh24,wtoxsnh49,wtoxsnh32,wtoxsnh64,wtoxsn32, # wtoxsn64,wtoxsnh26,wtoxsnsh64,wtoxssc,wtoxssw,wtoxsswc, # wtoxsswcfl external s14aaf,s21baf,x04cbf,x04abf,s09abf external g05caf,g05cbf,g05ccf,e04jaf external c05nbf,c02ajf,x02ajf,wtofsw,wtofsz,wtofst * if(osm.eq.'g'.and.oseed.eq.'y') then call g05ccf endif * ik= 0 do i=1,51 ifz(i)= 0 enddo idtl= 0 idt4= 0 idt7= 0 idt9= 0 * wm= wmi zm= zmi zg= zgi alsz= alszi * ars= rs pis= pi*pi em2= em*em amm2= rmm*rmm tm2= tm*tm rnm2= rnm*rnm uqm2= uqm*uqm dqm2= dqm*dqm bqm2= bqm*bqm cqm2= cqm*cqm sqm2= sqm*sqm bqm2= bqm*bqm dmy2= dmy*dmy hm= hmi * *-----alpha_s * nf= 5 if(iqcd.eq.2) then x1= 0.001d0 x2= 10.0d0 xacc= 1.d-12 nf= 5 qcdl= wtoqcdlam(nf,als,wm,x1,x2,xacc) else qcdl= 0.d0 endif als= wtoralphas(zm,wm,alsz,nf) * do i=1,4 rae(i)= ae(i)/rs enddo teq= rae(1)+rae(2)+rae(3)+rae(4) sct120= 2.d0*rae(1)*rae(2)*(1.d0-bfsa(1)) sct130= 2.d0*rae(1)*rae(3)*(1.d0-bfsa(2)) sct140= 2.d0*rae(1)*rae(4)*(1.d0-bfsa(3)) sct230= 2.d0*rae(2)*rae(3)*(1.d0-bfsa(4)) sct240= 2.d0*rae(2)*rae(4)*(1.d0-bfsa(5)) sct340= 2.d0*rae(3)*rae(4)*(1.d0-bfsa(6)) * *-----parameters * alpha= 1.d0/alphai alw= 1.d0/alwi pih= 0.5d0*pi s= rs*rs qgfopi= 0.25d0*gf/pi wm2= wm*wm wmc= wm2*wm zm2= zm*zm wg0= 1.5d0*gf*wmc/pi if(oqcd.eq.'n') then wg= wg0 else wg= wg0*(1.d0+2.d0/3.d0*als/pi) endif wg2= wg*wg if(ios.eq.1) then sth2= 0.5d0*pi*alw/gf/wm2 g2= 4.d0*pi*alw/sth2 else if(ios.eq.2) then sth2= 1.d0-wm2/zm2 g2= 8.d0*gf*wm2 else if(ios.eq.3) then scal= -zm2 call wtopself(scal,pggf) derl= 0.25d0*alpha/pi*pggf call wtohadr5(zm,derh,ederh) der= derl+derh alz= (1.d0-der)/alpha args= 1.d0-2.d0*pi/alz/gf/zm2 sth2= 0.5d0*(1.d0-sqrt(args)) g2= 8.d0*gf*zm2*(1.d0-sth2) endif g4= g2*g2 g8= g4*g4 cth2= 1.d0-sth2 tth2= sth2/cth2 hsth2= 0.5d0*sth2 tsth2= 2.d0*sth2 scth2= 0.125d0/cth2 asth2= 0.5d0-sth2 sth4= sth2*sth2 cth4= cth2*cth2 psg(1)= 1.d0 psg(2)= -1.d0 psg(3)= -1.d0 psg(4)= 1.d0 if(ofl.eq.'a'.and.oanom.eq.'r') then acg1g= 0.d0 aclg= rlg ackg= dkg acg4g= 0.d0 acktg= 0.d0 acltg= 0.d0 acg5g= 0.d0 acg1z= dg1z aclz= rlg ackz= dg1z-sth2/cth2*dkg acg4z= 0.d0 acktz= 0.d0 acltz= 0.d0 acg5z= 0.d0 endif * *-----mass&width parameters * rwm= wm/rs rwm2= wm2/s rwg= wg/rs swg= wg/wm swgs= swg*swg if(oww.eq.'r'.or.oww.eq.'f') then rwmg= rwm*rwg else if(oww.eq.'i') then rwmg= rwm*(1.d0-rwg*rwg/rwm/rwm)*rwg else if(oww.eq.'t') then rwmg= (1.d0+swgs)*rwm*rwg endif opswgs= 1.d0+swgs s0w= rwm2/opswgs rzm= zm/rs rzm2= zm2/s rzg= zg/rs if(ozz.eq.'r'.or.ozz.eq.'f') then rzmg= rzm*rzg else if(ozz.eq.'i') then rzmg= rzm*(1.d0-rzg*rzg/rzm/rzm)*rzg endif szg= zg/zm szgs= szg*szg opszgs= 1.d0+szgs rszw= zg/zm rszw2= rszw*rszw s0z= rzm2/opszgs * tqm= tqmi if(ofl.eq.'c') then if(rio.eq.'i') then tqm2= tqm*tqm else if(rio.eq.'a') then if(otop.eq.'f') then tqm2= tqm*tqm else if(otop.eq.'d') then * *-----Iteration for m_t starts * gmt= 150.d0 st(1)= gmt*gmt tol= sqrt(x02ajf()) ifail= 1 call c05nbf(wtofst,nt,st,fvect,tol,wat,lwat,ifail) tqm2= st(1) tqm= sqrt(tqm2) endif endif * *-----alpha(zm) * pzr= -zm2 pzi= 0.d0 call wtopoleg(pzr,pzi,pggfz,pgglqz,pggnpz) call wtopolez0(pggf0,fw0) xal= alphai-0.25d0*(pggfz(1)-pggf0+pggnpz)/pi yal= -0.25d0*(pggfz(2)+pgglqz(2))/pi alzir= xal alzii= yal alzr= xal/(xal*xal+yal*yal) alzi= -yal/(xal*xal+yal*yal) * *-----Re(1/g^2(zm)) * apis= 16.d0*pis p2zr= -zm2 p2zi= 0.d0 call wtopole(p2zr,p2zi,p3qz,fzz,fwz) bx= 1.d0/8.d0/gf/zm2+(fw0-fzz(1))/zm2/apis xih= -p3qz(2)/apis agc= 4.d0*pi*alzr bgc= -1.d0-8.d0*pi*xih*alzi cgc= -4.d0*pi*xih*xih*alzr+bx ifail= 0 call c02ajf(agc,bgc,cgc,zs,zl,ifail) g2z= zs(1) * *-----Iteration for sw starts * sw(1)= wm2 sw(2)= -wm*wg tol= sqrt(x02ajf()) ifail= 1 call c05nbf(wtofsw,nbf,sw,fvecw,tol,waw,lwa,ifail) wmu= sqrt(sw(1)) wlg= -sw(2)/wmu do l=1,2 rsw(l)= sw(l)/s enddo * *-----Iteration for sz starts * dszi= 1.d-1 szi0= zm*zg izt= 1 1005 sz(1)= zm2 sz(2)= -szi0*(1.d0-(izt-1)*dszi) tol= sqrt(x02ajf()) ifail= 1 call c05nbf(wtofsz,nbf,sz,fvecz,tol,waz,lwa,ifail) if(ifail.eq.0) then zmu= sqrt(sz(1)) zlg= -sz(2)/zmu do l=1,2 rsz(l)= sz(l)/s enddo else if(izt.lt.20) then izt= izt+1 go to 1005 else print*,' Failure in searching for Z pole' endif endif endif tqm2= tqm*tqm * if((otype.eq.'nc21'.or.otype.eq.'nc25'.or. # otype.eq.'nc33'.or.otype.eq.'nc50'.or. # otype.eq.'nc68').and.omssm.eq.'n') then call wtocorrqcd(hm,als,bqm2,cqm2,rbm2,rcm2,rsm2) if(oqcd.eq.'y') then nf= 5 * alsh= wtoralphas(wm,hm,als,nf) * hg= qgfopi*hm*(3.d0*(rbm2+rcm2)*(1.d0+alsh/pi*( * # 17.d0/3.d0+(35.94d0-1.36d0*nf)*alsh/pi))+tm*tm) * hgb= qgfopi*hm*(3.d0*rbm2*(1.d0+alsh/pi*( * # 17.d0/3.d0+(35.94d0-1.36d0*nf)*alsh/pi))) * hgr= qgfopi*hm*(3.d0*rcm2*(1.d0+alsh/pi*( * # 17.d0/3.d0+(35.94d0-1.36d0*nf)*alsh/pi))+tm*tm) *///////////////////////////////////////////////////////////////////////// alsh= wtoralphas(wm,hm,als,nf) scalh= -hm*hm call wtopself(scalh,pggfh) derlh= 0.25d0*alpha/pi*pggfh call wtohadr5(hm,derhh,ederhh) derh= derlh+derhh alh= alpha/(1.d0-derh) hm2= hm*hm hctotc= alsh/pi*(17.d0/3.d0+(35.94d0-1.36d0*nf)* # alsh/pi) hcorrb= -6.d0*rbm2/hm2+0.472*alh/pi+alsh/pi*(0.651d0* # alh/pi+5.667d0-40.d0*rbm2/hm2+alsh/pi*( # 29.14671d0-87.725d0*rbm2/hm2+41.75760d0*alsh/pi)) clnt= 2*log(hm/tqm) hcorrt= alsh*alsh/pis*(3.111d0-0.667*clnt+rbm2/hm2*(-10.d0+ # 4.d0*clnt+4.d0/3.d0*log(rbm2/hm2))+hm2/tqm2*( # 0.241d0-0.070d0*clnt)+alsh/pi*( # 50.474d0-clnt*(8.167d0+1.278*clnt)))+gf*tqm2/ # 8.d0/pis*(1+alsh/pi*(-4.913d0+alsh/pi*(-72.117d0- # 20.945d0*clnt))) hctotb= hcorrb+hcorrt hg= qgfopi*hm*(3.d0*rcm2*(1.d0+hctotc)+3.d0*rbm2* # (1.d0+hctotb)+tm*tm) print*,qgfopi*hm*3.d0*rcm2*(1.d0+hctotc), # qgfopi*hm*3.d0*rbm2*(1.d0+hctotb), # qgfopi*hm*tm*tm hgb= qgfopi*hm*3.d0*rbm2*(1.d0+hctotb) hgr= qgfopi*hm*(3.d0*rcm2*(1.d0+hctotc)+tm*tm) *///////////////////////////////////////////////////////////////////////// else if(oqcd.eq.'n') then hg= qgfopi*hm*(3.d0*(rbm2+rcm2)+tm*tm) hgb= qgfopi*hm*3.d0*rbm2 hgr= qgfopi*hm*(3.d0*rcm2+tm*tm) endif if(oglu.eq.'y') then nf= 5 alsh= wtoralphas(wm,hm,als,nf) hg= hg+gf/36.d0/pi**3*alsh*alsh*hm**3*(1.d0+ # (95.d0/4.d0-7.d0/6.d0*nf)*alsh/pi) print*,gf/36.d0/pi**3*alsh*alsh*hm**3*(1.d0+ # (95.d0/4.d0-7.d0/6.d0*nf)*alsh/pi) hgr= hgr+gf/36.d0/pi**3*alsh*alsh*hm**3*(1.d0+ # (95.d0/4.d0-7.d0/6.d0*nf)*alsh/pi) endif * hm2= hm*hm hm3= hm2*hm if(hm.gt.(2.d0*wm)) then xwh= wm2/hm2 btw= sqrt(1.d0-4.d0*xwh) hgw= gf*hm3/8.d0/pi*(1.d0-4.d0*xwh*(1.d0-3.d0*xwh))* # btw hg= hg+hgw hgr= hgr+hgw endif if(hm.gt.(2.d0*zm)) then xzh= zm2/hm2 btz= sqrt(1.d0-4.d0*xzh) hgz= gf*hm3/16.d0/pi*(1.d0-4.d0*xzh*(1.d0-3.d0*xzh))* # btz hg= hg+hgz hgr= hgr+hgz endif betab= sqrt(1.d0-4.d0*rbm2/hm/hm) if(omh.eq.'y') then bfact= betab**3*(1.d0+hgb/hgr)/(1.d0+betab**3*hgb/hgr) else bfact= 1 endif else hg= 0.d0 bfact= 1.d0 endif * inpts= ipm inrand= irm * rhm= hm/rs rhm2= hm2/s rhg= hg/rs rhmg= rhm*rhg shg= hg/hm shgs= shg*shg opshgs= 1.d0+shgs * if(omssm.eq.'y') then ams= am*am rmus= rmu*rmu rmuc= rmus*rmu rmuq= rmus*rmus scalms= scalm*scalm scalmq= scalms*scalms bats= bat*bat babs= bab*bab tbetas= tbeta*tbeta sbeta= tbeta/sqrt(1.d0+tbeta*tbeta) cbeta= 1.d0/sqrt(1.d0+tbeta*tbeta) sbetas= sbeta*sbeta cbetas= cbeta*cbeta vev= 0.5d0/sqrt(gf) vevs= vev*vev gcg1s= 8.d0*gf*wm2 gcg2s= sth2/cth2*gcg1s gcg3s= gcg2s-gcg1s nf= 5 alst= wtoralphas(wm,tqm,als,nf) ttqm= tqm/(1.d0+4.d0/3.d0/pi*alst) call wtocorrqcd(tqm,als,bqm2,cqm2,tbm2,tcm2,tsm2) rht= ttqm/vev/sbeta rhb= sqrt(tbm2)/vev/cbeta rhts= rht*rht rhtq= rhts*rhts rhbs= rhb*rhb rhbq= rhbs*rhbs if(am.gt.tqm) then tbeta= tbeta*(1.d0+3.d0/16.d0/pis*(rhts-rhbs)* # log(am/tqm)) tbetas= tbeta*tbeta sbeta= tbeta/sqrt(1.d0+tbeta*tbeta) cbeta= 1.d0/sqrt(1.d0+tbeta*tbeta) sbetas= sbeta*sbeta cbetas= cbeta*cbeta endif tvar= 2.d0*log(scalm/tqm) tvars= tvar*tvar bxt= 2.d0*bats/scalms*(1.d0-bats/12.d0/scalms) bxb= 2.d0*babs/scalms*(1.d0-babs/12.d0/scalms) batb= 1.d0/6.d0*(-6.d0*rmus/scalms-(rmus-bat*bab)* # (rmus-bat*bab)/scalmq+3.d0*(bat+bab)*(bat+bab)/ # scalms) * hcl1= 0.25d0*(gcg1s+gcg2s)*(1.d0-3.d0/8.d0/pis*rhbs*tvar)+ # 3.d0/8.d0/pis*rhbq*(tvar+0.5d0*bxb+1.d0/16.d0/pis*( # 3.d0/2.d0*rhbs+0.5d0*rhts-8.d0*gcg3s)*(bxb+tvar)*tvar)- # 3.d0/96.d0/pis*rhtq*rmuq/scalmq*(1.d0+1.d0/16.d0/pis*( # 9.d0*rhts-5.d0*rhbs-16.d0*gcg3s)*tvar) hcl2= 0.25d0*(gcg1s+gcg2s)*(1.d0-3.d0/8.d0/pis*rhts*tvar)+ # 3.d0/8.d0/pis*rhtq*(tvar+0.5d0*bxt+1.d0/16.d0/pis*( # 3.d0/2.d0*rhts+0.5d0*rhbs-8.d0*gcg3s)*(bxt+tvar)*tvar)- # 3.d0/96.d0/pis*rhbq*rmuq/scalmq*(1.d0+1.d0/16.d0/pis*( # 9.d0*rhbs-5.d0*rhts-16.d0*gcg3s)*tvar) hcl3= 0.25d0*(gcg1s-gcg2s)*(1.d0-3.d0/16.d0/pis*(rhbs+rhts)* # tvar)+ # 3.d0/8.d0/pis*rhbs*rhts*(tvar+0.5d0*batb+1.d0/16.d0/pis*( # rhbs+rhts-8.d0*gcg3s)*(batb+tvar)*tvar)+ # 3.d0/96.d0/pis*rhtq*(3.d0*rmus/scalms-rmus*bats/scalmq)* # (1.d0+1.d0/16.d0/pis*( # 6.d0*rhts-2.d0*rhbs-16.d0*gcg3s)*tvar)+ # 3.d0/96.d0/pis*rhbq*(3.d0*rmus/scalms-rmus*babs/scalmq)* # (1.d0+1.d0/16.d0/pis*( # 6.d0*rhbs-2.d0*rhts-16.d0*gcg3s)*tvar) hcl4= -0.5d0*gcg2s*(1.d0-3.d0/16.d0/pis*(rhbs+rhts)* # tvar)- # 3.d0/8.d0/pis*rhbs*rhts*(tvar+0.5d0*batb+1.d0/16.d0/pis*( # rhbs+rhts-8.d0*gcg3s)*(batb+tvar)*tvar)+ # 3.d0/96.d0/pis*rhtq*(3.d0*rmus/scalms-rmus*bats/scalmq)* # (1.d0+1.d0/16.d0/pis*( # 6.d0*rhts-2.d0*rhbs-16.d0*gcg3s)*tvar)+ # 3.d0/96.d0/pis*rhbq*(3.d0*rmus/scalms-rmus*babs/scalmq)* # (1.d0+1.d0/16.d0/pis*( # 6.d0*rhbs-2.d0*rhts-16.d0*gcg3s)*tvar) hcl5= -3.d0/96.d0/pis*rhtq*rmus*bats/scalmq*(1.d0-1.d0/16.d0/ # pis*(2.d0*rhbs-6.d0*rhts+16.d0*gcg3s)*tvar)- # 3.d0/96.d0/pis*rhbq*rmus*babs/scalmq*(1.d0-1.d0/16.d0/ # pis*(2.d0*rhts-6.d0*rhbs+16.d0*gcg3s)*tvar) hcl6= 3.d0/96.d0/pis*rhtq*rmuc*bat/scalmq*(1.d0-1.d0/16.d0/ # pis*(7.d0/2.d0*rhbs-15.d0/2.d0*rhts+16.d0*gcg3s)*tvar)+ # 3.d0/96.d0/pis*rhbq*rmu/scalms*bab*(babs/scalms-6.d0)* # (1.d0-1.d0/16.d0/pis*(0.5d0*rhts-9.d0/2.d0*rhbs+16.d0* # gcg3s)*tvar) hcl7= 3.d0/96.d0/pis*rhbq*rmuc*bab/scalmq*(1.d0-1.d0/16.d0/ # pis*(7.d0/2.d0*rhts-15.d0/2.d0*rhbs+16.d0*gcg3s)*tvar)+ # 3.d0/96.d0/pis*rhtq*rmu/scalms*bat*(bats/scalms-6.d0)* # (1.d0-1.d0/16.d0/pis*(0.5d0*rhbs-9.d0/2.d0*rhts+16.d0* # gcg3s)*tvar) * rm12s= 2.d0*vevs*(sbeta*cbeta*(hcl3+hcl4)+hcl6*cbetas+ # hcl7*sbetas)-ams*sbeta*cbeta rm11s= 2.d0*vevs*(hcl1*cbetas+2.d0*hcl6*cbeta*sbeta+ # hcl5*sbetas)+ams*sbetas rm22s= 2.d0*vevs*(hcl2*sbetas+2.d0*hcl7*cbeta*sbeta+ # hcl5*cbetas)+ams*cbetas trm2= rm11s+rm22s dtm2= rm11s*rm22s-rm12s*rm12s * shms= 0.5d0*(trm2-sqrt(trm2*trm2-4.d0*dtm2)) bhms= 0.5d0*(trm2+sqrt(trm2*trm2-4.d0*dtm2)) shm= sqrt(shms) bhm= sqrt(bhms) chms= ams+(hcl5-hcl4)*vevs chm= sqrt(chms) * c2alpha= (rm11s-rm22s)/sqrt(trm2*trm2-4.d0*dtm2) salphas= 0.5d0*(1.d0-c2alpha) salpha= sqrt(salphas) calphas= 0.5d0*(1.d0+c2alpha) calpha= sqrt(calphas) * sbma= sbeta*calpha-cbeta*salpha cbma= cbeta*calpha+sbeta*salpha * call wtocorrqcd(shm,als,bqm2,cqm2,shbm2,shcm2,shsm2) call wtocorrqcd(bhm,als,bqm2,cqm2,bhbm2,bhcm2,bhsm2) call wtocorrqcd(am,als,bqm2,cqm2,abm2,acm2,asm2) call wtocorrqcd(chm,als,bqm2,cqm2,chbm2,chcm2,chsm2) alssh= wtoralphas(wm,shm,als,nf) alsbh= wtoralphas(wm,bhm,als,nf) alsa= wtoralphas(wm,am,als,nf) alsch= wtoralphas(wm,chm,als,nf) if(oqcd.eq.'y') then shg= qgfopi*shm*(3.d0*(salphas/cbetas*shbm2+ # calphas/sbetas*shcm2)*(1.d0+alssh/pi* # 17.d0/3.d0)+salphas/cbetas*tm*tm) bhg= qgfopi*bhm*(3.d0*(calphas/cbetas*bhbm2+ # salphas/sbetas*bhcm2)*(1.d0+alsbh/pi* # 17.d0/3.d0)+calphas/cbetas*tm*tm) ag= qgfopi*am*(3.d0*(tbetas*abm2+ # acm2/tbetas)*(1.d0+alsa/pi* # 17.d0/3.d0)+tbetas*tm*tm) chg= qgfopi*chm*(3.d0*(chcm2/tbetas+chsm2*tbetas)* # (1.d0+17.d0/3.d0*alsch/pi)+tm*tm*tbetas) else if(oqcd.eq.'n') then shg= qgfopi*shm*(3.d0*(salphas/cbetas*shbm2+ # calphas/sbetas*shcm2)+salphas/cbetas*tm*tm) bhg= qgfopi*bhm*(3.d0*(calphas/cbetas*bhbm2+ # salphas/sbetas*bhcm2)+calphas/cbetas*tm*tm) ag= qgfopi*am*(3.d0*(tbetas*abm2+ # acm2/tbetas)+tbetas*tm*tm) chg= qgfopi*chm*(3.d0*(chcm2/tbetas+chsm2*tbetas)+ # tm*tm*tbetas) endif * rshm= shm/rs rshm2= shms/s rshg= shg/rs rshmg= rshm*rshg sshg= shg/shm sshgs= sshg*sshg opsshgs= 1.d0+sshgs rbhm= bhm/rs rbhm2= bhms/s rbhg= bhg/rs rbhmg= rbhm*rbhg sbhg= bhg/bhm sbhgs= sbhg*sbhg opsbhgs= 1.d0+sbhgs ram= am/rs ram2= ams/s rag= ag/rs ramg= ram*rag sag= ag/am sags= sag*sag opsags= 1.d0+sags rchm= chm/rs rchm2= chms/s rchg= chg/rs rchmg= rchm*rchg schg= chg/chm schgs= schg*schg opschgs= 1.d0+schgs endif if(omssm.eq.'n') then rbqm2= rbm2/s else if(omssm.eq.'y') then rbqm2= shbm2/s endif * chf2= chf*chf chfp2= chfp*chfp conc(1)= chf*chfp*sth4/16.d0 conc(2)= 0.25d0*chf*tth2 conc(3)= 0.25d0*chfp*tth2 conc(4)= 0.25d0*chf*chfp*tth2 conc(5)= 6.25d-2/cth4 conc(6)= chf*chfp2*sth4/16.d0 conc(7)= chf2*chfp*sth4/16.d0 conc(8)= 3.125d-2/cth2 conc(9)= chf*sth2/8.d0 conc(10)= chf*sth2/16.d0 * ve= -0.5d0+2.d0*sth2 vf= tif-2.d0*chf*sth2 vfp= tifp-2.d0*chfp*sth2 vel= ve-0.5d0 vfl= vf+tif vfpl= vfp+tifp ver= ve+0.5d0 velr= ver*vel vfr= vf-tif vfpr= vfp-tifp vel2= vel*vel ver2= ver*ver vfl2= vfl*vfl vfr2= vfr*vfr vfpl2= vfpl*vfpl vfpr2= vfpr*vfpr vdp= -0.5d0-2.d0*chdp*sth2 vdpl= vdp-0.5d0 vdpr= vdp+0.5d0 vup= 0.5d0-2.d0*chup*sth2 vupl= vup+0.5d0 vupr= vup-0.5d0 if(otype.eq.'nc48'.or.otype.eq.'nc50') then vqr= vfr vql= vfl ver3= ver2*ver vel3= vel2*vel vqr2= vqr*vqr vql2= vql*vql * hbo(1)= ver*vql hbo(2)= velr hbo(3)= ver2*vel*vql hbo(4)= vel*vqr hbo(5)= velr hbo(6)= ver*vel2*vqr hbo(7)= ver*vqr hbo(8)= velr hbo(9)= ver2*vel*vqr hbo(10)= vel*vql hbo(11)= velr hbo(12)= ver*vel2*vql hbo(13)= ver*vqr hbo(14)= ver2 hbo(15)= ver3*vqr hbo(16)= vel*vql hbo(17)= vel2 hbo(18)= vel3*vql hbo(19)= ver*vql hbo(20)= ver2 hbo(21)= ver3*vql hbo(22)= vel*vqr hbo(23)= vel2 hbo(24)= vel3*vqr * hbe(1)= vel*vql hbe(2)= velr hbe(3)= ver*vel2*vql hbe(4)= ver*vqr hbe(5)= velr hbe(6)= ver2*vel*vqr hbe(7)= vel*vqr hbe(8)= velr hbe(9)= ver*vel2*vqr hbe(10)= ver*vql hbe(11)= velr hbe(12)= ver2*vel*vql hbe(13)= ver*vqr hbe(14)= ver2 hbe(15)= ver3*vqr hbe(16)= vel*vql hbe(17)= vel2 hbe(18)= vel3*vql hbe(19)= ver*vql hbe(20)= ver2 hbe(21)= ver3*vql hbe(22)= vel*vqr hbe(23)= vel2 hbe(24)= vel3*vqr * hmp(1)= vel*vql hmp(2)= ver*vql hmp(3)= velr*vql2 hmp(4)= ver*vqr hmp(5)= vel*vqr hmp(6)= velr*vqr2 hmp(7)= vel*vqr hmp(8)= ver*vqr hmp(9)= velr*vqr2 hmp(10)= ver*vql hmp(11)= vel*vql hmp(12)= velr*vql2 hmp(13)= ver*vqr hmp(14)= ver*vqr hmp(15)= ver2*vqr2 hmp(16)= vel*vql hmp(17)= vel*vql hmp(18)= vel2*vql2 hmp(19)= ver*vql hmp(20)= ver*vql hmp(21)= ver2*vql2 hmp(22)= vel*vqr hmp(23)= vel*vqr hmp(24)= vel2*vqr2 * else vqr= 0.d0 vql= 0.d0 do i=1,24 hbe(i)= 0.d0 hbo(i)= 0.d0 hmp(i)= 0.d0 enddo endif * hch(1)= vfl*ver/16.d0 hch(2)= vfr*vel/16.d0 hch(3)= vfr*ver/16.d0 hch(4)= vfl*vel/16.d0 hch(5)= vfpl*ver/16.d0 hch(6)= vfpr*vel/16.d0 hch(7)= vfpr*ver/16.d0 hch(8)= vfpl*vel/16.d0 hch(9)= ver2*vfl*vfpl/16.d0 hch(10)= vel2*vfr*vfpr/16.d0 hch(11)= ver2*vfl*vfpr/16.d0 hch(12)= vel2*vfr*vfpl/16.d0 hch(13)= ver2*vfr*vfpl/16.d0 hch(14)= vel2*vfl*vfpr/16.d0 hch(15)= ver2*vfr*vfpr/16.d0 hch(16)= vel2*vfl*vfpl/16.d0 hch(17)= vfr*vfpl/16.d0 hch(18)= vfl*vfpr/16.d0 hch(19)= vfl*vfpl/16.d0 hch(20)= vfr*vfpr/16.d0 hch(21)= ver*vfr*vfpl2/16.d0 hch(22)= vel*vfl*vfpr2/16.d0 hch(23)= ver*vfl*vfpl2/16.d0 hch(24)= vel*vfr*vfpr2/16.d0 hch(25)= ver*vfr*vfpr2/16.d0 hch(26)= vel*vfl*vfpl2/16.d0 hch(27)= vel*vfr*vfpl2/16.d0 hch(28)= ver*vfl*vfpr2/16.d0 hch(29)= ver*vfr2*vfpl/16.d0 hch(30)= vel*vfl2*vfpr/16.d0 hch(31)= ver*vfl2*vfpl/16.d0 hch(32)= vel*vfr2*vfpr/16.d0 hch(33)= ver*vfr2*vfpr/16.d0 hch(34)= vel*vfl2*vfpl/16.d0 hch(35)= ver*vfl2*vfpr/16.d0 hch(36)= vel*vfr2*vfpl/16.d0 * hchup= tsth2*chup hchu= tsth2*chu hchdp= tsth2*chdp hchd= tsth2*chd * *-----CC * g2f= 8.d0*gf*wm2 cth= sqrt(cth2) gca= sqrt(g2) gcf= sqrt(g2f) gvea= 0.25d0*gca/cth*(4.d0*sth2-1.d0) gaea= -0.25d0*gca/cth gwfa= 0.5d0*gca/sqrt(2.d0) g3va= -gca*cth gvef= 0.25d0*gcf/cth*(4.d0*sth2-1.d0) gaef= -0.25d0*gcf/cth gwff= 0.5d0*gcf/sqrt(2.d0) g3vf= -gcf*cth * *-----e.m. quantities are initialized * aexp= alpha/pi rln= log(s/em2) beta= 2.d0*aexp*(rln-1.d0) isf= iosf if(iosf.eq.1) then eta= beta eob= 1.d0 else if(iosf.gt.1) then eta= 2.d0*aexp*rln eob= rln/(rln-1.d0) endif feta= eta/16.d0 hbet= beta/2.d0 heta= eta/2.d0 hbeti= 2.d0/beta omhb= 1.d0-hbet ophb= 1.d0+hbet if(iosf.lt.3) then arge= hbet*(3.d0/4.d0-ge) else if(iosf.eq.3) then arge= -hbet*ge+3.d0/4.d0*heta endif ifg= 1 gamb= s14aaf(ophb,ifg) if(ifg.ne.0) then stop else d0gl= exp(arge)/gamb endif * *-----splitting&transformationsh * pin= 256.d0*pi**7 if(oxcm.eq.'c') then tfact= g8/128.d0/pin*fcuc*fcdc*cfct else if(oxcm.eq.'n') then if((otype.eq.'nc64'.or.otype.eq.'nc68'). # and.ofs.eq.'qq') then tfact= g8/128.d0/pin*cfct else tfact= g8/128.d0/pin*fcun*fcdn*cfct endif else if(oxcm.eq.'m') then tfact= g8/128.d0/pin*cfct endif if(otype.eq.'nc68') then tfacth= g8/128.d0/pin*fcun*fcdn*cfct else tfacth= 0.d0 endif * do i=1,6 rl(i)= aim(i)*aim(i) rr(i)= bim(i)*bim(i) srl(i)= sqrt(rl(i)) enddo dsm= rl(1) usm= rr(1) dsp= rl(6) usp= rr(6) * sdsm= sqrt(dsm) sdsp= sqrt(dsp) * sdsm= sqrt(dsm) sdsp= sqrt(dsp) * *-----limits for u from cuts on IM * vvl1= (sdsm+sdsp)*(sdsm+sdsp) vvl2= (srl(2)+srl(5))*(srl(2)+srl(5)) vvl3= (srl(3)+srl(4))*(srl(3)+srl(4)) ul= dmax1(vvl1,vvl2,vvl3) omul= 1.d0-ul suml= 0.d0 do j=1,4 suml= suml+rae(j) enddo * if(oxcm.eq.'c') then if(omdist.eq.'y') then ndim= 2 else if(itc.eq.11) then if(isf.eq.0) then Print*,' Born not allowed ' stop else ndim= 7 endif else if(itc.ge.7) then if(isf.eq.0) then ndim= 6 else ndim= 8 endif else if(isf.eq.0) then ndim= 7 else ndim= 9 endif endif endif else if(oxcm.eq.'n'.or.oxcm.eq.'m') then if(itc.gt.0) then if(isf.eq.0) then ndim= 6 else ndim= 8 endif else if(isf.eq.0) then ndim= 7 else ndim= 9 endif endif else if(isf.eq.0) then ndim= 7 else ndim= 9 endif endif * if(ndim.eq.2) then if(inpts.le.6) then jnpts= inpts else if(inpts.eq.7) then vk(1)= 1.d0 vk(2)= 0.141603d6 endif else if(ndim.eq.6.or.ndim.eq.8) then if(inpts.le.6) then jnpts= inpts else if(inpts.eq.8) then jnpts= 10193*101 if(ndim.eq.6) then vk(1)= 1.d0 vk(2)= 0.672849d+06 vk(3)= 0.53093d+05 vk(4)= 0.164857d+06 vk(5)= 0.114815d+06 vk(6)= 0.3215d+04 else if(ndim.eq.8) then vk(1)= 1.d0 vk(2)= 0.757894d+06 vk(3)= 0.784365d+06 vk(4)= 0.236855d+06 vk(5)= 0.347946d+06 vk(6)= 0.524281d+06 vk(7)= 0.128976d+06 vk(8)= 0.805687d+06 endif else if(inpts.eq.9) then jnpts= 22807*151 if(ndim.eq.6) then vk(1)= 1.d0 vk(2)= 0.2903683d+07 vk(3)= 0.278237d+06 vk(4)= 0.413956d+06 vk(5)= 0.1366666d+07 vk(6)= 0.1522064d+07 else if(ndim.eq.8) then vk(1)= 1.d0 vk(2)= 0.1464524d+07 vk(3)= 0.329469d+07 vk(4)= 0.2417287d+07 vk(5)= 0.33812d+05 vk(6)= 0.2709542d+07 vk(7)= 0.1615901d+07 vk(8)= 0.249863d+06 endif else print*,' NPTS choice is Not allowed ' stop endif else if(inpts.le.6) then jnpts= inpts else if(inpts.eq.7) then jnpts= 2011*101 if(ndim.eq.7) then vk(1)= 0.1d+1 vk(2)= 0.72725d+5 vk(3)= 0.118296d+6 vk(4)= 0.107084d+6 vk(5)= 0.1938d+4 vk(6)= 0.185127d+6 vk(7)= 0.14844d+6 else if(ndim.eq.9) then vk(1)= 0.1d+1 vk(2)= 0.82871d+5 vk(3)= 0.13509d+5 vk(4)= 0.159618d+6 vk(5)= 0.99403d+5 vk(6)= 0.53186d+5 vk(7)= 0.68306d+5 vk(8)= 0.86067d+5 vk(9)= 0.12481d+5 endif else if(inpts.eq.8) then jnpts= 10193*101 if(ndim.eq.7) then vk(1)= 0.1d+1 vk(2)= 0.1019431d+7 vk(3)= 0.35353d+6 vk(4)= 0.708948d+6 vk(5)= 0.951714d+6 vk(6)= 0.197618d+6 vk(7)= 0.54816d+6 else if(ndim.eq.9) then vk(1)= 0.1d+01 vk(2)= 0.755242d+6 vk(3)= 0.911407d+6 vk(4)= 0.442285d+6 vk(5)= 0.850204d+6 vk(6)= 0.572366d+6 vk(7)= 0.1026802d+7 vk(8)= 0.892453d+6 vk(9)= 0.685582d+6 endif else if(inpts.eq.9) then jnpts= 22807*151 if(ndim.eq.7) then vk(1)= 0.1d+01 vk(2)= 0.1619459d+7 vk(3)= 0.226133d+7 vk(4)= 0.256381d+7 vk(5)= 0.230245d+7 vk(6)= 0.855081d+6 vk(7)= 0.609193d+6 else if(ndim.eq.9) then vk(1)= 0.1d+01 vk(2)= 0.2055793d+7 vk(3)= 0.767734d+6 vk(4)= 0.3183104d+7 vk(5)= 0.2813063d+7 vk(6)= 0.2463708d+7 vk(7)= 0.58258d+5 vk(8)= 0.2817562d+7 vk(9)= 0.1276513d+7 endif else if(inpts.eq.10) then jnpts= 32771*181 if(ndim.eq.7) then vk(1)= 0.1d+01 vk(2)= 0.4960373d+7 vk(3)= 0.4851623d+7 vk(4)= 0.3262017d+7 vk(5)= 0.722217d+6 vk(6)= 0.4644124d+7 vk(7)= 0.3212165d+7 else if(ndim.eq.9) then vk(1)= 0.1d+01 vk(2)= 0.2498015d+7 vk(3)= 0.4246511d+7 vk(4)= 0.4724489d+7 vk(5)= 0.3448063d+7 vk(6)= 0.1119927d+7 vk(7)= 0.2141959d+7 vk(8)= 0.115857d+7 vk(9)= 0.287463d+7 endif else if(inpts.eq.11) then jnpts= 32771*379 if(ndim.eq.7) then vk(1)= 0.1d+01 vk(2)= 0.9646626d+07 vk(3)= 0.8128723d+07 vk(4)= 0.9521278d+07 vk(5)= 0.4720279d+07 vk(6)= 0.7036407d+07 vk(7)= 0.1868554d+07 else if(ndim.eq.9) then vk(1)= 0.1d+01 vk(2)= 0.4365962d+07 vk(3)= 0.650771d+07 vk(4)= 0.6559665d+07 vk(5)= 0.9838408d+07 vk(6)= 0.2561851d+07 vk(7)= 0.11842175d+08 vk(8)= 0.208211d+06 vk(9)= 0.6217272d+07 endif endif endif * ik= 0 do i=1,51 ifz(i)= 0 enddo * *-----Not available for changing * otrans= 'n' * *-----4(8) xs and 4(8) max are computed * *-----begin generation-evaluation branch * if(om.eq.'g') then if(osm.eq.'n') then * *-----all CC but CC12, XS is SIGAl * if(oxcm.eq.'c'.and.otype.ne.'cc12') then do i0=1,2 do j0=1,2 ix0= i0 it0= j0 ip= 2*(i0-1)+j0 pa0= 0.d0 pb0= 1.d0 np0= 9 call g05faf(pa0,pb0,np0,prx0) do i=1,npos xscmx(i)= 0.d0 enddo npts= jnpts itrans= 0 nrand= inrand ifail= 0 call g05cbf(0) call d01gcf(ndim,wtoxsc,wtoregion,npts,vk,nrand, # itrans,sig,esig,ifail) siga(ip)= sig do i=1,npos xscmx0(i)= xscmx(i) enddo * *-----Improving upon the maximum * ibound= 0 do i=1,npos call wtomxregion(i,bl,bu) istop= 0 do j=1,9 if((xmx(i,j).lt.bl(j)).or. # (xmx(i,j).gt.bu(j))) then istop= istop+1 endif enddo if(istop.gt.0.or.ostop.eq.'s') then xsmx(i,ip)= xscmx0(i) else irst= 0 if(xscmx0(i).eq.0.d0) then xscmx0(i)= 1.d0 do j=1,9 xmx0(j)= 0.5d0*(bl(j)+bu(j)) enddo else do j=1,9 xmx0(j)= xmx(i,j) enddo endif 1200 ifail= 1 ireg= i call e04jaf(ndim,ibound,bl,bu,xmx0,fmx, # iw,liw,wmx,lw,ifail) if(ifail.eq.2) then if(irst.eq.0) then irst= 1 go to 1200 else tmx1= xscmx0(i) tmx2= -fmx*xscmx0(i) xsmx(i,ip)= dmax1(tmx1,tmx2) endif else if(ifail.gt.5) then xsmx(i,ip)= xscmx0(i) else xsmx(i,ip)= -fmx*xscmx0(i) endif endif enddo * if(ip.eq.1) then print 1999,ip,siga(ip) open(nin,file='dummy',status='new') do i=1,9 write(nin,fmt=*) prx0(i) enddo do i=1,npos write(nin,fmt=*) siga(ip),xsmx(i,ip) enddo close(nin) else print 1999,ip,siga(ip) open(nin,file='dummy',status='old',access='append') do i=1,npos write(nin,fmt=*) siga(ip),xsmx(i,ip) enddo close(nin) endif enddo enddo else if(oxcm.eq.'n'.or.(oxcm.eq.'c'.and. # otype.eq.'cc12').or.oxcm.eq.'m') then pa0= 0.d0 pb0= 1.d0 np0= 9 call g05faf(pa0,pb0,np0,prx0) * *----- NC19/NC24/NC48/NC64, XS is SIGAN * if(otype.eq.'nc24'.or.otype.eq.'nc19'.or. # otype.eq.'nc48'.or.otype.eq.'nc64') then do i0=1,2 do j0=1,2 ix0= i0 it0= j0 ip= 2*(i0-1)+j0 do i=1,npos xshmx(i)= 0.d0 enddo npts= jnpts itrans= 0 nrand= inrand ifail= 0 if(otype.eq.'nc64') then do j=1,npos ifl(j)= 1 enddo call g05cbf(0) iint= 0 call d01gcf(ndim,wtoxsn64,wtoregion,npts,vk,nrand, # itrans,sig,esig,ifail) else call g05cbf(0) call d01gcf(ndim,wtoxsn,wtoregion,npts,vk,nrand, # itrans,sig,esig,ifail) endif sigan(ip)= sig * if(ostop.eq.'s') then do j=1,npos do i=1,ifl(j) axstr(i)= stry(j,i) enddo ifail= 0 mm1= 1 mm2= ifl(j) ord= 'a' call m01daf(axstr,mm1,mm2,ord,ir,ifail) do i=1,ifl(j) stry(j,ir(i))= axstr(i) enddo if(ifl(j).eq.2) then xshmx0(j)= stry(j,ifl(j)) else if(ifl(j).eq.1) then xshmx0(j)= 0.d0 else ifmx= ifl(j) ifm= ifl(j)-1 1133 if(ifm.eq.1) then xshmx0(j)= stry(j,ifmx) else avr= 0.d0 avr2= 0.d0 do i=1,ifm avr= avr+stry(j,i) avr2= avr2+stry(j,i)*stry(j,i) enddo avr= avr/ifm avr2= avr2/ifm std2= ifm/(ifm-1.d0)*(avr2-avr*avr) std= sqrt(std2) dev= (stry(j,ifmx)-avr)/std if(dev.lt.6.d0) then isort= 1 else isort= 0 endif if(isort.eq.0) then ifm= ifm-1 ifmx= ifmx-1 go to 1133 else xshmx0(j)= stry(j,ifm) endif endif endif enddo else do i=1,npos xshmx0(i)= xshmx(i) enddo endif * ibound= 0 do i=1,npos call wtomxregion(i,bl,bu) istop= 0 if(ostop.eq.'i') then do j=1,9 if((xmxh(i,j).lt.bl(j)).or. # (xmxh(i,j).gt.bu(j))) then istop= istop+1 endif enddo endif if(istop.gt.0.or.ostop.eq.'s') then xsmxn(i,ip)= xshmx0(i) else irst= 0 if(xshmx0(i).eq.0.d0) then xshmx0(i)= 1.d0 do j=1,9 xmx0(j)= 0.5d0*(bl(j)+bu(j)) enddo else do j=1,9 xmx0(j)= xmxh(i,j) enddo endif 1930 ifail= 1 ireg= i call e04jaf(ndim,ibound,bl,bu,xmx0,fmx, # iw,liw,wmx,lw,ifail) if(ifail.eq.2) then if(irst.eq.0) then irst= 1 go to 1930 else tmx1= xshmx0(i) tmx2= -fmx*xshmx0(i) xsmxn(i,ip)= dmax1(tmx1,tmx2) endif else if(ifail.gt.5) then xsmxn(i,ip)= xshmx0(i) else xsmxn(i,ip)= -fmx*xshmx0(i) endif endif enddo * if(ip.eq.1) then print 1999,ip,sigan(ip) open(nin,file='dummy',status='new') do i=1,9 write(nin,fmt=*) prx0(i) enddo do i=1,npos write(nin,fmt=*) sigan(ip),xsmxn(i,ip) enddo close(nin) else print 1999,ip,sigan(ip) open(nin,file='dummy',status='old',access='append') do i=1,npos write(nin,fmt=*) sigan(ip),xsmxn(i,ip) enddo close(nin) endif enddo enddo * *---- NC68, XS is SIGAHB * else if(otype.eq.'nc68') then * if(omssm.eq.'n') then k0mx= 5 else if(omssm.eq.'y') then k0mx= 9 endif do i0=1,2 do j0=1,2 do k0=1,k0mx ix0= i0 it0= j0 ip0= k0 ip= k0mx*(2*(i0-1)+j0-1)+k0 do i=1,npos xshmx(i)= 0.d0 enddo npts= jnpts itrans= 0 nrand= inrand ifail= 0 do j=1,npos ifl(j)= 1 enddo call g05cbf(0) if(omssm.eq.'n') then call d01gcf(ndim,wtoxsnh64,wtoregion,npts,vk,nrand, # itrans,sig,esig,ifail) else if(omssm.eq.'y') then call d01gcf(ndim,wtoxsnsh64,wtoregion,npts,vk,nrand, # itrans,sig,esig,ifail) endif sigahb(ip)= sig * if(ostop.eq.'s') then do j=1,npos do i=1,ifl(j) axstr(i)= stry(j,i) enddo ifail= 0 mm1= 1 mm2= ifl(j) ord= 'a' call m01daf(axstr,mm1,mm2,ord,ir,ifail) do i=1,ifl(j) stry(j,ir(i))= axstr(i) enddo if(ifl(j).eq.2) then xshmx0(j)= stry(j,ifl(j)) else if(ifl(j).eq.1) then xshmx0(j)= 0.d0 else ifmx= ifl(j) ifm= ifl(j)-1 1233 if(ifm.eq.1) then xshmx0(j)= stry(j,ifmx) else avr= 0.d0 avr2= 0.d0 do i=1,ifm avr= avr+stry(j,i) avr2= avr2+stry(j,i)*stry(j,i) enddo avr= avr/ifm avr2= avr2/ifm std2= ifm/(ifm-1.d0)*(avr2-avr*avr) std= sqrt(std2) dev= (stry(j,ifmx)-avr)/std if(dev.lt.6.d0) then isort= 1 else isort= 0 endif if(isort.eq.0) then ifm= ifm-1 ifmx= ifmx-1 go to 1233 else xshmx0(j)= stry(j,ifm) endif endif endif enddo else do i=1,npos xshmx0(i)= xshmx(i) enddo endif * ibound= 0 do i=1,npos call wtomxregion(i,bl,bu) istop= 0 if(ostop.eq.'i') then do j=1,9 if((xmxh(i,j).lt.bl(j)).or. # (xmxh(i,j).gt.bu(j))) then istop= istop+1 endif enddo endif if(istop.gt.0.or.ostop.eq.'s') then xsmxhb(i,ip)= xshmx0(i) else irst= 0 if(xshmx0(i).eq.0.d0) then xshmx0(i)= 1.d0 do j=1,9 xmxh0(j)= 0.5d0*(bl(j)+bu(j)) enddo else do j=1,9 xmxh0(j)= xmxh(i,j) enddo endif 6901 ifail= 1 ireg= i call e04jaf(ndim,ibound,bl,bu,xmxh0,fmxh, # iw,liw,wmx,lw,ifail) if(ifail.eq.2) then if(irst.eq.0) then irst= 1 go to 6901 else tmx1= xshmx0(i) tmx2= -fmxh*xshmx0(i) xsmxhb(i,ip)= dmax1(tmx1,tmx2) endif else if(ifail.gt.5) then xsmxhb(i,ip)= xshmx0(i) else xsmxhb(i,ip)= -fmxh*xshmx0(i) endif endif enddo if(ip.eq.1) then print 1999,ip,sigahb(ip) open(nin,file='dummy',status='new') do i=1,9 write(nin,fmt=*) prx0(i) enddo do i=1,npos write(nin,fmt=*) sigahb(ip),xsmxhb(i,ip) enddo close(nin) else print 1999,ip,sigahb(ip) open(nin,file='dummy',status='old',access='append') do i=1,npos write(nin,fmt=*) sigahb(ip),xsmxhb(i,ip) enddo close(nin) endif enddo enddo enddo * *---- NC26/NC33, XS is SIGAHQ5 * else if(otype.eq.'nc33'.or.otype.eq.'nc26') then * do i0=1,2 do j0=1,2 do k0=1,3 ix0= i0 it0= j0 ip0= k0 ip= 3*(2*(i0-1)+j0-1)+k0 do i=1,npos xshmx(i)= 0.d0 enddo npts= jnpts itrans= 0 nrand= inrand ifail= 0 call g05cbf(0) if(otype.eq.'nc33') then call d01gcf(ndim,wtoxsnh32,wtoregion,npts,vk,nrand, # itrans,sig,esig,ifail) else if(otype.eq.'nc26') then call d01gcf(ndim,wtoxsnh26,wtoregion,npts,vk,nrand, # itrans,sig,esig,ifail) endif sigahq(ip)= sig do i=1,npos xshmx0(i)= xshmx(i) enddo ibound= 0 do i=1,npos call wtomxregion(i,bl,bu) istop= 0 do j=1,9 if((xmxh(i,j).lt.bl(j)).or. # (xmxh(i,j).gt.bu(j))) then istop= istop+1 endif enddo if(istop.gt.0.or.ostop.eq.'s') then xsmxhq(i,ip)= xshmx0(i) else irst= 0 if(xshmx0(i).eq.0.d0) then xshmx0(i)= 1.d0 do j=1,9 xmxh0(j)= 0.5d0*(bl(j)+bu(j)) enddo else do j=1,9 xmxh0(j)= xmxh(i,j) enddo endif 1901 ifail= 1 ireg= i call e04jaf(ndim,ibound,bl,bu,xmxh0,fmxh, # iw,liw,wmx,lw,ifail) if(ifail.eq.2) then if(irst.eq.0) then irst= 1 go to 1901 else tmx1= xshmx0(i) tmx2= -fmxh*xshmx0(i) xsmxhq(i,ip)= dmax1(tmx1,tmx2) endif else if(ifail.gt.5) then xsmxhq(i,ip)= xshmx0(i) else xsmxhq(i,ip)= -fmxh*xshmx0(i) endif endif enddo if(ip.eq.1) then print 1999,ip,sigahq(ip) open(nin,file='dummy',status='new') do i=1,9 write(nin,fmt=*) prx0(i) enddo do i=1,npos write(nin,fmt=*) sigahq(ip),xsmxhq(i,ip) enddo close(nin) else print 1999,ip,sigahq(ip) open(nin,file='dummy',status='old',access='append') do i=1,npos write(nin,fmt=*) sigahq(ip),xsmxhq(i,ip) enddo close(nin) endif enddo enddo enddo * *---- NC50, XS is SIGAH * else if(otype.eq.'nc50') then * do i0=1,2 do j0=1,2 do k0=1,2 ix0= i0 it0= j0 ip0= k0 ip= 2*(2*(i0-1)+j0-1)+k0 do i=1,npos xshmx(i)= 0.d0 enddo npts= jnpts itrans= 0 nrand= inrand ifail= 0 call g05cbf(0) call d01gcf(ndim,wtoxsnh49,wtoregion,npts,vk,nrand, # itrans,sig,esig,ifail) sigah(ip)= sig do i=1,npos xshmx0(i)= xshmx(i) enddo ibound= 0 do i=1,npos call wtomxregion(i,bl,bu) istop= 0 do j=1,9 if((xmxh(i,j).lt.bl(j)).or. # (xmxh(i,j).gt.bu(j))) then istop= istop+1 endif enddo if(istop.gt.0.or.ostop.eq.'s') then xsmxh(i,ip)= xshmx0(i) else irst= 0 if(xshmx0(i).eq.0.d0) then xshmx0(i)= 1.d0 do j=1,9 xmxh0(j)= 0.5d0*(bl(j)+bu(j)) enddo else do j=1,9 xmxh0(j)= xmxh(i,j) enddo endif 1201 ifail= 1 ireg= i call e04jaf(ndim,ibound,bl,bu,xmxh0,fmxh, # iw,liw,wmx,lw,ifail) if(ifail.eq.2) then if(irst.eq.0) then irst= 1 go to 1201 else tmx1= xshmx0(i) tmx2= -fmxh*xshmx0(i) xsmxh(i,ip)= dmax1(tmx1,tmx2) endif else if(ifail.gt.5) then xsmxh(i,ip)= xshmx0(i) else xsmxh(i,ip)= -fmxh*xshmx0(i) endif endif enddo if(ip.eq.1) then print 1999,ip,sigah(ip) open(nin,file='dummy',status='new') do i=1,9 write(nin,fmt=*) prx0(i) enddo do i=1,npos write(nin,fmt=*) sigah(ip),xsmxh(i,ip) enddo close(nin) else print 1999,ip,sigah(ip) open(nin,file='dummy',status='old',access='append') do i=1,npos write(nin,fmt=*) sigah(ip),xsmxh(i,ip) enddo close(nin) endif enddo enddo enddo * *-----CC12/NC25/NC32/MX43, XS is SIGAH * else if(otype.eq.'nc25'.or.otype.eq.'nc32'.or. # otype.eq.'cc12'.or.otype.eq.'mx43') then * if(otype.eq.'mx43') then if(ofs.eq.'ll') then k0mx= 1 else if(ofs.eq.'qq'.and.opglu.eq.'n') then k0mx= 1 else if(ofs.eq.'qq'.and.opglu.eq.'y') then k0mx= 2 endif else if(otype.eq.'nc32') then if(opglu.eq.'n') then k0mx= 1 else if(opglu.eq.'y') then k0mx= 2 endif else k0mx= 2 endif * do i0=1,2 do j0=1,2 do k0=1,k0mx ix0= i0 it0= j0 ip0= k0 ip= k0mx*(2*(i0-1)+j0-1)+k0 do i=1,npos xshmx(i)= 0.d0 enddo npts= jnpts itrans= 0 nrand= inrand ifail= 0 if(otype.eq.'nc25') then call g05cbf(0) call d01gcf(ndim,wtoxsnh24,wtoregion,npts,vk,nrand, # itrans,sig,esig,ifail) else if(otype.eq.'nc32') then call g05cbf(0) call d01gcf(ndim,wtoxsn32,wtoregion,npts,vk,nrand, # itrans,sig,esig,ifail) else if(otype.eq.'cc12') then call g05cbf(0) call d01gcf(ndim,wtoxssc,wtoregion,npts,vk,nrand, # itrans,sig,esig,ifail) else if(otype.eq.'mx43') then iint= 0 call g05cbf(0) call d01gcf(ndim,wtoxsm43,wtoregion,npts,vk,nrand, # itrans,sig,esig,ifail) endif sigah(ip)= sig do i=1,npos xshmx0(i)= xshmx(i) enddo ibound= 0 * do i=1,npos call wtomxregion(i,bl,bu) istop= 0 do j=1,9 if((xmxh(i,j).lt.bl(j)).or. # (xmxh(i,j).gt.bu(j))) then istop= istop+1 endif enddo if(istop.gt.0.or.ostop.eq.'s') then xsmxh(i,ip)= xshmx0(i) else irst= 0 if(xshmx0(i).eq.0.d0) then xshmx0(i)= 1.d0 do j=1,9 xmxh0(j)= 0.5d0*(bl(j)+bu(j)) enddo else do j=1,9 xmxh0(j)= xmxh(i,j) enddo endif 1202 ifail= 1 ireg= i call e04jaf(ndim,ibound,bl,bu,xmxh0,fmxh, # iw,liw,wmx,lw,ifail) if(ifail.eq.2) then if(irst.eq.0) then irst= 1 go to 1202 else tmx1= xshmx0(i) tmx2= -fmxh*xshmx0(i) xsmxh(i,ip)= dmax1(tmx1,tmx2) endif else if(ifail.gt.5) then xsmxh(i,ip)= xshmx0(i) else xsmxh(i,ip)= -fmxh*xshmx0(i) endif endif enddo if(ip.eq.1) then print 1999,ip,sigah(ip) open(nin,file='dummy',status='new') do i=1,9 write(nin,fmt=*) prx0(i) enddo do i=1,npos write(nin,fmt=*) sigah(ip),xsmxh(i,ip) enddo close(nin) else print 1999,ip,sigah(ip) open(nin,file='dummy',status='old',access='append') do i=1,npos write(nin,fmt=*) sigah(ip),xsmxh(i,ip) enddo close(nin) endif enddo enddo enddo * *---- NC21, XS is SIGAH * else if(otype.eq.'nc21') then * do i0=1,2 do j0=1,2 do k0=1,2 ix0= i0 it0= j0 ip0= k0 ip= 2*(2*(i0-1)+j0-1)+k0 do i=1,npos xshmx(i)= 0.d0 enddo npts= jnpts itrans= 0 nrand= inrand ifail= 0 call g05cbf(0) call d01gcf(ndim,wtoxsnh19,wtoregion,npts,vk,nrand, # itrans,sig,esig,ifail) sigah(ip)= sig do i=1,npos xshmx0(i)= xshmx(i) enddo ibound= 0 do i=1,npos call wtomxregion(i,bl,bu) istop= 0 do j=1,9 if((xmxh(i,j).lt.bl(j)).or. # (xmxh(i,j).gt.bu(j))) then istop= istop+1 endif enddo if(istop.gt.0.or.ostop.eq.'s') then xsmxh(i,ip)= xshmx0(i) else irst= 0 if(xshmx0(i).eq.0.d0) then xshmx0(i)= 1.d0 do j=1,9 xmxh0(j)= 0.5d0*(bl(j)+bu(j)) enddo else do j=1,9 xmxh0(j)= xmxh(i,j) enddo endif 1203 ifail= 1 ireg= i call e04jaf(ndim,ibound,bl,bu,xmxh0,fmxh, # iw,liw,wmx,lw,ifail) if(ifail.eq.2) then if(irst.eq.0) then irst= 1 go to 1203 else tmx1= xshmx0(i) tmx2= -fmxh*xshmx0(i) xsmxh(i,ip)= dmax1(tmx1,tmx2) endif else if(ifail.gt.5) then xsmxh(i,ip)= xshmx0(i) else xsmxh(i,ip)= -fmxh*xshmx0(i) endif endif enddo if(ip.eq.1) then print 1999,ip,sigah(ip) open(nin,file='dummy',status='new') do i=1,9 write(nin,fmt=*) prx0(i) enddo do i=1,npos write(nin,fmt=*) sigah(ip),xsmxh(i,ip) enddo close(nin) else print 1999,ip,sigah(ip) open(nin,file='dummy',status='old',access='append') do i=1,npos write(nin,fmt=*) sigah(ip),xsmxh(i,ip) enddo close(nin) endif enddo enddo enddo endif endif * else if(osm.eq.'g') then if(itmx.gt.itmxmx) then print*,' Number of events exceedes the maximum allowed ' stop endif if(oxcm.eq.'c') then open(nin,file='dummy',status='old') do i=1,9 read(nin,fmt=*) prx0n(i) enddo do ip=1,4 do i=1,npos read(nin,fmt=*) siga(ip),xsmx(i,ip) enddo enddo close(nin) nint= 4 call wtovolume(npos,nint,xsmx,oxsmx,gbl,gbu,vol,volt) else if(oxcm.eq.'n'.or.oxcm.eq.'m') then if(otype.eq.'nc24'.or.otype.eq.'nc19'.or. # otype.eq.'nc48'.or.otype.eq.'nc64') then open(nin,file='dummy',status='old') do i=1,9 read(nin,fmt=*) prx0n(i) enddo do ip=1,4 do i=1,npos read(nin,fmt=*) sigan(ip),xsmxn(i,ip) enddo enddo close(nin) nint= 4 call wtovolume(npos,nint,xsmxn,oxsmxn,gbl,gbu,vol,volt) else if(otype.eq.'nc68') then open(nin,file='dummy',status='old') do i=1,9 read(nin,fmt=*) prx0n(i) enddo if(omssm.eq.'n') then ipmx= 20 else if(omssm.eq.'y') then ipmx= 36 endif do ip=1,ipmx do i=1,npos read(nin,fmt=*) sigahb(ip),xsmxhb(i,ip) enddo enddo close(nin) else if(otype.eq.'nc33'.or.otype.eq.'nc26') then open(nin,file='dummy',status='old') do i=1,9 read(nin,fmt=*) prx0n(i) enddo do ip=1,12 do i=1,npos read(nin,fmt=*) sigahq(ip),xsmxhq(i,ip) enddo enddo close(nin) else open(nin,file='dummy',status='old') do i=1,9 read(nin,fmt=*) prx0n(i) enddo if(otype.eq.'mx43') then if(ofs.eq.'ll') then ipmx= 4 else if(ofs.eq.'qq'.and.opglu.eq.'n') then ipmx= 4 else if(ofs.eq.'qq'.and.opglu.eq.'y') then ipmx= 8 endif else if(otype.eq.'nc32') then if(opglu.eq.'n') then ipmx= 4 else if(opglu.eq.'y') then ipmx= 8 endif else ipmx= 8 endif * do ip=1,ipmx do i=1,npos read(nin,fmt=*) sigah(ip),xsmxh(i,ip) enddo enddo close(nin) endif if(otype.eq.'nc68') then nint= ipmx call wtovolume(npos,nint,xsmxhb,oxsmxhb,gblhb,gbuhb, # volhb,volhtb) else if(otype.eq.'nc33'.or.otype.eq.'nc26') then nint= 12 call wtovolume(npos,nint,xsmxhq,oxsmxhq,gblhq,gbuhq, # volhq,volhtq) else nint= ipmx call wtovolume(npos,nint,xsmxh,oxsmxh,gblh,gbuh,volh, # volht) endif endif endif * *-----generation for CC starts * nlim= npos-1 do i= 1,itmx amps1(i)= 0.d0 amps2(i)= 0.d0 enddo if(oxcm.eq.'c') then * it= 1 sigt= siga(1)+siga(2)+siga(3)+siga(4) sigt2= siga(2)+siga(3)+siga(4) sigt3= siga(3)+siga(4) tsig= sigt 1504 wxs= g05caf(dmy) wxs0= siga(1)/sigt if(wxs.lt.wxs0) then ix0= 1 it0= 1 1500 jp= 1 volp= volt(1) 3000 test0= vol(jp,1)/volp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbu(jp,l,1)-gbl(jp,l,1))*aprx(l)+ # gbl(jp,l,1) enddo else if(jp.lt.nlim) then volp= volp-vol(jp,1) jp= jp+1 go to 3000 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbu(jp,l,1)-gbl(jp,l,1))*aprx(l)+ # gbl(jp,l,1) enddo endif endif pru= g05caf(pru) xt0= wtoxsc(ndim,prx) xt= xt0/siga(1) prh= oxsmx(jp,1)/siga(1) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xap(i) enddo amps(it)= xt0 it= it+1 else go to 1500 endif else wxs2= g05caf(dmy) wxs0= siga(2)/sigt2 if(wxs2.lt.wxs0) then ix0= 1 it0= 2 1501 jp= 1 volp= volt(2) 3001 test0= vol(jp,2)/volp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbu(jp,l,2)-gbl(jp,l,2))*aprx(l)+ # gbl(jp,l,2) enddo else if(jp.lt.nlim) then volp= volp-vol(jp,2) jp= jp+1 go to 3001 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbu(jp,l,2)-gbl(jp,l,2))*aprx(l)+ # gbl(jp,l,2) enddo endif endif pru= g05caf(pru) xt0= wtoxsc(ndim,prx) xt= xt0/siga(2) prh= oxsmx(jp,2)/siga(2) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xap(i) enddo amps(it)= xt0 it= it+1 else go to 1501 endif else wxs3= g05caf(dmy) wxs0= siga(3)/sigt3 if(wxs3.lt.wxs0) then ix0= 2 it0= 1 1502 jp= 1 volp= volt(3) 3002 test0= vol(jp,3)/volp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbu(jp,l,3)-gbl(jp,l,3))*aprx(l)+ # gbl(jp,l,3) enddo else if(jp.lt.nlim) then volp= volp-vol(jp,3) jp= jp+1 go to 3002 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbu(jp,l,3)-gbl(jp,l,3))*aprx(l)+ # gbl(jp,l,3) enddo endif endif pru= g05caf(pru) xt0= wtoxsc(ndim,prx) xt= xt0/siga(3) prh= oxsmx(jp,3)/siga(3) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xap(i) enddo amps(it)= xt0 it= it+1 else go to 1502 endif else ix0= 2 it0= 2 1503 jp= 1 volp= volt(4) 3003 test0= vol(jp,4)/volp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbu(jp,l,4)-gbl(jp,l,4))*aprx(l)+ # gbl(jp,l,4) enddo else if(jp.lt.nlim) then volp= volp-vol(jp,4) jp= jp+1 go to 3003 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbu(jp,l,4)-gbl(jp,l,4))*aprx(l)+ # gbl(jp,l,4) enddo endif endif pru= g05caf(pru) xt0= wtoxsc(ndim,prx) xt= xt0/siga(4) prh= oxsmx(jp,4)/siga(4) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xap(i) enddo amps(it)= xt0 it= it+1 else go to 1503 endif endif endif endif if(it.le.itmx) then go to 1504 endif * *-----generation for Higgs/NC/Mix starts * else if(oxcm.eq.'n'.or.oxcm.eq.'m') then * if(otype.eq.'nc24'.or.otype.eq.'nc19'.or. # otype.eq.'nc48'.or.otype.eq.'nc64') then if(otype.eq.'nc64'.and.ofs.eq.'qq'.and.opglu.eq.'n'.and. # oint.eq.'y') then iselsp= 1 else iselsp= 0 endif it= 1 sigt= sigan(1)+sigan(2)+sigan(3)+sigan(4) sigt2= sigan(2)+sigan(3)+sigan(4) sigt3= sigan(3)+sigan(4) tsig= sigt 1564 wxs= g05caf(dmy) wxs0= sigan(1)/sigt if(wxs.lt.wxs0) then ix0= 1 it0= 10 1560 jp= 1 volp= volt(1) 3060 test0= vol(jp,1)/volp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbu(jp,l,1)-gbl(jp,l,1))*aprx(l)+gbl(jp,l,1) enddo else if(jp.lt.nlim) then volp= volp-vol(jp,1) jp= jp+1 go to 3060 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbu(jp,l,1)-gbl(jp,l,1))*aprx(l)+gbl(jp,l,1) enddo endif endif pru= g05caf(pru) if(otype.eq.'nc64') then iint= 0 xt0= wtoxsn64(ndim,prx) if(iselsp.eq.1) then iint= 1 xt1= wtoxsn64(ndim,prx) iint= 2 xt2= wtoxsn64(ndim,prx) endif xt= xt0/sigan(1) else xt0= wtoxsn(ndim,prx) xt= xt0/sigan(1) endif prh= oxsmxn(jp,1)/sigan(1) prtst= xt-pru*prh if(prtst.ge.0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 if(iselsp.eq.1) then amps1(it)= xt1 amps2(it)= xt2 endif it= it+1 else go to 1560 endif else wxs2= g05caf(dmy) wxs0= sigan(2)/sigt2 if(wxs2.lt.wxs0) then ix0= 1 it0= 2 1561 jp= 1 volp= volt(2) 3061 test0= vol(jp,2)/volp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbu(jp,l,2)-gbl(jp,l,2))*aprx(l)+ # gbl(jp,l,2) enddo else if(jp.lt.nlim) then volp= volp-vol(jp,2) jp= jp+1 go to 3061 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbu(jp,l,2)-gbl(jp,l,2))*aprx(l)+ # gbl(jp,l,2) enddo endif endif pru= g05caf(pru) if(otype.eq.'nc64') then iint= 0 xt0= wtoxsn64(ndim,prx) if(iselsp.eq.1) then iint= 1 xt1= wtoxsn64(ndim,prx) iint= 2 xt2= wtoxsn64(ndim,prx) endif xt= xt0/sigan(2) else xt0= wtoxsn(ndim,prx) xt= xt0/sigan(2) endif prh= oxsmxn(jp,2)/sigan(2) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 if(iselsp.eq.1) then amps1(it)= xt1 amps2(it)= xt2 endif it= it+1 else go to 1561 endif else wxs3= g05caf(dmy) wxs0= sigan(3)/sigt3 if(wxs3.lt.wxs0) then ix0= 2 it0= 1 1562 jp= 1 volp= volt(3) 3062 test0= vol(jp,3)/volp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbu(jp,l,3)-gbl(jp,l,3))*aprx(l)+ # gbl(jp,l,3) enddo else if(jp.lt.nlim) then volp= volp-vol(jp,3) jp= jp+1 go to 3062 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbu(jp,l,3)-gbl(jp,l,3))*aprx(l)+ # gbl(jp,l,3) enddo endif endif pru= g05caf(pru) if(otype.eq.'nc64') then iint= 0 xt0= wtoxsn64(ndim,prx) if(iselsp.eq.1) then iint= 1 xt1= wtoxsn64(ndim,prx) iint= 2 xt2= wtoxsn64(ndim,prx) endif xt= xt0/sigan(3) else xt0= wtoxsn(ndim,prx) xt= xt0/sigan(3) endif prh= oxsmxn(jp,3)/sigan(3) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 if(iselsp.eq.1) then amps1(it)= xt1 amps2(it)= xt2 endif it= it+1 else go to 1562 endif else ix0= 2 it0= 2 1563 jp= 1 volp= volt(4) 3063 test0= vol(jp,4)/volp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbu(jp,l,4)-gbl(jp,l,4))*aprx(l)+ # gbl(jp,l,4) enddo else if(jp.lt.nlim) then volp= volp-vol(jp,4) jp= jp+1 go to 3063 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbu(jp,l,4)-gbl(jp,l,4))*aprx(l)+ # gbl(jp,l,4) enddo endif endif pru= g05caf(pru) if(otype.eq.'nc64') then iint= 0 xt0= wtoxsn64(ndim,prx) if(iselsp.eq.1) then iint= 1 xt1= wtoxsn64(ndim,prx) iint= 2 xt2= wtoxsn64(ndim,prx) endif xt= xt0/sigan(4) else xt0= wtoxsn(ndim,prx) xt= xt0/sigan(4) endif prh= oxsmxn(jp,4)/sigan(4) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 if(iselsp.eq.1) then amps1(it)= xt1 amps2(it)= xt2 endif it= it+1 else go to 1563 endif endif endif endif if(it.le.itmx) then go to 1564 endif * else if(otype.eq.'nc68') then it= 1 sigt= 0.d0 do i=1,ipmx sigt= sigt+sigahb(i) enddo tsig= sigt sigp= sigt 8517 do i0=1,2 do j0=1,2 do k0=1,k0mx ix0= i0 it0= j0 ip0= k0 ip= k0mx*(2*(i0-1)+j0-1)+k0 wxs= g05caf(dmy) wxs0= sigahb(ip)/sigp if(wxs.lt.wxs0) then 8505 jp= 1 volhp= volhtb(ip) 8005 test0= volhb(jp,ip)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhb(jp,l,ip)-gblhb(jp,l,ip))* # aprx(l)+gblhb(jp,l,ip) enddo else if(jp.lt.nlim) then volhp= volhp-volhb(jp,ip) jp= jp+1 go to 8005 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhb(jp,l,ip)-gblhb(jp,l,ip))* # aprx(l)+gblhb(jp,l,ip) enddo endif endif pru= g05caf(dmy) if(omssm.eq.'n') then xt0= wtoxsnh64(ndim,prx) xt= xt0/sigahb(ip) else if(omssm.eq.'y') then xt0= wtoxsnsh64(ndim,prx) xt= xt0/sigahb(ip) endif prh= oxsmxhb(jp,ip)/sigahb(ip) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 8505 endif endif if(ip.le.(ipmx-1)) then sigp= sigp-sigahb(ip) else 8506 jp= 1 volhp= volhtb(ip) 8006 test0= volhb(jp,ip)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhb(jp,l,ip)-gblhb(jp,l,ip))* # aprx(l)+gblhb(jp,l,ip) enddo else if(jp.lt.nlim) then volhp= volhp-volhb(jp,ip) jp= jp+1 go to 8006 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhb(jp,l,ip)-gblhb(jp,l,ip))* # aprx(l)+gblhb(jp,l,ip) enddo endif endif pru= g05caf(dmy) if(omssm.eq.'n') then xt0= wtoxsnh64(ndim,prx) xt= xt0/sigahb(ip) else if(omssm.eq.'y') then xt0= wtoxsnsh64(ndim,prx) xt= xt0/sigahb(ip) endif prh= oxsmxhb(jp,ip)/sigahb(ip) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 8506 endif endif enddo enddo enddo if(it.le.itmx) then go to 8517 endif * else if(otype.eq.'nc33'.or.otype.eq.'nc26') then it= 1 sigt= 0.d0 do i=1,12 sigt= sigt+sigahq(i) enddo tsig= sigt sigt2= sigt-sigahq(1) sigt3= sigt2-sigahq(2) sigt4= sigt3-sigahq(3) sigt5= sigt4-sigahq(4) sigt6= sigt5-sigahq(5) sigt7= sigt6-sigahq(6) sigt8= sigt7-sigahq(7) sigt9= sigt8-sigahq(8) sigt10= sigt9-sigahq(9) sigt11= sigt10-sigahq(10) 2517 wxs= g05caf(dmy) wxs0= sigahq(1)/sigt if(wxs.lt.wxs0) then ix0= 1 it0= 1 ip0= 1 2505 jp= 1 volhp= volhtq(1) 4005 test0= volhq(jp,1)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,1)-gblhq(jp,l,1))* # aprx(l)+gblhq(jp,l,1) enddo else if(jp.lt.nlim) then volhp= volhp-volhq(jp,1) jp= jp+1 go to 4005 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,1)-gblhq(jp,l,1))* # aprx(l)+gblhq(jp,l,1) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc33') then xt0= wtoxsnh32(ndim,prx) xt= xt0/sigahq(1) else if(otype.eq.'nc26') then xt0= wtoxsnh26(ndim,prx) xt= xt0/sigahq(1) endif prh= oxsmxhq(jp,1)/sigahq(1) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 2505 endif else wxs2= g05caf(dmy) wxs0= sigahq(2)/sigt2 if(wxs2.lt.wxs0) then ix0= 1 it0= 1 ip0= 2 2506 jp= 1 volhp= volhtq(2) 4006 test0= volhq(jp,2)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,2)-gblhq(jp,l,2))* # aprx(l)+gblhq(jp,l,2) enddo else if(jp.lt.nlim) then volhp= volhp-volhq(jp,2) jp= jp+1 go to 4006 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,2)-gblhq(jp,l,2))* # aprx(l)+gblhq(jp,l,2) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc33') then xt0= wtoxsnh32(ndim,prx) xt= xt0/sigahq(2) else if(otype.eq.'nc26') then xt0= wtoxsnh26(ndim,prx) xt= xt0/sigahq(2) endif prh= oxsmxhq(jp,2)/sigahq(2) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 2506 endif else wxs3= g05caf(dmy) wxs0= sigahq(3)/sigt3 if(wxs3.lt.wxs0) then ix0= 1 it0= 1 ip0= 3 2507 jp= 1 volhp= volhtq(3) 4007 test0= volhq(jp,3)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,3)-gblhq(jp,l,3))* # aprx(l)+gblhq(jp,l,3) enddo else if(jp.lt.nlim) then volhp= volhp-volhq(jp,3) jp= jp+1 go to 4007 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,3)-gblhq(jp,l,3))* # aprx(l)+gblhq(jp,l,3) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc33') then xt0= wtoxsnh32(ndim,prx) xt= xt0/sigahq(3) else if(otype.eq.'nc26') then xt0= wtoxsnh26(ndim,prx) xt= xt0/sigahq(3) endif prh= oxsmxhq(jp,3)/sigahq(3) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 2507 endif else wxs4= g05caf(dmy) wxs0= sigahq(4)/sigt4 if(wxs4.lt.wxs0) then ix0= 1 it0= 2 ip0= 1 2508 jp= 1 volhp= volhtq(4) 4008 test0= volhq(jp,4)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,4)-gblhq(jp,l,4))* # aprx(l)+gblhq(jp,l,4) enddo else if(jp.lt.nlim) then volhp= volhp-volhq(jp,4) jp= jp+1 go to 4008 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,4)-gblhq(jp,l,4))* # aprx(l)+gblhq(jp,l,4) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc33') then xt0= wtoxsnh32(ndim,prx) xt= xt0/sigahq(4) else if(otype.eq.'nc26') then xt0= wtoxsnh26(ndim,prx) xt= xt0/sigahq(4) endif prh= oxsmxhq(jp,4)/sigahq(4) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 2508 endif else wxs5= g05caf(dmy) wxs0= sigahq(5)/sigt5 if(wxs5.lt.wxs0) then ix0= 1 it0= 2 ip0= 2 2509 jp= 1 volhp= volhtq(5) 4009 test0= volhq(jp,5)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,5)-gblhq(jp,l,5))* # aprx(l)+gblhq(jp,l,5) enddo else if(jp.lt.nlim) then volhp= volhp-volhq(jp,5) jp= jp+1 go to 4009 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,5)-gblhq(jp,l,5))* # aprx(l)+gblhq(jp,l,5) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc33') then xt0= wtoxsnh32(ndim,prx) xt= xt0/sigahq(5) else if(otype.eq.'nc26') then xt0= wtoxsnh26(ndim,prx) xt= xt0/sigahq(5) endif prh= oxsmxhq(jp,5)/sigahq(5) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 2509 endif else wxs6= g05caf(dmy) wxs0= sigahq(6)/sigt6 if(wxs6.lt.wxs0) then ix0= 1 it0= 2 ip0= 3 2510 jp= 1 volhp= volhtq(6) 4010 test0= volhq(jp,6)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,6)-gblhq(jp,l,6))* # aprx(l)+gblhq(jp,l,6) enddo else if(jp.lt.nlim) then volhp= volhp-volhq(jp,6) jp= jp+1 go to 4010 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,6)-gblhq(jp,l,6))* # aprx(l)+gblhq(jp,l,6) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc33') then xt0= wtoxsnh32(ndim,prx) xt= xt0/sigahq(6) else if(otype.eq.'nc26') then xt0= wtoxsnh26(ndim,prx) xt= xt0/sigahq(6) endif prh= oxsmxhq(jp,6)/sigahq(6) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 2510 endif else wxs7= g05caf(dmy) wxs0= sigahq(7)/sigt7 if(wxs7.lt.wxs0) then ix0= 2 it0= 1 ip0= 1 2511 jp= 1 volhp= volhtq(7) 4011 test0= volhq(jp,7)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,7)-gblhq(jp,l,7))* # aprx(l)+gblhq(jp,l,7) enddo else if(jp.lt.nlim) then volhp= volhp-volhq(jp,7) jp= jp+1 go to 4011 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,7)-gblhq(jp,l,7))* # aprx(l)+gblhq(jp,l,7) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc33') then xt0= wtoxsnh32(ndim,prx) xt= xt0/sigahq(7) else if(otype.eq.'nc26') then xt0= wtoxsnh26(ndim,prx) xt= xt0/sigahq(7) endif prh= oxsmxhq(jp,7)/sigahq(7) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 2511 endif else wxs8= g05caf(dmy) wxs0= sigahq(8)/sigt8 if(wxs8.lt.wxs0) then ix0= 2 it0= 1 ip0= 2 2512 jp= 1 volhp= volhtq(8) 4012 test0= volhq(jp,8)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,8)-gblhq(jp,l,8))* # aprx(l)+gblhq(jp,l,8) enddo else if(jp.lt.nlim) then volhp= volhp-volhq(jp,8) jp= jp+1 go to 4012 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,8)-gblhq(jp,l,8))* # aprx(l)+gblhq(jp,l,8) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc33') then xt0= wtoxsnh32(ndim,prx) xt= xt0/sigahq(8) else if(otype.eq.'nc26') then xt0= wtoxsnh26(ndim,prx) xt= xt0/sigahq(8) endif prh= oxsmxhq(jp,8)/sigahq(8) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1h else go to 2512 endif else wxs9= g05caf(dmy) wxs0= sigahq(9)/sigt9 if(wxs9.lt.wxs0) then ix0= 2 it0= 1 ip0= 3 2513 jp= 1 volhp= volhtq(9) 4013 test0= volhq(jp,9)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,9)-gblhq(jp,l,9))* # aprx(l)+gblhq(jp,l,9) enddo else if(jp.lt.nlim) then volhp= volhp-volhq(jp,9) jp= jp+1 go to 4013 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,9)-gblhq(jp,l,9))* # aprx(l)+gblhq(jp,l,9) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc33') then xt0= wtoxsnh32(ndim,prx) xt= xt0/sigahq(9) else if(otype.eq.'nc26') then xt0= wtoxsnh26(ndim,prx) xt= xt0/sigahq(9) endif prh= oxsmxhq(jp,9)/sigahq(9) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 2513 endif else wxs10= g05caf(dmy) wxs0= sigahq(10)/sigt10 if(wxs10.lt.wxs0) then ix0= 2 it0= 2 ip0= 1 2514 jp= 1 volhp= volhtq(10) 4014 test0= volhq(jp,10)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,10)-gblhq(jp,l,10))* # aprx(l)+gblhq(jp,l,10) enddo else if(jp.lt.nlim) then volhp= volhp-volhq(jp,10) jp= jp+1 go to 4014 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,10)-gblhq(jp,l,10))* # aprx(l)+gblhq(jp,l,10) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc33') then xt0= wtoxsnh32(ndim,prx) xt= xt0/sigahq(10) else if(otype.eq.'nc26') then xt0= wtoxsnh26(ndim,prx) xt= xt0/sigahq(10) endif prh= oxsmxhq(jp,10)/sigahq(10) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 2514 endif else wxs11= g05caf(dmy) wxs0= sigahq(11)/sigt11 if(wxs11.lt.wxs0) then ix0= 2 it0= 2 ip0= 2 2515 jp= 1 volhp= volhtq(11) 4015 test0= volhq(jp,11)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,11)-gblhq(jp,l,11))* # aprx(l)+gblhq(jp,l,11) enddo else if(jp.lt.nlim) then volhp= volhp-volhq(jp,11) jp= jp+1 go to 4015 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,11)-gblhq(jp,l,11))* # aprx(l)+gblhq(jp,l,11) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc33') then xt0= wtoxsnh32(ndim,prx) xt= xt0/sigahq(11) else if(otype.eq.'nc26') then xt0= wtoxsnh26(ndim,prx) xt= xt0/sigahq(11) endif prh= oxsmxhq(jp,11)/sigahq(11) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 2515 endif else ix0= 2 it0= 2 ip0= 3 2516 jp= 1 volhp= volhtq(12) 4016 test0= volhq(jp,12)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,12)-gblhq(jp,l,12))* # aprx(l)+gblhq(jp,l,12) enddo else if(jp.lt.nlim) then volhp= volhp-volhq(jp,12) jp= jp+1 go to 4016 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuhq(jp,l,12)-gblhq(jp,l,12))* # aprx(l)+gblhq(jp,l,12) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc33') then xt0= wtoxsnh32(ndim,prx) xt= xt0/sigahq(12) else if(otype.eq.'nc26') then xt0= wtoxsnh26(ndim,prx) xt= xt0/sigahq(12) endif prh= oxsmxhq(jp,12)/sigahq(12) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 2516 endif endif endif endif endif endif endif endif endif endif endif endif if(it.le.itmx) then go to 2517 endif * else * if((otype.eq.'mx43'.and.ofs.eq.'ll').or. # (otype.eq.'mx43'.and.ofs.eq.'qq'.and.opglu.eq.'n').or. # (otype.eq.'nc32'.and.opglu.eq.'n')) then * if(otype.eq.'mx43'.and.ofs.eq.'qq'.and.opglu.eq.'n'.and. # oint.eq.'y') then iselsp= 1 else iselsp= 0 endif * it= 1 sigt= 0.d0 do i=1,4 sigt= sigt+sigah(i) enddo tsig= sigt sigt2= sigt-sigah(1) sigt3= sigt2-sigah(2) 9513 wxs= g05caf(dmy) wxs0= sigah(1)/sigt if(wxs.lt.wxs0) then ix0= 1 it0= 1 ip0= 1 9505 jp= 1 volhp= volht(1) 9005 test0= volh(jp,1)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,1)-gblh(jp,l,1))* # aprx(l)+gblh(jp,l,1) enddo else if(jp.lt.nlim) then volhp= volhp-volh(jp,1) jp= jp+1 go to 9005 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,1)-gblh(jp,l,1))* # aprx(l)+gblh(jp,l,1) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc32') then xt0= wtoxsn32(ndim,prx) xt= xt0/sigah(1) else if(otype.eq.'mx43') then iint= 0 xt0= wtoxsm43(ndim,prx) xt= xt0/sigah(1) endif prh= oxsmxh(jp,1)/sigah(1) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 if(iselsp.eq.1) then iint= 1 xt1= wtoxsm43(ndim,prx) iint= 2 xt2= wtoxsm43(ndim,prx) amps1(it)= xt1 amps2(it)= xt2 endif it= it+1 else go to 9505 endif else wxs2= g05caf(dmy) wxs0= sigah(2)/sigt2 if(wxs2.lt.wxs0) then ix0= 1 it0= 2 ip0= 1 9506 jp= 1 volhp= volht(2) 9006 test0= volh(jp,2)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,2)-gblh(jp,l,2))* # aprx(l)+gblh(jp,l,2) enddo else if(jp.lt.nlim) then volhp= volhp-volh(jp,2) jp= jp+1 go to 9006 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,2)-gblh(jp,l,2))* # aprx(l)+gblh(jp,l,2) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc32') then xt0= wtoxsn32(ndim,prx) xt= xt0/sigah(2) else if(otype.eq.'mx43') then xt0= wtoxsm43(ndim,prx) xt= xt0/sigah(2) endif prh= oxsmxh(jp,2)/sigah(2) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 if(iselsp.eq.1) then iint= 1 xt1= wtoxsm43(ndim,prx) iint= 2 xt2= wtoxsm43(ndim,prx) amps1(it)= xt1 amps2(it)= xt2 endif it= it+1 else go to 9506 endif else wxs3= g05caf(dmy) wxs0= sigah(3)/sigt3 if(wxs3.lt.wxs0) then ix0= 2 it0= 1 ip0= 1 9507 jp= 1 volhp= volht(3) 9007 test0= volh(jp,3)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,3)-gblh(jp,l,3))* # aprx(l)+gblh(jp,l,3) enddo else if(jp.lt.nlim) then volhp= volhp-volh(jp,3) jp= jp+1 go to 9007 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,3)-gblh(jp,l,3))* # aprx(l)+gblh(jp,l,3) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc32') then xt0= wtoxsn32(ndim,prx) xt= xt0/sigah(3) else if(otype.eq.'mx43') then iint= 0 xt0= wtoxsm43(ndim,prx) xt= xt0/sigah(3) endif prh= oxsmxh(jp,3)/sigah(3) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 if(iselsp.eq.1) then iint= 1 xt1= wtoxsm43(ndim,prx) iint= 2 xt2= wtoxsm43(ndim,prx) amps1(it)= xt1 amps2(it)= xt2 endif it= it+1 else go to 9507 endif else ix0= 2 it0= 2 ip0= 1 9512 jp= 1 volhp= volht(4) 9012 test0= volh(jp,4)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,4)-gblh(jp,l,4))* # aprx(l)+gblh(jp,l,4) enddo else if(jp.lt.nlim) then volhp= volhp-volh(jp,4) jp= jp+1 go to 9012 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,4)-gblh(jp,l,4))* # aprx(l)+gblh(jp,l,4) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc32') then xt0= wtoxsn32(ndim,prx) xt= xt0/sigah(4) else if(otype.eq.'mx43') then iint= 0 xt0= wtoxsm43(ndim,prx) xt= xt0/sigah(4) endif prh= oxsmxh(jp,4)/sigah(4) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 if(iselsp.eq.1) then iint= 1 xt1= wtoxsm43(ndim,prx) iint= 2 xt2= wtoxsm43(ndim,prx) amps1(it)= xt1 amps2(it)= xt2 endif it= it+1 else go to 9512 endif endif endif endif if(it.le.itmx) then go to 9513 endif * else * it= 1 sigt= 0.d0 do i=1,8 sigt= sigt+sigah(i) enddo tsig= sigt sigt2= sigt-sigah(1) sigt3= sigt2-sigah(5) sigt4= sigt3-sigah(3) sigt5= sigt4-sigah(2) sigt6= sigt5-sigah(7) sigt7= sigt6-sigah(6) call g05cbf(0) 1513 wxs= g05caf(dmy) wxs0= sigah(1)/sigt if(wxs.lt.wxs0) then ix0= 1 it0= 1 ip0= 1 1505 jp= 1 volhp= volht(1) 3005 test0= volh(jp,1)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,1)-gblh(jp,l,1))* # aprx(l)+gblh(jp,l,1) enddo else if(jp.lt.nlim) then volhp= volhp-volh(jp,1) jp= jp+1 go to 3005 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,1)-gblh(jp,l,1))* # aprx(l)+gblh(jp,l,1) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc50') then xt0= wtoxsnh49(ndim,prx) else if(otype.eq.'nc25') then xt0= wtoxsnh24(ndim,prx) else if(otype.eq.'nc21') then xt0= wtoxsnh19(ndim,prx) else if(otype.eq.'nc32') then xt0= wtoxsn32(ndim,prx) else if(otype.eq.'cc12') then xt0= wtoxssc(ndim,prx) else if(otype.eq.'mx43') then iint= 0 xt0= wtoxsm43(ndim,prx) endif xt= xt0/sigah(1) prh= oxsmxh(jp,1)/sigah(1) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 1505 endif else wxs2= g05caf(dmy) wxs0= sigah(5)/sigt2 if(wxs2.lt.wxs0) then ix0= 2 it0= 1 ip0= 1 1506 jp= 1 volhp= volht(5) 3006 test0= volh(jp,5)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,5)-gblh(jp,l,5))* # aprx(l)+gblh(jp,l,5) enddo else if(jp.lt.nlim) then volhp= volhp-volh(jp,5) jp= jp+1 go to 3006 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,5)-gblh(jp,l,5))* # aprx(l)+gblh(jp,l,5) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc50') then xt0= wtoxsnh49(ndim,prx) else if(otype.eq.'nc25') then xt0= wtoxsnh24(ndim,prx) else if(otype.eq.'nc21') then xt0= wtoxsnh19(ndim,prx) else if(otype.eq.'nc32') then xt0= wtoxsn32(ndim,prx) else if(otype.eq.'cc12') then xt0= wtoxssc(ndim,prx) else if(otype.eq.'mx43') then iint= 0 xt0= wtoxsm43(ndim,prx) endif xt= xt0/sigah(5) prh= oxsmxh(jp,5)/sigah(5) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 1506 endif else wxs3= g05caf(dmy) wxs0= sigah(3)/sigt3 if(wxs3.lt.wxs0) then ix0= 1 it0= 2 ip0= 1 1507 jp= 1 volhp= volht(3) 3007 test0= volh(jp,3)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,3)-gblh(jp,l,3))* # aprx(l)+gblh(jp,l,3) enddo else if(jp.lt.nlim) then volhp= volhp-volh(jp,3) jp= jp+1 go to 3007 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,3)-gblh(jp,l,3))* # aprx(l)+gblh(jp,l,3) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc50') then xt0= wtoxsnh49(ndim,prx) else if(otype.eq.'nc25') then xt0= wtoxsnh24(ndim,prx) else if(otype.eq.'nc21') then xt0= wtoxsnh19(ndim,prx) else if(otype.eq.'nc32') then xt0= wtoxsn32(ndim,prx) else if(otype.eq.'cc12') then xt0= wtoxssc(ndim,prx) else if(otype.eq.'mx43') then iint= 0 xt0= wtoxsm43(ndim,prx) endif xt= xt0/sigah(3) prh= oxsmxh(jp,3)/sigah(3) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 1507 endif else wxs4= g05caf(dmy) wxs0= sigah(2)/sigt4 if(wxs4.lt.wxs0) then ix0= 1 it0= 1 ip0= 2 1508 jp= 1 volhp= volht(2) 3008 test0= volh(jp,2)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,2)-gblh(jp,l,2))* # aprx(l)+gblh(jp,l,2) enddo else if(jp.lt.nlim) then volhp= volhp-volh(jp,2) jp= jp+1 go to 3008 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,2)-gblh(jp,l,2))* # aprx(l)+gblh(jp,l,2) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc50') then xt0= wtoxsnh49(ndim,prx) else if(otype.eq.'nc25') then xt0= wtoxsnh24(ndim,prx) else if(otype.eq.'nc21') then xt0= wtoxsnh19(ndim,prx) else if(otype.eq.'nc32') then xt0= wtoxsn32(ndim,prx) else if(otype.eq.'cc12') then xt0= wtoxssc(ndim,prx) else if(otype.eq.'mx43') then iint= 0 xt0= wtoxsm43(ndim,prx) endif xt= xt0/sigah(2) prh= oxsmxh(jp,2)/sigah(2) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 1508 endif else wxs5= g05caf(dmy) wxs0= sigah(7)/sigt5 if(wxs5.lt.wxs0) then ix0= 2 it0= 2 ip0= 1 1509 jp= 1 volhp= volht(7) 3009 test0= volh(jp,7)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,7)-gblh(jp,l,7))* # aprx(l)+gblh(jp,l,7) enddo else if(jp.lt.nlim) then volhp= volhp-volh(jp,7) jp= jp+1 go to 3009 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,7)-gblh(jp,l,7))* # aprx(l)+gblh(jp,l,7) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc50') then xt0= wtoxsnh49(ndim,prx) else if(otype.eq.'nc25') then xt0= wtoxsnh24(ndim,prx) else if(otype.eq.'nc21') then xt0= wtoxsnh19(ndim,prx) else if(otype.eq.'nc32') then xt0= wtoxsn32(ndim,prx) else if(otype.eq.'cc12') then xt0= wtoxssc(ndim,prx) else if(otype.eq.'mx43') then iint= 0 xt0= wtoxsm43(ndim,prx) endif xt= xt0/sigah(7) prh= oxsmxh(jp,7)/sigah(7) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 1509 endif else wxs6= g05caf(dmy) wxs0= sigah(6)/sigt6 if(wxs6.lt.wxs0) then ix0= 2 it0= 1 ip0= 2 1510 jp= 1 volhp= volht(6) 3010 test0= volh(jp,6)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,6)-gblh(jp,l,6))* # aprx(l)+gblh(jp,l,6) enddo else if(jp.lt.nlim) then volhp= volhp-volh(jp,6) jp= jp+1 go to 3010 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,6)-gblh(jp,l,6))* # aprx(l)+gblh(jp,l,6) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc50') then xt0= wtoxsnh49(ndim,prx) else if(otype.eq.'nc25') then xt0= wtoxsnh24(ndim,prx) else if(otype.eq.'nc21') then xt0= wtoxsnh19(ndim,prx) else if(otype.eq.'nc32') then xt0= wtoxsn32(ndim,prx) else if(otype.eq.'cc12') then xt0= wtoxssc(ndim,prx) else if(otype.eq.'mx43') then iint= 0 xt0= wtoxsm43(ndim,prx) endif xt= xt0/sigah(6) prh= oxsmxh(jp,6)/sigah(6) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 1510 endif else wxs7= g05caf(dmy) wxs0= sigah(4)/sigt7 if(wxs7.lt.wxs0) then ix0= 1 it0= 2 ip0= 2 1511 jp= 1 volhp= volht(4) 3011 test0= volh(jp,4)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,4)-gblh(jp,l,4))* # aprx(l)+gblh(jp,l,4) enddo else if(jp.lt.nlim) then volhp= volhp-volh(jp,4) jp= jp+1 go to 3011 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,4)-gblh(jp,l,4))* # aprx(l)+gblh(jp,l,4) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc50') then xt0= wtoxsnh49(ndim,prx) else if(otype.eq.'nc25') then xt0= wtoxsnh24(ndim,prx) else if(otype.eq.'nc21') then xt0= wtoxsnh19(ndim,prx) else if(otype.eq.'nc32') then xt0= wtoxsn32(ndim,prx) else if(otype.eq.'cc12') then xt0= wtoxssc(ndim,prx) else if(otype.eq.'mx43') then iint= 0 xt0= wtoxsm43(ndim,prx) endif xt= xt0/sigah(4) prh= oxsmxh(jp,4)/sigah(4) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 1511 endif else ix0= 2 it0= 2 ip0= 2 1512 jp= 1 volhp= volht(8) 3012 test0= volh(jp,8)/volhp pru0= g05caf(dmy) if(pru0.lt.test0) then do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,8)-gblh(jp,l,8))* # aprx(l)+gblh(jp,l,8) enddo else if(jp.lt.nlim) then volhp= volhp-volh(jp,8) jp= jp+1 go to 3012 else jp= jp+1 do l=1,9 aprx(l)= g05caf(dmy) prx(l)= (gbuh(jp,l,8)-gblh(jp,l,8))* # aprx(l)+gblh(jp,l,8) enddo endif endif pru= g05caf(dmy) if(otype.eq.'nc50') then xt0= wtoxsnh49(ndim,prx) else if(otype.eq.'nc25') then xt0= wtoxsnh24(ndim,prx) else if(otype.eq.'nc21') then xt0= wtoxsnh19(ndim,prx) else if(otype.eq.'nc32') then xt0= wtoxsn32(ndim,prx) else if(otype.eq.'cc12') then xt0= wtoxssc(ndim,prx) else if(otype.eq.'mx43') then iint= 0 xt0= wtoxsm43(ndim,prx) endif xt= xt0/sigah(8) prh= oxsmxh(jp,8)/sigah(8) prtst= xt-pru*prh if(prtst.ge.0.d0.and.xt.ne.0.d0) then do i=1,ninv prxa(i,it)= xaph(i) enddo amps(it)= xt0 it= it+1 else go to 1512 endif endif endif endif endif endif endif endif if(it.le.itmx) then go to 1513 endif * *-----CC or Higgs/NC generation ends * endif endif endif * *-----CC or Higgs storage starts * if(ostore.eq.'d') then if(oxcm.eq.'c') then xint= 10.d0 xintl= 25.d0 xintwp= 2.d0 xinte= 25.d0 xint4j= 40.d0 do i=1,im exi= 75.d0+(i-1.d0)*xint/im exf= 75.d0+i*xint/im exli= -15.d0+(i-1.d0)*xintl/im exlf= -15.d0+i*xintl/im exwpi= -1.d0+(i-1.d0)*xintwp/im exwpf= -1.d0+i*xintwp/im exei= (i-1.d0)*xinte/im exef= i*xinte/im ex4ji= 90.d0+(i-1.d0)*xint4j/im ex4jf= 75.d0+i*xint4j/im exc(i)= 0.d0 excl(i)= 0.d0 excwp(i)= 0.d0 exce(i)= 0.d0 exc4j(i)= 0.d0 do j=1,itmx cvv= prxa(1,j)*prxa(2,j) extst(j)= cvv*prxa(4,j) if(extst(j).gt.0.d0) then extst(j)= sqrt(extst(j))*rs if(exf.gt.extst(j).and.extst(j).gt.exi) then exc(i)= exc(i)+1.d0 endif he1= prxa(3,j)+prxa(6,j)+prxa(7,j) ht1= prxa(9,j) if((prxa(2,j)*(he1-ht1)).gt.0.d0) then qtrf(j)= log(prxa(2,j)*(he1-ht1)) if(exlf.gt.qtrf(j).and.qtrf(j).gt.exli) then excl(i)= excl(i)+1.d0 endif endif endif hewp= 0.5d0*(prxa(2,j)*(1.d0+prxa(4,j)- # prxa(3,j))-(prxa(2,j)-prxa(1,j))* # (1.d0-prxa(8,j))) hbwps= 1.d0-cvv*prxa(4,j)/hewp/hewp if(hbwps.gt.0.d0) then hbwp= sqrt(hbwps) cwpst(j)= (1.d0-prxa(1,j)*(1.d0-prxa(8,j))/ # hewp)/hbwp if(exwpf.gt.cwpst(j).and.cwpst(j).gt. # exwpi) then excwp(i)= excwp(i)+1.d0 endif endif exetst(j)= prxa(2,j)*(prxa(3,j)+prxa(6,j)+ # prxa(7,j))-(prxa(2,j)-prxa(1,j))* # prxa(9,j)*rs if(exef.gt.exetst(j).and.exetst(j).gt.exei) then exce(i)= exce(i)+1.d0 endif xiv1= sqrt(cvv*prxa(3,j))*rs xiv2= sqrt(cvv*prxa(4,j))*rs xiv3= sqrt(cvv*prxa(5,j))*rs xiv4= sqrt(cvv*prxa(6,j))*rs xiv5= sqrt(cvv*prxa(7,j))*rs xiv6= sqrt(cvv)*sqrt(1.d0-prxa(3,j)- # prxa(4,j)-prxa(5,j)-prxa(6,j)- # prxa(7,j))*rs div1= abs(xiv1-xiv2) div2= abs(xiv3-xiv4) div3= abs(xiv5-xiv6) divall= dmin1(div1,div2,div3) if(divall.eq.div1) then sum4j= xiv1+xiv2 else if(divall.eq.div2) then sum4j= xiv3+xiv4 else if(divall.eq.div3) then sum4j= xiv5+xiv6 endif if(ex4ji.lt.sum4j.and.sum4j.lt.ex4jf) then exc4j(i)= exc4j(i)+1 endif enddo enddo * ern= xint/im open(nout,file='mco',status='new') do i=1,im exb= 75.d0+(i-1.d0)*xint/im ecx= 75.d0+(i-0.5d0)*xint/im write(nout,fmt=2500) exb,ecx,exc(i)/itmx/ern*tsig enddo close(nout) ern4j= xint4j/im open(nout,file='mco4j',status='new') do i=1,im exb4j= 90.d0+(i-1.d0)*xint4j/im ecx4j= 90.d0+(i-0.5d0)*xint4j/im write(nout,fmt=2500) exb4j,ecx4j,exc4j(i)/itmx/ # ern4j*tsig enddo close(nout) ernl= xintl/im open(nout,file='mcol',status='new') do i=1,im exb= -15.d0+(i-1.d0)*xintl/im ecx= -15.d0+(i-0.5d0)*xintl/im write(nout,fmt=2500) exb,ecx,excl(i)/itmx/ernl*tsig enddo close(nout) ernwp= xintwp/im open(nout,file='mcowp',status='new') do i=1,im exb= -1.d0+(i-1.d0)*xintwp/im ecx= -1.d0+(i-0.5d0)*xintwp/im write(nout,fmt=2500) exb,ecx, # excwp(i)/itmx/ernwp*tsig enddo close(nout) erne= xinte/im open(nout,file='mcoe',status='new') do i=1,im exeb= (i-1.d0)*xinte/im ecex= (i-0.5d0)*xinte/im write(nout,fmt=2500) exeb,ecex,exce(i)/itmx/erne*tsig enddo close(nout) else if(oxcm.eq.'n') then xint= 60.d0 xintv= 100.d0 xintqt= 120.d0 do i=1,im exi= 60.d0+(i-1.d0)*xint/im exf= 60.d0+i*xint/im exiv= 30.d0+(i-1.d0)*xintv/im exfv= 30.d0+i*xintv/im if(otype.eq.'nc68') then excm(i)= 0.d0 do j=1,itmx bthr= 50.d0 hvv= prxa(1,j)*prxa(2,j) extst1(j)= sqrt(hvv*prxa(3,j))*rs extst2(j)= sqrt(hvv*prxa(4,j))*rs extst3(j)= sqrt(hvv*prxa(5,j))*rs extst4(j)= sqrt(hvv*prxa(6,j))*rs extst5(j)= sqrt(hvv*prxa(7,j))*rs if(oev.eq.'a1') then secr= 1.d0-prxa(3,j)-prxa(4,j)-prxa(5,j)- # prxa(6,j)-prxa(7,j) if(secr.gt.0.d0) then extst6(j)= sqrt(hvv*secr)*rs if(exf.gt.extst1(j).and.extst1(j).gt.exi. # and.extst2(j).gt.bthr) then excm(i)= excm(i)+1.d0 else if(exf.gt.extst2(j).and. # extst2(j).gt.exi.and. # extst1(j).gt.bthr) then excm(i)= excm(i)+1.d0 else if(exf.gt.extst3(j).and. # extst3(j).gt.exi.and. # extst4(j).gt.bthr) then excm(i)= excm(i)+1.d0 else if(exf.gt.extst4(j).and. # extst4(j).gt.exi.and. # extst3(j).gt.bthr) then excm(i)= excm(i)+1.d0 else if(exf.gt.extst5(j).and. # extst5(j).gt.exi.and. # extst6(j).gt.bthr) then excm(i)= excm(i)+1.d0 else if(exf.gt.extst6(j).and. # extst6(j).gt.exi.and. # extst5(j).gt.bthr) then excm(i)= excm(i)+1.d0 endif endif else if(oev.eq.'a2') then rim1= extst1(j) rim2= extst2(j) rim3= extst3(j) rim4= extst4(j) rim5= extst5(j) secr= 1.d0-prxa(3,j)-prxa(4,j)-prxa(5,j)- # prxa(6,j)-prxa(7,j) if(secr.gt.0.d0) then extst6(j)= sqrt(hvv*secr)*rs rim6= extst6(j) iev= ievs(exi,exf,rim1,rim2,rim3,rim4,rim5,rim6) if(iev.eq.1) then excm(i)= excm(i)+1.d0 endif endif endif enddo endif exc1(i)= 0.d0 exc2(i)= 0.d0 do j=1,itmx hvv= prxa(1,j)*prxa(2,j) extst1(j)= sqrt(hvv*prxa(4,j))*rs extst2(j)= sqrt(hvv*prxa(3,j))*rs if(exf.gt.extst1(j).and.extst1(j).gt.exi) then exc1(i)= exc1(i)+1.d0 endif if(exf.gt.extst2(j).and.extst2(j).gt.exi) then exc2(i)= exc2(i)+1.d0 endif enddo enddo ern= xint/im ern= tsig/itmx/ern if(otype.eq.'nc68') then open(nout,file='mcombb',status='new') do i=1,im exb= 60.d0+(i-1.d0)*xint/im ecx= 60.d0+(i-0.5d0)*xint/im write(nout,fmt=2500) exb,ecx,excm(i)*ern enddo close(nout) else open(nout,file='mcombb',status='new') do i=1,im exb= 60.d0+(i-1.d0)*xint/im ecx= 60.d0+(i-0.5d0)*xint/im write(nout,fmt=2500) exb,ecx,exc1(i)*ern, # exc2(i)*ern enddo close(nout) endif * if(ipr.eq.31.or.ipr.eq.34.or.ipr.eq.35.or. # ipr.eq.10.or.ipr.eq.13.or.ipr.eq.14) then icl= 1 else if(ipr.eq.30.or.ipr.eq.32.or.ipr.eq.33.or. # ipr.eq.7.or.ipr.eq.19.or.ipr.eq.25) then icl= 2 endif do i=1,im exi= -1.d0+(i-1.d0)*2.d0/im exf= -1.d0+i*2.d0/im exim= 60.d0+(i-1.d0)*xint/im exfm= 60.d0+i*xint/im exiv= 50.d0+(i-1.d0)*xintv/im exfv= 50.d0+i*xintv/im exiqt= (i-1.d0)*xintqt/im exfqt= i*xintqt/im exc1(i)= 0.d0 exc2(i)= 0.d0 exc3(i)= 0.d0 exc4(i)= 0.d0 exc5(i)= 0.d0 exc6(i)= 0.d0 exc7(i)= 0.d0 do j=1,itmx hvv= prxa(1,j)*prxa(2,j) he1= prxa(3,j)+prxa(6,j)+prxa(7,j) he2= 1.d0-prxa(4,j)-prxa(6,j)-prxa(7,j) he3= prxa(4,j)+prxa(5,j)+prxa(7,j) he4= 1.d0-prxa(3,j)-prxa(5,j)-prxa(7,j) ht1= prxa(9,j) ht2= prxa(8,j)-prxa(9,j) ht3= prxa(10,j) ht4= 1.d0-prxa(8,j)-prxa(10,j) hxp= prxa(2,j) hxm= prxa(1,j) hct1= 1.d0-2.d0*hxm*ht1/(hxp*he1-(hxp-hxm)*ht1) hct2= 1.d0-2.d0*hxm*ht2/(hxp*he2-(hxp-hxm)*ht2) hct3= 1.d0-2.d0*hxm*ht3/(hxp*he3-(hxp-hxm)*ht3) hct4= 1.d0-2.d0*hxm*ht4/(hxp*he4-(hxp-hxm)*ht4) hct12= 1.d0-0.5d0*hvv*prxa(3,j)/he1/he2 hct34= 1.d0-0.5d0*hvv*prxa(4,j)/he3/he4 he12= 1.d0+prxa(3,j)-prxa(4,j) he34= 1.d0+prxa(4,j)-prxa(3,j) hem= 0.5d0*(hxp*he12-(hxp-hxm)*prxa(8,j)) hbem= 1.d0-hvv*prxa(3,j)/hem/hem if(hbem.le.0.d0) then isrm= 0 else isrm= 1 hbem= sqrt(hbem) hctsumm= (1.d0-hxm*prxa(8,j)/hem)/hbem endif hep= 0.5d0*(hxp*he34-(hxp-hxm)*(1.d0-prxa(8,j))) hbep= 1.d0-hvv*prxa(4,j)/hep/hep if(hbep.le.0.d0) then isrp= 0 else isrp= 1 hbep= sqrt(hbep) hctsump= (1.d0-hxm*(1.d0-prxa(8,j))/hep)/hbep endif rmmm= s*(1.d0-2.d0*hem+hvv*prxa(3,j)) rmmp= s*(1.d0-2.d0*hep+hvv*prxa(4,j)) rmmm= sqrt(rmmm) rmmp= sqrt(rmmp) if(icl.eq.2) then tq2= -hvv*prxa(3,j)+hem*hem beqs= 1.d0-hvv*prxa(3,j)/hem/hem if(beqs.lt.0.d0.or.tq2.lt.0.d0) then iqt= 0 endif beq= sqrt(beqs) tq= sqrt(tq2) ctq= (1.d0-hxm*prxa(8,j)/hem)/beq ctqf= 1.d0-ctq*ctq if(ctqf.lt.0.d0) then iqt= 0 endif tqt= tq*sqrt(ctqf)*rs iqt= 1 endif if(icl.eq.1) then if((exf.gt.hct3.and.hct3.gt.exi).or. # (exf.gt.hct4.and.hct4.gt.exi)) then exc1(i)= exc1(i)+1.d0 endif if(exf.gt.hct34.and.hct34.gt.exi) then exc2(i)= exc2(i)+1.d0 endif if((isrp.eq.1).and.(exf.gt.hctsump.and. # hctsump.gt.exi)) then exc3(i)= exc3(i)+1.d0 endif if(exfm.gt.rmmm.and.rmmm.gt.exim) then exc4(i)= exc4(i)+1.d0 endif if((isrm.eq.1).and.(exfm.gt.rmmm.and. # rmmm.gt.exim).and.(0.8d0.gt.hctsumm.and. # hctsumm.gt.-0.8d0)) then exc5(i)= exc5(i)+1.d0 endif if(exfv.gt.(hep*rs).and.(hep*rs).gt.exiv) then exc6(i)= exc6(i)+1.d0 endif else if(icl.eq.2) then if((exf.gt.hct1.and.hct1.gt.exi).or. # (exf.gt.hct2.and.hct2.gt.exi)) then exc1(i)= exc1(i)+1.d0 endif if(exf.gt.hct12.and.hct12.gt.exi) then exc2(i)= exc2(i)+1.d0 endif if((isrm.eq.1).and.(exf.gt.hctsumm.and. # hctsumm.gt.exi)) then exc3(i)= exc3(i)+1.d0 endif if(exfm.gt.rmmp.and.rmmp.gt.exim) then exc4(i)= exc4(i)+1.d0 endif if((isrp.eq.1).and.(exfm.gt.rmmp.and. # rmmp.gt.exim).and.(0.8d0.gt.hctsump.and. # hctsump.gt.-0.8d0)) then exc5(i)= exc5(i)+1.d0 endif if(exfv.gt.(hem*rs).and.(hem*rs).gt.exiv) then exc6(i)= exc6(i)+1.d0 endif if((iqt.eq.1).and.(exfv.gt.tqt.and. # tqt.gt.exiv)) then exc7(i)= exc7(i)+1.d0 endif endif enddo enddo * if(otype.ne.'nc68') then ern= 2.d0/im ern= tsig/itmx/ern open(nout,file='mcocbb',status='new') do i=1,im exb= -1.d0+(i-1.d0)*2.d0/im ecx= -1.d0+(i-0.5d0)*2.d0/im write(nout,fmt=2501) exb,ecx,exc1(i)*ern, # exc2(i)*ern,exc3(i)*ern enddo close(nout) open(nout,file='mcomm',status='new') ern= xint/im ern= tsig/itmx/ern do i=1,im exb= 60.d0+(i-1.d0)*xint/im ecx= 60.d0+(i-0.5d0)*xint/im write(nout,fmt=2500) exb,ecx,exc4(i)*ern, # exc5(i)*ern enddo close(nout) ern= xintv/im ern= tsig/itmx/ern open(nout,file='mcoev',status='new') do i=1,im exb= 30.d0+(i-1.d0)*xintv/im ecx= 30.d0+(i-0.5d0)*xintv/im write(nout,fmt=2500) exb,ecx,exc6(i)*ern enddo close(nout) if(icl.eq.2) then ern= xintqt/im ern= tsig/itmx/ern open(nout,file='mcoqt',status='new') do i=1,im exb= (i-1.d0)*xintqt/im ecx= (i-0.5d0)*xintqt/im write(nout,fmt=2500) exb,ecx,exc7(i)*ern enddo close(nout) endif endif * endif print 2503,itmx else if(ostore.eq.'m') then open(nout,file='mom',status='new') if(iosf.eq.0) then nhep= 4 else nhep= 6 endif * if(iosf.gt.0) then idhep(5)= 22 idhep(6)= 22 endif if(ipr.eq.1) then idhep(1)= 16 idhep(2)= -15 idhep(3)= 13 idhep(4)= -14 else if(ipr.eq.2) then if(ockm.eq.'n') then idhep(1)= 2 idhep(2)= -1 idhep(3)= 13 idhep(4)= -14 else if(ickm.eq.1) then idhep(1)= 2 idhep(2)= -1 idhep(3)= 13 idhep(4)= -14 else if(ickm.eq.2) then idhep(1)= 2 idhep(2)= -3 idhep(3)= 13 idhep(4)= -14 else if(ickm.eq.3) then idhep(1)= 2 idhep(2)= -5 idhep(3)= 13 idhep(4)= -14 else if(ickm.eq.4) then idhep(1)= 4 idhep(2)= -1 idhep(3)= 13 idhep(4)= -14 else if(ickm.eq.5) then idhep(1)= 4 idhep(2)= -3 idhep(3)= 13 idhep(4)= -14 else if(ickm.eq.6) then idhep(1)= 4 idhep(2)= -5 idhep(3)= 13 idhep(4)= -14 endif endif else if(ipr.eq.3) then if(ockm.eq.'n') then idhep(1)= 4 idhep(2)= -3 idhep(3)= 1 idhep(4)= -2 else if(ickm.eq.1) then idhep(1)= 4 idhep(2)= -3 idhep(3)= 1 idhep(4)= -2 else if(ickm.eq.2) then idhep(1)= 4 idhep(2)= -1 idhep(3)= 1 idhep(4)= -2 else if(ickm.eq.3) then idhep(1)= 4 idhep(2)= -5 idhep(3)= 1 idhep(4)= -2 else if(ickm.eq.4) then idhep(1)= 4 idhep(2)= -3 idhep(3)= 3 idhep(4)= -2 else if(ickm.eq.5) then idhep(1)= 4 idhep(2)= -1 idhep(3)= 3 idhep(4)= -2 else if(ickm.eq.6) then idhep(1)= 4 idhep(2)= -5 idhep(3)= 3 idhep(4)= -2 else if(ickm.eq.7) then idhep(1)= 4 idhep(2)= -3 idhep(3)= 5 idhep(4)= -2 else if(ickm.eq.8) then idhep(1)= 4 idhep(2)= -1 idhep(3)= 5 idhep(4)= -2 else if(ickm.eq.9) then idhep(1)= 4 idhep(2)= -5 idhep(3)= 5 idhep(4)= -2 else if(ickm.eq.10) then idhep(1)= 2 idhep(2)= -3 idhep(3)= 1 idhep(4)= -2 else if(ickm.eq.11) then idhep(1)= 2 idhep(2)= -5 idhep(3)= 1 idhep(4)= -2 else if(ickm.eq.12) then idhep(1)= 2 idhep(2)= -5 idhep(3)= 3 idhep(4)= -2 else if(ickm.eq.13) then idhep(1)= 4 idhep(2)= -3 idhep(3)= 1 idhep(4)= -4 else if(ickm.eq.14) then idhep(1)= 4 idhep(2)= -5 idhep(3)= 1 idhep(4)= -4 else if(ickm.eq.15) then idhep(1)= 4 idhep(2)= -5 idhep(3)= 3 idhep(4)= -4 endif endif else if(ipr.eq.4) then idhep(1)= 14 idhep(2)= -13 idhep(3)= 11 idhep(4)= -12 else if(ipr.eq.5) then if(ockm.eq.'n') then idhep(1)= 2 idhep(2)= -1 idhep(3)= 11 idhep(4)= -12 else if(ickm.eq.1) then idhep(1)= 2 idhep(2)= -1 idhep(3)= 11 idhep(4)= -12 else if(ickm.eq.2) then idhep(1)= 2 idhep(2)= -3 idhep(3)= 11 idhep(4)= -12 else if(ickm.eq.3) then idhep(1)= 2 idhep(2)= -5 idhep(3)= 11 idhep(4)= -12 else if(ickm.eq.4) then idhep(1)= 4 idhep(2)= -1 idhep(3)= 11 idhep(4)= -12 else if(ickm.eq.5) then idhep(1)= 4 idhep(2)= -3 idhep(3)= 11 idhep(4)= -12 else if(ickm.eq.6) then idhep(1)= 4 idhep(2)= -5 idhep(3)= 11 idhep(4)= -12 endif endif else if(ipr.eq.6) then idhep(1)= 13 idhep(2)= -13 idhep(3)= 16 idhep(4)= -16 else if(ipr.eq.7) then idhep(1)= 1 idhep(2)= -1 idhep(3)= 14 idhep(4)= -14 else if(ipr.eq.8) then idhep(1)= 2 idhep(2)= -2 idhep(3)= 14 idhep(4)= -14 else if(ipr.eq.9) then idhep(1)= 13 idhep(2)= -13 idhep(3)= 15 idhep(4)= -15 else if(ipr.eq.10) then idhep(1)= 13 idhep(2)= -13 idhep(3)= 1 idhep(4)= -1 else if(ipr.eq.11) then idhep(1)= 13 idhep(2)= -13 idhep(3)= 2 idhep(4)= -2 else if(ipr.eq.12) then idhep(1)= 14 idhep(2)= -14 idhep(3)= 16 idhep(4)= -16 else if(ipr.eq.13) then idhep(1)= 3 idhep(2)= -3 idhep(3)= 2 idhep(4)= -2 else if(ipr.eq.14) then idhep(1)= 1 idhep(2)= -1 idhep(3)= 3 idhep(4)= -3 else if(ipr.eq.15) then idhep(1)= 2 idhep(2)= -2 idhep(3)= 4 idhep(4)= -4 else if(ipr.eq.16) then idhep(1)= 14 idhep(2)= -14 idhep(3)= 12 idhep(4)= -12 else if(ipr.eq.17) then idhep(1)= 13 idhep(2)= -13 idhep(3)= 12 idhep(4)= -12 else if(ipr.eq.18) then idhep(1)= 2 idhep(2)= -2 idhep(3)= 12 idhep(4)= -12 else if(ipr.eq.19) then idhep(1)= 1 idhep(2)= -1 idhep(3)= 12 idhep(4)= -12 else if(ipr.eq.20) then idhep(1)= 14 idhep(2)= -13 idhep(3)= 13 idhep(4)= -14 else if(ipr.eq.21) then if(ockm.eq.'n') then idhep(1)= 2 idhep(2)= -1 idhep(3)= 1 idhep(4)= -2 else if(ickm.eq.1) then idhep(1)= 2 idhep(2)= -1 idhep(3)= 1 idhep(4)= -2 else if(ickm.eq.2) then idhep(1)= 2 idhep(2)= -3 idhep(3)= 3 idhep(4)= -2 else if(ickm.eq.3) then idhep(1)= 2 idhep(2)= -5 idhep(3)= 5 idhep(4)= -2 else if(ickm.eq.4) then idhep(1)= 4 idhep(2)= -1 idhep(3)= 1 idhep(4)= -4 else if(ickm.eq.5) then idhep(1)= 4 idhep(2)= -3 idhep(3)= 3 idhep(4)= -4 else if(ickm.eq.6) then idhep(1)= 4 idhep(2)= -5 idhep(3)= 5 idhep(4)= -4 endif endif else if(ipr.eq.22) then idhep(1)= 13 idhep(2)= -13 idhep(3)= 11 idhep(4)= -11 else if(ipr.eq.23) then idhep(1)= 14 idhep(2)= -14 idhep(3)= 11 idhep(4)= -11 else if(ipr.eq.24) then idhep(1)= 2 idhep(2)= -2 idhep(3)= 11 idhep(4)= -11 else if(ipr.eq.25) then idhep(1)= 1 idhep(2)= -1 idhep(3)= 11 idhep(4)= -11 else if(ipr.eq.26) then idhep(1)= 13 idhep(2)= -13 idhep(3)= 13 idhep(4)= -13 else if(ipr.eq.27) then idhep(1)= 14 idhep(2)= -14 idhep(3)= 14 idhep(4)= -14 else if(ipr.eq.28) then idhep(1)= 2 idhep(2)= -2 idhep(3)= 2 idhep(4)= -2 else if(ipr.eq.29) then idhep(1)= 1 idhep(2)= -1 idhep(3)= 1 idhep(4)= -1 else if(ipr.eq.30) then idhep(1)= 5 idhep(2)= -5 idhep(3)= 14 idhep(4)= -14 else if(ipr.eq.31) then idhep(1)= 13 idhep(2)= -13 idhep(3)= 5 idhep(4)= -5 else if(ipr.eq.32) then idhep(1)= 5 idhep(2)= -5 idhep(3)= 12 idhep(4)= -12 else if(ipr.eq.33) then idhep(1)= 5 idhep(2)= -5 idhep(3)= 11 idhep(4)= -11 else if(ipr.eq.34) then idhep(1)= 2 idhep(2)= -2 idhep(3)= 5 idhep(4)= -5 else if(ipr.eq.35) then idhep(1)= 1 idhep(2)= -1 idhep(3)= 5 idhep(4)= -5 else if(ipr.eq.36) then idhep(1)= 5 idhep(2)= -5 idhep(3)= 5 idhep(4)= -5 else if(ipr.eq.37) then idhep(1)= 15 idhep(2)= -15 idhep(3)= 5 idhep(4)= -5 else if(ipr.eq.38) then idhep(1)= 15 idhep(2)= -16 idhep(3)= 4 idhep(4)= -3 endif * do i=1,itmx gxp= prxa(2,i) gxm= prxa(1,i) ge1= prxa(3,i)+prxa(6,i)+prxa(7,i) ge2= 1.d0+prxa(3,i)-prxa(4,i)-ge1 ge3= prxa(4,i)+prxa(5,i)+prxa(7,i) ge4= 1.d0+prxa(4,i)-prxa(3,i)-ge3 gt1= prxa(9,i) gt2= prxa(8,i)-prxa(9,i) gt3= prxa(10,i) gt4= 1.d0-prxa(8,i)-gt3 * phep(4,1)= gxp*ge1-(gxp-gxm)*gt1 phep(4,2)= gxp*ge2-(gxp-gxm)*gt2 phep(4,3)= gxp*ge3-(gxp-gxm)*gt3 phep(4,4)= gxp*ge4-(gxp-gxm)*gt4 phep(3,1)= gxp*ge1-(gxp+gxm)*gt1 phep(3,2)= gxp*ge2-(gxp+gxm)*gt2 phep(3,3)= gxp*ge3-(gxp+gxm)*gt3 phep(3,4)= gxp*ge4-(gxp+gxm)*gt4 phep(2,1)= 0.d0 phep(1,1)= phep(4,1)*phep(4,1)-phep(3,1)*phep(3,1) phep(1,1)= sqrt(phep(1,1)) phep(1,2)= (-phep(3,1)*phep(3,2)+phep(4,1)*phep(4,2)- # 2.d0*gxp*gxm*prxa(3,i))/phep(1,1) phep(2,2)= phep(4,2)*phep(4,2)-phep(3,2)*phep(3,2)- # phep(1,2)*phep(1,2) phep(2,2)= sqrt(phep(2,2)) phep(1,3)= (-phep(3,1)*phep(3,3)+phep(4,1)*phep(4,3)- # 2.d0*gxp*gxm*prxa(7,i))/phep(1,1) phep(2,3)= phep(4,3)*phep(4,3)-phep(3,3)*phep(3,3)- # phep(1,3)*phep(1,3) phep(2,3)= sqrt(phep(2,3)) gsgn= (-2.d0*gxp*gxm*prxa(5,i)-phep(1,2)*phep(1,3)- # phep(3,2)*phep(3,3)+phep(4,2)*phep(4,3))/phep(2,2)/ # phep(2,3) if(gsgn.gt.0.d0) then phep(2,3)= +phep(2,3) else phep(2,3)= -phep(2,3) endif phep(1,4)= -phep(1,1)-phep(1,2)-phep(1,3) phep(2,4)= -phep(2,1)-phep(2,2)-phep(2,3) phep(3,4)= -phep(3,1)-phep(3,2)-phep(3,3)+gxp-gxm phep(4,4)= -phep(4,1)-phep(4,2)-phep(3,4)+gxp+gxm do j=1,4 phep(5,j)= 0.d0 enddo if(osmass.eq.'y') then do k=1,4 sclm(k)= (phep(4,k)*phep(4,k)-fsm(k)*fsm(k))/ # (phep(1,k)*phep(1,k)+phep(2,k)*phep(2,k)+ # phep(3,k)*phep(3,k)) do kp=1,3 phep(kp,k)= sqrt(sclm(k))*phep(kp,k) enddo enddo endif * write(nout,fmt=*) nhep write(nout,fmt=*) (phep(j,1),j=1,5) write(nout,fmt=*) (phep(j,2),j=1,5) write(nout,fmt=*) (phep(j,3),j=1,5) write(nout,fmt=*) (phep(j,4),j=1,5) enddo close(nout) else if(ostore.eq.'i') then open(nout,file='inv',status='new') do i=1,itmx write(nout,fmt=*) (prxa(j,i),j=1,10) enddo close(nout) endif print 2503,itmx * *-----CC or Higgs storage ends * 2500 format(1x,4e20.7) 2501 format(1x,5e20.7) 2502 format(1x,3e20.7) 2503 format(1x,i6,' events successfully stored ') * *-----end generation branch * endif else * *-----start evaluation branch * npts= jnpts itrans= 0 nrand= inrand ifail= 0 if(oxcm.eq.'c') then if(otype.eq.'cc12') then call d01gcf(ndim,wtoxssc,wtoregion,npts,vk,nrand, # itrans,sig,esig,ifail) else if(otype.eq.'cc20'.and.osinw.eq.'e') then separ= -cos(separa*pi/180.d0) sigext= 0.d0 esigext= 0.d0 go to 9998 if((-asccsw).lt.separ) then sigext= 0.d0 esigext= 0.d0 else iac(1)= 1 do i=1,5 aim(i)= 0.d0 bim(i)= 1.d0 enddo aim(6)= sthqsw bim(6)= 1.d0 iac(2)= 0 ae(1)= 0.d0 ae(2)= 0.d0 ae(3)= 0.d0 ae(4)= 0.d0 iac(3)= 1 asa(1)= separ bsa(1)= +1.d0 do i=2,4 asa(i)= -1.d0 bsa(i)= +1.d0 enddo iac(4)= 0 do i=1,6 afsa(i)= -1.d0 bfsa(i)= +1.d0 enddo isaa= 0 isab= 0 do i=1,4 ombsa(i)= 1.d0-bsa(i) omasa(i)= 1.d0-asa(i) opbsa(i)= 1.d0+bsa(i) opasa(i)= 1.d0+asa(i) enddo do j=1,4 sgam(j)= 0.5d0*omasa(j) cgam(j)= 0.5d0*ombsa(j) enddo sg12= 0.5d0*(1.d0-afsa(1)) cg12= 0.5d0*(1.d0-bfsa(1)) sg13= 0.5d0*(1.d0-afsa(2)) cg13= 0.5d0*(1.d0-bfsa(2)) sg14= 0.5d0*(1.d0-afsa(3)) cg14= 0.5d0*(1.d0-bfsa(3)) sg23= 0.5d0*(1.d0-afsa(4)) cg23= 0.5d0*(1.d0-bfsa(4)) sg24= 0.5d0*(1.d0-afsa(5)) cg24= 0.5d0*(1.d0-bfsa(5)) sg34= 0.5d0*(1.d0-afsa(6)) cg34= 0.5d0*(1.d0-bfsa(6)) call g05cbf(0) call d01gcf(ndim,wtoxsc,wtoregion,npts,vk,nrand, # itrans,sig2,esig2,ifail) if(oprint.eq.'y') then print 652 print 20,ifail print 651,(ifz(i),i=1,51) ifzt= 0 do j=1,14 ifzt= ifzt+4*ifz(j) enddo do j=15,41 ifzt= ifzt+2*ifz(j) enddo do j=42,51 ifzt= ifzt+ifz(j) enddo if(ifzt.ne.0) then print 653,ifzt endif print 60,ik,4*ik print 65,1.d2*(1.d0-0.25d0*ifzt/ik) print 654 call wtooutput(ipr) print 2000,sig2,esig2 print 2001,1.d2*esig2/sig2 endif * iac(1)= 1 do i=1,5 aim(i)= 0.d0 bim(i)= 1.d0 enddo aim(6)= sthqsw bim(6)= 1.d0 iac(2)= 0 ae(1)= 0.d0 ae(2)= 0.d0 ae(3)= 0.d0 ae(4)= 0.d0 iac(3)= 1 asa(1)= -asccsw bsa(1)= +asccsw do i=2,4 asa(i)= -1.d0 bsa(i)= +1.d0 enddo iac(4)= 0 do i=1,6 afsa(i)= -1.d0 bfsa(i)= +1.d0 enddo isaa= 0 isab= 0 do i=1,4 ombsa(i)= 1.d0-bsa(i) omasa(i)= 1.d0-asa(i) opbsa(i)= 1.d0+bsa(i) opasa(i)= 1.d0+asa(i) enddo do j=1,4 sgam(j)= 0.5d0*omasa(j) cgam(j)= 0.5d0*ombsa(j) enddo sg12= 0.5d0*(1.d0-afsa(1)) cg12= 0.5d0*(1.d0-bfsa(1)) sg13= 0.5d0*(1.d0-afsa(2)) cg13= 0.5d0*(1.d0-bfsa(2)) sg14= 0.5d0*(1.d0-afsa(3)) cg14= 0.5d0*(1.d0-bfsa(3)) sg23= 0.5d0*(1.d0-afsa(4)) cg23= 0.5d0*(1.d0-bfsa(4)) sg24= 0.5d0*(1.d0-afsa(5)) cg24= 0.5d0*(1.d0-bfsa(5)) sg34= 0.5d0*(1.d0-afsa(6)) cg34= 0.5d0*(1.d0-bfsa(6)) ifail= 0 ik= 0 do i=1,51 ifz(i)= 0 enddo call g05cbf(0) call d01gcf(ndim,wtoxsc,wtoregion,npts,vk,nrand, # itrans,sig1,esig1,ifail) sigext= sig2-sig1 esigext= esig1*esig1+esig2*esig2 esigext= sqrt(esigext) if(oprint.eq.'y') then print 652 print 20,ifail print 651,(ifz(i),i=1,51) ifzt= 0 do j=1,14 ifzt= ifzt+4*ifz(j) enddo do j=15,41 ifzt= ifzt+2*ifz(j) enddo do j=42,51 ifzt= ifzt+ifz(j) enddo if(ifzt.ne.0) then print 653,ifzt endif print 60,ik,4*ik print 65,1.d2*(1.d0-0.25d0*ifzt/ik) print 654 call wtooutput(ipr) print 2000,sig1,esig1 print 2001,1.d2*esig1/sig1 endif endif * 9998 continue iac(1)= 1 do i=1,5 aim(i)= 0.d0 bim(i)= 1.d0 enddo aim(6)= sthqsw bim(6)= 1.d0 iac(2)= 0 ae(1)= 0.d0 ae(2)= 0.d0 ae(3)= 0.d0 ae(4)= 0.d0 iac(3)= 1 asa(1)= -1.d0 bsa(1)= separ do i=2,4 asa(i)= -1.d0 bsa(i)= +1.d0 enddo iac(4)= 0 do i=1,6 afsa(i)= -1.d0 bfsa(i)= +1.d0 enddo isaa= 0 isab= 0 do i=1,4 ombsa(i)= 1.d0-bsa(i) omasa(i)= 1.d0-asa(i) opbsa(i)= 1.d0+bsa(i) opasa(i)= 1.d0+asa(i) enddo do j=1,4 sgam(j)= 0.5d0*omasa(j) cgam(j)= 0.5d0*ombsa(j) enddo sg12= 0.5d0*(1.d0-afsa(1)) cg12= 0.5d0*(1.d0-bfsa(1)) sg13= 0.5d0*(1.d0-afsa(2)) cg13= 0.5d0*(1.d0-bfsa(2)) sg14= 0.5d0*(1.d0-afsa(3)) cg14= 0.5d0*(1.d0-bfsa(3)) sg23= 0.5d0*(1.d0-afsa(4)) cg23= 0.5d0*(1.d0-bfsa(4)) sg24= 0.5d0*(1.d0-afsa(5)) cg24= 0.5d0*(1.d0-bfsa(5)) sg34= 0.5d0*(1.d0-afsa(6)) cg34= 0.5d0*(1.d0-bfsa(6)) tc= separa/180.d0*pi tcs= tc*tc s0cut= aim(6)*aim(6) if(ipr.eq.4) then fcolsw= 1.d0 else if(ipr.eq.5) then fcolsw= 3.d0 endif tfactsw= g8*sth4/2048.d0/pi**5*fcolsw*cfct if(iosf.eq.1) then ndim= 6 else if(iosf.eq.0) then ndim= 4 endif if(inpts.le.6) then jnpts= inpts else if(inpts.eq.7) then jnpts= 2011*101 if(ndim.eq.4) then vk(1)= 0.1d+01 vk(2)= 0.12294d+05 vk(3)= 0.27852d+05 vk(4)= 0.170453d+06 else if(ndim.eq.6) then vk(1)= 0.1d+01 vk(2)= 0.9463d+05 vk(3)= 0.79132d+05 vk(4)= 0.167923d+06 vk(5)= 0.164405d+06 vk(6)= 0.154994d+06 endif else if(inpts.eq.8) then jnpts= 10193*101 if(ndim.eq.4) then vk(1)= 0.1d+01 vk(2)= 0.9601d+06 vk(3)= 0.449688d+06 vk(4)= 0.792432d+06 endif else if(inpts.eq.9) then jnpts= 22807*151 if(ndim.eq.4) then vk(1)= 0.1d+01 vk(2)= 0.2670673d+07 vk(3)= 0.124894d+07 vk(4)= 0.521697d+06 endif endif npts= jnpts itrans= 0 nrand= inrand ifail= 0 ik= 0 do i=1,51 ifz(i)= 0 enddo call g05cbf(0) call d01gcf(ndim,wtoxssw,wtoregion,npts,vk,nrand, # itrans,sigint,esigint,ifail) sig= sigint+sigext esig= esigint*esigint+esigext*esigext esig= sqrt(esig) if(oprint.eq.'y') then print 652 print 20,ifail print 651,(ifz(i),i=1,7) ifzt= 0 do j=1,6 ifzt= ifzt+5*ifz(j) enddo ifzt= ifzt+ifz(7) if(ifzt.ne.0) then print 653,ifzt endif print 60,ik,5*ik print 65,1.d2*(1.d0-ifzt/(ik*5.d0)) print 654 call wtooutput(ipr) print 2000,sigint,esigint print 2001,1.d2*esigint/sigint endif else if(otype.eq.'cc20'.and.osinw.eq.'c') then separ= -cos(separa*pi/180.d0) sigext= 0.d0 esigext= 0.d0 go to 9999 if((-asccsw).lt.separ) then sigext= 0.d0 esigext= 0.d0 else iac(1)= 0 do i=1,5 aim(i)= 0.d0 bim(i)= 1.d0 enddo aim(6)= sthqsw bim(6)= 1.d0 iac(2)= 1 ae(1)= 0.d0 ae(2)= 0.d0 ae(3)= eusw ae(4)= edsw iac(3)= 1 asa(1)= separ bsa(1)= +1.d0 do i=2,3 asa(i)= -1.d0 bsa(i)= +1.d0 enddo asa(4)= -alsw bsa(4)= +alsw iac(4)= 0 do i=1,6 afsa(i)= -1.d0 bfsa(i)= +1.d0 enddo isaa= 0 isab= 0 do i=1,4 ombsa(i)= 1.d0-bsa(i) omasa(i)= 1.d0-asa(i) opbsa(i)= 1.d0+bsa(i) opasa(i)= 1.d0+asa(i) enddo do j=1,4 sgam(j)= 0.5d0*omasa(j) cgam(j)= 0.5d0*ombsa(j) enddo sg12= 0.5d0*(1.d0-afsa(1)) cg12= 0.5d0*(1.d0-bfsa(1)) sg13= 0.5d0*(1.d0-afsa(2)) cg13= 0.5d0*(1.d0-bfsa(2)) sg14= 0.5d0*(1.d0-afsa(3)) cg14= 0.5d0*(1.d0-bfsa(3)) sg23= 0.5d0*(1.d0-afsa(4)) cg23= 0.5d0*(1.d0-bfsa(4)) sg24= 0.5d0*(1.d0-afsa(5)) cg24= 0.5d0*(1.d0-bfsa(5)) sg34= 0.5d0*(1.d0-afsa(6)) cg34= 0.5d0*(1.d0-bfsa(6)) call g05cbf(0) call d01gcf(ndim,wtoxsc,wtoregion,npts,vk,nrand, # itrans,sig2,esig2,ifail) if(oprint.eq.'y') then print 652 print 20,ifail print 651,(ifz(i),i=1,51) ifzt= 0 do j=1,14 ifzt= ifzt+4*ifz(j) enddo do j=15,41 ifzt= ifzt+2*ifz(j) enddo do j=42,51 ifzt= ifzt+ifz(j) enddo if(ifzt.ne.0) then print 653,ifzt endif print 60,ik,4*ik print 65,1.d2*(1.d0-0.25d0*ifzt/ik) print 654 call wtooutput(ipr) print 2000,sig2,esig2 print 2001,1.d2*esig2/sig2 endif * iac(1)= 0 do i=1,5 aim(i)= 0.d0 bim(i)= 1.d0 enddo aim(6)= sthqsw bim(6)= 1.d0 iac(2)= 1 ae(1)= 0.d0 ae(2)= 0.d0 ae(3)= eusw ae(4)= edsw iac(3)= 1 asa(1)= -asccsw bsa(1)= +asccsw do i=2,3 asa(i)= -1.d0 bsa(i)= +1.d0 enddo asa(4)= -alsw bsa(4)= +alsw iac(4)= 0 do i=1,6 afsa(i)= -1.d0 bfsa(i)= +1.d0 enddo isaa= 0 isab= 0 do i=1,4 ombsa(i)= 1.d0-bsa(i) omasa(i)= 1.d0-asa(i) opbsa(i)= 1.d0+bsa(i) opasa(i)= 1.d0+asa(i) enddo do j=1,4 sgam(j)= 0.5d0*omasa(j) cgam(j)= 0.5d0*ombsa(j) enddo sg12= 0.5d0*(1.d0-afsa(1)) cg12= 0.5d0*(1.d0-bfsa(1)) sg13= 0.5d0*(1.d0-afsa(2)) cg13= 0.5d0*(1.d0-bfsa(2)) sg14= 0.5d0*(1.d0-afsa(3)) cg14= 0.5d0*(1.d0-bfsa(3)) sg23= 0.5d0*(1.d0-afsa(4)) cg23= 0.5d0*(1.d0-bfsa(4)) sg24= 0.5d0*(1.d0-afsa(5)) cg24= 0.5d0*(1.d0-bfsa(5)) sg34= 0.5d0*(1.d0-afsa(6)) cg34= 0.5d0*(1.d0-bfsa(6)) ifail= 0 ik= 0 do i=1,51 ifz(i)= 0 enddo call g05cbf(0) call d01gcf(ndim,wtoxsc,wtoregion,npts,vk,nrand, # itrans,sig1,esig1,ifail) sigext= sig2-sig1 esigext= esig1*esig1+esig2*esig2 esigext= sqrt(esigext) if(oprint.eq.'y') then print 652 print 20,ifail print 651,(ifz(i),i=1,51) ifzt= 0 do j=1,14 ifzt= ifzt+4*ifz(j) enddo do j=15,41 ifzt= ifzt+2*ifz(j) enddo do j=42,51 ifzt= ifzt+ifz(j) enddo if(ifzt.ne.0) then print 653,ifzt endif print 60,ik,4*ik print 65,1.d2*(1.d0-0.25d0*ifzt/ik) print 654 call wtooutput(ipr) print 2000,sig1,esig1 print 2001,1.d2*esig1/sig1 endif endif 9999 continue iac(1)= 0 do i=1,5 aim(i)= 0.d0 bim(i)= 1.d0 enddo aim(6)= sthqsw bim(6)= 1.d0 iac(2)= 1 ae(1)= 0.d0 ae(2)= 0.d0 ae(3)= eusw ae(4)= edsw iac(3)= 1 asa(1)= -1.d0 bsa(1)= separ do i=2,4 asa(i)= -1.d0 bsa(i)= +1.d0 enddo do i=2,3 asa(i)= -1.d0 bsa(i)= +1.d0 enddo asa(4)= -alsw bsa(4)= +alsw iac(4)= 0 do i=1,6 afsa(i)= -1.d0 bfsa(i)= +1.d0 enddo isaa= 0 isab= 0 do i=1,4 ombsa(i)= 1.d0-bsa(i) omasa(i)= 1.d0-asa(i) opbsa(i)= 1.d0+bsa(i) opasa(i)= 1.d0+asa(i) enddo do j=1,4 sgam(j)= 0.5d0*omasa(j) cgam(j)= 0.5d0*ombsa(j) enddo sg12= 0.5d0*(1.d0-afsa(1)) cg12= 0.5d0*(1.d0-bfsa(1)) sg13= 0.5d0*(1.d0-afsa(2)) cg13= 0.5d0*(1.d0-bfsa(2)) sg14= 0.5d0*(1.d0-afsa(3)) cg14= 0.5d0*(1.d0-bfsa(3)) sg23= 0.5d0*(1.d0-afsa(4)) cg23= 0.5d0*(1.d0-bfsa(4)) sg24= 0.5d0*(1.d0-afsa(5)) cg24= 0.5d0*(1.d0-bfsa(5)) sg34= 0.5d0*(1.d0-afsa(6)) cg34= 0.5d0*(1.d0-bfsa(6)) * tc= separa/180.d0*pi tcs= tc*tc s0cut= aim(6)*aim(6) if(ipr.eq.4) then fcolsw= 1.d0 else if(ipr.eq.5) then fcolsw= 3.d0 endif tfactsw= g8*sth4/8192.d0/pi**6*fcolsw*cfct do ii=1,3 do jj=1,6 inm= ii jnm= jj npts= jnpts itrans= 0 nrand= inrand ifail= 0 ik= 0 do i=1,51 ifz(i)= 0 enddo call g05cbf(0) call d01gcf(ndim,wtoxsswc,wtoregion,npts,vk,nrand, # itrans,sigp,esigp,ifail) sigpart(ii,jj)= sigp esigpart(ii,jj)= esigp if(oprint.eq.'y') then print 652 print 20,ifail print 651,(ifz(i),i=1,47) ifzt= 0 do j=1,32 ifzt= ifzt+2*ifz(j) enddo do j=33,47 ifzt= ifzt+ifz(j) enddo if(ifzt.ne.0) then print 653,ifzt endif print 60,ik,2*ik print 65,1.d2*(1.d0-ifzt/(ik*2.d0)) print 654 call wtooutput(ipr) endif enddo enddo sigint= 0.d0 esigint= 0.d0 do ii=1,3 do jj=1,6 sigint= sigint+sigpart(ii,jj) esigint= esigint+esigpart(ii,jj)*esigpart(ii,jj) enddo enddo esigint= sqrt(esigint) sig= sigint+sigext esig= esigint*esigint+esigext*esigext esig= sqrt(esig) print 2000,sig,esig * *----------------------------------------------------------------------- * else if(otype.eq.'cc20'.and.osinw.eq.'f') then separ= -cos(separa*pi/180.d0) * sigext= 0.d0 * esigext= 0.d0 go to 9997 iac(1)= 0 do i=1,5 aim(i)= 0.d0 bim(i)= 1.d0 enddo aim(6)= sthqsw bim(6)= 1.d0 iac(2)= 1 ae(1)= 0.d0 ae(2)= 0.d0 ae(3)= eusw ae(4)= edsw iac(3)= 1 asa(1)= -1.d0 bsa(1)= +1.d0 do i=2,3 asa(i)= -1.d0 bsa(i)= +1.d0 enddo asa(4)= -alsw bsa(4)= +alsw iac(4)= 0 do i=1,6 afsa(i)= -1.d0 bfsa(i)= +1.d0 enddo isaa= 0 isab= 0 do i=1,4 ombsa(i)= 1.d0-bsa(i) omasa(i)= 1.d0-asa(i) opbsa(i)= 1.d0+bsa(i) opasa(i)= 1.d0+asa(i) enddo do j=1,4 sgam(j)= 0.5d0*omasa(j) cgam(j)= 0.5d0*ombsa(j) enddo sg12= 0.5d0*(1.d0-afsa(1)) cg12= 0.5d0*(1.d0-bfsa(1)) sg13= 0.5d0*(1.d0-afsa(2)) cg13= 0.5d0*(1.d0-bfsa(2)) sg14= 0.5d0*(1.d0-afsa(3)) cg14= 0.5d0*(1.d0-bfsa(3)) sg23= 0.5d0*(1.d0-afsa(4)) cg23= 0.5d0*(1.d0-bfsa(4)) sg24= 0.5d0*(1.d0-afsa(5)) cg24= 0.5d0*(1.d0-bfsa(5)) sg34= 0.5d0*(1.d0-afsa(6)) cg34= 0.5d0*(1.d0-bfsa(6)) call g05cbf(0) call d01gcf(ndim,wtoxsc,wtoregion,npts,vk,nrand, # itrans,sig2,esig2,ifail) if(oprint.eq.'y') then print 652 print 20,ifail print 651,(ifz(i),i=1,51) ifzt= 0 do j=1,14 ifzt= ifzt+4*ifz(j) enddo do j=15,41 ifzt= ifzt+2*ifz(j) enddo do j=42,51 ifzt= ifzt+ifz(j) enddo if(ifzt.ne.0) then print 653,ifzt endif print 60,ik,4*ik print 65,1.d2*(1.d0-0.25d0*ifzt/ik) print 654 call wtooutput(ipr) print 2000,sig2,esig2 print 2001,1.d2*esig2/sig2 endif * * 9997 sig2= 0.d0 * esig2= 0.d0 iac(1)= 0 do i=1,5 aim(i)= 0.d0 bim(i)= 1.d0 enddo aim(6)= sthqsw bim(6)= 1.d0 iac(2)= 1 ae(1)= 0.d0 ae(2)= 0.d0 ae(3)= eusw ae(4)= edsw iac(3)= 1 asa(1)= -asccsw bsa(1)= +1.d0 do i=2,3 asa(i)= -1.d0 bsa(i)= +1.d0 enddo asa(4)= -alsw bsa(4)= +alsw iac(4)= 0 do i=1,6 afsa(i)= -1.d0 bfsa(i)= +1.d0 enddo isaa= 0 isab= 0 do i=1,4 ombsa(i)= 1.d0-bsa(i) omasa(i)= 1.d0-asa(i) opbsa(i)= 1.d0+bsa(i) opasa(i)= 1.d0+asa(i) enddo do j=1,4 sgam(j)= 0.5d0*omasa(j) cgam(j)= 0.5d0*ombsa(j) enddo sg12= 0.5d0*(1.d0-afsa(1)) cg12= 0.5d0*(1.d0-bfsa(1)) sg13= 0.5d0*(1.d0-afsa(2)) cg13= 0.5d0*(1.d0-bfsa(2)) sg14= 0.5d0*(1.d0-afsa(3)) cg14= 0.5d0*(1.d0-bfsa(3)) sg23= 0.5d0*(1.d0-afsa(4)) cg23= 0.5d0*(1.d0-bfsa(4)) sg24= 0.5d0*(1.d0-afsa(5)) cg24= 0.5d0*(1.d0-bfsa(5)) sg34= 0.5d0*(1.d0-afsa(6)) cg34= 0.5d0*(1.d0-bfsa(6)) ifail= 0 ik= 0 do i=1,51 ifz(i)= 0 enddo call g05cbf(0) call d01gcf(ndim,wtoxsc,wtoregion,npts,vk,nrand, # itrans,sig1,esig1,ifail) sigext= sig2-sig1 esigext= esig1*esig1+esig2*esig2 esigext= sqrt(esigext) if(oprint.eq.'y') then print 652 print 20,ifail print 651,(ifz(i),i=1,51) ifzt= 0 do j=1,14 ifzt= ifzt+4*ifz(j) enddo do j=15,41 ifzt= ifzt+2*ifz(j) enddo do j=42,51 ifzt= ifzt+ifz(j) enddo if(ifzt.ne.0) then print 653,ifzt endif print 60,ik,4*ik print 65,1.d2*(1.d0-0.25d0*ifzt/ik) print 654 if(ofl.eq.'c') then print 2002,1.d0*idtl,1.d0*idt4,1.d0*idt7,1.d0*idt9 endif call wtooutput(ipr) print 2000,sig1,esig1 print 2001,1.d2*esig1/sig1 endif * 9997 continue iac(1)= 0 do i=1,5 aim(i)= 0.d0 bim(i)= 1.d0 enddo aim(6)= sthqsw bim(6)= 1.d0 iac(2)= 1 ae(1)= 0.d0 ae(2)= 0.d0 ae(3)= eusw ae(4)= edsw iac(3)= 1 asa(1)= -1.d0 bsa(1)= separ do i=2,4 asa(i)= -1.d0 bsa(i)= +1.d0 enddo do i=2,3 asa(i)= -1.d0 bsa(i)= +1.d0 enddo asa(4)= -alsw bsa(4)= +alsw iac(4)= 0 do i=1,6 afsa(i)= -1.d0 bfsa(i)= +1.d0 enddo isaa= 0 isab= 0 do i=1,4 ombsa(i)= 1.d0-bsa(i) omasa(i)= 1.d0-asa(i) opbsa(i)= 1.d0+bsa(i) opasa(i)= 1.d0+asa(i) enddo do j=1,4 sgam(j)= 0.5d0*omasa(j) cgam(j)= 0.5d0*ombsa(j) enddo sg12= 0.5d0*(1.d0-afsa(1)) cg12= 0.5d0*(1.d0-bfsa(1)) sg13= 0.5d0*(1.d0-afsa(2)) cg13= 0.5d0*(1.d0-bfsa(2)) sg14= 0.5d0*(1.d0-afsa(3)) cg14= 0.5d0*(1.d0-bfsa(3)) sg23= 0.5d0*(1.d0-afsa(4)) cg23= 0.5d0*(1.d0-bfsa(4)) sg24= 0.5d0*(1.d0-afsa(5)) cg24= 0.5d0*(1.d0-bfsa(5)) sg34= 0.5d0*(1.d0-afsa(6)) cg34= 0.5d0*(1.d0-bfsa(6)) * tc= separa/180.d0*pi tcs= tc*tc s0cut= aim(6)*aim(6) if(ipr.eq.4) then fcolsw= 1.d0 else if(ipr.eq.5) then fcolsw= 3.d0 endif tfactsw= fcolsw*cfct/8192.d0/pi**6 if(iosf.eq.0) then iimn= 2 else iimn= 1 endif do ii=iimn,4 do jj=1,6 inm= ii jnm= jj npts= jnpts itrans= 0 nrand= inrand ifail= 0 ik= 0 do i=1,51 ifz(i)= 0 enddo call g05cbf(0) call d01gcf(ndim,wtoxsswcfl,wtoregion,npts,vk,nrand, # itrans,sigp,esigp,ifail) print 5000,sigp,esigp 5000 format(1x,2e20.5) sigfl(ii,jj)= sigp esigfl(ii,jj)= esigp if(oprint.eq.'y') then print 652 print 20,ifail print 651,(ifz(i),i=1,40) ifzt= 0 do j=1,27 ifzt= ifzt+2*ifz(j) enddo do j=28,40 ifzt= ifzt+ifz(j) enddo if(ifzt.ne.0) then print 653,ifzt endif print 60,ik,2*ik print 65,1.d2*(1.d0-ifzt/(ik*2.d0)) print 654 call wtooutput(ipr) endif enddo enddo sigint= 0.d0 esigint= 0.d0 do ii=iimn,4 do jj=1,6 sigint= sigint+sigfl(ii,jj) esigint= esigint+esigfl(ii,jj)*esigfl(ii,jj) enddo enddo esigint= sqrt(esigint) sig= sigint+sigext esig= esigint*esigint+esigext*esigext esig= sqrt(esig) print 2000,sig,esig * *-------------------------------------------------------------------------- * else call d01gcf(ndim,wtoxsc,wtoregion,npts,vk,nrand, # itrans,sig,esig,ifail) endif endif else if(oxcm.eq.'n') then if(otype.eq.'nc68') then if(omssm.eq.'n') then call d01gcf(ndim,wtoxsnh64,wtoregion,npts,vk, # nrand,itrans,sig,esig,ifail) else if(omssm.eq.'y') then call d01gcf(ndim,wtoxsnsh64,wtoregion,npts,vk, # nrand,itrans,sig,esig,ifail) endif else if(otype.eq.'nc50') then call d01gcf(ndim,wtoxsnh49,wtoregion,npts,vk, # nrand,itrans,sig,esig,ifail) else if(otype.eq.'nc33') then call d01gcf(ndim,wtoxsnh32,wtoregion,npts,vk, # nrand,itrans,sig,esig,ifail) else if(otype.eq.'nc25') then call d01gcf(ndim,wtoxsnh24,wtoregion,npts,vk, # nrand,itrans,sig,esig,ifail) else if(otype.eq.'nc26') then call d01gcf(ndim,wtoxsnh26,wtoregion,npts,vk, # nrand,itrans,sig,esig,ifail) else if(otype.eq.'nc21') then call d01gcf(ndim,wtoxsnh19,wtoregion,npts,vk, # nrand,itrans,sig,esig,ifail) else if(otype.eq.'nc32') then call d01gcf(ndim,wtoxsn32,wtoregion,npts,vk, # nrand,itrans,sig,esig,ifail) else if(otype.eq.'nc64') then call d01gcf(ndim,wtoxsn64,wtoregion,npts,vk, # nrand,itrans,sig,esig,ifail) else call d01gcf(ndim,wtoxsn,wtoregion,npts,vk, # nrand,itrans,sig,esig,ifail) endif else if(oxcm.eq.'m') then call d01gcf(ndim,wtoxsm43,wtoregion,npts,vk, # nrand,itrans,sig,esig,ifail) endif * *-----end of generation-evaluation branch * endif * if(om.eq.'e') then if(otype.eq.'cc20'.and.osinw.ne.'n') then continue else if(oprint.eq.'y') then print 652 print 20,ifail print 651,(ifz(i),i=1,51) 652 format(/' Call DUMP -----------------------------------',/) 20 format(/' On exit IFAIL = ',i1,/) 651 format(1x,5i10) * ifzt= 0 if(omdist.eq.'y') then do j=1,51 ifzt= ifzt+ifz(j) enddo else do j=1,14 ifzt= ifzt+4*ifz(j) enddo do j=15,41 ifzt= ifzt+2*ifz(j) enddo do j=42,51 ifzt= ifzt+ifz(j) enddo endif if(ifzt.ne.0) then print 653,ifzt endif if(omdist.eq.'y') then print 60,ik,ik print 65,1.d2*(1.d0-1.d0*ifzt/ik) else print 60,ik,4*ik print 65,1.d2*(1.d0-0.25d0*ifzt/ik) endif print 654 if(ofl.eq.'c') then print 2002,1.d0*idtl,1.d0*idt4,1.d0*idt7,1.d0*idt9 endif call wtooutput(ipr) endif endif 60 format(/' # of calls = ',i10,1x, # ' # of integrands = ',i10) 65 format(/' 1-nc/int (%) = ',f10.3) 653 format(/' Total # of null calls = ',i10) 654 format(/' End of Call DUMP-----------------------------',/) endif 1999 format(/' Partial cross-section ',i2,' = ',e20.5) 2000 format(' sigma (nb) =',e20.5,3x,'+-',3x,e20.5) 2001 format(/' Rel. error of ',f10.3,' %',///) 2002 format(/' # of 0-4-7-8 Gram det ',/4e20.4) * return end