geant_gep/src/gustep.f - annotate

Return to gustep.f CVS log

Up to [HallC] / geant_gep / src

1 jones 1.1 subroutine gustep 2 c 3 c this subroutine was written by jpsullivan april 21-22, 1993 4 c the tracking realated part is relatively simple -- if the 5 c particle leave the volume called 'targ', throw it away. 6 c it also makes a bunch of histograms 7 c 8 c keep,gctrak. 9 -- author : 10 common/gctrak/vect(7),getot,gekin,vout(7),nmec,lmec(30),namec(30) 11 + ,nstep ,maxnst,destep,destel,safety,sleng 12 + ,step ,snext ,sfield 13 + ,tofg ,gekrat,upwght,ignext,inwvol,istop ,idecad,iekbin 14 + , ilosl, imull,ingoto,nldown,nlevin,nlvsav,istory 15 c 16 integer nmec,lmec,namec,nstep ,maxnst,ignext,inwvol,istop 17 + ,idecad,iekbin,ilosl, imull,ingoto,nldown,nlevin 18 + ,nlvsav,istory 19 real vect,getot,gekin,vout,destep,destel,safety,sleng ,step 20 + ,snext,sfield,tofg ,gekrat,upwght 21 c end gctrak 22 jones 1.1 * keep,gcvolu. 23 -- author : 24 common/gcvolu/nlevel,names(15),number(15), 25 + lvolum(15),lindex(15),infrom,nlevmx,nldev(15),linmx(15), 26 + gtran(3,15),grmat(10,15),gonly(15),glx(3) 27 c 28 integer nlevel,number,lvolum,lindex,infrom,nlevmx, 29 + nldev,linmx 30 character4 names 31 real gtran,grmat,gonly,glx 32 c end gcvolu 33 c 34 * keep,gcbank. 35 -- author : 36 integer iq,lq,nzebra,ixstor,ixdiv,ixcons,lmain,lr1 37 integer kwbank,kwwork,iws 38 real gversn,zversn,fendq,ws,q 39 c 40 parameter (kwbank=69000,kwwork=5200) 41 common/gcbank/nzebra,gversn,zversn,ixstor,ixdiv,ixcons,fendq(16) 42 + ,lmain,lr1,ws(kwbank) 43 jones 1.1 dimension iq(2),q(2),lq(8000),iws(2) 44 equivalence (q(1),iq(1),lq(9)),(lq(1),lmain),(iws(1),ws(1)) 45 common/gclink/jdigi ,jdraw ,jhead ,jhits ,jkine ,jmate ,jpart 46 + ,jrotm ,jrung ,jset ,jstak ,jgstat,jtmed ,jtrack,jvertx 47 + ,jvolum,jxyz ,jgpar ,jgpar2,jsklt 48 c 49 integer jdigi ,jdraw ,jhead ,jhits ,jkine ,jmate ,jpart 50 + ,jrotm ,jrung ,jset ,jstak ,jgstat,jtmed ,jtrack,jvertx 51 + ,jvolum,jxyz ,jgpar,jgpar2 ,jsklt 52 c 53 keep,gcking. 54 -- author : 55 common/gcking/kcase,ngkine,gkin(5,100),tofd(100),iflgk(100) 56 integer kcase,ngkine ,iflgk 57 real gkin,tofd 58 c end gcking 59 c 60 keep,gckine. 61 -- author : 62 -- author : 63 integer ikine,itra,istak,ivert,ipart,itrtyp,napart,ipaold 64 jones 1.1 real pkine,amass,charge,tlife,vert,pvert 65 common/gckine/ikine,pkine(10),itra,istak,ivert,ipart,itrtyp 66 + ,napart(5),amass,charge,tlife,vert(3),pvert(4),ipaold 67 c end gckine 68 c 69 integer ihset,ihdet,iset,idet,idtype,nvname,numbv 70 common/gcsets/ihset,ihdet,iset,idet,idtype,nvname,numbv(20) 71 real x1,y1,z1,lpar,v1,v2,v3,newdist,x1new,y1new,z1new 72 real xstr,ystr,zstr,xstrnew,ystrnew,zstrnew 73 real rotmat2,rotmat3,rotmat4,rotmat1
74 brash 1.6 c
75 jones 1.1 common/geomstep/rotmat1(3,3),rotmat2(3,3),rotmat3(3,3), 76 & rotmat4(3,3) 77 c 78 include 'fpp_local.h' 79 include 'geant_local.h' 80 c include 'parameter.h' 81 c include 'espace_type.h' 82 c include 'detector.h' 83 c include 'transport.h' 84 c include 'option.h' 85 c 86 c 87 integer i,make_hist,ioff,ihit,ieffcheck 88 real pt_pi,ppar_pi,arg1,arg2,rapid_pi,pchmb,hits(6) 89 real a,b,c,beta,z,y,straw,ypath,rndm(3),ycompare
90 brash 1.2 c write(6,)'entering gustep' 91 c write(6,)'inwvol =',inwvol 92 c write(6,)'position =',vect(1),vect(2),vect(3) 93 c write(6,)'names =',names(nlevel)
94 jones 1.1 c 95 c 96 c 97 if ( ngkine.gt.0. ) then 98 do i=1,nmec 99 if ( lmec(i).eq.12 ) then 100 c write ( 6,* ) ' gustep: hadronic interaction' 101 c write ( 6,* ) ' nevent=',nevent 102 end if 103 end do 104 mylast = min(100,ngkine) 105 do i=1,mylast 106 iflgk(i) = 1 107 if ( gkin(5,i).eq.4 ) iflgk(i) = 0 108 if ( gkin(4,i).gt.0.001 ) iflgk(i)=0 109 end do 110 c n_2nd = n_2nd + ngkine 111 endif 112 c 113 c the following call makes sure all of the secondary particles 114 c get tracked too (provided the flag iflgk(i) for that particle 115 jones 1.1 c was set in the loop above -- this point is not correctly or 116 c clearly documented in the version of the geant manual i have). 117 c 118 c if(sectrack) then 119 c call gsking ( 0 ) 120 c endif 121 c
122 brash 1.6 c write(,)'In GUSTEP: ... names(nlevel) = ',names(nlevel) 123
124 jones 1.1 make_hist = 0 125 if ( inwvol.eq.1 .and. names(nlevel).eq."hall" ) then 126 make_hist=1 127 c 128 c if we get here, the tracking is done for this track 129 c (istop=1) but make a bunch of histograms before exitting 130 c note that ipart=8 means a pi+ and 9 is a pi- 131 c 132 c istop = 1 133 else if ( istop.ne.0.and. names(nlevel).eq."aira" ) then 134 make_hist=0 135 else if ( istop.ne.0.and. names(nlevel).eq."airb" ) then 136 make_hist=0 137 else if ( istop.ne.0.and. names(nlevel).eq."airc" ) then 138 make_hist=0 139 else if ( istop.ne.0.and. names(nlevel).eq."aird" ) then 140 make_hist=0
141 brash 1.2 else if ( istop.ne.0.and. names(nlevel).eq."aire" ) then 142 make_hist=0 143 else if ( istop.ne.0.and. names(nlevel).eq."airf" ) then 144 make_hist=0
145 brash 1.6 else if ( names(nlevel).eq."airg" ) then 146 c write(,)'In airg ... inwvol = ',inwvol 147 c write(,)'Z-value = ',vect(3) 148 if (istop.ne.0) then 149 make_hist=0 150 endif 151 else if ( names(nlevel).eq."HALL" ) then 152 c write(,)'In HALL ... inwvol = ',inwvol 153 c write(,)'Z-value = ',vect(3) 154 if (istop.ne.0) then 155 make_hist=0 156 endif
157 brash 1.2 else if ( istop.ne.0.and. names(nlevel).eq."airh" ) then 158 make_hist=0 159 else if ( istop.ne.0.and. names(nlevel).eq."hch1" ) then 160 make_hist=0 161 else if ( istop.ne.0.and. names(nlevel).eq."hch2" ) then 162 make_hist=0
163 brash 1.4 else if ( names(nlevel).eq."fch1" ) then
164 brash 1.6 c write(,)'In fch1 ... inwvol = ',inwvol
165 brash 1.4 if(inwvol.eq.1) then
166 brash 1.5 c write(6,)'Coordinates at fch1' 167 c write(6,)'x=',vect(1),' y=',vect(2),' z=',vect(3)
168 brash 1.6 x1a=vect(1) 169 y1a=vect(2) 170 z1a=vect(3)
171 brash 1.4 endif
172 brash 1.6 if(inwvol.eq.2) then 173 x1b=vect(1) 174 y1b=vect(2) 175 z1b=vect(3) 176 177 call get_wire_numbers(x1a,y1a,z1a,x1b,y1b,z1b,n1u,n1x,n1v) 178
179 brash 1.7 call get_drift_distance(x1a,y1a,z1a,x1b,y1b,z1b, 180 & n1u,n1x,n1v,d1u,d1x,d1v) 181
182 brash 1.6 write(,)'Hit in first chamber ...' 183 write(,)'Wire Numbers: ',n1u,n1x,n1v 184 185 endif 186 187 if ( istop.ne.0 ) then 188 make_hist=0 189 endif 190 else if ( names(nlevel).eq."fch2" ) then 191 c write(,)'In fch2 ... inwvol = ',inwvol 192 c write(,)'Z-value = ',vect(3) 193 if(inwvol.eq.1) then 194 x2a=vect(1) 195 y2a=vect(2) 196 z2a=vect(3) 197 endif 198 if(inwvol.eq.2) then 199 x2b=vect(1) 200 y2b=vect(2) 201 z2b=vect(3) 202 203 brash 1.6 call get_wire_numbers(x2a,y2a,z2a,x2b,y2b,z2b,n2u,n2x,n2v) 204
205 brash 1.8 call get_drift_distance(x2a,y2a,z2a,x2b,y2b,z2b, 206 & n2u,n2x,n2v,d2u,d2x,d2v) 207
208 brash 1.6 write(,)'Hit in second chamber ...' 209 write(,)'Wire Numbers: ',n2u,n2x,n2v 210 211 endif 212 c 213 if ( istop.ne.0 ) then 214 make_hist=0 215 endif 216 else if ( names(nlevel).eq."fch3" ) then 217 if(inwvol.eq.1) then 218 x3a=vect(1) 219 y3a=vect(2) 220 z3a=vect(3) 221 endif 222 if(inwvol.eq.2) then 223 x3b=vect(1) 224 y3b=vect(2) 225 z3b=vect(3) 226 227 call get_wire_numbers(x3a,y3a,z3a,x3b,y3b,z3b,n3u,n3x,n3v) 228
229 brash 1.8 call get_drift_distance(x3a,y3a,z3a,x3b,y3b,z3b, 230 & n3u,n3x,n3v,d3u,d3x,d3v) 231
232 brash 1.6 write(,)'Hit in third chamber ...' 233 write(,)'Wire Numbers: ',n3u,n3x,n3v 234 235 endif 236 c 237 if ( istop.ne.0 ) then 238 make_hist=0 239 endif 240 else if ( names(nlevel).eq."fch4" ) then 241 if(inwvol.eq.1) then 242 x4a=vect(1) 243 y4a=vect(2) 244 z4a=vect(3) 245 endif 246 if(inwvol.eq.2) then 247 x4b=vect(1) 248 y4b=vect(2) 249 z4b=vect(3) 250 251 call get_wire_numbers(x4a,y4a,z4a,x4b,y4b,z4b,n4u,n4x,n4v) 252
253 brash 1.8 call get_drift_distance(x4a,y4a,z4a,x4b,y4b,z4b, 254 & n4u,n4x,n4v,d4u,d4x,d4v) 255
256 brash 1.6 write(,)'Hit in fourth chamber ...' 257 write(,)'Wire Numbers: ',n4u,n4x,n4v 258 259 endif 260 c
261 brash 1.4 if ( istop.ne.0 ) then 262 make_hist=0 263 endif
264 brash 1.2 else if ( istop.ne.0.and. names(nlevel).eq."sci1" ) then 265 make_hist=0 266 else if ( names(nlevel).eq."anl1" ) then
267 jones 1.1 if(inwvol.eq.1) then
268 brash 1.5 c write(6,)'We have a hit in the fist analyzer at' 269 c write(6,)'x=',vect(1),' y=',vect(2),' z=',vect(3)
270 brash 1.4 xdet=vect(1) 271 ydet=vect(2) 272 zdet=vect(3)
273 jones 1.1 endif 274 if ( istop.ne.0 ) then 275 make_hist=0 276 endif 277 end if 278 c 279 c 280 c store current track parameters (including position ) in jxyz structure. 281 c 282 call gsxyz 283 c 284 c moved histograming stuff to gulast 285 c 286 c write(6,*)'done in gustep' 287 9999 return
288 brash 1.6 end 289 290 subroutine get_wire_numbers(xa,ya,za,xb,yb,zb,nu,nx,nv) 291
292 brash 1.7 implicit none 293
294 brash 1.6 real*8 xa,ya,za,xb,yb,zb
295 brash 1.7 integer*4 nu,nx,nv
296 brash 1.6 297 real*8 zc,zu,zx,zv,zt,xp,yp,xu,xx,xv,yu,yx,yv,uw,xw,vw
298 brash 1.7 real*8 anu,anx,anv
299 brash 1.6 300 c We have the (x,y,z) coordinates of the entrance (a) and exit (b) points 301 c of the track. We can use this information to calculate the wire numbers that 302 c were hit in each plane. 303 c 304 zc=(zb-za)/2.0+za 305 zu=zc-1.60 306 zx=zc 307 zv=zc+1.60 308 zt=(zb-za) 309 xp=(xb-xa)/zt 310 yp=(yb-ya)/zt
311 brash 1.7 xu=xa+xp(zu-za) 312 yu=ya+yp(zu-za) 313 xx=xa+xp(zx-za) 314 yx=ya+yp(zx-za) 315 xv=xa+xp(zv-za) 316 yv=ya+yp(zv-za)
317 brash 1.6 c 318 uw=(xu+yu)/sqrt(2.0) 319 xw=xx 320 vw=(-xv+yv)/sqrt(2.0)
321 brash 1.7 c
322 brash 1.8 c write(,)'*******************' 323 c write(,)xa,ya,za 324 c write(,)xb,yb,zb 325 c write(,)xu,yu,zu 326 c write(,)xx,yx,zx 327 c write(,)xv,yv,zv 328 c write(,)uw,xw,vw 329 c write(,)'*******************'
330 brash 1.7 c 331 anu=(-uw-3.592+104.0)/2.0 332 anx=(-xw-5.080+84.0)/2.0 333 anv=(vw-3.592+104.0)/2.0 334 c
335 brash 1.8 c write(,)anu,anx,anv
336 brash 1.7 nu=anu 337 nv=anv 338 nx=anx 339 if((anu-nu).ge.0.500)nu=nu+1 340 if((anx-nx).ge.0.500)nx=nx+1 341 if((anv-nv).ge.0.500)nv=nv+1 342 c 343 return 344 end 345 346 subroutine get_drift_distance(xa,ya,za,xb,yb,zb, 347 & nu,nx,nv,du,dx,dv) 348 349 implicit none 350 351 real8 xa,ya,za,xb,yb,zb,du,dx,dv 352 integer4 nu,nx,nv 353 real8 tphi,ttheta 354 real8 uw,xw,vw 355 real8 c11,c12,c21,c22,d1,d2,a,zt,xt,yt 356 real8 xu,yu,zu
357 brash 1.8 real8 xv,yv,zv 358 real8 xx,yx,zx
359 brash 1.7 360 tphi=(xb-xa)/(zb-za) 361 ttheta=(yb-ya)/(zb-za) 362 363 uw=-2.0nu-3.592+104.0 364 xw=-2.0nx-5.080+84.0 365 vw=2.0*nv+3.592-104.0 366 zu=(zb+za)/2.0-1.60
367 brash 1.8 zv=(zb+za)/2.0+1.60 368 zx=(zb+za)/2.0
369 brash 1.7
370 brash 1.8 c write(,)uw,xw,vw
371 brash 1.7 372 c11=tphi-ttheta 373 c12=sqrt(2.0) 374 d1=-xa+ya 375 c21=tphitphi+tthetattheta+1.0 376 c22=-ttheta/sqrt(2.0)+tphi/sqrt(2.0) 377 d2=uw(ttheta+tphi)/sqrt(2.0)-xatphi-yattheta+zu-za 378 379 a=(d1c21-d2c11)/(c21c12-c22c11) 380 zt=(d1-c12a)/c11
381 brash 1.8 c write(,)zt,a
382 brash 1.7 383 xt=xa+zttphi 384 yt=ya+ztttheta 385 xu=(uw-a)/sqrt(2.0) 386 yu=(uw+a)/sqrt(2.0) 387 388 zt=zt+za 389
390 brash 1.8 c write(,)'Calculating Drift Distance ...' 391 c write(,)'A: ',xa,ya,za 392 c write(,)'B: ',xb,yb,zb 393 c write(,)'T: ',xt,yt,zt 394 c write(,)'U: ',xu,yu,zu
395 brash 1.7 396 du=sqrt((xt-xu)2+(yt-yu)2+(zt-zu)**2)
397 brash 1.8 398 c write(,)'Drift distance = ',du 399 400 c11=tphi+ttheta 401 c12=-sqrt(2.0) 402 d1=-xa-ya 403 c21=tphitphi+tthetattheta+1.0 404 c22=-ttheta/sqrt(2.0)-tphi/sqrt(2.0) 405 d2=vw(ttheta+tphi)/sqrt(2.0)-xatphi-yattheta+zv-za 406 407 a=(d1c21-d2c11)/(c21c12-c22c11) 408 zt=(d1-c12a)/c11 409 c write(,)zt,a 410 411 xt=xa+zttphi 412 yt=ya+ztttheta 413 xv=-(vw-a)/sqrt(2.0) 414 yv=(vw+a)/sqrt(2.0) 415 416 zt=zt+za 417 418 brash 1.8 c write(,)'Calculating Drift Distance ...' 419 c write(,)'A: ',xa,ya,za 420 c write(,)'B: ',xb,yb,zb 421 c write(,)'T: ',xt,yt,zt 422 c write(,)'U: ',xv,yv,zv 423 424 dv=sqrt((xt-xv)2+(yt-yv)2+(zt-zv)*2) 425 426 c write(,)'Drift distance = ',dv 427 428 c11=tphi 429 c12=-1.0 430 d1=-ya 431 c21=tphitphi+tthetattheta+1.0 432 c22=-ttheta 433 d2=xwtphi-xatphi-yattheta+zx-za 434 435 a=(d1c21-d2c11)/(c21c12-c22c11) 436 zt=(d1-c12a)/c11 437 c write(,)zt,a 438 439 brash 1.8 xt=xa+zttphi 440 yt=ya+ztttheta 441 xx=xw 442 yx=a 443 444 zt=zt+za 445 446 c write(,)'Calculating Drift Distance ...' 447 c write(,)'A: ',xa,ya,za 448 c write(,)'B: ',xb,yb,zb 449 c write(,)'T: ',xt,yt,zt 450 c write(,)'U: ',xx,yx,zx 451 452 dx=sqrt((xt-xx)2+(yt-yx)2+(zt-zx)*2)
453 brash 1.7
454 brash 1.8 c write(,)'Drift distance = ',dx
455 brash 1.6 456 return
457 jones 1.1 end
458 brash 1.7 459 460 461 462 463 464 465
466 jones 1.1

Analyzer/Replay: Mark Jones, Documents: Stephen Wood