geant_gep/src/gukine.f - annotate

Return to gukine.f CVS log

Up to [HallC] / geant_gep / src

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 brash 1.4 real r2,r3
24 jones 1.1 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.13 integer*4 junk1,ikinsetting,isectrack
29 jones 1.1 real8 xfp,yfp,tthfp,tphfp,pfp,junk2,junk3,thfp,phfp 30 real8 e0(10),eang(10),hang(10),targ_thick(10)
31 brash 1.2 real*8 xtgt,ytgt,thtgt,phtgt,ptgt,dptgt
32 brash 1.15 real yf,ff,thetaf
33 jones 1.1 character*80 junkline 34 c
35 brash 1.4 c The following variables are the SIMC variables 36 c 37 real8 x,y,z !(cm) 38 real8 dpp !delta p/p (%) 39 real8 dxdz,dydz !X,Y slope in spectrometer 40 real8 x_fp,y_fp,dx_fp,dy_fp !Focal plane values to return 41 real8 p_spec,th_spec !spectrometer setting 42 real8 fry !vertical position@tgt (+y=down) 43 real8 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 real8 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 brash 1.13 c 61 nhu1=0 62 nhx1=0 63 nhv1=0 64 nhu2=0 65 nhx2=0 66 nhv2=0 67 nhu3=0 68 nhx3=0 69 nhv3=0 70 nhu4=0 71 nhx4=0 72 nhv4=0
73 brash 1.15 theta_front=999.0
74 brash 1.13 c
75 brash 1.2 call grndm ( rndm , 3 )
76 jones 1.1 c
77 brash 1.2 111 format(a80) 78 iremain=nevent/1000 79 remain=nevent/1000.0 80 c if(iremain.eq.remain) write(,)'nevent =',nevent
81 jones 1.1 write(,)'nevent =',nevent 82 83 if(nevent.eq.0) then
84 brash 1.2 c write(,)'Getting kinematics setting ...'
85 jones 1.1 open(unit=1,file='geant_kinematics.dat',type='UNKNOWN') 86 read(1,*)ikinsetting
87 brash 1.13 read(1,*)isectrack
88 jones 1.1 close(unit=1)
89 brash 1.13 if(isectrack.eq.1) then 90 sectrack=.true. 91 else 92 sectrack=.false. 93 endif
94 jones 1.1 open(unit=1,file='hdr_gep.dat',status='old') 95 do i=1,10 96 read(1,)e0(i),eang(i),hang(i),targ_thick(i) 97 enddo 98 einc=e0(ikinsetting) 99 hrse_ang=eang(ikinsetting) 100 hrsh_ang=-1.0hang(ikinsetting)
101 brash 1.2 trg_thk=targ_thick(ikinsetting)
102 jones 1.1 close(unit=1) 103 endif 104
105 brash 1.2 c 106 c Now we have the incident electron energy and the angles of 107 c the two spectrometers, as well as the target thickness. 108 c We need to use a) kinematics routines and b) some knowledge 109 c of the extended target acceptances to choose the x,y,theta,phi, 110 c and dpmom at the target. 111 c 112 call kincalc(einc,hrse_ang,hrsh_ang,trg_thk,xtgt,ytgt, 113 $ thtgt,phtgt,ptgt,dptgt) 114 115 c
116 brash 1.15 c At this point, after using polynomials from Lerose, we should 117 c have x,y,theta,phi, and dpmom at the focal plane. 118 c 119 c write(,)'Back from kincalc ...' 120 c write(,)xtgt,ytgt,thtgt,phtgt,dptgt 121 c x(1)=tan(phtgt) 122 c x(2)=ytgt 123 c x(3)=tan(thtgt) 124 c x(4)=dptgt 125 y_target=ytgt100.0 126 c 127 c Here is where we call the SIMC routines to calculate the focal plane 128 c quantities from the target quantities. 129 c 130 c xfp=xfinal(x,4)100.0 131 c phfp=atan(thetaf(x,4))/3.14159265180.0 132 c yfp=yf(x,4)100.0 133 c thfp=atan(ff(x,4))/3.14159265*180.0
134 brash 1.14 p_spec=ptgt
135 brash 1.15 th_spec=hrsh_ang 136 dpp=dptgt100.00 137 x=xtgt 138 y=ytgt 139 z=0.0 140 dxdz=tan(phtgt3.14159265/180.0) 141 dydz=tan(thtgt*3.14159265/180.0)
142 brash 1.9 if(nevent.le.2)write(,)'****** mom = ',p_spec,' ****'
143 brash 1.15 c write(,)x,y,z 144 c write(,)dxdz,dydz
145 brash 1.4 c
146 brash 1.15 c Call to SIMC routine
147 brash 1.4 c
148 brash 1.15 decay_flag=.false. 149 wcs_flag=.false. 150 ms_flag=.false. 151 m2=(938.2796*2) 152 call mc_hms (p_spec, th_spec, dpp, x, y, z, dxdz, dydz, 153 > x_fp, dx_fp, y_fp, dy_fp, m2, 154 > ms_flag, wcs_flag, decay_flag, resmult, fry, 155 > ok_spec, pathlen) 156 c write(,)'Momentum = ',p_spec,' Angle = ',th_spec 157 c write(,)x,y,z 158 c write(,)dxdz,dydz 159 c write(6,)x_fp,dx_fp,y_fp,dy_fp,ok_spec 160 c write(,)m2,ms_flag,wcs_flag,decay_flag 161 c write(,)resmult,fry,ok_spec,pathlen 162 163 xfp=x_fp 164 phfp=atan(dx_fp)/3.14159265180.0 165 yfp=y_fp100.0 166 thfp=atan(dy_fp)/3.14159265180.0100.0
167 brash 1.8 c write(6,)'x,phi,y,theta,ok_spec =' 168 c write(6,)xfp,phfp,yfp,thfp,ok_spec
169 brash 1.4 c 170 do ii=1,20 171 ntuple_array(ii)=0.0 172 enddo 173 ntuple_array(11)=real(xfp) 174 ntuple_array(12)=real(yfp) 175 ntuple_array(13)=real(thfp) 176 ntuple_array(14)=real(phfp) 177 c
178 brash 1.3 pcentral=ptgt
179 brash 1.15 pfp=dptgt
180 brash 1.4 c
181 brash 1.15 c 182 c write(,)'Focal plane quantities ...' 183 c write(,)xfp,yfp,thfp,phfp,pcentral,pfp 184 185 c xfp=xtgt ! in cm 186 c yfp=ytgt ! in cm 187 c thfp=thtgt/3.14159265180.0 ! in degrees 188 c phfp=phtgt/3.14159265180.0 ! in degrees 189 c pcentral=ptgt ! in MeV/c 190 c pfp=dptgt ! fraction of central momentum 191 c
192 brash 1.2 c Finally, we convert this information over to a format that GEANT likes. 193 c
194 jones 1.1
195 brash 1.2 1135 dpmom=pfp 196 pfp=pcentral(1.0+pfp) 197 pmom=pfp 198 ea=sqrt(pfp2+938.2796*2) 199 tinit=ea-938.2796 200 vert(1)=xfp 201 vert(2)=yfp
202 brash 1.15 vert(3)=0.0
203 brash 1.2 ntbeam = 0.0 204 nttarg = 0.0 205 ubuf(1) = 0.0 206 207 1000 continue 208 c 209 ipart = 14 ! geant pid (8=pi+,9,pi-, 14 =p) 210 ptot=pfp/1000.
211 jones 1.1 call grndm ( rndm , 3 )
212 brash 1.2 angr=thfp3.14159265/180.0 213 phir=phfp3.14159265/180.0 214 psir=datan(dtan(phir)*dcos(angr)) 215 216 ang=thfp 217 phi=phfp
218 jones 1.1 c
219 brash 1.2 c these are the actual parameters for the track.
220 jones 1.1 c
221 brash 1.2 xinit=vert(1) 222 yinit=vert(2) 223 thinit=angr 224 phiinit=phir 225 226 c sptransport.l.particle.fp_h.ph=dtan(angr) 227 c sptransport.l.particle.fp_h.th=dtan(phir) 228 c 229 c sptransport.l.particle.fp_h.x=vert(1)/100.0 230 c 231 c sptransport.l.particle.fp_h.y=vert(2)/100.0
232 jones 1.1 c 233 c
234 brash 1.2 c next we misalign the track to simulate misalignment of the 235 c entire space frame with respect to the vdc's
236 jones 1.1 c
237 brash 1.2 c just include the translational offsets to start
238 jones 1.1 c
239 brash 1.2 c write(,)vert(1),vert(2),angr,phir,psir 240 241 xofff=0.0 242 yofff=0.0 243 thofff=0.0 244 phofff=0.0 245 psofff=0.0 246 c 247 c Define the inverse Euler rotation for (thofff,phiofff,psiofff) 248 c 249 c 250 rotmat(1,1)=dcos(psofff)dcos(thofff)+dsin(psofff)dsin(thofff) 251 $ dsin(phofff) 252 rotmat(1,2)=-dcos(phofff)dsin(thofff) 253 rotmat(1,3)=-dsin(psofff)dcos(thofff)+dcos(psofff)dsin(thofff) 254 $ dcos(phofff) 255 rotmat(2,1)=dcos(psofff)dsin(thofff)-dsin(psofff)dcos(thofff) 256 $ dsin(phofff) 257 rotmat(2,2)=dcos(phofff)dcos(thofff) 258 rotmat(2,3)=-dsin(psofff)dsin(thofff)-dcos(psofff)dcos(thofff) 259 $ dsin(phofff) 260 brash 1.2 rotmat(3,1)=dsin(psofff)dcos(phofff) 261 rotmat(3,2)=dsin(phofff) 262 rotmat(3,3)=dcos(psofff)dcos(phofff) 263 c 264 c 265 vert(1)=vert(1)-xofff 266 vert(2)=vert(2)-yofff 267 c 268 xyz(1)=dsin(psir) 269 xyz(2)=dcos(psir)dsin(angr) 270 xyz(3)=dcos(psir)dcos(angr) 271 c 272 do ic=1,3 273 xyznew(ic)=xyz(1)rotmat(ic,1)+ 274 & xyz(2)rotmat(ic,2)+xyz(3)rotmat(ic,3) 275 enddo 276 277 angr=datan(xyznew(2)/xyznew(3)) 278 termang=xyznew(3)/dcos(angr) 279 if(termang.gt.1.0) termang=1.0 280 if(termang.lt.-1.0) termang=-1.0 281 brash 1.2 psir=dacos(termang)abs(xyznew(1))/xyznew(1) 282 phir=datan(dtan(psir)/dcos(angr)) 283 284 c write(,)vert(1),vert(2),angr,phir,psir
285 jones 1.1 286 call gsvert ( vert,ntbeam,nttarg,ubuf,0,nvtx )
287 brash 1.2 etot = ea - 938.2796
288 jones 1.1 plab(1) = ptotdsin(psir) 289 plab(2) = ptotdsin(angr)dcos(psir) 290 plab(3) = ptotdcos(angr)*dcos(psir)
291 brash 1.2
292 jones 1.1 call hfill ( 100, etot, 0., 1. )
293 brash 1.2
294 jones 1.1 call gskine ( plab,ipart,nvtx,ubuf,0,nt )
295 brash 1.2
296 jones 1.1 if ( nt.le.0 ) then 297 write ( 6,* ) ' gukine: error defining track' 298 write ( 6,* ) ' i=',i,' nt=',nt 299 stop 300 end if
301 brash 1.2 302 nevent=nevent+1 303
304 jones 1.1 return
305 brash 1.2 306 end 307 308 subroutine kincalc(e0,eang,hang,trg,xtgt,ytgt, 309 $ thtgt,phtgt,ptgt,dptgt) 310 311 implicit none 312 313 real8 e0,eang,hang,trg,xtgt,ytgt,thtgt,phtgt,ptgt,dptgt 314 real8 fg,gf 315 integer i,j,k 316 real8 mt,mtg,mr,mpi,mn,mp,me,pi,mhe,alpha 317 real8 escat,pscat,pcentral,thetae,phie,thetap,phip 318 real rndm(3) 319 320 fg=3.14159265/180.0 321 gf=1.0/fg 322 323 me = 0.511 324 mpi = 139.57 325 mp = 938.2796 326 brash 1.2 mn = 939.5731 327 mhe = 2808.41 328 mt = mp 329 mtg = mt/1.e3 330 331 alpha=1./137. 332 333 334 1432 call grndm ( rndm , 3 ) 335 escat=mp/(1.0+mp/e0-cos(fgeang)) 336 pscat=sqrt(e02+escat2- 337 $ 2.0e0escatcos(fg*eang)) 338 pcentral=pscat 339 340 c First, we assume that whichever are is the more backward will 341 c determine the acceptance. We then randomly choose the theta 342 c and phi for the arm that determines the acceptance within the 343 c usual full HRS acceptance and then determine from kinematics 344 c what the corresponding theta and phi are for the other arm. 345 346 if (abs(eang).gt.abs(hang)) then
347 brash 1.7 c phie=-.025+rndm(1).050 348 c thetae=-.009+rndm(2).018 349 if(pcentral.ge.5000) then
350 brash 1.9 c phie=-.00+rndm(1).00 351 c thetae=-.00+rndm(2).00
352 brash 1.7 phie=-.130+rndm(1).260 353 thetae=-.065+rndm(2).130 354 else if(pcentral.ge.3000.and.pcentral.lt.5000) then 355 phie=-.067+rndm(1).135 356 thetae=-.034+rndm(2).067 357 else 358 phie=-.087+rndm(1).174 359 thetae=-.044+rndm(2).087 360 endif 361 c phie=-.00+rndm(1).00 362 c thetae=-.00+rndm(2).00
363 brash 1.2 escat=mp/(1.0+mp/e0-cos(fgeang+thetae)cos(phie)) 364 pscat=sqrt(e02+escat2- 365 $ 2.0e0escatcos(fgeang+thetae)cos(phie)) 366 phip=dasin(escatdsin(phie)/pscat) 367 thetap=dasin((escatsin(fgeang+thetae)cos(phie))/ 368 $ (pscatcos(phip)))-fg*hang 369 else 370 c Hadron arm defining acceptance 371 1221 call grndm ( rndm, 3)
372 brash 1.7 phie=-.080+rndm(1).160 373 thetae=-.030+rndm(2).060 374 c phie=-.00+rndm(1).00 375 c thetae=-.00+rndm(2).00
376 brash 1.2 escat=mp/(1.0+mp/e0-cos(fgeang+thetae)cos(phie)) 377 pscat=sqrt(e02+escat2- 378 $ 2.0e0escatcos(fgeang+thetae)cos(phie)) 379 phip=-1.0dasin(escatsin(phie)/pscat) 380 thetap=dasin((escatsin(fgeang+thetae)cos(phie))/ 381 $ (pscatcos(phip)))-fghang 382 if(abs(thetap).gt.0.030.or.abs(phip).gt.0.065) goto 1221 383 endif 384 385 ptgt=pscat
386 brash 1.5 c write(,)escat,eang,pscat,hang
387 brash 1.2 dptgt=(pcentral-pscat)/pcentral
388 jones 1.1 c
389 brash 1.2 c Following statement to fill just the high and low dp bins 390 c which have very low statistics. This is normally commented 391 c out. 392 c 393 c if(abs(dptgt).le.0.030) goto 1432 394 c 395 thtgt=thetap 396 phtgt=phip 397 xtgt=0.0 398 ytgt=0.0 399 400 return
401 jones 1.1 end

Analyzer/Replay: Mark Jones, Documents: Stephen Wood