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

Analyzer/Replay: Mark Jones, Documents: Stephen Wood