1 gaskelld 1.1 	subroutine rho_decay(orig,p_spec,epsilon,success)
  2              
  3              C+_____________________________________________________________________________
  4              !  rho_decay:  This routine takes the electroproduced rho0 and
  5              !   produces the decaying pi+ (or pi-) that we actually detect in the 
  6              !   spectrometer.
  7              !
  8              !   This routine should be called before mc_hms or mc_sos.
  9              !
 10              !   The outgoing pion will be thrown flat in phi and according to
 11              !   sin(theta)**2+epsilon*R*cos(theta)**2 (normalized) in the rho
 12              !   rest frame. 
 13              !
 14              !   Note that this routine depends on R = sigma_L/sigma_T (and epsilon).
 15              C-_____________________________________________________________________________
 16              
 17              	implicit none
 18              
 19              	include 'simulate.inc'
 20              
 21              	type(event):: orig
 22 gaskelld 1.1 C Local declarations.
 23              	real*8 p_spec,ph
 24              	real*8 rph,rth1,rth
 25              	real*8 beta,gamma,dlen
 26              	real*8 ef,pf,pxf,pyf,pzf,pxr,pyr,pzr
 27              	real*8 bx,by,bz,er,pr
 28              
 29              	real*8 th_oop,th_inp,cos_th_inp
 30              	real*8 pzprime
 31              
 32              	real*8 epsilon,R_rho
 33              
 34              	real*8 ctaurho !do not use "real" ctau (i.e. for pion) since
 35              	               !we're doing rho decay here
 36              	real*8 norm,dist
 37              	real*8 grnd
 38              
 39              	logical success
 40              
 41              	parameter (ctaurho=1.31467d-13)
 42              C ============================= Executable Code ================================
 43 gaskelld 1.1 
 44              C Calulate R = sigmaL/sigmaT
 45              C Warning!  Note that if you change the parameterization of R in
 46              C rho-physics, do it here as well for self consistency!
 47              C This parameterization is taken from HERMES data...
 48              
 49              	R_rho = 0.33*(orig%Q2/mrho2)**(0.61)
 50              	
 51              
 52              	ph = orig%p%P
 53              	beta = ph/sqrt(ph**2+Mh2)
 54              	gamma = 1./sqrt(1.-beta*beta)
 55              	dlen=ctaurho*beta*gamma	!1/e decay length (beta*c*tau*gamma)
 56              
 57              C Generate center/mass decay angles and momenta.
 58              	rph = grnd()*2.*pi
 59               100	rth1 = grnd()*2.-1.
 60              	rth = acos(rth1)
 61              
 62              

 63 gaskelld 1.2 	norm = (1.0+2.0*epsilon*R_rho)*grnd()
 64              	dist = sin(rth)**2+2.0*epsilon*R_rho*cos(rth)**2

 65 gaskelld 1.1 	if(dist.lt.norm) goto 100
 66              c       write(6,*) 'now decaying the rho',ctaurho,mh2
 67              	ntup%rhotheta=rth
 68              c	er = 384.65
 69              c	pr = 358.4353
 70              	er = ntup%rhomass/2.0

 71 gaskelld 1.3 	if(er.lt.Mpi) then
 72              	   success=.false.
 73              	   return
 74              	endif

 75 gaskelld 1.1 	pr = sqrt(er**2-Mpi2)
 76              	pxr = pr*sin(rth)*cos(rph)
 77              	pyr = pr*sin(rth)*sin(rph)
 78              	pzr = pr*cos(rth)
 79              
 80              C Boost to Lab frame, calculate new angles and momentum, finish drift
 81              
 82              	bx = -beta * orig%up%x
 83              	by = -beta * orig%up%y
 84              	bz = -beta * orig%up%z
 85              	call loren(gamma,bx,by,bz,er,pxr,pyr,pzr,ef,pxf,pyf,pzf,pf)
 86              
 87              
 88              C DJG check if pion is at least going toward to spectrometer!
 89              C DJG If not, I'm outta here.
 90              C DJG Need to rotate coordinate sytem counterclockwise (clockwise) by angle
 91              C DJG theta_spec about x-axis (points DOWN) for HMS (SOS).
 92              
 93              	if(hadron_arm.eq.1.or.hadron_arm.eq.3) then
 94              	   pzprime = pzf*cos(spec%p%theta)-pyf*sin(spec%p%theta)
 95              	elseif (hadron_arm.eq.2.or.hadron_arm.eq.4) then
 96 gaskelld 1.1 	   pzprime = pzf*cos(spec%p%theta)+pyf*sin(spec%p%theta)
 97              	else
 98              	   write(6,*) 'Unknown spectrometer setup dude!'
 99              	   stop
100              	endif
101              
102              	if(pzprime.lt.0.0) then
103              c	   write(6,*) 'Pion from rho decay moving away from spectrometer'
104              c	   write(6,*) 'I guess this event is a no-go!'
105              	   success = .false.
106              	   return
107              	endif
108              	   
109              C DJG If pion's heading for spectrometer, try to calculate xptar and yptar
110              
111              	th_oop=asin(pxf/pf)
112                      cos_th_inp = pzf/pf/cos(th_oop)
113              	th_inp = acos(cos_th_inp)
114              	if(hadron_arm.eq.1.or.hadron_arm.eq.3) then
115              	   if(pyf.lt.0.0) then
116              	      th_inp = spec%p%theta-th_inp 
117 gaskelld 1.1 	   else
118              	      th_inp = spec%p%theta+th_inp
119              	   endif
120              	elseif(hadron_arm.eq.2.or.hadron_arm.eq.4) then
121              	   if(pyf.gt.0.0) then
122              	      th_inp = th_inp - spec%p%theta
123              	   else
124              	      th_inp = th_inp + spec%p%theta
125              	   endif
126              	else
127              	   write(6,*) 'Rho decay not set up for your spectrometer!'
128              	   stop
129              	endif
130              
131              C DJG This should never happen!!!
132              
133              	if((th_inp.gt.pi/2.).or.(th_oop.gt.pi/2.)) then
134              	   write(6,*) 'Decay pion going backwards - fahgettaboutit!'
135              	   success=.false.
136              	   return
137              	endif
138 gaskelld 1.1 
139              	orig%up%x = pxf/pf
140              	orig%up%y = pyf/pf
141              	orig%up%z = pzf/pf
142              
143              	orig%p%xptar = tan(th_oop)
144              	orig%p%yptar = tan(th_inp)
145              	orig%p%delta = 100.*(pf/p_spec-1.)
146              	
147              	orig%p%P = pf
148              
149              	Mh = Mpi
150              	Mh2 = Mpi2
151              

152 gaskelld 1.4 	Mh2_final = Mh2

153 gaskelld 1.1 
154              	orig%p%E = sqrt(orig%p%P**2 + Mh2)
155              

156 gaskelld 1.3 C Calculate "physics" angles
157              	call physics_angles(spec%p%theta,spec%p%phi,
158                   &     orig%p%xptar,orig%p%yptar,orig%p%theta,orig%p%phi)

159 gaskelld 1.1 
160              	return
161              	end
162              
163              
164