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 brash 1.13 11 implicit none 12 13 integer4 mylast 14 integer4 nhu,nhv,nhx 15
16 jones 1.1 common/gctrak/vect(7),getot,gekin,vout(7),nmec,lmec(30),namec(30) 17 + ,nstep ,maxnst,destep,destel,safety,sleng 18 + ,step ,snext ,sfield 19 + ,tofg ,gekrat,upwght,ignext,inwvol,istop ,idecad,iekbin 20 + , ilosl, imull,ingoto,nldown,nlevin,nlvsav,istory 21 c 22 integer nmec,lmec,namec,nstep ,maxnst,ignext,inwvol,istop 23 + ,idecad,iekbin,ilosl, imull,ingoto,nldown,nlevin 24 + ,nlvsav,istory 25 real vect,getot,gekin,vout,destep,destel,safety,sleng ,step 26 + ,snext,sfield,tofg ,gekrat,upwght 27 c end gctrak 28 * keep,gcvolu. 29 -- author : 30 common/gcvolu/nlevel,names(15),number(15), 31 + lvolum(15),lindex(15),infrom,nlevmx,nldev(15),linmx(15), 32 + gtran(3,15),grmat(10,15),gonly(15),glx(3) 33 c 34 integer nlevel,number,lvolum,lindex,infrom,nlevmx, 35 + nldev,linmx 36 character4 names 37 jones 1.1 real gtran,grmat,gonly,glx 38 c end gcvolu 39 c 40 * keep,gcbank. 41 -- author : 42 integer iq,lq,nzebra,ixstor,ixdiv,ixcons,lmain,lr1 43 integer kwbank,kwwork,iws 44 real gversn,zversn,fendq,ws,q 45 c 46 parameter (kwbank=69000,kwwork=5200) 47 common/gcbank/nzebra,gversn,zversn,ixstor,ixdiv,ixcons,fendq(16) 48 + ,lmain,lr1,ws(kwbank) 49 dimension iq(2),q(2),lq(8000),iws(2) 50 equivalence (q(1),iq(1),lq(9)),(lq(1),lmain),(iws(1),ws(1)) 51 common/gclink/jdigi ,jdraw ,jhead ,jhits ,jkine ,jmate ,jpart 52 + ,jrotm ,jrung ,jset ,jstak ,jgstat,jtmed ,jtrack,jvertx 53 + ,jvolum,jxyz ,jgpar ,jgpar2,jsklt 54 c 55 integer jdigi ,jdraw ,jhead ,jhits ,jkine ,jmate ,jpart 56 + ,jrotm ,jrung ,jset ,jstak ,jgstat,jtmed ,jtrack,jvertx 57 + ,jvolum,jxyz ,jgpar,jgpar2 ,jsklt 58 jones 1.1 c 59 keep,gcking. 60 *-- author : 61 common/gcking/kcase,ngkine,gkin(5,100),tofd(100),iflgk(100) 62 integer kcase,ngkine ,iflgk 63 real gkin,tofd
64 brash 1.16 character*4 chcase
65 jones 1.1 c end gcking 66 c 67 * keep,gckine. 68 -- author : 69 -- author : 70 integer ikine,itra,istak,ivert,ipart,itrtyp,napart,ipaold 71 real pkine,amass,charge,tlife,vert,pvert 72 common/gckine/ikine,pkine(10),itra,istak,ivert,ipart,itrtyp 73 + ,napart(5),amass,charge,tlife,vert(3),pvert(4),ipaold 74 c end gckine 75 c 76 integer ihset,ihdet,iset,idet,idtype,nvname,numbv 77 common/gcsets/ihset,ihdet,iset,idet,idtype,nvname,numbv(20) 78 real x1,y1,z1,lpar,v1,v2,v3,newdist,x1new,y1new,z1new 79 real xstr,ystr,zstr,xstrnew,ystrnew,zstrnew 80 real rotmat2,rotmat3,rotmat4,rotmat1
81 brash 1.6 c
82 jones 1.1 common/geomstep/rotmat1(3,3),rotmat2(3,3),rotmat3(3,3), 83 & rotmat4(3,3) 84 c 85 include 'fpp_local.h' 86 include 'geant_local.h' 87 c include 'parameter.h' 88 c include 'espace_type.h' 89 c include 'detector.h' 90 c include 'transport.h' 91 c include 'option.h' 92 c 93 c 94 integer i,make_hist,ioff,ihit,ieffcheck 95 real pt_pi,ppar_pi,arg1,arg2,rapid_pi,pchmb,hits(6) 96 real a,b,c,beta,z,y,straw,ypath,rndm(3),ycompare
97 brash 1.14 logical idflag 98 real*8 d1uetemp,d1xetemp,d1vetemp
99 brash 1.2 c write(6,)'entering gustep' 100 c write(6,)'inwvol =',inwvol 101 c write(6,)'position =',vect(1),vect(2),vect(3) 102 c write(6,)'names =',names(nlevel)
103 jones 1.1 c 104 c 105 c 106 if ( ngkine.gt.0. ) then 107 do i=1,nmec 108 if ( lmec(i).eq.12 ) then 109 c write ( 6,* ) ' gustep: hadronic interaction' 110 c write ( 6,* ) ' nevent=',nevent 111 end if 112 end do 113 mylast = min(100,ngkine) 114 do i=1,mylast
115 brash 1.16 c write(,)'Secondaries Loop ...',i,' particle = ',gkin(5,i) 116 c if(gkin(5,i).eq.14)then 117 c write(,)'Total E =',gkin(4,i) 118 c write(,)'KE = ',gkin(4,i)-.93827 119 c elseif(gkin(5,i).eq.8.or.gkin(5,i).eq.9)then 120 c write(,)'Total E =',gkin(4,i) 121 c write(,)'KE = ',gkin(4,i)-.1395 122 c endif 123 if(gkin(5,i).eq.14.or.gkin(5,i).eq.8.or. 124 & gkin(5,i).eq.9)then 125 iflgk(i)=1 126 else 127 iflgk(i)=-1 128 endif 129 c iflgk(i) = 1 130 c if ( gkin(5,i).eq.4 ) iflgk(i) = 0 131 c if ( gkin(4,i).gt.0.001 ) iflgk(i)=0
132 jones 1.1 end do 133 c n_2nd = n_2nd + ngkine 134 endif 135 c 136 c the following call makes sure all of the secondary particles 137 c get tracked too (provided the flag iflgk(i) for that particle 138 c was set in the loop above -- this point is not correctly or 139 c clearly documented in the version of the geant manual i have). 140 c
141 brash 1.16 if(sectrack)then 142 c write(,)'SECONDARY TRACKING ...' 143 c write(,)'Number of secondaries = ',ngkine 144 c call uhtoc(kcase,4,chcase,4) 145 c write(,)'Source of interaction = ',chcase 146 c do i=1,ngkine 147 c write(,)'Secondary ',i,' ID =',gkin(5,i) 148 c if(gkin(5,i).eq.14.or.gkin(5,i).eq.8.or. 149 c & gkin(5,i).eq.9)then 150 c iflgk(i)=1 151 c else 152 c iflgk(i)=-1 153 c endif 154 c enddo 155 156 call gsking ( 0 ) 157 endif
158 jones 1.1 c
159 brash 1.6 c write(,)'In GUSTEP: ... names(nlevel) = ',names(nlevel)
160 brash 1.11 c write(,)'In GUSTEP: ... inwvol = ',inwvol
161 brash 1.6
162 jones 1.1 make_hist = 0 163 if ( inwvol.eq.1 .and. names(nlevel).eq."hall" ) then 164 make_hist=1 165 c 166 c if we get here, the tracking is done for this track 167 c (istop=1) but make a bunch of histograms before exitting 168 c note that ipart=8 means a pi+ and 9 is a pi- 169 c 170 c istop = 1 171 else if ( istop.ne.0.and. names(nlevel).eq."aira" ) then 172 make_hist=0 173 else if ( istop.ne.0.and. names(nlevel).eq."airb" ) then 174 make_hist=0 175 else if ( istop.ne.0.and. names(nlevel).eq."airc" ) then 176 make_hist=0 177 else if ( istop.ne.0.and. names(nlevel).eq."aird" ) then 178 make_hist=0
179 brash 1.2 else if ( istop.ne.0.and. names(nlevel).eq."aire" ) then 180 make_hist=0 181 else if ( istop.ne.0.and. names(nlevel).eq."airf" ) then 182 make_hist=0
183 brash 1.6 else if ( names(nlevel).eq."airg" ) then 184 c write(,)'In airg ... inwvol = ',inwvol 185 c write(,)'Z-value = ',vect(3) 186 if (istop.ne.0) then 187 make_hist=0 188 endif 189 else if ( names(nlevel).eq."HALL" ) then 190 c write(,)'In HALL ... inwvol = ',inwvol 191 c write(,)'Z-value = ',vect(3) 192 if (istop.ne.0) then 193 make_hist=0 194 endif
195 brash 1.2 else if ( istop.ne.0.and. names(nlevel).eq."airh" ) then 196 make_hist=0 197 else if ( istop.ne.0.and. names(nlevel).eq."hch1" ) then 198 make_hist=0
199 brash 1.11 else if ( names(nlevel).eq."hch2" ) then 200 c write(,)'In hch2' 201 if(inwvol.eq.1) then 202 c write(6,)'Coordinates at hch2' 203 c write(6,)'x=',vect(1),' y=',vect(2),' z=',vect(3) 204 xahch2=vect(1) 205 yahch2=vect(2) 206 zahch2=vect(3) 207 endif 208 if(inwvol.eq.2) then 209 xbhch2=vect(1) 210 ybhch2=vect(2) 211 zbhch2=vect(3) 212 endif 213 214 if ( istop.ne.0 ) then 215 make_hist=0 216 endif
217 brash 1.4 else if ( names(nlevel).eq."fch1" ) then
218 brash 1.6 c write(,)'In fch1 ... inwvol = ',inwvol
219 brash 1.4 if(inwvol.eq.1) then
220 brash 1.5 c write(6,)'Coordinates at fch1' 221 c write(6,)'x=',vect(1),' y=',vect(2),' z=',vect(3)
222 brash 1.6 x1a=vect(1) 223 y1a=vect(2) 224 z1a=vect(3)
225 brash 1.4 endif
226 brash 1.16 if(inwvol.eq.2.and.(ipart.eq.8.or.ipart.eq.9. 227 & or.ipart.eq.14)) then
228 brash 1.6 x1b=vect(1) 229 y1b=vect(2) 230 z1b=vect(3) 231
232 brash 1.13 call get_wire_numbers(x1a,y1a,z1a,x1b,y1b,z1b,nhu,nhx,nhv,
233 brash 1.12 & n1ua,n1xa,n1va,n1ub,n1xb,n1vb)
234 brash 1.6
235 brash 1.16 nhu1=nhu1+1 236 nhx1=nhx1+1 237 nhv1=nhv1+1
238 brash 1.13 239 if(nhu.eq.1.and.nhv.eq.1.and.nhx.eq.1) then
240 brash 1.14 idflag=.false.
241 brash 1.12 call get_drift_distance_ejb(x1a,y1a,z1a,x1b,y1b,z1b,
242 brash 1.14 & n1ua,n1xa,n1va,d1ue,d1xe,d1ve,idflag)
243 brash 1.16 c call get_drift_distance(x1a,y1a,z1a,x1b,y1b,z1b, 244 c & n1ua,n1xa,n1va,d1u,d1x,d1v)
245 brash 1.12 else 246 call get_drift_distance_ejb(x1a,y1a,z1a,x1b,y1b,z1b,
247 brash 1.14 & n1ua,n1xa,n1va,d1ue,d1xe,d1ve,idflag) 248 d1uetemp=d1ue 249 d1xetemp=d1xe 250 d1vetemp=d1ve
251 brash 1.12 call get_drift_distance_ejb(x1a,y1a,z1a,x1b,y1b,z1b,
252 brash 1.14 & n1ub,n1xb,n1vb,d1ue,d1xe,d1ve,idflag) 253 if(idflag) 254 & write(,)'Drift Distance 1a: ', 255 & d1uetemp,d1xetemp,d1vetemp 256 if(idflag) 257 & write(,)'Drift Distance 1b: ',d1ue,d1xe,d1ve
258 brash 1.16 c call get_drift_distance(x1a,y1a,z1a,x1b,y1b,z1b, 259 c & n1ua,n1xa,n1va,d1u,d1x,d1v) 260 c if(idflag) 261 c & write(,)'Drift Distance 1c: ',d1u,d1x,d1v 262 c call get_drift_distance(x1a,y1a,z1a,x1b,y1b,z1b, 263 c & n1ub,n1xb,n1vb,d1u,d1x,d1v) 264 c if(idflag) 265 c & write(,)'Drift Distance 1d: ',d1u,d1x,d1v
266 brash 1.12 endif
267 brash 1.11
268 brash 1.16 write(,)'Hit in first chamber ...' 269 write(,)'Number of Hits: ',nhu1,nhx1,nhv1
270 brash 1.12
271 brash 1.6 endif 272 273 if ( istop.ne.0 ) then 274 make_hist=0 275 endif 276 else if ( names(nlevel).eq."fch2" ) then 277 c write(,)'In fch2 ... inwvol = ',inwvol 278 c write(,)'Z-value = ',vect(3) 279 if(inwvol.eq.1) then 280 x2a=vect(1) 281 y2a=vect(2) 282 z2a=vect(3) 283 endif
284 brash 1.16 if(inwvol.eq.2.and.(ipart.eq.8.or.ipart.eq.9 285 & .or.ipart.eq.14)) then
286 brash 1.6 x2b=vect(1) 287 y2b=vect(2) 288 z2b=vect(3) 289
290 brash 1.13 call get_wire_numbers(x2a,y2a,z2a,x2b,y2b,z2b,nhu,nhx,nhv,
291 brash 1.12 & n2ua,n2xa,n2va,n2ub,n2xb,n2vb)
292 brash 1.11
293 brash 1.16 nhu2=nhu2+1 294 nhx2=nhx2+1 295 nhv2=nhv2+1
296 brash 1.13 297 if(nhu.eq.1.and.nhv.eq.1.and.nhx.eq.1) then
298 brash 1.12 call get_drift_distance_ejb(x2a,y2a,z2a,x2b,y2b,z2b,
299 brash 1.14 & n2ua,n2xa,n2va,d2ue,d2xe,d2ve,idflag)
300 brash 1.16 c call get_drift_distance(x2a,y2a,z2a,x2b,y2b,z2b, 301 c & n2ua,n2xa,n2va,d2u,d2x,d2v)
302 brash 1.12 else 303 call get_drift_distance_ejb(x2a,y2a,z2a,x2b,y2b,z2b,
304 brash 1.14 & n2ua,n2xa,n2va,d2ue,d2xe,d2ve,idflag) 305 c write(,)'Drift Distance 2a: ',d2ue,d2xe,d2ve
306 brash 1.12 call get_drift_distance_ejb(x2a,y2a,z2a,x2b,y2b,z2b,
307 brash 1.14 & n2ub,n2xb,n2vb,d2ue,d2xe,d2ve,idflag) 308 c write(,)'Drift Distance 2b: ',d2ue,d2xe,d2ve
309 brash 1.16 c call get_drift_distance(x2a,y2a,z2a,x2b,y2b,z2b, 310 c & n2ua,n2xa,n2va,d2u,d2x,d2v)
311 brash 1.14 c write(,)'Drift Distance 2c: ',d2u,d2x,d2v
312 brash 1.16 c call get_drift_distance(x2a,y2a,z2a,x2b,y2b,z2b, 313 c & n2ub,n2xb,n2vb,d2u,d2x,d2v)
314 brash 1.14 c write(,)'Drift Distance 2c: ',d2u,d2x,d2v
315 brash 1.12 endif
316 brash 1.15 317 call calc_theta_phi(xahch2,yahch2,zahch2, 318 & xbhch2,ybhch2,zbhch2, 319 & x1a,y1a,z1a,x2b,y2b,z2b, 320 & theta_front,phi_front)
321 brash 1.8
322 brash 1.16 write(,)'Hit in second chamber ...' 323 write(,)'Number of hits: ',nhu2,nhx2,nhv2
324 brash 1.6 325 endif 326 c 327 if ( istop.ne.0 ) then 328 make_hist=0 329 endif 330 else if ( names(nlevel).eq."fch3" ) then 331 if(inwvol.eq.1) then 332 x3a=vect(1) 333 y3a=vect(2) 334 z3a=vect(3) 335 endif
336 brash 1.16 if(inwvol.eq.2.and.(ipart.eq.8.or.ipart.eq.9 337 & .or.ipart.eq.14)) then
338 brash 1.6 x3b=vect(1) 339 y3b=vect(2) 340 z3b=vect(3) 341
342 brash 1.13 call get_wire_numbers(x3a,y3a,z3a,x3b,y3b,z3b,nhu,nhx,nhv,
343 brash 1.12 & n3ua,n3xa,n3va,n3ub,n3xb,n3vb)
344 brash 1.11
345 brash 1.16 nhu3=nhu3+1 346 nhx3=nhx3+1 347 nhv3=nhv3+1
348 brash 1.13 349 if(nhu.eq.1.and.nhv.eq.1.and.nhx.eq.1) then
350 brash 1.12 call get_drift_distance_ejb(x3a,y3a,z3a,x3b,y3b,z3b,
351 brash 1.14 & n3ua,n3xa,n3va,d3ue,d3xe,d3ve,idflag)
352 brash 1.16 c call get_drift_distance(x3a,y3a,z3a,x3b,y3b,z3b, 353 c & n3ua,n3xa,n3va,d3u,d3x,d3v)
354 brash 1.12 else 355 call get_drift_distance_ejb(x3a,y3a,z3a,x3b,y3b,z3b,
356 brash 1.14 & n3ua,n3xa,n3va,d3ue,d3xe,d3ve,idflag) 357 c write(,)'Drift Distance 3a: ',d3ue,d3xe,d3ve
358 brash 1.12 call get_drift_distance_ejb(x3a,y3a,z3a,x3b,y3b,z3b,
359 brash 1.14 & n3ub,n3xb,n3vb,d3ue,d3xe,d3ve,idflag) 360 c write(,)'Drift Distance 3b: ',d3ue,d3xe,d3ve
361 brash 1.16 c call get_drift_distance(x3a,y3a,z3a,x3b,y3b,z3b, 362 c & n3ua,n3xa,n3va,d3u,d3x,d3v)
363 brash 1.14 c write(,)'Drift Distance 3c: ',d3u,d3x,d3v
364 brash 1.12 endif 365
366 brash 1.8
367 brash 1.16 write(,)'Hit in third chamber ...' 368 write(,)'Number of Hits: ',nhu3,nhx3,nhv3
369 brash 1.6 370 endif 371 c 372 if ( istop.ne.0 ) then 373 make_hist=0 374 endif 375 else if ( names(nlevel).eq."fch4" ) then 376 if(inwvol.eq.1) then 377 x4a=vect(1) 378 y4a=vect(2) 379 z4a=vect(3) 380 endif
381 brash 1.16 if(inwvol.eq.2.and.(ipart.eq.8.or.ipart.eq.9 382 & .or.ipart.eq.14)) then
383 brash 1.6 x4b=vect(1) 384 y4b=vect(2) 385 z4b=vect(3) 386
387 brash 1.13 call get_wire_numbers(x4a,y4a,z4a,x4b,y4b,z4b,nhu,nhx,nhv,
388 brash 1.12 & n4ua,n4xa,n4va,n4ub,n4xb,n4vb)
389 brash 1.11
390 brash 1.16 nhu4=nhu4+1 391 nhx4=nhx4+1 392 nhv4=nhv4+1
393 brash 1.13 394 if(nhu.eq.1.and.nhv.eq.1.and.nhx.eq.1) then
395 brash 1.12 call get_drift_distance_ejb(x4a,y4a,z4a,x4b,y4b,z4b,
396 brash 1.14 & n4ua,n4xa,n4va,d4ue,d4xe,d4ve,idflag)
397 brash 1.16 c call get_drift_distance(x4a,y4a,z4a,x4b,y4b,z4b, 398 c & n4ua,n4xa,n4va,d4u,d4x,d4v)
399 brash 1.12 else 400 call get_drift_distance_ejb(x4a,y4a,z4a,x4b,y4b,z4b,
401 brash 1.14 & n4ua,n4xa,n4va,d4ue,d4xe,d4ve,idflag) 402 c write(,)'Drift Distance 4a: ',d4ue,d4xe,d4ve
403 brash 1.12 call get_drift_distance_ejb(x4a,y4a,z4a,x4b,y4b,z4b,
404 brash 1.14 & n4ub,n4xb,n4vb,d4ue,d4xe,d4ve,idflag) 405 c write(,)'Drift Distance 4b: ',d4ue,d4xe,d4ve
406 brash 1.16 c call get_drift_distance(x4a,y4a,z4a,x4b,y4b,z4b, 407 c & n4ua,n4xa,n4va,d4u,d4x,d4v)
408 brash 1.14 c write(,)'Drift Distance 4c: ',d4u,d4x,d4v
409 brash 1.12 endif 410
411 brash 1.8
412 brash 1.16 write(,)'Hit in fourth chamber ...' 413 write(,)'Wire Numbers: ',nhu4,nhx4,nhv4 414 c write(,)'Particle ID: ',ipart
415 brash 1.6 416 endif 417 c
418 brash 1.4 if ( istop.ne.0 ) then 419 make_hist=0 420 endif
421 brash 1.2 else if ( istop.ne.0.and. names(nlevel).eq."sci1" ) then 422 make_hist=0 423 else if ( names(nlevel).eq."anl1" ) then
424 jones 1.1 if(inwvol.eq.1) then
425 brash 1.5 c write(6,)'We have a hit in the fist analyzer at' 426 c write(6,)'x=',vect(1),' y=',vect(2),' z=',vect(3)
427 brash 1.4 xdet=vect(1) 428 ydet=vect(2) 429 zdet=vect(3)
430 jones 1.1 endif 431 if ( istop.ne.0 ) then 432 make_hist=0 433 endif 434 end if 435 c 436 c 437 c store current track parameters (including position ) in jxyz structure. 438 c 439 call gsxyz 440 c 441 c moved histograming stuff to gulast 442 c 443 c write(6,*)'done in gustep' 444 9999 return
445 brash 1.6 end 446
447 brash 1.13 subroutine get_wire_numbers(xa,ya,za,xb,yb,zb,nhu,nhx,nhv,
448 brash 1.12 % nu1,nx1,nv1,nu2,nx2,nv2)
449 brash 1.6
450 brash 1.7 implicit none 451
452 brash 1.6 real*8 xa,ya,za,xb,yb,zb
453 brash 1.12 integer4 nu1,nx1,nv1 454 integer4 nu2,nx2,nv2
455 brash 1.13 integer*4 nhu,nhx,nhv
456 brash 1.6 457 real*8 zc,zu,zx,zv,zt,xp,yp,xu,xx,xv,yu,yx,yv,uw,xw,vw
458 brash 1.7 real*8 anu,anx,anv
459 brash 1.6
460 brash 1.12 nu2=0 461 nx2=0 462 nv2=0 463 c
464 brash 1.6 c We have the (x,y,z) coordinates of the entrance (a) and exit (b) points 465 c of the track. We can use this information to calculate the wire numbers that 466 c were hit in each plane. 467 c 468 zc=(zb-za)/2.0+za 469 zu=zc-1.60 470 zx=zc 471 zv=zc+1.60 472 zt=(zb-za) 473 xp=(xb-xa)/zt 474 yp=(yb-ya)/zt
475 brash 1.12 c
476 brash 1.14 c Project to the FRONT of the "cell" associated with each plane. 477 c 478 xu=xa+xp(zu-za-0.8) 479 yu=ya+yp(zu-za-0.8) 480 xx=xa+xp(zx-za-0.8) 481 yx=ya+yp(zx-za-0.8) 482 xv=xa+xp(zv-za-0.8) 483 yv=ya+yp(zv-za-0.8) 484 c 485 c xu=xa 486 c yu=ya 487 c xx=xa 488 c yx=ya 489 c xv=xa 490 c yv=ya
491 brash 1.6 c 492 uw=(xu+yu)/sqrt(2.0) 493 xw=xx 494 vw=(-xv+yv)/sqrt(2.0)
495 brash 1.7 c
496 brash 1.8 c write(,)'********************'
497 brash 1.10 c write(,)'A: ',xa,ya,za 498 c write(,)'B: ',xb,yb,zb 499 c write(,)'U: ',xu,yu,zu 500 c write(,)'X: ',xx,yx,zx 501 c write(,)'V: ',xv,yv,zv 502 c write(,)'W: ',uw,xw,vw
503 brash 1.14 c write(,)'********************'
504 brash 1.7 c 505 anu=(-uw-3.592+104.0)/2.0 506 anx=(-xw-5.080+84.0)/2.0 507 anv=(vw-3.592+104.0)/2.0 508 c
509 brash 1.10 c write(,)'Wires: ',anu,anx,anv 510 c write(,)'********************'
511 brash 1.12 nu1=anu 512 nv1=anv 513 nx1=anx 514 if((anu-nu1).ge.0.500)nu1=nu1+1 515 if((anx-nx1).ge.0.500)nx1=nx1+1 516 if((anv-nv1).ge.0.500)nv1=nv1+1 517 c write(,)'using front of chamber: ',nu1,nx1,nv1 518 c
519 brash 1.14 c Now project to the BACK of the "cell" associated with each plane.
520 brash 1.12 c
521 brash 1.14 xu=xa+xp(zu-za+0.8) 522 yu=ya+yp(zu-za+0.8) 523 xx=xa+xp(zx-za+0.8) 524 yx=ya+yp(zx-za+0.8) 525 xv=xa+xp(zv-za+0.8) 526 yv=ya+yp(zv-za+0.8) 527 c 528 c xu=xb 529 c yu=yb 530 c xx=xb 531 c yx=yb 532 c xv=xb 533 c yv=yb
534 brash 1.12 c 535 uw=(xu+yu)/sqrt(2.0) 536 xw=xx 537 vw=(-xv+yv)/sqrt(2.0)
538 brash 1.7 c
539 brash 1.12 c write(,)'*******************' 540 c write(,)'A: ',xa,ya,za 541 c write(,)'B: ',xb,yb,zb 542 c write(,)'U: ',xu,yu,zu 543 c write(,)'X: ',xx,yx,zx 544 c write(,*)'V: ',xv,yv,zv
545 brash 1.10 c write(,)'W: ',uw,xw,vw 546 c write(,)'********************'
547 brash 1.12 c 548 anu=(-uw-3.592+104.0)/2.0 549 anx=(-xw-5.080+84.0)/2.0 550 anv=(vw-3.592+104.0)/2.0 551 c 552 c write(,)'Wires: ',anu,anx,anv 553 c write(,)'********************' 554 nu2=anu 555 nv2=anv 556 nx2=anx 557 if((anu-nu2).ge.0.500)nu2=nu2+1 558 if((anx-nx2).ge.0.500)nx2=nx2+1 559 if((anv-nv2).ge.0.500)nv2=nv2+1
560 brash 1.13 nhu=1 561 nhx=1 562 nhv=1 563 if(nu1.ne.nu2)nhu=2 564 if(nx1.ne.nx2)nhx=2 565 if(nv1.ne.nv2)nhv=2
566 brash 1.7
567 brash 1.12 return 568 end
569 brash 1.8
570 brash 1.12 c subroutine get_drift_distance_ejb(xa,ya,za,xb,yb,zb, 571 c & nu,nx,nv,du,dx,dv) 572 c 573 c implicit none 574 c 575 c real8 xa,ya,za,xb,yb,zb,du,dx,dv 576 c integer4 nu,nx,nv 577 c real8 tphi,ttheta 578 c real8 uw,xw,vw 579 c real8 c11,c12,c21,c22,d1,d2,a,zt,xt,yt 580 c real8 xu,yu,zu 581 c real8 xv,yv,zv 582 c real8 xx,yx,zx 583 c 584 c tphi=(xb-xa)/(zb-za) 585 c ttheta=(yb-ya)/(zb-za) 586 c 587 c uw=-2.0nu-3.592+104.0 588 c xw=-2.0nx-5.080+84.0 589 c vw=2.0nv+3.592-104.0 590 c zu=(zb+za)/2.0-1.60 591 brash 1.12 c zv=(zb+za)/2.0+1.60 592 c zx=(zb+za)/2.0 593 c 594 cc write(,)'W: ',uw,xw,vw 595 cc write(,)'******************' 596 c 597 c c11=tphi-ttheta 598 c c12=sqrt(2.0) 599 c d1=-xa+ya 600 c c21=tphitphi+tthetattheta+1.0 601 c c22=-ttheta/sqrt(2.0)+tphi/sqrt(2.0) 602 c d2=uw(ttheta+tphi)/sqrt(2.0)-xatphi-yattheta+zu-za 603 c 604 c a=(d1c21-d2c11)/(c21c12-c22c11) 605 c zt=(d1-c12a)/c11 606 cc write(,)'U Plane zt,a: ',zt,a 607 c 608 c xt=xa+zttphi 609 c yt=ya+ztttheta 610 c xu=(uw-a)/sqrt(2.0) 611 c yu=(uw+a)/sqrt(2.0) 612 brash 1.12 c 613 c zt=zt+za 614 c 615 c du=sqrt((xt-xu)2+(yt-yu)2+(zt-zu)2) 616 c 617 c c11=tphi+ttheta 618 c c12=-sqrt(2.0) 619 c d1=-xa-ya 620 c c21=tphitphi+tthetattheta+1.0 621 c c22=-ttheta/sqrt(2.0)-tphi/sqrt(2.0) 622 c d2=vw(ttheta+tphi)/sqrt(2.0)-xatphi-yattheta+zv-za 623 c 624 c a=(d1c21-d2c11)/(c21c12-c22c11) 625 c zt=(d1-c12a)/c11 626 cc write(,)'V Plane zt,a: ',zt,a 627 c 628 c xt=xa+zttphi 629 c yt=ya+ztttheta 630 c xv=-(vw-a)/sqrt(2.0) 631 c yv=(vw+a)/sqrt(2.0) 632 c 633 brash 1.12 c zt=zt+za 634 c 635 c dv=sqrt((xt-xv)2+(yt-yv)2+(zt-zv)2) 636 c 637 c c11=tphi 638 c c12=-1.0 639 c d1=-ya 640 c c21=tphitphi+tthetattheta+1.0 641 c c22=-ttheta 642 c d2=xwtphi-xatphi-yattheta+zx-za 643 c 644 c a=(d1c21-d2c11)/(c21c12-c22c11) 645 c zt=(d1-c12a)/c11 646 cc write(,)'X Plane zt,a: ',zt,a 647 c 648 c xt=xa+zttphi 649 c yt=ya+ztttheta 650 c xx=xw 651 c yx=a 652 c 653 c zt=zt+za 654 brash 1.12 c 655 c dx=sqrt((xt-xx)2+(yt-yx)2+(zt-zx)*2) 656 c 657 cc if((du.gt.1.00)) then
658 brash 1.11 c write(,)'Calculating Drift Distance ...' 659 c write(,)'A: ',xa,ya,za 660 c write(,)'B: ',xb,yb,zb
661 brash 1.12 cc write(,)'U: ',xu,yu,zu 662 cc write(,)'V: ',xv,yv,zv 663 cc write(,)'X: ',xx,yx,zx
664 brash 1.11 c write(,)'U Drift distance = ',du 665 c write(,)'V Drift distance = ',dv 666 c write(,)'X Drift distance = ',dx
667 brash 1.12 cc endif 668 c 669 c return 670 c end
671 brash 1.11 672 673 subroutine get_drift_distance(xa,ya,za,xb,yb,zb, 674 & nu,nx,nv,distu,distx,distv) 675 c 676 c This subroutine uses the entrance (a) and exit (b) 677 c points of a chamber to define the line of the track. 678 c It uses the wire number and wire direction to define 679 c the wire line. It then uses these to build parallel 680 c planes by computing a normal vector to both of the lines. 681 c Finally, it calculates the distance between these two 682 c planes which is the distance of shortest approach. 683 c 684 implicit none 685 c 686 real8 xa,ya,za,xb,yb,zb 687 real8 distu,distx,distv,xp,yp,zc 688 integer4 nu,nx,nv 689 real8 uw,xw,vw 690 real8 zt,xt,yt 691 real8 xu,yu,zu 692 brash 1.11 real8 xv,yv,zv 693 real8 xx,yx,zx 694 c 695 c the direction of each of the lines 696 c vect1=track vectu=u wire
697 brash 1.14 real8 vect1(1:3), vectu(1:3) 698 real8 vectx(1:3), vectv(1:3)
699 brash 1.11 c 700 c the normal vector to both lines
701 brash 1.14 real8 normu(1:3) 702 real8 normx(1:3) 703 real*8 normv(1:3)
704 brash 1.11 c 705 c the coefficients of the plane
706 brash 1.14 real8 pvectu(1:4) 707 real8 pvectx(1:4) 708 real*8 pvectv(1:4)
709 brash 1.11 c 710 vect1(1)=xb-xa 711 vect1(2)=yb-ya 712 vect1(3)=zb-za 713 c 714 zu=(zb+za)/2.0-1.60 715 zv=(zb+za)/2.0+1.60 716 zx=(zb+za)/2.0 717 c 718 c use line number to calculate distance relative to 719 c wire plane, then convert to x and y 720 c write(,)'nu nx nv',nu,nx,nv 721 uw=-2.0nu-3.592+104.0 722 xw=-2.0nx-5.080+84.0 723 vw=2.0*nv+3.592-104.0 724 xu=uw/sqrt(2.0) 725 yu=uw/sqrt(2.0) 726 xx=xw 727 yx=0 728 xv=-vw/sqrt(2.0) 729 yv=vw/sqrt(2.0)
730 brash 1.12 c write(,)'uw xu yu zu',uw,xu,yu,zu
731 brash 1.11 c write(,)'xw xx yx zx',xw,xx,yx,zx 732 c write(,)'vw xv yv zv',vw,xv,yv,zv 733 c 734 c define direction vector for wires, will be the same 735 c for each wire in a given plane, and is known 736 c for each plane 737 vectu(1)=1.0/sqrt(2.0) 738 vectu(2)=1.0/sqrt(2.0) 739 vectu(3)=0.0 740 vectx(1)=0.0 741 vectx(2)=1.0 742 vectx(3)=0.0 743 vectv(1)=-1.0/sqrt(2.0) 744 vectv(2)=1.0/sqrt(2.0) 745 vectv(3)=0.0 746 c
747 brash 1.12 c write(,)'distance calculations .....' 748 c write(,)xa,ya,za 749 c write(,)vect1(1),vect1(2),vect1(3) 750 c write(,)xu,yu,zu 751 c write(,)vectu(1),vectu(2),vectu(3) 752 c write(,)xx,yx,zx 753 c write(,)vectx(1),vectx(2),vectx(3) 754 c write(,)xv,yv,zv 755 c write(,)vectv(1),vectv(2),vectv(3) 756 c write(,)'distance calculations .....' 757 c
758 brash 1.11 c cross product 759 normu(1)=vect1(2)vectu(3)-vect1(3)vectu(2) 760 normu(2)=vect1(1)vectu(3)-vect1(3)vectu(1) 761 normu(3)=vect1(1)vectu(2)-vect1(2)vectu(1) 762 normx(1)=vect1(2)vectx(3)-vect1(3)vectx(2) 763 normx(2)=vect1(1)vectx(3)-vect1(3)vectx(1) 764 normx(3)=vect1(1)vectx(2)-vect1(2)vectx(1) 765 normv(1)=vect1(2)vectv(3)-vect1(3)vectv(2) 766 normv(2)=vect1(1)vectv(3)-vect1(3)vectv(1) 767 normv(3)=vect1(1)vectv(2)-vect1(2)vectv(1) 768 c 769 pvectu(1)=normu(1) 770 pvectu(2)=normu(2) 771 pvectu(3)=normu(3) 772 pvectu(4)=normu(1)(-xu)+normu(2)(-yu)+normu(3)(-zu) 773 pvectx(1)=normx(1) 774 pvectx(2)=normx(2) 775 pvectx(3)=normx(3) 776 pvectx(4)=normx(1)(-xx)+normx(2)(-yx)+normx(3)(-zx) 777 pvectv(1)=normv(1) 778 pvectv(2)=normv(2) 779 brash 1.11 pvectv(3)=normv(3) 780 pvectv(4)=normv(1)(-xv)+normv(2)(-yv)+normv(3)(-zv) 781 c 782 c distance formula 783 distu=(pvectu(1)xa+pvectu(2)ya+pvectu(3)za+pvectu(4)) 784 & /sqrt(normu(1)2+normu(2)2+normu(3)*2) 785 distx=(pvectx(1)xa+pvectx(2)ya+pvectx(3)za+pvectx(4)) 786 & /sqrt(normx(1)2+normx(2)2+normx(3)*2) 787 distv=(pvectv(1)xa+pvectv(2)ya+pvectv(3)za+pvectv(4)) 788 & /sqrt(normv(1)2+normv(2)2+normv(3)*2) 789 c write(,*)'Drift distance: ',distu, distx, distv
790 brash 1.14 c 791 c write(,)'Brads routine ....' 792 c write(,)xa,ya,za 793 c write(,)xb,yb,zb 794 c write(,)nu,nx,nv 795 c write(,)uw,xw,vw 796 c write(,)'Drift distance: ',distu, distx, distv 797
798 brash 1.11 return 799 end 800
801 brash 1.12 subroutine get_drift_distance_ejb(xa,ya,za,xb,yb,zb,
802 brash 1.14 & nu,nx,nv,distu,distx,distv,idflag)
803 brash 1.12 c 804 c Author: Ed Brash - December 15th, 2005 805 c Yet another attempt at a full drift distance calculation 806 c 807 implicit none 808 c 809 real8 xa,ya,za,xb,yb,zb 810 real8 distu,distx,distv,xp,yp,zc 811 integer4 nu,nx,nv 812 real8 uw,xw,vw 813 real8 zt,xt,yt 814 real8 xu,yu,zu 815 real8 xv,yv,zv 816 real8 xx,yx,zx
817 brash 1.14 logical idflag
818 brash 1.12 c 819 c the direction of each of the lines 820 c vect1=track vectu=u wire 821 real8 vect1(1:3), vectu(1:3) 822 real8 vectx(1:3), vectv(1:3) 823 824 c the difference vector between the defining points 825 real8 du(1:3),dx(1:3),dv(1:3) 826 c 827 c the normal vector to both lines 828 real8 normu(1:3) 829 real8 normx(1:3) 830 real8 normv(1:3) 831 real8 normumag,normxmag,normvmag 832 c 833 c the coefficients of the distance vector 834 real8 dvectu(1:4) 835 real8 dvectx(1:4) 836 real8 dvectv(1:4) 837 c
838 brash 1.14 idflag=.false.
839 brash 1.12 vect1(1)=xb-xa 840 vect1(2)=yb-ya 841 vect1(3)=zb-za 842 c 843 zu=(zb+za)/2.0-1.60 844 zv=(zb+za)/2.0+1.60 845 zx=(zb+za)/2.0 846 c 847 c use line number to calculate distance relative to 848 c wire plane, then convert to x and y 849 c write(,)'nu nx nv',nu,nx,nv 850 uw=-2.0nu-3.592+104.0 851 xw=-2.0nx-5.080+84.0 852 vw=2.0nv+3.592-104.0 853 xu=uw/sqrt(2.0) 854 yu=uw/sqrt(2.0) 855 xx=xw 856 yx=0 857 xv=-vw/sqrt(2.0) 858 yv=vw/sqrt(2.0) 859 c write(,)'uw xu yu zu',uw,xu,yu,zu 860 brash 1.12 c write(,)'xw xx yx zx',xw,xx,yx,zx 861 c write(,)'vw xv yv zv',vw,xv,yv,zv 862 c 863 c define direction vector for wires, will be the same 864 c for each wire in a given plane, and is known 865 c for each plane 866 vectu(1)=1.0/sqrt(2.0) 867 vectu(2)=-1.0/sqrt(2.0) 868 vectu(3)=0.0 869 vectx(1)=0.0 870 vectx(2)=1.0 871 vectx(3)=0.0 872 vectv(1)=1.0/sqrt(2.0) 873 vectv(2)=1.0/sqrt(2.0) 874 vectv(3)=0.0 875 c 876 c write(,)'distance calculations .....' 877 c write(,)xa,ya,za 878 c write(,)vect1(1),vect1(2),vect1(3) 879 c write(,)xu,yu,zu 880 c write(,)vectu(1),vectu(2),vectu(3) 881 brash 1.12 c write(,)xx,yx,zx 882 c write(,)vectx(1),vectx(2),vectx(3) 883 c write(,)xv,yv,zv 884 c write(,)vectv(1),vectv(2),vectv(3) 885 c write(,)'distance calculations .....' 886 c 887 c cross product 888 normu(1)=vect1(2)vectu(3)-vect1(3)vectu(2) 889 normu(2)=vect1(3)vectu(1)-vect1(1)vectu(3) 890 normu(3)=vect1(1)vectu(2)-vect1(2)vectu(1) 891 normx(1)=vect1(2)vectx(3)-vect1(3)vectx(2) 892 normx(2)=vect1(3)vectx(1)-vect1(1)vectx(3) 893 normx(3)=vect1(1)vectx(2)-vect1(2)vectx(1) 894 normv(1)=vect1(2)vectv(3)-vect1(3)vectv(2) 895 normv(2)=vect1(3)vectv(1)-vect1(1)vectv(3) 896 normv(3)=vect1(1)vectv(2)-vect1(2)vectv(1) 897 c write(,)normu(1),normu(2),normu(3) 898 c 899 normumag=sqrt(normu(1)2+normu(2)2+normu(3)2) 900 normxmag=sqrt(normx(1)2+normx(2)2+normx(3)2) 901 normvmag=sqrt(normv(1)2+normv(2)2+normv(3)*2) 902 brash 1.12 normu(1)=normu(1)/normumag 903 normu(2)=normu(2)/normumag 904 normu(3)=normu(3)/normumag
905 brash 1.14 normx(1)=normx(1)/normxmag 906 normx(2)=normx(2)/normxmag 907 normx(3)=normx(3)/normxmag 908 normv(1)=normv(1)/normvmag 909 normv(2)=normv(2)/normvmag 910 normv(3)=normv(3)/normvmag
911 brash 1.12 c write(,)normumag 912 c 913 du(1)=xa-xu 914 du(2)=ya-yu 915 du(3)=za-zu 916 dx(1)=xa-xx 917 dx(2)=ya-yx 918 dx(3)=za-zx 919 dv(1)=xa-xv 920 dv(2)=ya-yv 921 dv(3)=za-zv 922 c 923 c 924 c distance formula 925 distu=du(1)normu(1)+du(2)normu(2)+du(3)normu(3) 926 distx=dx(1)normx(1)+dx(2)normx(2)+dx(3)normx(3) 927 distv=dv(1)normv(1)+dv(2)normv(2)+dv(3)*normv(3) 928 c
929 brash 1.14 if(distu.gt.1.0.or.distx.gt.1.0.or.distv.gt.1.0) 930 & idflag=.true. 931 if(distu.gt.1.28.or.distx.gt.1.28.or.distv.gt.1.28)then 932 write(,)'Problem Child !!!' 933 write(,)'distance calculations .....' 934 write(,)xa,ya,za 935 c write(,)vect1(1),vect1(2),vect1(3) 936 write(,)xu,yu,zu 937 c write(,)vectu(1),vectu(2),vectu(3) 938 write(,)xx,yx,zx 939 c write(,)vectx(1),vectx(2),vectx(3) 940 write(,)xv,yv,zv 941 c write(,)vectv(1),vectv(2),vectv(3) 942 write(,)'normalization factors' 943 write(,)normu(1),normu(2),normu(3) 944 write(,)normumag 945 write(,)normx(1),normx(2),normx(3) 946 write(,)normxmag 947 write(,)normv(1),normv(2),normv(3) 948 write(,)normvmag 949 write(,)'Drift distance: ',distu, distx, distv 950 brash 1.14 endif
951 brash 1.12 c 952 return 953 end 954
955 brash 1.11 956 subroutine calc_theta_phi(xin1,yin1,zin1,xin2,yin2,zin2, 957 & xsc1,ysc1,zsc1,xsc2,ysc2,zsc2,theta,phi) 958 c 959 implicit none
960 brash 1.12 include 'fpp_local.h' 961 include 'geant_local.h'
962 brash 1.11 c 963 real8 xin1,yin1,zin1,xin2,yin2,zin2 964 real8 xsc1,ysc1,zsc1,xsc2,ysc2,zsc2,theta,phi 965 real*8 ftheta, fphi, fpsi
966 brash 1.12 real*8 lin,lout,theta_ejb,phi_ejb
967 brash 1.11 c 968 real invect(1:3) 969 real scvect(1:3) 970 real scvect2(1:3)
971 brash 1.12 real in(1:3) 972 real out(1:3) 973 real scat(1:3)
974 brash 1.11 c 975 invect(1)=xin2-xin1 976 invect(2)=yin2-yin1 977 invect(3)=zin2-zin1 978 scvect(1)=xsc2-xsc1 979 scvect(2)=ysc2-ysc1 980 scvect(3)=zsc2-zsc1
981 brash 1.12 c write(,)'INCOMING: ',invect(1),invect(2),invect(3) 982 c write(,)'SCATTERED: ',scvect(1),scvect(2),scvect(3) 983 c 984 c EJB calculation of theta and phi 985 c 986 in(1)=invect(1)/invect(3) 987 in(2)=invect(2)/invect(3) 988 in(3)=invect(3)/invect(3) 989 out(1)=scvect(1)/scvect(3) 990 out(2)=scvect(2)/scvect(3) 991 out(3)=scvect(3)/scvect(3) 992 lin=sqrt(in(1)2+in(2)2+in(3)2) 993 lout=sqrt(out(1)2+out(2)2+out(3)2) 994 scat(1)=out(1)-in(1) 995 scat(2)=out(2)-in(2) 996 scat(3)=out(3) 997 x_ejb=scat(1) 998 y_ejb=scat(2) 999 z_ejb=scat(3) 1000 if(scat(1).ge.0.0.and.scat(2).gt.0.0)then 1001 phi_ejb=atan(scat(1)/scat(2))57.2957795 1002 brash 1.12 else if(scat(1).ge.0.0.and.scat(2).lt.0.0)then 1003 phi_ejb=atan(scat(1)/scat(2))57.2957795+180.00 1004 else if(scat(1).le.0.0.and.scat(2).lt.0.0)then 1005 phi_ejb=atan(scat(1)/scat(2))57.2957795+180.00 1006 else if(scat(1).le.0.0.and.scat(2).gt.0.0)then 1007 phi_ejb=atan(scat(1)/scat(2))57.2957795+360.00 1008 endif 1009 c 1010 theta_ejb=acos((in(1)out(1)+in(2)out(2)+in(3)out(3))/ 1011 & (linlout))57.2957795 1012 c write(,)'EJB Incoming Vector = ',in(1),in(2),in(3) 1013 c write(,)'EJB Outgoing Vector = ',out(1),out(2),out(3) 1014 c write(,)'EJB Scattered Vector = ',scat(1),scat(2),scat(3) 1015 c write(,*)'EJB Thetas = ',theta_ejb,phi_ejb
1016 brash 1.11 c
1017 brash 1.12 c end EJB calculation 1018 c 1019 1020
1021 brash 1.11 ftheta=acos(invect(3)/sqrt(invect(1)2+invect(3)2)) 1022 fphi=acos(invect(3)/sqrt(invect(2)2+invect(3)2)) 1023 fpsi=acos(sqrt(invect(2)2+invect(3)2)/sqrt(invect(1)2 1024 & +invect(2)2+invect(3)**2)) 1025 c
1026 brash 1.12 c write(,)'ftheta, fphi, fpsi', 1027 c & ftheta57.296,fphi57.296,fpsi*57.296
1028 brash 1.11 scvect2(1)=scvect(1)cos(fpsi)-sin(fpsi)(scvect(2)sin(fphi) 1029 & +scvect(3)cos(fphi)) 1030 scvect2(2)=scvect(2)cos(fphi)-scvect(3)sin(fphi) 1031 scvect2(3)=scvect(1)sin(fpsi)+cos(fpsi)(scvect(2)sin(fphi) 1032 & +scvect(3)cos(fphi)) 1033 c
1034 brash 1.12 c write(,)'SCATTERED 2: ',scvect2(1),scvect2(2),scvect2(3)
1035 brash 1.11 theta=atan(sqrt(scvect2(1)2+scvect2(2)2)/scvect2(3))57.2957795 1036 phi=atan(scvect2(1)/scvect2(2))57.2957795 1037 if (scvect2(1).lt.0.0.and.scvect2(2).gt.0.0) 1038 & phi=phi+360.00 1039 if (scvect2(1).lt.0.0.and.scvect2(2).lt.0.0) 1040 & phi=phi+180.00 1041 if (scvect2(1).gt.0.0.and.scvect2(2).lt.0.0) 1042 & phi=phi+180.00
1043 brash 1.12 1044 c write(,)'Theta,phi =',theta,phi
1045 brash 1.15 theta=theta_ejb 1046 phi=phi_ejb
1047 brash 1.11 c 1048 return 1049 end 1050
1051 brash 1.7 1052 1053 1054 1055 1056 1057 1058
1059 jones 1.1

Analyzer/Replay: Mark Jones, Documents: Stephen Wood