1 jones 1.1 subroutine gukine
2 c
|
3 brash 1.2 c
4 c This is a subroutine to generate (x,y,theta,phi,delta-p) at
5 c the focal plane for proton from the H(e,e'p) reaction.
6 c It makes use of a simple kinematics routine, some knowledge
7 c of the extended target acceptances, and polynomials for the
8 c transport through the spectrometer from John Lerose.
|
9 jones 1.1 c
10 implicit none
11 integer ikine,itra,istak,ivert,ipart,itrtyp,napart,ipaold
12 real pkine,amass,charge,tlife,vert,pvert
13 common/gckine/ikine,pkine(10),itra,istak,ivert,ipart,itrtyp
14 + ,napart(5),amass,charge,tlife,vert(3),pvert(4),ipaold
15 integer nubuf
16 parameter (nubuf=10)
17 real ubuf(nubuf),rndm(3)
18 integer ntbeam,nttarg,nvtx
19 save nvtx
20 real plab(3),ptot,etot
21 real*8 ang,angr,phi,phir,psir
|
22 brash 1.2 real mp,ea,pa,remain
|
23 jones 1.1 real r2,r3
24 integer nt,ierri,i0,i1,i2,i3
|
25 brash 1.2 integer i,j,k,ichoice,ichoice2,iremain
|
26 jones 1.1 real*8 rotmat,xyz(3),xyznew(3),termang
27 integer ic,ii,jj
|
28 brash 1.2 integer*4 junk1,ikinsetting
|
29 jones 1.1 real*8 xfp,yfp,tthfp,tphfp,pfp,junk2,junk3,thfp,phfp
30 real*8 e0(10),eang(10),hang(10),targ_thick(10)
|
31 brash 1.2 real*8 xtgt,ytgt,thtgt,phtgt,ptgt,dptgt
32 real yf,ff,thetaf
|
33 jones 1.1 character*80 junkline
34 c
|
35 brash 1.2 c The following variables are the SIMC variables
36 c
37 real*8 x,y,z !(cm)
38 real*8 dpp !delta p/p (%)
39 real*8 dxdz,dydz !X,Y slope in spectrometer
40 real*8 x_fp,y_fp,dx_fp,dy_fp !Focal plane values to return
41 real*8 p_spec,th_spec !spectrometer setting
42 real*8 fry !vertical position@tgt (+y=down)
43 real*8 pathlen
44 logical ms_flag !mult. scattering flag
45 logical wcs_flag !wire chamber smearing flag
46 logical decay_flag !check for particle decay
47 logical ok_spec !true if particle makes it
48 real*8 resmult,m2
49 c
|
50 jones 1.1 include 'fpp_local.h'
51 include 'geant_local.h'
52 c
53 c include 'parameter.h'
54 c include 'espace_type.h'
55 c include 'detector.h'
56 c include 'transport.h'
57 c include 'option.h'
58 c
|
59 brash 1.2 common/kincom/rotmat(3,3)
60 call grndm ( rndm , 3 )
|
61 jones 1.1 c
|
62 brash 1.2 111 format(a80)
63 iremain=nevent/1000
64 remain=nevent/1000.0
65 c if(iremain.eq.remain) write(*,*)'nevent =',nevent
|
66 jones 1.1 write(*,*)'nevent =',nevent
67
68 if(nevent.eq.0) then
|
69 brash 1.2 c write(*,*)'Getting kinematics setting ...'
|
70 jones 1.1 open(unit=1,file='geant_kinematics.dat',type='UNKNOWN')
71 read(1,*)ikinsetting
72 close(unit=1)
73 open(unit=1,file='hdr_gep.dat',status='old')
74 do i=1,10
75 read(1,*)e0(i),eang(i),hang(i),targ_thick(i)
76 enddo
77 einc=e0(ikinsetting)
78 hrse_ang=eang(ikinsetting)
79 hrsh_ang=-1.0*hang(ikinsetting)
|
80 brash 1.2 trg_thk=targ_thick(ikinsetting)
|
81 jones 1.1 close(unit=1)
82 endif
83
|
84 brash 1.2 c
85 c Now we have the incident electron energy and the angles of
86 c the two spectrometers, as well as the target thickness.
87 c We need to use a) kinematics routines and b) some knowledge
88 c of the extended target acceptances to choose the x,y,theta,phi,
89 c and dpmom at the target.
90 c
91 call kincalc(einc,hrse_ang,hrsh_ang,trg_thk,xtgt,ytgt,
92 $ thtgt,phtgt,ptgt,dptgt)
93
94 c
95 c At this point, after using polynomials from Lerose, we should
96 c have x,y,theta,phi, and dpmom at the focal plane.
97 c
98 c write(*,*)'Back from kincalc ...'
99 c write(*,*)xtgt,ytgt,thtgt,phtgt,dptgt
100 c x(1)=tan(phtgt)
101 c x(2)=ytgt
102 c x(3)=tan(thtgt)
103 c x(4)=dptgt
104 y_target=ytgt*100.0
105 brash 1.2 c
106 c Here is where we call the SIMC routines to calculate the focal plane
107 c quantities from the target quantities.
108 c
109 c xfp=xfinal(x,4)*100.0
110 c phfp=atan(thetaf(x,4))/3.14159265*180.0
111 c yfp=yf(x,4)*100.0
112 c thfp=atan(ff(x,4))/3.14159265*180.0
113 p_spec=ptgt
114 th_spec=hrsh_ang
115 dpp=dptgt*100.00
116 x=xtgt
117 y=ytgt
118 z=0.0
119 dxdz=tan(phtgt*3.14159265/180.0)
120 dydz=tan(thtgt*3.14159265/180.0)
121 write(*,*)p_spec,th_spec,dpp
122 write(*,*)x,y,z
123 write(*,*)dxdz,dydz
124 c
125 c Call to SIMC routine
126 brash 1.2 c
127 call mc_hms (p_spec, th_spec, dpp, x, y, z, dxdz, dydz,
128 > x_fp, dx_fp, y_fp, dy_fp, m2,
129 > ms_flag, wcs_flag, decay_flag, resmult, fry,
130 > ok_spec, pathlen)
131 write(*,*)p_spec,th_spec,dpp
132 write(*,*)x,y,z
133 write(*,*)dxdz,dydz
134 write(6,*)x_fp,dx_fp,y_fp,dy_fp
135 write(*,*)m2,ms_flag,wcs_flag,decay_flag
136 write(*,*)resmult,fry,ok_spec,pathlen
137
138 xfp=x_fp
139 phfp=atan(dx_fp)/3.14159265*180.0
140 yfp=y_fp
141 thfp=atan(dy_fp)/3.14159265*180.0
142
143 c
144 pcentral=ptgt
145 pfp=dptgt
146 c
147 brash 1.2 c
148 c write(*,*)'Focal plane quantities ...'
149 c write(*,*)xfp,yfp,thfp,phfp,pcentral,pfp
150
151 c xfp=xtgt ! in cm
152 c yfp=ytgt ! in cm
153 c thfp=thtgt/3.14159265*180.0 ! in degrees
154 c phfp=phtgt/3.14159265*180.0 ! in degrees
155 c pcentral=ptgt ! in MeV/c
156 c pfp=dptgt ! fraction of central momentum
157 c
158 c Finally, we convert this information over to a format that GEANT likes.
159 c
|
160 jones 1.1
|
161 brash 1.2 1135 dpmom=pfp
162 pfp=pcentral*(1.0+pfp)
163 pmom=pfp
164 ea=sqrt(pfp**2+938.2796**2)
165 tinit=ea-938.2796
166 vert(1)=xfp
167 vert(2)=yfp
168 vert(3)=0.0
169 ntbeam = 0.0
170 nttarg = 0.0
171 ubuf(1) = 0.0
172
173 1000 continue
174 c
175 ipart = 14 ! geant pid (8=pi+,9,pi-, 14 =p)
176 ptot=pfp/1000.
|
177 jones 1.1 call grndm ( rndm , 3 )
|
178 brash 1.2 angr=thfp*3.14159265/180.0
179 phir=phfp*3.14159265/180.0
180 psir=datan(dtan(phir)*dcos(angr))
181
182 ang=thfp
183 phi=phfp
|
184 jones 1.1 c
|
185 brash 1.2 c these are the actual parameters for the track.
|
186 jones 1.1 c
|
187 brash 1.2 xinit=vert(1)
188 yinit=vert(2)
189 thinit=angr
190 phiinit=phir
191
192 c sptransport.l.particle.fp_h.ph=dtan(angr)
193 c sptransport.l.particle.fp_h.th=dtan(phir)
194 c
195 c sptransport.l.particle.fp_h.x=vert(1)/100.0
196 c
197 c sptransport.l.particle.fp_h.y=vert(2)/100.0
|
198 jones 1.1 c
199 c
|
200 brash 1.2 c next we misalign the track to simulate misalignment of the
201 c entire space frame with respect to the vdc's
|
202 jones 1.1 c
|
203 brash 1.2 c just include the translational offsets to start
|
204 jones 1.1 c
|
205 brash 1.2 c write(*,*)vert(1),vert(2),angr,phir,psir
206
207 xofff=0.0
208 yofff=0.0
209 thofff=0.0
210 phofff=0.0
211 psofff=0.0
212 c
213 c Define the inverse Euler rotation for (thofff,phiofff,psiofff)
214 c
215 c
216 rotmat(1,1)=dcos(psofff)*dcos(thofff)+dsin(psofff)*dsin(thofff)
217 $ *dsin(phofff)
218 rotmat(1,2)=-dcos(phofff)*dsin(thofff)
219 rotmat(1,3)=-dsin(psofff)*dcos(thofff)+dcos(psofff)*dsin(thofff)
220 $ *dcos(phofff)
221 rotmat(2,1)=dcos(psofff)*dsin(thofff)-dsin(psofff)*dcos(thofff)
222 $ *dsin(phofff)
223 rotmat(2,2)=dcos(phofff)*dcos(thofff)
224 rotmat(2,3)=-dsin(psofff)*dsin(thofff)-dcos(psofff)*dcos(thofff)
225 $ *dsin(phofff)
226 brash 1.2 rotmat(3,1)=dsin(psofff)*dcos(phofff)
227 rotmat(3,2)=dsin(phofff)
228 rotmat(3,3)=dcos(psofff)*dcos(phofff)
229 c
230 c
231 vert(1)=vert(1)-xofff
232 vert(2)=vert(2)-yofff
233 c
234 xyz(1)=dsin(psir)
235 xyz(2)=dcos(psir)*dsin(angr)
236 xyz(3)=dcos(psir)*dcos(angr)
237 c
238 do ic=1,3
239 xyznew(ic)=xyz(1)*rotmat(ic,1)+
240 & xyz(2)*rotmat(ic,2)+xyz(3)*rotmat(ic,3)
241 enddo
242
243 angr=datan(xyznew(2)/xyznew(3))
244 termang=xyznew(3)/dcos(angr)
245 if(termang.gt.1.0) termang=1.0
246 if(termang.lt.-1.0) termang=-1.0
247 brash 1.2 psir=dacos(termang)*abs(xyznew(1))/xyznew(1)
248 phir=datan(dtan(psir)/dcos(angr))
249
250 c write(*,*)vert(1),vert(2),angr,phir,psir
|
251 jones 1.1
252 call gsvert ( vert,ntbeam,nttarg,ubuf,0,nvtx )
|
253 brash 1.2 etot = ea - 938.2796
|
254 jones 1.1 plab(1) = ptot*dsin(psir)
255 plab(2) = ptot*dsin(angr)*dcos(psir)
256 plab(3) = ptot*dcos(angr)*dcos(psir)
|
257 brash 1.2
|
258 jones 1.1 call hfill ( 100, etot, 0., 1. )
|
259 brash 1.2
|
260 jones 1.1 call gskine ( plab,ipart,nvtx,ubuf,0,nt )
|
261 brash 1.2
|
262 jones 1.1 if ( nt.le.0 ) then
263 write ( 6,* ) ' gukine: error defining track'
264 write ( 6,* ) ' i=',i,' nt=',nt
265 stop
266 end if
|
267 brash 1.2
268 nevent=nevent+1
269
|
270 jones 1.1 return
|
271 brash 1.2
272 end
273
274 subroutine kincalc(e0,eang,hang,trg,xtgt,ytgt,
275 $ thtgt,phtgt,ptgt,dptgt)
276
277 implicit none
278
279 real*8 e0,eang,hang,trg,xtgt,ytgt,thtgt,phtgt,ptgt,dptgt
280 real*8 fg,gf
281 integer i,j,k
282 real*8 mt,mtg,mr,mpi,mn,mp,me,pi,mhe,alpha
283 real*8 escat,pscat,pcentral,thetae,phie,thetap,phip
284 real rndm(3)
285
286 fg=3.14159265/180.0
287 gf=1.0/fg
288
289 me = 0.511
290 mpi = 139.57
291 mp = 938.2796
292 brash 1.2 mn = 939.5731
293 mhe = 2808.41
294 mt = mp
295 mtg = mt/1.e3
296
297 alpha=1./137.
298
299
300 1432 call grndm ( rndm , 3 )
301 escat=mp/(1.0+mp/e0-cos(fg*eang))
302 pscat=sqrt(e0**2+escat**2-
303 $ 2.0*e0*escat*cos(fg*eang))
304 pcentral=pscat
305
306 c First, we assume that whichever are is the more backward will
307 c determine the acceptance. We then randomly choose the theta
308 c and phi for the arm that determines the acceptance within the
309 c usual full HRS acceptance and then determine from kinematics
310 c what the corresponding theta and phi are for the other arm.
311
312 if (abs(eang).gt.abs(hang)) then
313 brash 1.2 phie=-.065+rndm(1)*.130
314 thetae=-.030+rndm(2)*.060
315 escat=mp/(1.0+mp/e0-cos(fg*eang+thetae)*cos(phie))
316 pscat=sqrt(e0**2+escat**2-
317 $ 2.0*e0*escat*cos(fg*eang+thetae)*cos(phie))
318 phip=dasin(escat*dsin(phie)/pscat)
319 thetap=dasin((escat*sin(fg*eang+thetae)*cos(phie))/
320 $ (pscat*cos(phip)))-fg*hang
321 else
322 c Hadron arm defining acceptance
323 1221 call grndm ( rndm, 3)
324 phie=-.065+rndm(1)*.130
325 thetae=-.030+rndm(2)*.060
326 escat=mp/(1.0+mp/e0-cos(fg*eang+thetae)*cos(phie))
327 pscat=sqrt(e0**2+escat**2-
328 $ 2.0*e0*escat*cos(fg*eang+thetae)*cos(phie))
329 phip=-1.0*dasin(escat*sin(phie)/pscat)
330 thetap=dasin((escat*sin(fg*eang+thetae)*cos(phie))/
331 $ (pscat*cos(phip)))-fg*hang
332 if(abs(thetap).gt.0.030.or.abs(phip).gt.0.065) goto 1221
333 endif
334 brash 1.2
335 ptgt=pscat
336 dptgt=(pcentral-pscat)/pcentral
|
337 jones 1.1 c
|
338 brash 1.2 c Following statement to fill just the high and low dp bins
339 c which have very low statistics. This is normally commented
340 c out.
341 c
342 c if(abs(dptgt).le.0.030) goto 1432
343 c
344 thtgt=thetap
345 phtgt=phip
346 xtgt=0.0
347 ytgt=0.0
348
349 return
|
350 jones 1.1 end
|