program pion_kin implicit none real*8 pi real*8 Ebeam,estep,W,Q2,theta_piq,xbj real*8 Eprime,qabs,qpar,qperp,omega,theta_e,theta_pi,Ppi,theta_pq real*8 mp,mn,mtar,mrec,mpi real*8 a,b,c,QA,QB,QC,radical,Epi,Kpi,t real*8 epsilon real*8 betacm,gammacm,ppicm,qcm,ppar,pperp,ppiparcm real*8 costhcm,thcm real*8 jacobian real*8 gammaflux real*8 sig3,xsec integer piplus_flag,i external sig3 parameter(pi=3.141592654) parameter(mp = 0.93827231) parameter(mn = 0.93956563) parameter(mpi = 0.13956995) write(6,*) 'Enter beam energy in GeV' read(5,*) Ebeam c write(6,*) 'Enter xbj (GeV)' c read(5,*) xbj write(6,*) 'Enter W (GeV)' read(5,*) W write(6,*) 'Enter Q2 (GeV^2)' read(5,*) Q2 write(6,*) 'Pion angle relative to q (degrees)' read(5,*) theta_pq theta_pq = pi/180. * theta_pq write(6,*) 'pi+ or pi- (1 = pi+)' read(5,*) piplus_flag c w = sqrt(-Q2 + mp**2 + Q2/xbj) xbj = Q2/(w**2+Q2-mp**2) write(6,*) 'x is',xbj write(6,*) 'W is',W if(piplus_flag.eq.1) then mtar = mp mrec = mn else mtar = mn mrec = mp endif Eprime = Ebeam - (W**2+Q2 - mtar**2)/2./mtar c Eprime = Ebeam - Q2/2./mtar/xbj omega = Ebeam - Eprime theta_e = 2.*asin(sqrt(Q2/4./Eprime/Ebeam)) qpar = Ebeam - Eprime*cos(theta_e) qperp = -Eprime*sin(theta_e) qabs = sqrt(qpar**2+qperp**2) theta_pi = acos(qpar/qabs) + theta_pq a = -1.*qabs*cos(theta_pq) b = qabs**2 c = Ebeam - Eprime + mtar QA = 4.*(a**2 - c**2) QB = 4.*c*(c**2-b+mpi**2-mrec**2) QC = -4.*a**2*mpi**2 - (c**2-b+mpi**2-mrec**2)**2 radical = QB**2 - 4.*QA*QC if(radical.lt.0.) then write(6,*) 'Kinematics are impossible! You lose!' else Epi = (-QB-sqrt(radical))/2./QA endif Ppi = sqrt(Epi**2-mpi**2) kpi = sqrt(Ppi**2 + qabs**2 -2.*Ppi*qabs*cos(theta_pq) ) t = -Q2 + mpi**2 - 2.*(Ebeam-Eprime)*Epi + > 2.*qabs*Ppi*cos(theta_pq) epsilon=1./(1.+2.*(Q2+(Ebeam-Eprime)**2)/ > Q2*(tan(theta_e/2.))**2) betacm = qabs/(omega+mtar) gammacm = (omega+mtar)/w qcm = gammacm*(qabs-omega*betacm) ppar=ppi*cos(theta_pq) pperp=ppi*sin(theta_pq) ppiparcm= gammacm*(ppar-epi*betacm) ppicm = sqrt(ppiparcm**2+pperp**2) costhcm = ppiparcm/ppicm thcm = acos(costhcm) jacobian = mtar*Ppi*qabs/(ppicm*qcm) * > 1./(mtar+omega-Epi/Ppi*qabs*cos(theta_pq)) gammaflux = 1./137.*1./(2.*3.14159**2)*Eprime/Ebeam* > (w**2-mtar**2)/2./mtar/(1-epsilon)/q2 xsec=sig3(w,t,q2,ebeam,eprime,epsilon,qabs,omega, > betacm,gammacm,ppi,ppicm,theta_pq,0.d0) write(6,*) 'Scattered electron energy =',Eprime,' GeV' write(6,*) 'omega (nu, whatever..) =',omega,' GeV' write(6,*) 'q is =',qabs,' GeV' write(6,*) 'Electron (lab) scattering angle=',theta_e*180/pi, > ' degrees' write(6,*) 'Pion (lab) scattering angle =',theta_pi*180/pi, > ' degrees' write(6,*) 'Pion momentum =',Ppi,' GeV' write(6,*) 'Pion energy =',Epi,' GeV' write(6,*) 'Kpi (or Precoil if you like..) =',Kpi,' GeV' write(6,*) 'virtual pi energy =',sqrt(t+kpi**2), > ' GeV' write(6,*) 't =',t,' GeV^2' write(6,*) 'epsilon =',epsilon write(6,*) 'Ppi cm =',ppicm,' GeV' write(6,*) 'q cm =',qcm,' GeV' write(6,*) 'Theta-cm =',thcm*180.0/pi,' deg' write(6,*) 'Jacobian (dW_cm -> dW_lab) =',jacobian write(6,*) 'Gammaflux =',gammaflux write(6,*) 'Xsec =',xsec 2 format(7(f7.3,3x)) 3 format(7(a7,3x)) end function sig3(w,t,q2,ee,eepr,epsilon,pgam,egam,betcm,gamcm, > ppi,ppicm,thpi,phix) ! this function calculates cross-sections d3sigma/dEe'dOmegae'Omegapi in ! units of ub/GeV/sr^2 implicit none real*8 ee,eepr,mp,mpi,q2,t,alpha,pi,epsilon,pgam real*8 ppi,ppicm,pgamcm,a,mn,phix,egam real*8 sigl,sigt,siglt,sigtt,gtpr,sig3,mrho,fpi,fpi07 real*8 gamcm,betcm,w,fpi2,fpi072,r,epi,thpi,cmtolab mpi = 0.13957 mp = 0.93828 mn = 0.93957 alpha = 1./137.036 pi = 3.14156 epi=sqrt(mpi**2+ppi**2) ! parameterization from Brauel and Bebek data. sigl=30.6*exp(-14*abs(t))+1.15*exp(-0.6*abs(t)) sigt=2.82*exp(-2*abs(t)) siglt=-0.14 sigtt=-2.23 ! now scale sigl by Q2*pion_formfactor mrho=0.712 !determined by empirical fit to formfactor if (q2.lt.2.19) then fpi=1./(1.+q2/mrho**2) else fpi=1./(1.+2.19/mrho**2)/q2 endif fpi2=fpi**2 fpi07=1./(1.+0.70/mrho**2) fpi072=fpi07**2 sigl=sigl*(fpi2*q2)/(fpi072*0.7) ! Brauel scaled all his data to W=2.19 GeV, so rescale by 1/(w**2-mp**2)**2 sigl=sigl*(2.19**2-mp**2)**2/(w**2-mp**2)**2 sigt=sigt*(2.19**2-mp**2)**2/(w**2-mp**2)**2 ! outside the range of validity of Brauel's paper, use Bebek's Q2=10 formfactor ! paper (PRD 17 (1978) 1693). if (q2.ge.2. .or. abs(t).ge.0.2) then r=q2/4.2 sigl=(sigl+sigt)/(r+1) sigt=sigl*r siglt=0. sigtt=0. endif c if (sigl .lt. (1.*sigt)) then c sig3 = 1.e-10 c return c endif gtpr=alpha/2./(pi**2)*eepr/ee*(w**2-mp**2)/2./mp/q2/(1.-epsilon) pgamcm = (pgam-betcm*egam)/gamcm/(1.-betcm**2) a = pgamcm*ppicm/pi write(6,*) 'dtdphi->dOmega_cm',a sig3= gtpr*( a*sigt + epsilon*a*sigl > + epsilon*a*cos(2.*phix)*sigtt > + sqrt(.5*epsilon*(1.+epsilon))*2.*a*cos(phix)*siglt ) cmtolab = ppi*ppi*w/ppicm/ > abs( (mp+egam)*ppi - epi*pgam*cos(thpi) ) write(6,*) 'cmtolab is',cmtolab sig3 = sig3*cmtolab return end