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.2 integer*4 junk1,ikinsetting
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.10 real yf,ff,thetaf,dptgt2
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 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 brash 1.10 p_spec=ptgt
96 brash 1.2 c
97 brash 1.9 if(nevent.le.2)write(,)'****** mom = ',p_spec,' ****'
98 brash 1.4 c
99 brash 1.10 call grndm(rndm,3) 100 dptgt2=rndm(1)0.200-0.100 101 xfp=rndm(2)80.00-40.00 102 phfp=rndm(3)5.0-2.5 103 c 104 call grndm(rndm,3) 105 yfp=rndm(1)60.00-30.00 106 thfp=rndm(2)*5.0-2.5
107 brash 1.4 c
108 brash 1.8 c write(6,)'x,phi,y,theta,ok_spec =' 109 c write(6,)xfp,phfp,yfp,thfp,ok_spec
110 brash 1.4 c 111 do ii=1,20 112 ntuple_array(ii)=0.0 113 enddo 114 ntuple_array(11)=real(xfp) 115 ntuple_array(12)=real(yfp) 116 ntuple_array(13)=real(thfp) 117 ntuple_array(14)=real(phfp) 118 c
119 brash 1.3 pcentral=ptgt
120 brash 1.10 pfp=dptgt2
121 brash 1.4 c
122 brash 1.2 c Finally, we convert this information over to a format that GEANT likes. 123 c
124 jones 1.1
125 brash 1.2 1135 dpmom=pfp 126 pfp=pcentral(1.0+pfp) 127 pmom=pfp 128 ea=sqrt(pfp2+938.2796*2) 129 tinit=ea-938.2796 130 vert(1)=xfp 131 vert(2)=yfp 132 vert(3)=0.0 133 ntbeam = 0.0 134 nttarg = 0.0 135 ubuf(1) = 0.0 136 137 1000 continue 138 c 139 ipart = 14 ! geant pid (8=pi+,9,pi-, 14 =p) 140 ptot=pfp/1000.
141 jones 1.1 call grndm ( rndm , 3 )
142 brash 1.2 angr=thfp3.14159265/180.0 143 phir=phfp3.14159265/180.0 144 psir=datan(dtan(phir)*dcos(angr)) 145 146 ang=thfp 147 phi=phfp
148 jones 1.1 c
149 brash 1.2 c these are the actual parameters for the track.
150 jones 1.1 c
151 brash 1.2 xinit=vert(1) 152 yinit=vert(2) 153 thinit=angr 154 phiinit=phir 155 156 c sptransport.l.particle.fp_h.ph=dtan(angr) 157 c sptransport.l.particle.fp_h.th=dtan(phir) 158 c 159 c sptransport.l.particle.fp_h.x=vert(1)/100.0 160 c 161 c sptransport.l.particle.fp_h.y=vert(2)/100.0
162 jones 1.1 c 163 c
164 brash 1.2 c next we misalign the track to simulate misalignment of the 165 c entire space frame with respect to the vdc's
166 jones 1.1 c
167 brash 1.2 c just include the translational offsets to start
168 jones 1.1 c
169 brash 1.2 c write(,)vert(1),vert(2),angr,phir,psir 170 171 xofff=0.0 172 yofff=0.0 173 thofff=0.0 174 phofff=0.0 175 psofff=0.0 176 c 177 c Define the inverse Euler rotation for (thofff,phiofff,psiofff) 178 c 179 c 180 rotmat(1,1)=dcos(psofff)dcos(thofff)+dsin(psofff)dsin(thofff) 181 $ dsin(phofff) 182 rotmat(1,2)=-dcos(phofff)dsin(thofff) 183 rotmat(1,3)=-dsin(psofff)dcos(thofff)+dcos(psofff)dsin(thofff) 184 $ dcos(phofff) 185 rotmat(2,1)=dcos(psofff)dsin(thofff)-dsin(psofff)dcos(thofff) 186 $ dsin(phofff) 187 rotmat(2,2)=dcos(phofff)dcos(thofff) 188 rotmat(2,3)=-dsin(psofff)dsin(thofff)-dcos(psofff)dcos(thofff) 189 $ dsin(phofff) 190 brash 1.2 rotmat(3,1)=dsin(psofff)dcos(phofff) 191 rotmat(3,2)=dsin(phofff) 192 rotmat(3,3)=dcos(psofff)dcos(phofff) 193 c 194 c 195 vert(1)=vert(1)-xofff 196 vert(2)=vert(2)-yofff 197 c 198 xyz(1)=dsin(psir) 199 xyz(2)=dcos(psir)dsin(angr) 200 xyz(3)=dcos(psir)dcos(angr) 201 c 202 do ic=1,3 203 xyznew(ic)=xyz(1)rotmat(ic,1)+ 204 & xyz(2)rotmat(ic,2)+xyz(3)rotmat(ic,3) 205 enddo 206 207 angr=datan(xyznew(2)/xyznew(3)) 208 termang=xyznew(3)/dcos(angr) 209 if(termang.gt.1.0) termang=1.0 210 if(termang.lt.-1.0) termang=-1.0 211 brash 1.2 psir=dacos(termang)abs(xyznew(1))/xyznew(1) 212 phir=datan(dtan(psir)/dcos(angr)) 213 214 c write(,)vert(1),vert(2),angr,phir,psir
215 jones 1.1 216 call gsvert ( vert,ntbeam,nttarg,ubuf,0,nvtx )
217 brash 1.2 etot = ea - 938.2796
218 jones 1.1 plab(1) = ptotdsin(psir) 219 plab(2) = ptotdsin(angr)dcos(psir) 220 plab(3) = ptotdcos(angr)*dcos(psir)
221 brash 1.2
222 jones 1.1 call hfill ( 100, etot, 0., 1. )
223 brash 1.2
224 jones 1.1 call gskine ( plab,ipart,nvtx,ubuf,0,nt )
225 brash 1.2
226 jones 1.1 if ( nt.le.0 ) then 227 write ( 6,* ) ' gukine: error defining track' 228 write ( 6,* ) ' i=',i,' nt=',nt 229 stop 230 end if
231 brash 1.2 232 nevent=nevent+1 233
234 jones 1.1 return
235 brash 1.2 236 end 237 238 subroutine kincalc(e0,eang,hang,trg,xtgt,ytgt, 239 $ thtgt,phtgt,ptgt,dptgt) 240 241 implicit none 242 243 real8 e0,eang,hang,trg,xtgt,ytgt,thtgt,phtgt,ptgt,dptgt 244 real8 fg,gf 245 integer i,j,k 246 real8 mt,mtg,mr,mpi,mn,mp,me,pi,mhe,alpha 247 real8 escat,pscat,pcentral,thetae,phie,thetap,phip 248 real rndm(3) 249 250 fg=3.14159265/180.0 251 gf=1.0/fg 252 253 me = 0.511 254 mpi = 139.57 255 mp = 938.2796 256 brash 1.2 mn = 939.5731 257 mhe = 2808.41 258 mt = mp 259 mtg = mt/1.e3 260 261 alpha=1./137. 262 263 264 1432 call grndm ( rndm , 3 ) 265 escat=mp/(1.0+mp/e0-cos(fgeang)) 266 pscat=sqrt(e02+escat2- 267 $ 2.0e0escatcos(fg*eang)) 268 pcentral=pscat 269 270 c First, we assume that whichever are is the more backward will 271 c determine the acceptance. We then randomly choose the theta 272 c and phi for the arm that determines the acceptance within the 273 c usual full HRS acceptance and then determine from kinematics 274 c what the corresponding theta and phi are for the other arm. 275 276 if (abs(eang).gt.abs(hang)) then
277 brash 1.7 c phie=-.025+rndm(1).050 278 c thetae=-.009+rndm(2).018 279 if(pcentral.ge.5000) then
280 brash 1.9 c phie=-.00+rndm(1).00 281 c thetae=-.00+rndm(2).00
282 brash 1.7 phie=-.130+rndm(1).260 283 thetae=-.065+rndm(2).130 284 else if(pcentral.ge.3000.and.pcentral.lt.5000) then 285 phie=-.067+rndm(1).135 286 thetae=-.034+rndm(2).067 287 else 288 phie=-.087+rndm(1).174 289 thetae=-.044+rndm(2).087 290 endif 291 c phie=-.00+rndm(1).00 292 c thetae=-.00+rndm(2).00
293 brash 1.2 escat=mp/(1.0+mp/e0-cos(fgeang+thetae)cos(phie)) 294 pscat=sqrt(e02+escat2- 295 $ 2.0e0escatcos(fgeang+thetae)cos(phie)) 296 phip=dasin(escatdsin(phie)/pscat) 297 thetap=dasin((escatsin(fgeang+thetae)cos(phie))/ 298 $ (pscatcos(phip)))-fg*hang 299 else 300 c Hadron arm defining acceptance 301 1221 call grndm ( rndm, 3)
302 brash 1.7 phie=-.080+rndm(1).160 303 thetae=-.030+rndm(2).060 304 c phie=-.00+rndm(1).00 305 c thetae=-.00+rndm(2).00
306 brash 1.2 escat=mp/(1.0+mp/e0-cos(fgeang+thetae)cos(phie)) 307 pscat=sqrt(e02+escat2- 308 $ 2.0e0escatcos(fgeang+thetae)cos(phie)) 309 phip=-1.0dasin(escatsin(phie)/pscat) 310 thetap=dasin((escatsin(fgeang+thetae)cos(phie))/ 311 $ (pscatcos(phip)))-fghang 312 if(abs(thetap).gt.0.030.or.abs(phip).gt.0.065) goto 1221 313 endif 314 315 ptgt=pscat
316 brash 1.5 c write(,)escat,eang,pscat,hang
317 brash 1.2 dptgt=(pcentral-pscat)/pcentral
318 jones 1.1 c
319 brash 1.2 c Following statement to fill just the high and low dp bins 320 c which have very low statistics. This is normally commented 321 c out. 322 c 323 c if(abs(dptgt).le.0.030) goto 1432 324 c 325 thtgt=thetap 326 phtgt=phip 327 xtgt=0.0 328 ytgt=0.0 329 330 return
331 jones 1.1 end

Analyzer/Replay: Mark Jones, Documents: Stephen Wood