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

Analyzer/Replay: Mark Jones, Documents: Stephen Wood