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

Analyzer/Replay: Mark Jones, Documents: Stephen Wood