! ! Note that several changes were made to the code below. These include: ! ! 1. Wrapped main program as function "genevents" which is now called from ! C. ! ! 2. Added tgtvars and initialstate common blocks and moved the rhob ! common block all into the cppvars.h file that is #included in ! several functions below. ! ! 3. Changed names of members in rhob common block and then used ! EQUIVALENCE statements to reconnect them with their original ! names inside functions. ! ! 4. Commented out data lines setting values in common blocks that are ! now set in the C-program. ! ! 5. Write output in stdhep ascii format ! ! n.b. Comments associtated with these changes use '!' to help make ! them eaiser to identify. ! c c To compile: ‘gfortran -ffixed-line-length-132 2pi_event_gen_v7.f -o event_gen.exe’ c To run: ‘./event_gen.exe’ c real*8 FUNCTION genevents() implicit none ! Here we include the common block definitions and then declare some of the members ! as equivalent to some local variables. This is done to preserve the original naming ! of the local variables, while avoiding conflicts in other functions below that ! define some (but not all) of the variables. #include "cppvars.h" EQUIVALENCE (z_target,ztgt),(a_target,Atgt),(target_len,tgtlen),(radiation_len,radlen) EQUIVALENCE (rho_norm,norm),(rho_0,rho0),(c_density,c_den),(a_density,a_den) EQUIVALENCE (Ebeam,Egamma),(polarization,pol),(gamma_beam_rate,gamma_rate) EQUIVALENCE (pion_or_muon,pion),(rho_background,rho) EQUIVALENCE (num_events,nevent),(num_test_events,ntest),(rand_iseed,iseed) real*8 Egamma,Wpipi,thetapipi, & kp(2,3),k(2,3),klab(2,3),kproton(3), & thetapi,phipi,dsigma,dsigmamax,dsigmatest,dsigma_prim,dsigma_rho, & costheta,costheta_min,Emin,Emax,beta,kpicm,Epicm,phipipi, & pi,mtgt,mpi,energy1,energy2,pipimass,sigma_sum,cross_section,ztgt, & tgtlen,radlen,Atgt,gamma_rate,rate,q2,theta_lab(2),k_lab(2),theta_max,deltaE,delta_theta, & array(200,20,2),norm,pol,t,vx,vy,vz,th1,th2,t_test, & q2_array(0:100),theta_array(0:100),q2_max,xmin,xmax,x,costhetapipi,theta1,theta2,E1,E2,theta_cut, & theta,c_den,a_den,rho0,hbarc,norm_vec(2,3),drift,xhit(2),yhit(2),deltax,deltay, & scale,delta_sum,cosphi logical pion,rho c c c Make sure that you, c A. select the target c B. select pion or muon final state c C. turn rho on (true) or off (false) c C. select appropriate rho photoproduction cross section c D. select rho width c E. change iseed c double precision xlo,xhi, ZBQLUAB,ZBQLNOR c integer iseed,ievent,nevent,ntest,j,jj c ! data iseed /3457/ ! data nevent,ntest /100000,1000000/ ! data Egamma,pol /5.5,1./ !beam energy and polarization data Emin,Emax /.280,.500/ !min and max pipi cm energy data theta_max,q2_max /2.0,2.0/ !maximum pi0 angle in degrees, Q2 cut, if no cut make >1 data hbarc /.197327/ data xlo,xhi /0.,1./ c A. Select the target c data ztgt,Atgt,tgtlen,radlen /1.,1., .0338,63.04/ !based on 30 cm target c data ztgt,Atgt,tgtlen,radlen /6.,12. ,.05, 42.70/ c data ztgt,Atgt,tgtlen,radlen,c_den,a_den /14.,28., .05, 21.82, 3.14, .537/ !units of c_den and a_den are fm c data ztgt,Atgt,tgtlen,radlen,c_den,a_den /20.,40., .05, 16.12, 3.51, .563/ !target length is in units of RL ! data ztgt,Atgt,tgtlen,radlen,c_den,a_den /26.,56., .05, 13.84, 3.971, .5935/ !RL is g/cm^2 c data ztgt,Atgt,tgtlen,radlen,c_den,a_den /50.,116.,.05, 8.82, 5.416, .552/ c data ztgt,Atgt,tgtlen,radlen /82.,208.,.05,6.37/ !90% cut is at 0.6 degree ! data gamma_rate /1.e7/ c B. Select pion or muon ! data pion /.true./ c C. Select rho0 background on (true) or off (false) ! data rho /.false./ c ! common /rhob/ norm,rho0,c_den,a_den c c BLOCK DATA ZBQLBD01 * * Initializes seed array etc. for random number generator. * The values below have themselves been generated using the * NAG generator. * COMMON /ZBQL0001/ ZBQLIX,B,C DOUBLE PRECISION ZBQLIX(43),B,C INTEGER I DATA (ZBQLIX(I),I=1,43) /8.001441D7,5.5321801D8, +1.69570999D8,2.88589940D8,2.91581871D8,1.03842493D8, +7.9952507D7,3.81202335D8,3.11575334D8,4.02878631D8, +2.49757109D8,1.15192595D8,2.10629619D8,3.99952890D8, +4.12280521D8,1.33873288D8,7.1345525D7,2.23467704D8, +2.82934796D8,9.9756750D7,1.68564303D8,2.86817366D8, +1.14310713D8,3.47045253D8,9.3762426D7 ,1.09670477D8, +3.20029657D8,3.26369301D8,9.441177D6,3.53244738D8, +2.44771580D8,1.59804337D8,2.07319904D8,3.37342907D8, +3.75423178D8,7.0893571D7 ,4.26059785D8,3.95854390D8, +2.0081010D7,5.9250059D7,1.62176640D8,3.20429173D8, +2.63576576D8/ DATA B / 4.294967291D9 / DATA C / 0.0D0 / c pi=acos(-1.) c c_den=c_den/hbarc a_den=a_den/hbarc c if (pion) then mpi=0.139570 !do pion case else mpi=0.105658 !do muon case end if c call ZBQLINI(ISEED) if (rho) call init_rho(norm) call density_init(rho0,ztgt,c_den,a_den) !for initializing the nuclear form factor (rho0 is density const.) c do i=0,100 q2_array(i)=0. theta_array(i)=0. end do c costheta_min=cos(theta_max*pi/180.) c xmax=log(1.-costheta_min) !change of variable to get better resolution near cos theta =1 xmin=log(1.-(1.-.5*( .01*Emin**2/(2.*Egamma**2) )**2)) !integrate down to 1/100 of the Primakoff peak angle c mtgt=Atgt*.931494 c c Find the maximum cross section over the allowed phase space c dsigmamax=0. sigma_sum=0. do i=1,ntest Wpipi=Emin+(Emax-Emin)*ZBQLUAB(xlo,xhi) c costheta=costheta_min+(1.-costheta_min)*ZBQLUAB(xlo,xhi) x=xmin+(xmax-xmin)*ZBQLUAB(xlo,xhi) costhetapipi=1.-exp(x) c thetapipi=acos(costhetapipi) phipipi=2.*pi*ZBQLUAB(xlo,xhi) c costheta=-1.+2.*ZBQLUAB(xlo,xhi) thetapi=acos(costheta) phipi= 2.*pi*ZBQLUAB(xlo,xhi) c dsigma= dsigma_prim & (Egamma,Wpipi,thetapipi,thetapi,phipipi,phipi,kpicm,beta,Epicm,q2,ztgt,mtgt,pol,pion,mpi) if (rho) dsigma=dsigma+ & dsigma_rho( Egamma,Wpipi,thetapipi,thetapi,phipipi,phipi,kpicm,beta,Epicm,t,ztgt,mtgt,pol) c dsigma = dsigma*(1.-costhetapipi) if (q2.gt.q2_max) dsigma=0. !Q2 cut, acceptance cut c if (dsigma.gt.dsigmamax) dsigmamax=dsigma sigma_sum=sigma_sum+dsigma end do c cross_section=sigma_sum/float(ntest)*4.*pi*(Emax-Emin)*2.*pi*(1.-costheta_min) cross_section=sigma_sum/float(ntest)*4.*pi*(Emax-Emin)*2.*pi*(xmax-xmin) rate=cross_section*1.e-9*100.*1.e-26*tgtlen*radlen/Atgt*6.022e23*gamma_rate*3600. print *, ' maximum cross scection nb/sr^2/GeV', dsigmamax print *, ' total cross section nb', cross_section print *, ' rate per hour', rate c open(unit=2,status='new',file='2pi.stdhep') c c Start event generation c do ievent=1,nevent 100 continue Wpipi=Emin+(Emax-Emin)*ZBQLUAB(xlo,xhi) c costheta=costheta_min+(1.-costheta_min)*ZBQLUAB(xlo,xhi) x=xmin+(xmax-xmin)*ZBQLUAB(xlo,xhi) costhetapipi=1.-exp(x) c thetapipi=acos(costhetapipi) phipipi=2.*pi*ZBQLUAB(xlo,xhi) c costheta=-1.+2.*ZBQLUAB(xlo,xhi) thetapi=acos(costheta) phipi=2.*pi*ZBQLUAB(xlo,xhi) c dsigma= dsigma_prim & (Egamma,Wpipi,thetapipi,thetapi,phipipi,phipi,kpicm,beta,Epicm,q2,ztgt,mtgt,pol,pion,mpi) if (rho) dsigma=dsigma + & dsigma_rho( Egamma,Wpipi,thetapipi,thetapi,phipipi,phipi,kpicm,beta,Epicm,t,ztgt,mtgt,pol) dsigma = dsigma*(1.-costhetapipi) if (q2.gt.q2_max) dsigma=0. !Q2 cut c dsigmatest=dsigmamax*ZBQLUAB(xlo,xhi) if (dsigma.lt.dsigmatest) go to 100 c c Find lab frame momenta. Boost in direction of CM frame kp(1,3)=( kpicm*cos(thetapi)+beta*Epicm)/sqrt(1.-beta**2) kp(2,3)=(-kpicm*cos(thetapi)+beta*Epicm)/sqrt(1.-beta**2) kp(1,1)=kpicm*sin(thetapi)*cos(phipi) kp(1,2)=kpicm*sin(thetapi)*sin(phipi) kp(2,1)=-kp(1,1) kp(2,2)=-kp(1,2) c First rotate around x-axis c find direction vector for pipi system in lab frame vz=cos(thetapipi) vx=sin(thetapipi)*cos(phipipi) vy=sin(thetapipi)*sin(phipipi) c do rotation around x axis th1=asin(vy) do i=1,2 k(i,1)=kp(i,1) k(i,2)=kp(i,2)*cos(th1)+kp(i,3)*sin(th1) k(i,3)=kp(i,3)*cos(th1)-kp(i,2)*sin(th1) end do c Next, rotate around y axis th2=asin(vx/sqrt(vx**2+vz**2)) do i=1,2 klab(i,2)=k(i,2) klab(i,1)=k(i,1)*COS(th2)+k(i,3)*sin(th2) klab(i,3)=k(i,3)*cos(th2)-k(i,1)*sin(th2) end do c energy1=sqrt(klab(1,1)**2+klab(1,2)**2+klab(1,3)**2+mpi**2) energy2=sqrt(klab(2,1)**2+klab(2,2)**2+klab(2,3)**2+mpi**2) t=2.*mtgt*(energy1+energy2-Egamma) t_test=(Egamma-energy1-energy2)**2-(klab(1,1)+klab(2,1))**2-(klab(1,2)+klab(2,2))**2 & -(Egamma-klab(1,3)-klab(2,3))**2 c i=int(1.e4*q2) j=int(thetapipi*180./pi*50.) if (i.gt.100) i=100 if (j.gt.100) j=100 ! max of 2 degrees q2_array(i)=q2_array(i)+1. theta_array(j)=theta_array(j)+1. theta1=atan(sqrt(klab(1,1)**2+klab(1,2)**2)/klab(1,3))*180./pi theta2=atan(sqrt(klab(2,1)**2+klab(2,2)**2)/klab(2,3))*180./pi E1=sqrt(klab(1,1)**2+klab(1,2)**2+klab(1,3)**2+mpi**2) E2=sqrt(klab(2,1)**2+klab(2,2)**2+klab(2,3)**2+mpi**2) c c theta_cut=0.8 c drift=6.2 !meters to FCAL c if (theta1.gt.theta_cut.and.theta2.gt.theta_cut) then c write(2,10) E1,E2 c c This analysis is for finding separation of hits on the TOF c do j=1,2 c do jj=1,3 c norm_vec(j,jj)=klab(j,jj)/sqrt(klab(j,1)**2+klab(j,2)**2+klab(j,3)**2) c end do c scale=drift/norm_vec(j,3) c xhit(j)=scale*norm_vec(j,1) c yhit(j)=scale*norm_vec(j,2) c end do c deltax=(xhit(1)-xhit(2))*100. c deltay=(yhit(1)-yhit(2))*100. c delta_sum=sqrt(deltax**2+deltay**2) c cosphi= ( norm_vec(1,1)*norm_vec(2,1) + norm_vec(1,2)*norm_vec(2,2) ) c & /sqrt(norm_vec(1,1)**2+norm_vec(1,2)**2)/sqrt(norm_vec(2,1)**2+norm_vec(2,2)**2) c write(2,10) cosphi,delta_sum c c end if c write(2,*) '2 ',ievent,' 2' write(2,*) '1 pi+ ',mpi,' 1 ',klab(1,1),klab(1,2),klab(1,3),E1 write(2,*) '2 pi- ',mpi,' -1 ',klab(2,1),klab(2,2),klab(2,3),E2 ! write(2,10) klab(1,1),klab(1,2),klab(1,3),klab(2,1),klab(2,2),klab(2,3) 10 format(2x,10(f13.5,1x)) end do close(2) c do i=0, 100 q2=float(i)*1.e-4 thetapi=float(i)/50. print *, q2, q2_array(i), thetapi, theta_array(i) end do c c stop return end c real*8 function dsigma_prim & (Egamma,Wpipi,thetapipi,thetapi,phipipi,phipi,kpicm,beta,Epicm,q2,ztgt,mtgt,pol,pion,mpi) implicit none real*8 alphaem,ztgt,pi,dsigma_pipi,Wpipi,thetapi,phipipi,phipi,alpha,beta,Egamma,Q2,thetapipi,kpicm, & a,b,c,kpipi,Epicm,mpi,mtgt,FF2,Epipi,pol,dNdE,gamma logical pion data alphaem /7.297e-3/ c pi=acos(-1.) alpha=Egamma*Mtgt+Wpipi**2/2. beta=Egamma*cos(thetapipi) a=(Egamma+Mtgt)**2-beta**2 b=-2.*alpha*beta c=(Egamma+Mtgt)**2*Wpipi**2-alpha**2 kpipi=(-b+sqrt(b**2-4.*a*c))/(2.*a) !3-momentum of pi-pi system Epipi=sqrt(kpipi**2+Wpipi**2) c Q2=-2.*mtgt*(Epipi-Egamma) !this comes from expression for t Q2=(kpipi**2+Egamma**2-2.*kpipi*Egamma*cos(thetapipi)) & -(Egamma-sqrt(Wpipi**2+kpipi**2))**2 !Q2 transfer to nucleus c kpicm=sqrt(Wpipi**2/4.-mpi**2) !pi momentum in CM frame Epicm=Wpipi/2. !pi energy in CM frame beta=kpipi/Epipi !lab velocity of CM frame c c dsigma_prim=2.*alphaem*Ztgt**2/pi**2*dsigma_pipi(Wpipi,thetapi,phipi,pion,mpi)*beta**2*Egamma**4 & /(Wpipi*Q2**2)*sin(thetapipi)**2*FF2(Q2,Ztgt) !units of nb/sr^2/GeV if (pion) then dsigma_prim=dsigma_prim*(1.+pol*cos(2.*phipipi)) !pipi , units of nb/sr^2/GeV else dsigma_prim=dsigma_prim*(1.-pol*cos(2.*phipipi)) !mumu, units of nb/sr^2/GeV end if c c print *, ' dsigma_prim', dsigma_prim, ' Q2', q2, ' Epipi', Epipi return end c real*8 function dsigma_rho(Egamma,Wpipi,thetapipi,thetapi,phipipi,phipi,kpicm,beta,Epicm,t,ztgt,mtgt,pol) implicit none #include "cppvars.h" EQUIVALENCE (rho_norm,norm),(rho_0,rho0),(c_density,c_den),(a_density,a_den),(sigma_rho_0,sigma_0),(A_t,At) real*8 alphaem,ztgt,pi,dsigma_pipi,Wpipi,thetapi,phipipi,phipi,alpha,beta,Egamma,thetapipi,kpicm, & a,b,c,kpipi,Epicm,mpi,mtgt,FF2,Epipi,t,domega_dt,rho,sigma_0,At,norm,pol, & rho0,c_den,a_den data alphaem,pi,mpi /7.297e-3,3.14159,.139570/ c c C. Select rho photoproduction cross section c data sigma_0,At /109.e3,6.0/ !Ballam on proton nb/GeV^2, Gev^-2 c data sigma_0,At /130.e3,7.9/ !ABBHHM on proton nb/GeV^2 ! data sigma_0,At /647.e6,290./ c data sigma_0,At /37.9e6,100./ !Alvensleben on silicon, extrapolated c data sigma_0,At /67.0e6,102./ !Alvensleben on 40Ca, extrapolated , should be 102 c data sigma_0,At /312.e6,196./ !Alvensleben on 116Sn, extrapolated , ! common /rhob/ norm,rho0,c_den,a_den c alpha=Egamma*Mtgt+Wpipi**2/2. beta=Egamma*cos(thetapipi) a=(Egamma+Mtgt)**2-beta**2 b=-2.*alpha*beta c=(Egamma+Mtgt)**2*Wpipi**2-alpha**2 kpipi=(-b+sqrt(b**2-4.*a*c))/(2.*a) !lab 3-momentum of pi-pi system Epipi=sqrt(kpipi**2+Wpipi**2) !lab energy of pi-pi system t=2.*mtgt*(Epipi-Egamma) c kpicm=sqrt(Wpipi**2/4.-mpi**2) !pi momentum in CM frame Epicm=Wpipi/2. !pi energy in CM frame beta=kpipi/Epipi !lab velocity of CM frame c domega_dt=(1./Egamma+1./mtgt)/(2.*kpipi)-( (t-Wpipi**2)/(2.*Egamma)+Epipi )/kpipi**2*(Epipi/(kpipi*2.*mtgt)) domega_dt=domega_dt*2.*pi dsigma_rho=sigma_0*rho(wpipi)*exp(At*t)/domega_dt/norm* & 3./(8.*pi)*sin(thetapi)**2*(1.+pol*cos(2.*phipi)) !units of nb/sr^2/GeV, think this is right! return end c real*8 function dsigma_pipi(Wpipi,thetapi,phipi,pion,mpi) !units of nb/sr implicit none real*8 Wpipi,thetapi,phipi,mpi,pi,muon_xsctn logical pion c pi=acos(-1.) c if (pion) then if (Wpipi.lt..300) then dsigma_pipi=300./(0.6*4.*pi)*(Wpipi-2.*mpi)/(.300-2.*mpi) !linear rise, isotropic, CB data only 60% coverage else dsigma_pipi=300./(0.6*4.*pi)*(.300/Wpipi)**1.29 !power fall off, isotropic end if else !then muon dsigma_pipi=muon_xsctn(Wpipi,thetapi,phipi,mpi) end if return end c real*8 function FF2(Q2,Z) implicit none #include "cppvars.h" EQUIVALENCE (rho_norm,norm),(rho_0,rho0),(c_density,c),(a_density,a_den) real*8 Q2, q02,rms,hbarc,Q,gamma,r(12),A(12),FF,Z,rho0,c,a_den,pi,norm integer*4 i c data q02 /0.71/ data rms,hbarc /5.51,.197/ ! common /rhob/ norm,rho0,c,a_den c data r(1),A(1) /0.1,0.003845/ data r(2),A(2) /0.7,0.009724/ data r(3),A(3) /1.6,0.033093/ data r(4),A(4) /2.1,0.000120/ data r(5),A(5) /2.7,0.083107/ data r(6),A(6) /3.5,0.080869/ data r(7),A(7) /4.2,0.139957/ data r(8),A(8) /5.1,0.260892/ data r(9),A(9) /6.0,0.336013/ data r(10),A(10) /6.6,0.033637/ data r(11),A(11) /7.6,0.018729/ data r(12),A(12) /8.7,0.000020/ data gamma /1.388/ c c Select the Primakoff FF^2 pi=acos(-1.) Q=sqrt(Q2) c if (z.eq.1.) then !proton FF2=(1./(1.+Q2/q02)**2)**2 c else if (z.eq.82.) then !lead FF=0. do i=1,12 FF=FF + A(i)*(gamma**2*cos(Q*r(i)/hbarc)+2.*r(i)*hbarc/Q*sin(Q*r(i)/hbarc))/(gamma**2+2.*r(i)**2)* & exp(-Q2/4.*gamma**2/hbarc**2) end do FF2=FF**2 c else !for everything else use 2-parameter fermi model FF=4.*pi**2*rho0*a_den**3/( (q*a_den)**2*(sinh(pi*q*a_den))**2 )* & ( pi*q*a_den*cosh(pi*q*a_den)*sin(q*c)-q*c*cos(q*c)*sinh(pi*q*a_den) ) do i=1,10 FF=FF+8.*pi*rho0*a_den**3*(-1.)**(i-1)*float(i)*exp(-float(i)*c/a_den)/(float(i)**2+(q*a_den)**2)**2 end do FF2=FF**2 c end if c return end c subroutine density_init(rho0,ztgt,c,a_den) implicit none real*8 pi,c,a_den,rho0,w,ztgt integer*4 i c rho0=0. if (ztgt.eq.82..or.ztgt.eq.6..or.ztgt.eq.1.) return pi=acos(-1.) w=4.*pi*c/3.*( (pi*a_den)**2+c**2 ) do i=1,10 w=w+8.*pi*a_den**3*(-1.)**(i-1)*exp(-float(i)*c/a_den)/float(i)**3 end do rho0=1./w return end c c real*8 function rho(mpipi) implicit none #include "cppvars.h" EQUIVALENCE (rho_width_model, rho_model) real*8 mpipi,gamma0,gamma,mrho,kpicm,kcm,norm,mpi integer rho_model data gamma0,mrho,mpi /0.145,0.770,.139570/ c c D. Select rho width c kpicm=sqrt(mpipi**2/4.-mpi**2) !pi momentum in CM frame kcm =sqrt(mrho**2/4.-mpi**2) !pi momentum in CM frame at rho peak c Zeus parameterization #1 if ( rho_model.eq.1 ) then gamma=gamma0*(kpicm/kcm)**3*mrho/mpipi rho=mpipi*mrho*gamma/( (mpipi**2-mrho**2)**2+mrho**2*gamma**2 ) elseif ( rho_model.eq.2 ) then c Bulos (nuclear targets) gamma= gamma0*mrho/mpipi*( (mpipi**2-4.*mpi**2)/(mrho**2-4.*mpi**2) )**1.5 rho=mpipi*gamma/((mrho**2-mpipi**2)**2 + mrho**2*gamma**2 ) elseif ( rho_model.eq.3 ) then c Alvensleben (nuclear targets) gamma= gamma0*mrho/mpipi*( (mpipi**2-4.*mpi**2)/(mrho**2-4.*mpi**2) )**1.5 rho=mrho*gamma/((mrho**2-mpipi**2)**2 + mrho**2*gamma**2 ) else print*, 'Invalid value for rho_width_model: ',rho_width_model,' should be 1,2, or 3' stop -1 endif c return end c subroutine init_rho(norm) implicit none real*8 norm,mpi,m_min,m_max,sigma_sum,mpipi,rho double precision ZBQLUAB,zlo,zhi integer*4 i,nevent c data m_min,m_max,mpi,zlo,zhi /2.,12.,.139570,0.,1./ data nevent /10000/ c m_min=m_min*mpi m_max=m_max*mpi sigma_sum=0. do i=1,nevent mpipi=m_min+(m_max-m_min)*ZBQLUAB(zlo,zhi) sigma_sum=sigma_sum+rho(mpipi) end do norm=sigma_sum/float(nevent)*(m_max-m_min) return end c real*8 function muon_xsctn(M_cm,theta,phi,M_mu) implicit none real*8 E_plus,M_mu,P_plus,beta,K_cm,K1x,K1y,E1,E2,cos_theta_lab,theta_lab, & jacobian,cos_opening,K_mu_cm,plus_minus_to_gamma_gamma, & M_cm,theta_cm,BD,cos_cm,sin_cm,theta,phi c cos_cm=cos(theta) sin_cm=sqrt(1.-cos_cm**2) E_plus=(M_cm**2-2.*M_mu**2)/(2.*M_mu) P_plus=sqrt(E_plus**2-M_mu**2) beta=P_plus/(E_plus+M_mu) K_cm=M_cm/2. K_mu_cm=sqrt(K_cm**2-M_mu**2) K1x=(K_cm*cos_cm+beta*K_cm)/sqrt(1.-beta**2) K1y=K_cm*sin_cm E1=sqrt(K1x**2+K1y**2) E2=E_plus+M_mu-E1 cos_theta_lab=K1x/sqrt(K1x**2+K1y**2) theta_lab=acos(cos_theta_lab) jacobian=cos(theta_lab)**2*sqrt(1.-beta**2)* & ( cos_cm/(cos_cm+beta)+sin_cm**2/(cos_cm+beta)**2 ) if (sin_cm.ne.0.) then jacobian=jacobian*sin(theta_lab)/sin_cm end if cos_opening=1.-M_cm**2/(2.*E1*E2) plus_minus_to_gamma_gamma=jacobian*BD(M_mu,E_plus,P_plus,theta_lab,cos_opening) muon_xsctn=2.*plus_minus_to_gamma_gamma*(K_mu_cm/K_cm)**2 return end c c real*8 function BD(M_mu,E_plus,P_plus,theta_lab,cos_opening) !units of nb/sr implicit none real*8 alpha,M_mu,E_plus,P_plus,theta_lab,cos_opening,hbarc data alpha, hbarc /7.297e-3,.197/ BD=alpha**2*(M_mu+E_plus)/(8.*P_plus*(M_mu+E_plus-P_plus*cos(theta_lab)))**2 & *( 4.*(E_plus-P_plus*cos(theta_lab))/M_mu+4.*M_mu/(E_plus-P_plus*cos(theta_lab))+8. & -4.*(1.+cos_opening**2) )*hbarc**2/100.*1.e9 return end ******************************************************************* ******** FILE: randgen.f *********** ******** AUTHORS: Richard Chandler *********** ******** (richard@stats.ucl.ac.uk) *********** ******** Paul Northrop *********** ******** (northrop@stats.ox.ac.uk) *********** ******** LAST MODIFIED: 26/8/03 *********** ******** See file randgen.txt for details *********** ******************************************************************* ****************************************************************** ****************************************************************** SUBROUTINE ZBQLINI(SEED) ****************************************************************** * To initialize the random number generator - either * repeatably or nonrepeatably. Need double precision * variables because integer storage can't handle the * numbers involved ****************************************************************** * ARGUMENTS * ========= * SEED (integer, input). User-input number which generates * elements of the array ZBQLIX, which is subsequently used * in the random number generation algorithm. If SEED=0, * the array is seeded using the system clock if the * FORTRAN implementation allows it. ****************************************************************** * PARAMETERS * ========== * LFLNO (integer). Number of lowest file handle to try when * opening a temporary file to copy the system clock into. * Default is 80 to keep out of the way of any existing * open files (although the program keeps searching till * it finds an available handle). If this causes problems, * (which will only happen if handles 80 through 99 are * already in use), decrease the default value. ****************************************************************** INTEGER LFLNO PARAMETER (LFLNO=80) ****************************************************************** * VARIABLES * ========= * SEED See above * ZBQLIX Seed array for the random number generator. Defined * in ZBQLBD01 * B,C Used in congruential initialisation of ZBQLIX * SS,MM,} System clock secs, mins, hours and days * HH,DD } * FILNO File handle used for temporary file * INIT Indicates whether generator has already been initialised * INTEGER SEED,SS,MM,HH,DD,FILNO,I INTEGER INIT DOUBLE PRECISION ZBQLIX(43),B,C DOUBLE PRECISION TMPVAR1,DSS,DMM,DHH,DDD COMMON /ZBQL0001/ ZBQLIX,B,C SAVE INIT * * Ensure we don't call this more than once in a program * IF (INIT.GE.1) THEN IF(INIT.EQ.1) THEN WRITE(*,1) INIT = 2 ENDIF RETURN ELSE INIT = 1 ENDIF * * If SEED = 0, cat the contents of the clock into a file * and transform to obtain ZQBLIX(1), then use a congr. * algorithm to set remaining elements. Otherwise take * specified value of SEED. * *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> *>>>>>>> NB FOR SYSTEMS WHICH DO NOT SUPPORT THE >>>>>>> *>>>>>>> (NON-STANDARD) 'CALL SYSTEM' COMMAND, >>>>>>> *>>>>>>> THIS WILL NOT WORK, AND THE FIRST CLAUSE >>>>>>> *>>>>>>> OF THE FOLLOWING IF BLOCK SHOULD BE >>>>>>> *>>>>>>> COMMENTED OUT. >>>>>>> *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> IF (SEED.EQ.0) THEN *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> *>>>>>>> COMMENT OUT FROM HERE IF YOU DON'T HAVE >>>>>>> *>>>>>>> 'CALL SYSTEM' CAPABILITY ... >>>>>>> *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> CALL SYSTEM(' date +%S%M%H%j > zbql1234.tmp') * * Try all file numbers for LFLNO to 999 * FILNO = LFLNO 10 OPEN(FILNO,FILE='zbql1234.tmp',ERR=11) GOTO 12 11 FILNO = FILNO + 1 IF (FILNO.GT.999) THEN WRITE(*,2) RETURN ENDIF GOTO 10 12 READ(FILNO,'(3(I2),I3)') SS,MM,HH,DD CLOSE(FILNO) CALL SYSTEM('rm zbql1234.tmp') DSS = DINT((DBLE(SS)/6.0D1) * B) DMM = DINT((DBLE(MM)/6.0D1) * B) DHH = DINT((DBLE(HH)/2.4D1) * B) DDD = DINT((DBLE(DD)/3.65D2) * B) TMPVAR1 = DMOD(DSS+DMM+DHH+DDD,B) *<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< *<<<<<<<< ... TO HERE (END OF COMMENTING OUT FOR <<<<<<< *<<<<<<<< USERS WITHOUT 'CALL SYSTEM' CAPABILITY <<<<<<< *<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ELSE TMPVAR1 = DMOD(DBLE(SEED),B) ENDIF ZBQLIX(1) = TMPVAR1 DO 100 I = 2,43 TMPVAR1 = ZBQLIX(I-1)*3.0269D4 TMPVAR1 = DMOD(TMPVAR1,B) ZBQLIX(I) = TMPVAR1 100 CONTINUE 1 FORMAT(//5X,'****WARNING**** You have called routine ZBQLINI ', +'more than',/5X,'once. I''m ignoring any subsequent calls.',//) 2 FORMAT(//5X,'**** ERROR **** In routine ZBQLINI, I couldn''t', +' find an',/5X, +'available file number. To rectify the problem, decrease the ', +'value of',/5X, +'the parameter LFLNO at the start of this routine (in file ', +'randgen.f)',/5X, +'and recompile. Any number less than 100 should work.') END ****************************************************************** FUNCTION ZBQLU01(DUMMY) * * Returns a uniform random number between 0 & 1, using * a Marsaglia-Zaman type subtract-with-borrow generator. * Uses double precision, rather than integer, arithmetic * throughout because MZ's integer constants overflow * 32-bit integer storage (which goes from -2^31 to 2^31). * Ideally, we would explicitly truncate all integer * quantities at each stage to ensure that the double * precision representations do not accumulate approximation * error; however, on some machines the use of DNINT to * accomplish this is *seriously* slow (run-time increased * by a factor of about 3). This double precision version * has been tested against an integer implementation that * uses long integers (non-standard and, again, slow) - * the output was identical up to the 16th decimal place * after 10^10 calls, so we're probably OK ... * DOUBLE PRECISION ZBQLU01,DUMMY,B,C,ZBQLIX(43),X,B2,BINV INTEGER CURPOS,ID22,ID43 COMMON /ZBQL0001/ ZBQLIX,B,C SAVE /ZBQL0001/ SAVE CURPOS,ID22,ID43 DATA CURPOS,ID22,ID43 /1,22,43/ B2 = B BINV = 1.0D0/B 5 X = ZBQLIX(ID22) - ZBQLIX(ID43) - C IF (X.LT.0.0D0) THEN X = X + B C = 1.0D0 ELSE C = 0.0D0 ENDIF ZBQLIX(ID43) = X * * Update array pointers. Do explicit check for bounds of each to * avoid expense of modular arithmetic. If one of them is 0 the others * won't be * CURPOS = CURPOS - 1 ID22 = ID22 - 1 ID43 = ID43 - 1 IF (CURPOS.EQ.0) THEN CURPOS=43 ELSEIF (ID22.EQ.0) THEN ID22 = 43 ELSEIF (ID43.EQ.0) THEN ID43 = 43 ENDIF * * The integer arithmetic there can yield X=0, which can cause * problems in subsequent routines (e.g. ZBQLEXP). The problem * is simply that X is discrete whereas U is supposed to * be continuous - hence if X is 0, go back and generate another * X and return X/B^2 (etc.), which will be uniform on (0,1/B). * IF (X.LT.BINV) THEN B2 = B2*B GOTO 5 ENDIF ZBQLU01 = X/B2 END ****************************************************************** FUNCTION ZBQLUAB(A,B) * * Returns a random number uniformly distributed on (A,B) * DOUBLE PRECISION A,B,ZBQLU01,ZBQLUAB * * Even if A > B, this will work as B-A will then be -ve * IF (A.NE.B) THEN ZBQLUAB = A + ( (B-A)*ZBQLU01(0.0D0) ) ELSE ZBQLUAB = A WRITE(*,1) ENDIF 1 FORMAT(/5X,'****WARNING**** (function ZBQLUAB) Upper and ', +'lower limits on uniform',/5X,'distribution are identical',/) END ****************************************************************** FUNCTION ZBQLEXP(MU) * * Returns a random number exponentially distributed with * mean MU * DOUBLE PRECISION MU,ZBQLEXP,ZBQLU01 ZBQLEXP = 0.0D0 IF (MU.LT.0.0D0) THEN WRITE(*,1) RETURN ENDIF ZBQLEXP = -DLOG(ZBQLU01(0.0D0))*MU 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLEXP',/) END ****************************************************************** FUNCTION ZBQLNOR(MU,SIGMA) * * Returns a random number Normally distributed with mean * MU and standard deviation |SIGMA|, using the Box-Muller * algorithm * DOUBLE PRECISION THETA,R,ZBQLNOR,ZBQLU01,PI,MU,SIGMA DOUBLE PRECISION SPARE INTEGER STATUS SAVE STATUS,SPARE,PI DATA STATUS /-1/ IF (STATUS.EQ.-1) PI = 4.0D0*DATAN(1.0D0) IF (STATUS.LE.0) THEN THETA = 2.0D0*PI*ZBQLU01(0.0D0) R = DSQRT( -2.0D0*DLOG(ZBQLU01(0.0D0)) ) ZBQLNOR = (R*DCOS(THETA)) SPARE = (R*DSIN(THETA)) STATUS = 1 ELSE ZBQLNOR = SPARE STATUS = 0 ENDIF ZBQLNOR = MU + (SIGMA*ZBQLNOR) END ****************************************************************** FUNCTION ZBQLBIN(N,P) * * Returns a random number binomially distributed (N,P) * DOUBLE PRECISION P,ZBQLBET1 DOUBLE PRECISION PP,PPP,G,Y,TINY INTEGER N,ZBQLBIN,ZBQLGEO,IZ,NN TINY = 1.0D-8 ZBQLBIN = 0 IF (.NOT.( (P.GE.0.0D0).AND.(P.LE.1.0D0) )) THEN WRITE(*,1) RETURN ELSEIF(N.LE.0) THEN WRITE(*,1) RETURN ENDIF * * First step: if NP > 10, say, things will be expensive, and * we can get into the right ballpark by guessing a value for * ZBQLBIN (IZ, say), and simulating Y from a Beta distribution * with parameters IZ and NN-IZ+1 (NN starts off equal to N). * If Y is less than PP (which starts off as P) then the IZth order * statistic from NN U(0,1) variates is less than P, and we know * that there are at least IZ successes. In this case we focus on * the remaining (NN-IZ) order statistics and count how many are * less than PP, which is binomial (NN-IZ,(PP-Y)/(1-Y)). * Otherwise, if Y is greater than PP there must be less * than IZ successes, so we can count the number of order statistics * under PP, which is binomial (IZ-1,P/Y). When we've got NN*PP * small enough, we go to the next stage of the algorithm and * generate the final bits directly. * NN = N PP = P 10 IZ = INT(DBLE(NN)*PP) + 1 IF ( (IZ.GT.10).AND.(IZ.LT.NN-10) ) THEN Y = ZBQLBET1(DBLE(IZ),DBLE(NN-IZ+1)) IF (Y.LT.PP) THEN ZBQLBIN = ZBQLBIN + IZ NN = NN - IZ PP = (PP-Y) / (1.0D0-Y) ELSE NN = IZ-1 PP = PP/Y ENDIF GOTO 10 ENDIF * * PP is the probability of the binomial we're currently * simulating from. For the final part, we simulate either number * of failures or number of success, depending which is cheaper. * 20 IF (PP.GT.0.5) THEN PPP = 1.0D0-PP ELSE PPP = PP ENDIF G = 0 IZ = 0 * * ZBQLGEO falls over for miniscule values of PPP, so ignore these * (tiny probability of any successes in this case, anyway) * IF (PPP.GT.TINY) THEN 30 G = G + ZBQLGEO(PPP) IF (G.LE.NN) THEN IZ = IZ + 1 GOTO 30 ENDIF ENDIF IF (PP.GT.0.5) IZ = NN - IZ ZBQLBIN = ZBQLBIN + IZ 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLBIN',/) END ****************************************************************** FUNCTION ZBQLGEO(P) * * Returns a random number geometrically distributed with * parameter P ie. mean 1/P * DOUBLE PRECISION P,ZBQLU01,U,TINY INTEGER ZBQLGEO TINY = 1.0D-12 ZBQLGEO = 0 IF (.NOT.( (P.GE.0.0D0).AND.(P.LE.1.0D0) )) THEN WRITE(*,1) RETURN ENDIF IF (P.GT.0.9D0) THEN 10 ZBQLGEO = ZBQLGEO + 1 U = ZBQLU01(0.0D0) IF (U.GT.P) GOTO 10 ELSE U = ZBQLU01(0.0D0) * * For tiny P, 1-p will be stored inaccurately and log(1-p) may * be zero. In this case approximate log(1-p) by -p * IF (P.GT.TINY) THEN ZBQLGEO = 1 + INT( DLOG(U)/DLOG(1.0D0-P) ) ELSE ZBQLGEO = 1 + INT(-DLOG(U)/P) ENDIF ENDIF 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLGEO',/) END ****************************************************************** FUNCTION ZBQLPOI(MU) * * Returns a random number Poisson distributed with mean MU * DOUBLE PRECISION ZBQLU01,X,Y,MU,PI DOUBLE PRECISION ZBQLLG,ZBQLGAM,MU1,TMP1,TMP2,T INTEGER ZBQLPOI,ZBQLBIN,K,INIT SAVE INIT,PI DATA INIT /0/ IF (INIT.EQ.0) THEN PI = 4.0D0*DATAN(1.0D0) INIT = 1 ENDIF ZBQLPOI = 0 IF (MU.LT.0.0D0) THEN WRITE(*,1) RETURN ENDIF * * For small MU, generate exponentials till their sum exceeds 1 * (equivalently, uniforms till their product falls below e^-MU) * IF (MU.LE.1.0D3) THEN MU1 = MU * * For values of MU less than 1000, use order statistics - the Kth * event in a Poisson process of rate MU has a Gamma distribution * with parameters (MU,K); if it's greater than 1 we know that there * are less than K events in (0,1) (and the exact number is binomial) * and otherwise the remaining number is another Poisson. Choose K so * that we'll get pretty close to 1 in the first go but are unlikely * to overshoot it. * 19 IF (MU1.GT.1.0D1) THEN K = INT(MU1-DSQRT(MU1)) Y = ZBQLGAM(DBLE(K),MU1) IF (Y.GT.1.0D0) THEN ZBQLPOI = ZBQLPOI + ZBQLBIN(K-1,(1.0D0/Y)) RETURN ENDIF ZBQLPOI = ZBQLPOI + K MU1 = MU1 * (1.0D0-Y) GOTO 19 ENDIF Y = DEXP(-MU1) X = 1.0D0 20 X = X*ZBQLU01(0.0D0) IF (X.GT.Y) THEN ZBQLPOI = ZBQLPOI + 1 GOTO 20 ENDIF * * For really huge values of MU, use rejection sampling as in * Press et al (1992) - large numbers mean some accuracy may be * lost, but it caps the execution time. * ELSE TMP1 = DSQRT(2.0D0*MU) TMP2 = ZBQLLG(MU+1.0D0)-(MU*DLOG(MU)) 30 Y = DTAN(PI*ZBQLU01(0.0D0)) ZBQLPOI = INT(MU + (TMP1*Y)) IF (ZBQLPOI.LT.0) GOTO 30 X = DBLE(ZBQLPOI) T = (X*DLOG(MU)-ZBQLLG(X+1.0D0)) + TMP2 IF (DABS(T).LT.1.0D2) THEN T = 0.9D0*(1.0D0+(Y*Y))*DEXP(T) IF (ZBQLU01(0.0D0).GT.T) GOTO 30 ELSE T = DLOG(0.9D0*(1.0D0+(Y*Y))) + T IF (DLOG(ZBQLU01(0.0D0)).GT.T) GOTO 30 ENDIF ENDIF 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLPOI',/) END ****************************************************************** FUNCTION ZBQLGAM(G,H) * * Returns a random number with a gamma distribution with mean * G/H and variance G/(H^2). (ie. shape parameter G & scale * parameter H) * DOUBLE PRECISION C,D,R,ZBQLGAM,ZBQLU01,G,H,A,z1,z2,B1,B2,M DOUBLE PRECISION U1,U2,U,V,TEST,X double precision c1,c2,c3,c4,c5,w ZBQLGAM = 0.0D0 IF ( (G.LE.0.0D0).OR.(H.LT.0.0D0) ) THEN WRITE(*,1) RETURN ENDIF IF (G.LT.1.0D0) THEN 889 u=ZBQLU01(0.0d0) v=ZBQLU01(0.0d0) if (u.gt.exp(1.0d0)/(g+exp(1.0d0))) goto 891 ZBQLGAM=((g+exp(1.0d0))*u/exp(1.0d0))**(1.0d0/g) if (v.gt.exp(-ZBQLGAM)) then goto 889 else goto 892 endif 891 ZBQLGAM=-log((g+exp(1.0d0))*(1.0d0-u)/(g*exp(1.0d0))) if (v.gt.ZBQLGAM**(g-1.0)) goto 889 892 ZBQLGAM=ZBQLGAM/h RETURN ELSEIF (G.LT.2.0D0) THEN M = 0.0D0 elseif (g.gt.10.0d0) then c1=g-1.0d0 c2=(g-1.0d0/(6.0d0*g))/c1 c3=2.0d0/c1 c4=c3+2.0d0 c5=1.0d0/sqrt(g) 777 u=ZBQLU01(0.0d0) v=ZBQLU01(0.0d0) if (g.gt.2.50d0) then u=v+c5*(1.0d0-1.860d0*u) endif if (u.le.0.0d0.or.u.ge.1.0d0) goto 777 w=c2*v/u if (c3*u+w+1.0d0/w.le.c4) goto 778 if (c3*log(u)-log(w)+w.ge.1.0d0) goto 777 778 ZBQLGAM=c1*w/h return ELSE M = -(G-2.0D0) ENDIF R = 0.50D0 a = ((g-1.0d0)/exp(1.0d0))**((g-1.0d0)/(r+1.0d0)) C = (R*(M+G)+1.0D0)/(2.0D0*R) D = M*(R+1.0D0)/R z1 = C-DSQRT(C*C-D) * * On some systems (e.g. g77 0.5.24 on Linux 2.4.24), C-DSQRT(C*C) * is not exactly zero - this needs trapping if negative. * IF ((Z1-M.LT.0.0D0).AND.(Z1-M.GT.-1.0D-12)) Z1 = M z2 = C+DSQRT(C*C-D) B1=(z1*(z1-M)**(R*(G-1.0D0)/(R+1.0D0)))*DEXP(-R*(z1-M)/(R+1.0D0)) B2=(z2*(z2-M)**(R*(G-1.0D0)/(R+1.0D0)))*DEXP(-R*(z2-M)/(R+1.0D0)) 50 U1=ZBQLU01(0.0D0) U2=ZBQLU01(0.0D0) U=A*U1 V=B1+(B2-B1)*U2 X=V/(U**R) IF (X.LE.M) GOTO 50 TEST = ((X-M)**((G-1)/(R+1)))*EXP(-(X-M)/(R+1.0D0)) IF (U.LE.TEST) THEN ZBQLGAM = (X-M)/H ELSE GOTO 50 ENDIF 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLGAM',/5X, '(both parameters must be positive)',/) END *************************************************************** FUNCTION ZBQLBET1(NU1,NU2) * * Returns a random number, beta distributed with degrees * of freedom NU1 and NU2. * DOUBLE PRECISION NU1,NU2,ZBQLGAM,ZBQLBET1,ZBQLU01,X1,X2 ZBQLBET1 = 0.0D0 IF ( (NU1.LE.0.0).OR.(NU2.LE.0.0) ) THEN WRITE(*,1) RETURN ENDIF * * If parameters are too small, gamma subroutine tends to return zero * as all the probability goes to the origin and we get rounding * errors, even with double precision. In this case, we use Johnk's * method, suitably scaled to avoid rounding errors as much as possible. * IF ( (NU1.LT.0.9D0).AND.(NU2.LT.0.9D0) ) THEN 10 X1 = ZBQLU01(0.0D0) X2 = ZBQLU01(0.0D0) IF ( (X1**(1.0D0/NU1))+(X2**(1.0D0/NU2)).GT.1.0D0) GOTO 10 X1 = (DLOG(X2)/NU2) - (DLOG(X1)/NU1) ZBQLBET1 = (1.0D0 + DEXP(X1))**(-1) IF (ZBQLBET1.GT.1.0D0) GOTO 10 ELSE X1 = ZBQLGAM(NU1,1.0D0) X2 = ZBQLGAM(NU2,1.0D0) ZBQLBET1 = X1/(X1+X2) ENDIF RETURN 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLBET1',/5X, +'(both degrees of freedom must be positive)',/) END *************************************************************** FUNCTION ZBQLWEI(A,B) * * Returns a random number, Weibull distributed with shape parameter * A and location parameter B, i.e. density is * f(x) = ( A/(B**A) ) * x**(A-1) * EXP( -(x/B)**A ) * DOUBLE PRECISION A,B,ZBQLU01,ZBQLWEI,U ZBQLWEI = 0.0D0 IF ( (A.LE.0.0).OR.(B.LE.0.0) ) THEN WRITE(*,1) RETURN ENDIF U = ZBQLU01(0.0D0) ZBQLWEI = B * ( (-DLOG(U))**(1.0D0/A) ) 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLWEI',/5X, +'(both parameters must be positive)',/) END *************************************************************** FUNCTION ZBQLNB(R,P) * * Returns a pseudo-random number according to a Negative * Binomial distribution with parameters (R,P). NB these are * both DOUBLE - it copes with non-integer R as well. The * form of the distribution is *not* the no. of trials to * the Rth success - see documentation for full spec. * DOUBLE PRECISION R,P,ZBQLGAM,Y INTEGER ZBQLNB,ZBQLPOI ZBQLNB = 0 IF ( (R.LE.0.0D0).OR.(P.LE.0.0D0).OR.(P.GE.1.0D0) ) THEN WRITE(*,1) RETURN ENDIF Y = ZBQLGAM(R,1.0D0) Y = Y*P/(1.0D0-P) ZBQLNB = ZBQLPOI(Y) 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLNB') END *************************************************************** FUNCTION ZBQLPAR(A,B) * * Returns a random number, Pareto distributed with parameters * A and B. The density is A*(B**A) / (B+X)**(A+1) for X > 0. * (this is slightly nonstandard - see documentation in * randgen.txt). The algorithm is straightforward - it uses the * inverse CDF method. * DOUBLE PRECISION A,B,ZBQLPAR,ZBQLU01,U ZBQLPAR = 0.0D0 IF ( (A.LE.0.0D0).OR.(B.LE.0.0D0) ) THEN WRITE(*,1) RETURN ENDIF U = ZBQLU01(0.0D0) ZBQLPAR = B * (U**(-1.0D0/A)-1.0D0) 1 FORMAT(/5X,'****ERROR**** Illegal parameter value in ', +' ZBQLPAR',/5X, +'(both parameters must be positive)',/) END *************************************************************** FUNCTION ZBQLLG(X) * * Returns log(G(X)) where G is the Gamma function. The algorithm is * that given in Press et al (1992), Section 6.1, although this * version also allows for arguments less than 1. * DOUBLE PRECISION X,Z,Z2,ZBQLLG,PI,RLN2P,C(0:6),TMP,SUM INTEGER INIT,I SAVE INIT,C,RLN2P,PI DATA INIT /0/ DATA (C(I),I=0,6) / + 1.000000000190015D0,76.18009172947146D0, + -86.50532032941677D0,24.01409824083091D0, + -1.231739572450155D0,0.1208650973866179D-2, + -0.5395239384953D-5/ IF (INIT.EQ.0) THEN PI = 4.0D0*DATAN(1.0D0) RLN2P = 0.5D0*DLOG(2.0D0*PI) INIT = 1 ENDIF * * Compute for x > 1, then use transformation if necessary. Z is * our working argument. * IF (X.GE.1.0D0) THEN Z = X ELSE Z = 2.0D0-X Z2 = 1.0D0-X ENDIF IF (DABS(Z-1.0D0).LT.1.0D-12) THEN ZBQLLG = 0.0D0 RETURN ENDIF TMP = Z + 4.5D0 TMP = ( (Z-0.5D0)*DLOG(TMP) ) - TMP + RLN2P SUM = C(0) DO 50 I=1,6 SUM = SUM + (C(I)/(Z+DBLE(I-1))) 50 CONTINUE ZBQLLG = TMP + DLOG(SUM) * * Transformation required if X<1 * IF (X.LT.1.0D0) THEN TMP = PI*Z2 ZBQLLG = DLOG(TMP/DSIN(TMP)) - ZBQLLG ENDIF END