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.22 call cone_test(xahch2,yahch2,zahch2, 366 & xbhch2,ybhch2,zbhch2,theta_front,phi_front, 367 & zclose_front,z2b) 368
369 brash 1.17 endif 370 371 c write(,)'Hit in second chamber ...' 372 c write(,)'Number of hits: ',nhu2,nhx2,nhv2
373 brash 1.6 374 endif 375 c 376 if ( istop.ne.0 ) then 377 make_hist=0 378 endif 379 else if ( names(nlevel).eq."fch3" ) then 380 if(inwvol.eq.1) then 381 x3a=vect(1) 382 y3a=vect(2) 383 z3a=vect(3) 384 endif
385 brash 1.16 if(inwvol.eq.2.and.(ipart.eq.8.or.ipart.eq.9 386 & .or.ipart.eq.14)) then
387 brash 1.6 x3b=vect(1) 388 y3b=vect(2) 389 z3b=vect(3) 390
391 brash 1.13 call get_wire_numbers(x3a,y3a,z3a,x3b,y3b,z3b,nhu,nhx,nhv,
392 brash 1.12 & n3ua,n3xa,n3va,n3ub,n3xb,n3vb)
393 brash 1.11
394 brash 1.16 nhu3=nhu3+1 395 nhx3=nhx3+1 396 nhv3=nhv3+1
397 brash 1.13 398 if(nhu.eq.1.and.nhv.eq.1.and.nhx.eq.1) then
399 brash 1.12 call get_drift_distance_ejb(x3a,y3a,z3a,x3b,y3b,z3b,
400 brash 1.14 & n3ua,n3xa,n3va,d3ue,d3xe,d3ve,idflag)
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.12 else 404 call get_drift_distance_ejb(x3a,y3a,z3a,x3b,y3b,z3b,
405 brash 1.14 & n3ua,n3xa,n3va,d3ue,d3xe,d3ve,idflag) 406 c write(,)'Drift Distance 3a: ',d3ue,d3xe,d3ve
407 brash 1.12 call get_drift_distance_ejb(x3a,y3a,z3a,x3b,y3b,z3b,
408 brash 1.14 & n3ub,n3xb,n3vb,d3ue,d3xe,d3ve,idflag) 409 c write(,)'Drift Distance 3b: ',d3ue,d3xe,d3ve
410 brash 1.16 c call get_drift_distance(x3a,y3a,z3a,x3b,y3b,z3b, 411 c & n3ua,n3xa,n3va,d3u,d3x,d3v)
412 brash 1.14 c write(,)'Drift Distance 3c: ',d3u,d3x,d3v
413 brash 1.12 endif 414
415 brash 1.8
416 brash 1.17 c write(,)'Hit in third chamber ...' 417 c write(,)'Number of Hits: ',nhu3,nhx3,nhv3
418 brash 1.6 419 endif 420 c 421 if ( istop.ne.0 ) then 422 make_hist=0 423 endif 424 else if ( names(nlevel).eq."fch4" ) then 425 if(inwvol.eq.1) then 426 x4a=vect(1) 427 y4a=vect(2) 428 z4a=vect(3) 429 endif
430 brash 1.16 if(inwvol.eq.2.and.(ipart.eq.8.or.ipart.eq.9 431 & .or.ipart.eq.14)) then
432 brash 1.6 x4b=vect(1) 433 y4b=vect(2) 434 z4b=vect(3) 435
436 brash 1.13 call get_wire_numbers(x4a,y4a,z4a,x4b,y4b,z4b,nhu,nhx,nhv,
437 brash 1.12 & n4ua,n4xa,n4va,n4ub,n4xb,n4vb)
438 brash 1.11
439 brash 1.16 nhu4=nhu4+1 440 nhx4=nhx4+1 441 nhv4=nhv4+1
442 brash 1.13 443 if(nhu.eq.1.and.nhv.eq.1.and.nhx.eq.1) then
444 brash 1.12 call get_drift_distance_ejb(x4a,y4a,z4a,x4b,y4b,z4b,
445 brash 1.14 & n4ua,n4xa,n4va,d4ue,d4xe,d4ve,idflag)
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.12 else 449 call get_drift_distance_ejb(x4a,y4a,z4a,x4b,y4b,z4b,
450 brash 1.14 & n4ua,n4xa,n4va,d4ue,d4xe,d4ve,idflag) 451 c write(,)'Drift Distance 4a: ',d4ue,d4xe,d4ve
452 brash 1.12 call get_drift_distance_ejb(x4a,y4a,z4a,x4b,y4b,z4b,
453 brash 1.14 & n4ub,n4xb,n4vb,d4ue,d4xe,d4ve,idflag) 454 c write(,)'Drift Distance 4b: ',d4ue,d4xe,d4ve
455 brash 1.16 c call get_drift_distance(x4a,y4a,z4a,x4b,y4b,z4b, 456 c & n4ua,n4xa,n4va,d4u,d4x,d4v)
457 brash 1.14 c write(,)'Drift Distance 4c: ',d4u,d4x,d4v
458 brash 1.12 endif 459
460 brash 1.18 call calc_theta_phi(x1a,y1a,z1a,x2b,y2b,z2b, 461 & x3a,y3a,z3a,x4b,y4b,z4b, 462 & theta_temp,phi_temp) 463
464 brash 1.19 c write(,)'Raw Variables' 465 c write(,)x1a,y1a,z1a 466 c write(,)x2b,y2b,z2b 467 c write(,)x3a,y3a,z3a 468 c write(,)x4b,y4b,z4b 469 470 tphif=(x2b-x1a)/(z2b-z1a) 471 tthetaf=(y2b-y1a)/(z2b-z1a) 472 tphir=(x4b-x3a)/(z4b-z3a) 473 tthetar=(y4b-y3a)/(z4b-z3a) 474 zmid=(z3a+z2b)/2.0 475 x0f=x2b+(zmid-z2b)tphif 476 y0f=y2b+(zmid-z2b)tthetaf 477 x0r=x3a+(zmid-z3a)tphir 478 y0r=y3a+(zmid-z3a)tthetar 479 480 c write(,)'Rear Scattering:' 481 c write(,)x0f,y0f,tphif,tthetaf 482 c write(,)x0r,y0r,tphir,tthetar 483
484 brash 1.21 call calc_zclose_sclose(zmid,z4b, 485 & x0f,y0f,tphif,tthetaf, 486 & x0r,y0r,tphir,tthetar, 487 & theta_temp,phi_temp, 488 & zclose_temp,sclose_temp,cone_temp) 489
490 brash 1.19 c write(,)'zclose = ',zclose_temp 491 c write(,)'sclose = ',sclose_temp 492 493 c write(,)'Old scattering angle = ',theta_rear 494 c write(,)'New scattering angle = ',theta_temp 495 c write(,)'**************************'
496 brash 1.18 if(theta_temp.lt.theta_rear)then 497 theta_rear=theta_temp 498 phi_rear=phi_temp
499 brash 1.19 zclose_rear=zclose_temp 500 sclose_rear=sclose_temp
501 brash 1.21 icone_rear=cone_temp
502 brash 1.18 endif 503
504 brash 1.8
505 brash 1.17 c write(,)'Hit in fourth chamber ...' 506 c write(,)'Wire Numbers: ',nhu4,nhx4,nhv4
507 brash 1.16 c write(,)'Particle ID: ',ipart
508 brash 1.6 509 endif 510 c
511 brash 1.4 if ( istop.ne.0 ) then 512 make_hist=0 513 endif
514 brash 1.2 else if ( istop.ne.0.and. names(nlevel).eq."sci1" ) then 515 make_hist=0 516 else if ( names(nlevel).eq."anl1" ) then
517 jones 1.1 if(inwvol.eq.1) then
518 brash 1.5 c write(6,)'We have a hit in the fist analyzer at' 519 c write(6,)'x=',vect(1),' y=',vect(2),' z=',vect(3)
520 brash 1.4 xdet=vect(1) 521 ydet=vect(2) 522 zdet=vect(3)
523 jones 1.1 endif 524 if ( istop.ne.0 ) then 525 make_hist=0 526 endif 527 end if 528 c 529 c 530 c store current track parameters (including position ) in jxyz structure. 531 c 532 call gsxyz 533 c 534 c moved histograming stuff to gulast 535 c 536 c write(6,*)'done in gustep' 537 9999 return
538 brash 1.6 end 539
540 brash 1.13 subroutine get_wire_numbers(xa,ya,za,xb,yb,zb,nhu,nhx,nhv,
541 brash 1.12 % nu1,nx1,nv1,nu2,nx2,nv2)
542 brash 1.6
543 brash 1.7 implicit none 544
545 brash 1.6 real*8 xa,ya,za,xb,yb,zb
546 brash 1.12 integer4 nu1,nx1,nv1 547 integer4 nu2,nx2,nv2
548 brash 1.13 integer*4 nhu,nhx,nhv
549 brash 1.6 550 real*8 zc,zu,zx,zv,zt,xp,yp,xu,xx,xv,yu,yx,yv,uw,xw,vw
551 brash 1.7 real*8 anu,anx,anv
552 brash 1.6
553 brash 1.12 nu2=0 554 nx2=0 555 nv2=0 556 c
557 brash 1.6 c We have the (x,y,z) coordinates of the entrance (a) and exit (b) points 558 c of the track. We can use this information to calculate the wire numbers that 559 c were hit in each plane. 560 c 561 zc=(zb-za)/2.0+za 562 zu=zc-1.60 563 zx=zc 564 zv=zc+1.60 565 zt=(zb-za) 566 xp=(xb-xa)/zt 567 yp=(yb-ya)/zt
568 brash 1.12 c
569 brash 1.14 c Project to the FRONT of the "cell" associated with each plane. 570 c 571 xu=xa+xp(zu-za-0.8) 572 yu=ya+yp(zu-za-0.8) 573 xx=xa+xp(zx-za-0.8) 574 yx=ya+yp(zx-za-0.8) 575 xv=xa+xp(zv-za-0.8) 576 yv=ya+yp(zv-za-0.8) 577 c 578 c xu=xa 579 c yu=ya 580 c xx=xa 581 c yx=ya 582 c xv=xa 583 c yv=ya
584 brash 1.6 c 585 uw=(xu+yu)/sqrt(2.0) 586 xw=xx 587 vw=(-xv+yv)/sqrt(2.0)
588 brash 1.7 c
589 brash 1.8 c write(,)'********************'
590 brash 1.10 c write(,)'A: ',xa,ya,za 591 c write(,)'B: ',xb,yb,zb 592 c write(,)'U: ',xu,yu,zu 593 c write(,)'X: ',xx,yx,zx 594 c write(,)'V: ',xv,yv,zv 595 c write(,)'W: ',uw,xw,vw
596 brash 1.14 c write(,)'********************'
597 brash 1.7 c 598 anu=(-uw-3.592+104.0)/2.0 599 anx=(-xw-5.080+84.0)/2.0 600 anv=(vw-3.592+104.0)/2.0 601 c
602 brash 1.10 c write(,)'Wires: ',anu,anx,anv 603 c write(,)'********************'
604 brash 1.12 nu1=anu 605 nv1=anv 606 nx1=anx 607 if((anu-nu1).ge.0.500)nu1=nu1+1 608 if((anx-nx1).ge.0.500)nx1=nx1+1 609 if((anv-nv1).ge.0.500)nv1=nv1+1 610 c write(,)'using front of chamber: ',nu1,nx1,nv1 611 c
612 brash 1.14 c Now project to the BACK of the "cell" associated with each plane.
613 brash 1.12 c
614 brash 1.14 xu=xa+xp(zu-za+0.8) 615 yu=ya+yp(zu-za+0.8) 616 xx=xa+xp(zx-za+0.8) 617 yx=ya+yp(zx-za+0.8) 618 xv=xa+xp(zv-za+0.8) 619 yv=ya+yp(zv-za+0.8) 620 c 621 c xu=xb 622 c yu=yb 623 c xx=xb 624 c yx=yb 625 c xv=xb 626 c yv=yb
627 brash 1.12 c 628 uw=(xu+yu)/sqrt(2.0) 629 xw=xx 630 vw=(-xv+yv)/sqrt(2.0)
631 brash 1.7 c
632 brash 1.12 c write(,)'*******************' 633 c write(,)'A: ',xa,ya,za 634 c write(,)'B: ',xb,yb,zb 635 c write(,)'U: ',xu,yu,zu 636 c write(,)'X: ',xx,yx,zx 637 c write(,*)'V: ',xv,yv,zv
638 brash 1.10 c write(,)'W: ',uw,xw,vw 639 c write(,)'********************'
640 brash 1.12 c 641 anu=(-uw-3.592+104.0)/2.0 642 anx=(-xw-5.080+84.0)/2.0 643 anv=(vw-3.592+104.0)/2.0 644 c 645 c write(,)'Wires: ',anu,anx,anv 646 c write(,)'********************' 647 nu2=anu 648 nv2=anv 649 nx2=anx 650 if((anu-nu2).ge.0.500)nu2=nu2+1 651 if((anx-nx2).ge.0.500)nx2=nx2+1 652 if((anv-nv2).ge.0.500)nv2=nv2+1
653 brash 1.13 nhu=1 654 nhx=1 655 nhv=1 656 if(nu1.ne.nu2)nhu=2 657 if(nx1.ne.nx2)nhx=2 658 if(nv1.ne.nv2)nhv=2
659 brash 1.7
660 brash 1.12 return 661 end
662 brash 1.8
663 brash 1.12 c subroutine get_drift_distance_ejb(xa,ya,za,xb,yb,zb, 664 c & nu,nx,nv,du,dx,dv) 665 c 666 c implicit none 667 c 668 c real8 xa,ya,za,xb,yb,zb,du,dx,dv 669 c integer4 nu,nx,nv 670 c real8 tphi,ttheta 671 c real8 uw,xw,vw 672 c real8 c11,c12,c21,c22,d1,d2,a,zt,xt,yt 673 c real8 xu,yu,zu 674 c real8 xv,yv,zv 675 c real8 xx,yx,zx 676 c 677 c tphi=(xb-xa)/(zb-za) 678 c ttheta=(yb-ya)/(zb-za) 679 c 680 c uw=-2.0nu-3.592+104.0 681 c xw=-2.0nx-5.080+84.0 682 c vw=2.0nv+3.592-104.0 683 c zu=(zb+za)/2.0-1.60 684 brash 1.12 c zv=(zb+za)/2.0+1.60 685 c zx=(zb+za)/2.0 686 c 687 cc write(,)'W: ',uw,xw,vw 688 cc write(,)'******************' 689 c 690 c c11=tphi-ttheta 691 c c12=sqrt(2.0) 692 c d1=-xa+ya 693 c c21=tphitphi+tthetattheta+1.0 694 c c22=-ttheta/sqrt(2.0)+tphi/sqrt(2.0) 695 c d2=uw(ttheta+tphi)/sqrt(2.0)-xatphi-yattheta+zu-za 696 c 697 c a=(d1c21-d2c11)/(c21c12-c22c11) 698 c zt=(d1-c12a)/c11 699 cc write(,)'U Plane zt,a: ',zt,a 700 c 701 c xt=xa+zttphi 702 c yt=ya+ztttheta 703 c xu=(uw-a)/sqrt(2.0) 704 c yu=(uw+a)/sqrt(2.0) 705 brash 1.12 c 706 c zt=zt+za 707 c 708 c du=sqrt((xt-xu)2+(yt-yu)2+(zt-zu)2) 709 c 710 c c11=tphi+ttheta 711 c c12=-sqrt(2.0) 712 c d1=-xa-ya 713 c c21=tphitphi+tthetattheta+1.0 714 c c22=-ttheta/sqrt(2.0)-tphi/sqrt(2.0) 715 c d2=vw(ttheta+tphi)/sqrt(2.0)-xatphi-yattheta+zv-za 716 c 717 c a=(d1c21-d2c11)/(c21c12-c22c11) 718 c zt=(d1-c12a)/c11 719 cc write(,)'V Plane zt,a: ',zt,a 720 c 721 c xt=xa+zttphi 722 c yt=ya+ztttheta 723 c xv=-(vw-a)/sqrt(2.0) 724 c yv=(vw+a)/sqrt(2.0) 725 c 726 brash 1.12 c zt=zt+za 727 c 728 c dv=sqrt((xt-xv)2+(yt-yv)2+(zt-zv)2) 729 c 730 c c11=tphi 731 c c12=-1.0 732 c d1=-ya 733 c c21=tphitphi+tthetattheta+1.0 734 c c22=-ttheta 735 c d2=xwtphi-xatphi-yattheta+zx-za 736 c 737 c a=(d1c21-d2c11)/(c21c12-c22c11) 738 c zt=(d1-c12a)/c11 739 cc write(,)'X Plane zt,a: ',zt,a 740 c 741 c xt=xa+zttphi 742 c yt=ya+ztttheta 743 c xx=xw 744 c yx=a 745 c 746 c zt=zt+za 747 brash 1.12 c 748 c dx=sqrt((xt-xx)2+(yt-yx)2+(zt-zx)*2) 749 c 750 cc if((du.gt.1.00)) then
751 brash 1.11 c write(,)'Calculating Drift Distance ...' 752 c write(,)'A: ',xa,ya,za 753 c write(,)'B: ',xb,yb,zb
754 brash 1.12 cc write(,)'U: ',xu,yu,zu 755 cc write(,)'V: ',xv,yv,zv 756 cc write(,)'X: ',xx,yx,zx
757 brash 1.11 c write(,)'U Drift distance = ',du 758 c write(,)'V Drift distance = ',dv 759 c write(,)'X Drift distance = ',dx
760 brash 1.12 cc endif 761 c 762 c return 763 c end
764 brash 1.11 765 766 subroutine get_drift_distance(xa,ya,za,xb,yb,zb, 767 & nu,nx,nv,distu,distx,distv) 768 c 769 c This subroutine uses the entrance (a) and exit (b) 770 c points of a chamber to define the line of the track. 771 c It uses the wire number and wire direction to define 772 c the wire line. It then uses these to build parallel 773 c planes by computing a normal vector to both of the lines. 774 c Finally, it calculates the distance between these two 775 c planes which is the distance of shortest approach. 776 c 777 implicit none 778 c 779 real8 xa,ya,za,xb,yb,zb 780 real8 distu,distx,distv,xp,yp,zc 781 integer4 nu,nx,nv 782 real8 uw,xw,vw 783 real8 zt,xt,yt 784 real8 xu,yu,zu 785 brash 1.11 real8 xv,yv,zv 786 real8 xx,yx,zx 787 c 788 c the direction of each of the lines 789 c vect1=track vectu=u wire
790 brash 1.14 real8 vect1(1:3), vectu(1:3) 791 real8 vectx(1:3), vectv(1:3)
792 brash 1.11 c 793 c the normal vector to both lines
794 brash 1.14 real8 normu(1:3) 795 real8 normx(1:3) 796 real*8 normv(1:3)
797 brash 1.11 c 798 c the coefficients of the plane
799 brash 1.14 real8 pvectu(1:4) 800 real8 pvectx(1:4) 801 real*8 pvectv(1:4)
802 brash 1.11 c 803 vect1(1)=xb-xa 804 vect1(2)=yb-ya 805 vect1(3)=zb-za 806 c 807 zu=(zb+za)/2.0-1.60 808 zv=(zb+za)/2.0+1.60 809 zx=(zb+za)/2.0 810 c 811 c use line number to calculate distance relative to 812 c wire plane, then convert to x and y 813 c write(,)'nu nx nv',nu,nx,nv 814 uw=-2.0nu-3.592+104.0 815 xw=-2.0nx-5.080+84.0 816 vw=2.0*nv+3.592-104.0 817 xu=uw/sqrt(2.0) 818 yu=uw/sqrt(2.0) 819 xx=xw 820 yx=0 821 xv=-vw/sqrt(2.0) 822 yv=vw/sqrt(2.0)
823 brash 1.12 c write(,)'uw xu yu zu',uw,xu,yu,zu
824 brash 1.11 c write(,)'xw xx yx zx',xw,xx,yx,zx 825 c write(,)'vw xv yv zv',vw,xv,yv,zv 826 c 827 c define direction vector for wires, will be the same 828 c for each wire in a given plane, and is known 829 c for each plane 830 vectu(1)=1.0/sqrt(2.0) 831 vectu(2)=1.0/sqrt(2.0) 832 vectu(3)=0.0 833 vectx(1)=0.0 834 vectx(2)=1.0 835 vectx(3)=0.0 836 vectv(1)=-1.0/sqrt(2.0) 837 vectv(2)=1.0/sqrt(2.0) 838 vectv(3)=0.0 839 c
840 brash 1.12 c write(,)'distance calculations .....' 841 c write(,)xa,ya,za 842 c write(,)vect1(1),vect1(2),vect1(3) 843 c write(,)xu,yu,zu 844 c write(,)vectu(1),vectu(2),vectu(3) 845 c write(,)xx,yx,zx 846 c write(,)vectx(1),vectx(2),vectx(3) 847 c write(,)xv,yv,zv 848 c write(,)vectv(1),vectv(2),vectv(3) 849 c write(,)'distance calculations .....' 850 c
851 brash 1.11 c cross product 852 normu(1)=vect1(2)vectu(3)-vect1(3)vectu(2) 853 normu(2)=vect1(1)vectu(3)-vect1(3)vectu(1) 854 normu(3)=vect1(1)vectu(2)-vect1(2)vectu(1) 855 normx(1)=vect1(2)vectx(3)-vect1(3)vectx(2) 856 normx(2)=vect1(1)vectx(3)-vect1(3)vectx(1) 857 normx(3)=vect1(1)vectx(2)-vect1(2)vectx(1) 858 normv(1)=vect1(2)vectv(3)-vect1(3)vectv(2) 859 normv(2)=vect1(1)vectv(3)-vect1(3)vectv(1) 860 normv(3)=vect1(1)vectv(2)-vect1(2)vectv(1) 861 c 862 pvectu(1)=normu(1) 863 pvectu(2)=normu(2) 864 pvectu(3)=normu(3) 865 pvectu(4)=normu(1)(-xu)+normu(2)(-yu)+normu(3)(-zu) 866 pvectx(1)=normx(1) 867 pvectx(2)=normx(2) 868 pvectx(3)=normx(3) 869 pvectx(4)=normx(1)(-xx)+normx(2)(-yx)+normx(3)(-zx) 870 pvectv(1)=normv(1) 871 pvectv(2)=normv(2) 872 brash 1.11 pvectv(3)=normv(3) 873 pvectv(4)=normv(1)(-xv)+normv(2)(-yv)+normv(3)(-zv) 874 c 875 c distance formula 876 distu=(pvectu(1)xa+pvectu(2)ya+pvectu(3)za+pvectu(4)) 877 & /sqrt(normu(1)2+normu(2)2+normu(3)*2) 878 distx=(pvectx(1)xa+pvectx(2)ya+pvectx(3)za+pvectx(4)) 879 & /sqrt(normx(1)2+normx(2)2+normx(3)*2) 880 distv=(pvectv(1)xa+pvectv(2)ya+pvectv(3)za+pvectv(4)) 881 & /sqrt(normv(1)2+normv(2)2+normv(3)*2) 882 c write(,*)'Drift distance: ',distu, distx, distv
883 brash 1.14 c 884 c write(,)'Brads routine ....' 885 c write(,)xa,ya,za 886 c write(,)xb,yb,zb 887 c write(,)nu,nx,nv 888 c write(,)uw,xw,vw 889 c write(,)'Drift distance: ',distu, distx, distv 890
891 brash 1.11 return 892 end 893
894 brash 1.12 subroutine get_drift_distance_ejb(xa,ya,za,xb,yb,zb,
895 brash 1.14 & nu,nx,nv,distu,distx,distv,idflag)
896 brash 1.12 c 897 c Author: Ed Brash - December 15th, 2005 898 c Yet another attempt at a full drift distance calculation 899 c 900 implicit none 901 c 902 real8 xa,ya,za,xb,yb,zb 903 real8 distu,distx,distv,xp,yp,zc 904 integer4 nu,nx,nv 905 real8 uw,xw,vw 906 real8 zt,xt,yt 907 real8 xu,yu,zu 908 real8 xv,yv,zv 909 real8 xx,yx,zx
910 brash 1.14 logical idflag
911 brash 1.12 c 912 c the direction of each of the lines 913 c vect1=track vectu=u wire 914 real8 vect1(1:3), vectu(1:3) 915 real8 vectx(1:3), vectv(1:3) 916 917 c the difference vector between the defining points 918 real8 du(1:3),dx(1:3),dv(1:3) 919 c 920 c the normal vector to both lines 921 real8 normu(1:3) 922 real8 normx(1:3) 923 real8 normv(1:3) 924 real8 normumag,normxmag,normvmag 925 c 926 c the coefficients of the distance vector 927 real8 dvectu(1:4) 928 real8 dvectx(1:4) 929 real8 dvectv(1:4) 930 c
931 brash 1.14 idflag=.false.
932 brash 1.12 vect1(1)=xb-xa 933 vect1(2)=yb-ya 934 vect1(3)=zb-za 935 c 936 zu=(zb+za)/2.0-1.60 937 zv=(zb+za)/2.0+1.60 938 zx=(zb+za)/2.0 939 c 940 c use line number to calculate distance relative to 941 c wire plane, then convert to x and y 942 c write(,)'nu nx nv',nu,nx,nv 943 uw=-2.0nu-3.592+104.0 944 xw=-2.0nx-5.080+84.0 945 vw=2.0nv+3.592-104.0 946 xu=uw/sqrt(2.0) 947 yu=uw/sqrt(2.0) 948 xx=xw 949 yx=0 950 xv=-vw/sqrt(2.0) 951 yv=vw/sqrt(2.0) 952 c write(,)'uw xu yu zu',uw,xu,yu,zu 953 brash 1.12 c write(,)'xw xx yx zx',xw,xx,yx,zx 954 c write(,)'vw xv yv zv',vw,xv,yv,zv 955 c 956 c define direction vector for wires, will be the same 957 c for each wire in a given plane, and is known 958 c for each plane 959 vectu(1)=1.0/sqrt(2.0) 960 vectu(2)=-1.0/sqrt(2.0) 961 vectu(3)=0.0 962 vectx(1)=0.0 963 vectx(2)=1.0 964 vectx(3)=0.0 965 vectv(1)=1.0/sqrt(2.0) 966 vectv(2)=1.0/sqrt(2.0) 967 vectv(3)=0.0 968 c 969 c write(,)'distance calculations .....' 970 c write(,)xa,ya,za 971 c write(,)vect1(1),vect1(2),vect1(3) 972 c write(,)xu,yu,zu 973 c write(,)vectu(1),vectu(2),vectu(3) 974 brash 1.12 c write(,)xx,yx,zx 975 c write(,)vectx(1),vectx(2),vectx(3) 976 c write(,)xv,yv,zv 977 c write(,)vectv(1),vectv(2),vectv(3) 978 c write(,)'distance calculations .....' 979 c 980 c cross product 981 normu(1)=vect1(2)vectu(3)-vect1(3)vectu(2) 982 normu(2)=vect1(3)vectu(1)-vect1(1)vectu(3) 983 normu(3)=vect1(1)vectu(2)-vect1(2)vectu(1) 984 normx(1)=vect1(2)vectx(3)-vect1(3)vectx(2) 985 normx(2)=vect1(3)vectx(1)-vect1(1)vectx(3) 986 normx(3)=vect1(1)vectx(2)-vect1(2)vectx(1) 987 normv(1)=vect1(2)vectv(3)-vect1(3)vectv(2) 988 normv(2)=vect1(3)vectv(1)-vect1(1)vectv(3) 989 normv(3)=vect1(1)vectv(2)-vect1(2)vectv(1) 990 c write(,)normu(1),normu(2),normu(3) 991 c 992 normumag=sqrt(normu(1)2+normu(2)2+normu(3)2) 993 normxmag=sqrt(normx(1)2+normx(2)2+normx(3)2) 994 normvmag=sqrt(normv(1)2+normv(2)2+normv(3)*2) 995 brash 1.12 normu(1)=normu(1)/normumag 996 normu(2)=normu(2)/normumag 997 normu(3)=normu(3)/normumag
998 brash 1.14 normx(1)=normx(1)/normxmag 999 normx(2)=normx(2)/normxmag 1000 normx(3)=normx(3)/normxmag 1001 normv(1)=normv(1)/normvmag 1002 normv(2)=normv(2)/normvmag 1003 normv(3)=normv(3)/normvmag
1004 brash 1.12 c write(,)normumag 1005 c 1006 du(1)=xa-xu 1007 du(2)=ya-yu 1008 du(3)=za-zu 1009 dx(1)=xa-xx 1010 dx(2)=ya-yx 1011 dx(3)=za-zx 1012 dv(1)=xa-xv 1013 dv(2)=ya-yv 1014 dv(3)=za-zv 1015 c 1016 c 1017 c distance formula 1018 distu=du(1)normu(1)+du(2)normu(2)+du(3)normu(3) 1019 distx=dx(1)normx(1)+dx(2)normx(2)+dx(3)normx(3) 1020 distv=dv(1)normv(1)+dv(2)normv(2)+dv(3)*normv(3) 1021 c
1022 brash 1.14 if(distu.gt.1.0.or.distx.gt.1.0.or.distv.gt.1.0) 1023 & idflag=.true. 1024 if(distu.gt.1.28.or.distx.gt.1.28.or.distv.gt.1.28)then 1025 write(,)'Problem Child !!!' 1026 write(,)'distance calculations .....' 1027 write(,)xa,ya,za 1028 c write(,)vect1(1),vect1(2),vect1(3) 1029 write(,)xu,yu,zu 1030 c write(,)vectu(1),vectu(2),vectu(3) 1031 write(,)xx,yx,zx 1032 c write(,)vectx(1),vectx(2),vectx(3) 1033 write(,)xv,yv,zv 1034 c write(,)vectv(1),vectv(2),vectv(3)
1035 brash 1.19 c write(,)'normalization factors' 1036 c write(,)normu(1),normu(2),normu(3) 1037 c write(,)normumag 1038 c write(,)normx(1),normx(2),normx(3) 1039 c write(,)normxmag 1040 c write(,)normv(1),normv(2),normv(3) 1041 c write(,)normvmag 1042 c write(,)'Drift distance: ',distu, distx, distv
1043 brash 1.14 endif
1044 brash 1.12 c 1045 return 1046 end 1047
1048 brash 1.11 1049 subroutine calc_theta_phi(xin1,yin1,zin1,xin2,yin2,zin2, 1050 & xsc1,ysc1,zsc1,xsc2,ysc2,zsc2,theta,phi) 1051 c 1052 implicit none
1053 brash 1.12 include 'fpp_local.h' 1054 include 'geant_local.h'
1055 brash 1.11 c 1056 real8 xin1,yin1,zin1,xin2,yin2,zin2 1057 real8 xsc1,ysc1,zsc1,xsc2,ysc2,zsc2,theta,phi 1058 real*8 ftheta, fphi, fpsi
1059 brash 1.12 real*8 lin,lout,theta_ejb,phi_ejb
1060 brash 1.11 c 1061 real invect(1:3) 1062 real scvect(1:3) 1063 real scvect2(1:3)
1064 brash 1.12 real in(1:3) 1065 real out(1:3) 1066 real scat(1:3)
1067 brash 1.11 c 1068 invect(1)=xin2-xin1 1069 invect(2)=yin2-yin1 1070 invect(3)=zin2-zin1 1071 scvect(1)=xsc2-xsc1 1072 scvect(2)=ysc2-ysc1 1073 scvect(3)=zsc2-zsc1
1074 brash 1.12 c write(,)'INCOMING: ',invect(1),invect(2),invect(3) 1075 c write(,)'SCATTERED: ',scvect(1),scvect(2),scvect(3) 1076 c 1077 c EJB calculation of theta and phi 1078 c 1079 in(1)=invect(1)/invect(3) 1080 in(2)=invect(2)/invect(3) 1081 in(3)=invect(3)/invect(3) 1082 out(1)=scvect(1)/scvect(3) 1083 out(2)=scvect(2)/scvect(3) 1084 out(3)=scvect(3)/scvect(3) 1085 lin=sqrt(in(1)2+in(2)2+in(3)2) 1086 lout=sqrt(out(1)2+out(2)2+out(3)2) 1087 scat(1)=out(1)-in(1) 1088 scat(2)=out(2)-in(2) 1089 scat(3)=out(3) 1090 x_ejb=scat(1) 1091 y_ejb=scat(2) 1092 z_ejb=scat(3) 1093 if(scat(1).ge.0.0.and.scat(2).gt.0.0)then 1094 phi_ejb=atan(scat(1)/scat(2))57.2957795 1095 brash 1.12 else if(scat(1).ge.0.0.and.scat(2).lt.0.0)then 1096 phi_ejb=atan(scat(1)/scat(2))57.2957795+180.00 1097 else if(scat(1).le.0.0.and.scat(2).lt.0.0)then 1098 phi_ejb=atan(scat(1)/scat(2))57.2957795+180.00 1099 else if(scat(1).le.0.0.and.scat(2).gt.0.0)then 1100 phi_ejb=atan(scat(1)/scat(2))57.2957795+360.00 1101 endif 1102 c 1103 theta_ejb=acos((in(1)out(1)+in(2)out(2)+in(3)out(3))/ 1104 & (linlout))57.2957795 1105 c write(,)'EJB Incoming Vector = ',in(1),in(2),in(3) 1106 c write(,)'EJB Outgoing Vector = ',out(1),out(2),out(3) 1107 c write(,)'EJB Scattered Vector = ',scat(1),scat(2),scat(3) 1108 c write(,*)'EJB Thetas = ',theta_ejb,phi_ejb
1109 brash 1.11 c
1110 brash 1.12 c end EJB calculation 1111 c 1112 1113
1114 brash 1.11 ftheta=acos(invect(3)/sqrt(invect(1)2+invect(3)2)) 1115 fphi=acos(invect(3)/sqrt(invect(2)2+invect(3)2)) 1116 fpsi=acos(sqrt(invect(2)2+invect(3)2)/sqrt(invect(1)2 1117 & +invect(2)2+invect(3)**2)) 1118 c
1119 brash 1.12 c write(,)'ftheta, fphi, fpsi', 1120 c & ftheta57.296,fphi57.296,fpsi*57.296
1121 brash 1.11 scvect2(1)=scvect(1)cos(fpsi)-sin(fpsi)(scvect(2)sin(fphi) 1122 & +scvect(3)cos(fphi)) 1123 scvect2(2)=scvect(2)cos(fphi)-scvect(3)sin(fphi) 1124 scvect2(3)=scvect(1)sin(fpsi)+cos(fpsi)(scvect(2)sin(fphi) 1125 & +scvect(3)cos(fphi)) 1126 c
1127 brash 1.12 c write(,)'SCATTERED 2: ',scvect2(1),scvect2(2),scvect2(3)
1128 brash 1.11 theta=atan(sqrt(scvect2(1)2+scvect2(2)2)/scvect2(3))57.2957795 1129 phi=atan(scvect2(1)/scvect2(2))57.2957795 1130 if (scvect2(1).lt.0.0.and.scvect2(2).gt.0.0) 1131 & phi=phi+360.00 1132 if (scvect2(1).lt.0.0.and.scvect2(2).lt.0.0) 1133 & phi=phi+180.00 1134 if (scvect2(1).gt.0.0.and.scvect2(2).lt.0.0) 1135 & phi=phi+180.00
1136 brash 1.12 1137 c write(,)'Theta,phi =',theta,phi
1138 brash 1.15 theta=theta_ejb 1139 phi=phi_ejb
1140 brash 1.11 c 1141 return 1142 end 1143
1144 brash 1.21 subroutine calc_zclose_sclose(zmid,zback, 1145 & x0f,y0f,tphif,tthetaf, 1146 & x0r,y0r,tphir,tthetar, 1147 & theta,phi, 1148 & zclose,sclose,conetest) 1149 1150 implicit none 1151 1152 include 'geant_local.h'
1153 brash 1.19 1154 real8 x0f,y0f,tphif,tthetaf 1155 real8 x0r,y0r,tphir,tthetar 1156 real8 zclose,sclose 1157 real8 term1,term2,term3,term4,term5,term6
1158 brash 1.21 real8 rbig,zmid,zback 1159 real8 xback,xfront,yback,yfront,radius 1160 real8 theta,ttheta,phi,fg,r1x,r1y,r2x,r2y 1161 real8 xpt1,xpt2,xpt3,xpt4 1162 real8 ypt1,ypt2,ypt3,ypt4 1163 integer4 conetest
1164 brash 1.19 1165 rbig=1.0e15
1166 brash 1.21 conetest=1 1167 1168 xback=x0r+(zback-zmid)tphir 1169 yback=y0r+(zback-zmid)tthetar 1170 xfront=x0f+(zback-zmid)tphif 1171 yfront=y0f+(zback-zmid)tthetaf 1172 1173 c write(,)' ' 1174 c write(,)'*******************' 1175 c write(,)'Back Projection = ',xback,yback 1176 c write(,)'Front Projection = ',xfront,yfront 1177 1178 radius=sqrt((xback-xfront)2+(yback-yfront)2) 1179 1180 c write(,*)'Radius = ',radius
1181 brash 1.19 1182 term1=(x0r-x0f)(tphir-tphif) 1183 term2=(y0r-y0f)(tthetar-tthetaf) 1184 term3=(tphir-tphif)2 1185 term4=(tthetar-tthetaf)2 1186 if((term3+term4).ne.0) then 1187 zclose=-(term1+term2)/(term3+term4) 1188 else 1189 zclose=rbig 1190 endif 1191 term5=(x0r-x0f+(tphir-tphif)zclose) 1192 term6=(y0r-y0f+(tthetar-tthetaf)zclose) 1193 sclose=sqrt(term52+term62)
1194 brash 1.21 zclose=zclose+zmid 1195 1196 c write(,)'Scattering angles = ',theta,phi 1197 c write(,)'zclose = ',zclose 1198 c write(,)'zmid = ',zmid 1199 1200 fg=3.14159265/180.0 1201 ttheta=tan(thetafg) 1202 1203 r1x = (zback-zclose)(tphif + (ttheta-tphif)/(1.0+tthetatphif)) 1204 r2x = (zback-zclose)((ttheta+tphif)/(1.0-tthetatphif) - tphif) 1205 r1y = (zback-zclose)(tthetaf + 1206 & (ttheta-tthetaf)/(1.0+tthetatthetaf)) 1207 r2y = (zback-zclose)((ttheta+tthetaf)/(1.0-tthetatthetaf) 1208 & - tthetaf) 1209 1210 c write(,)'Back radii = ',r1x,r2x,r1y,r2y 1211 xpt1=xback-abs(r1x) 1212 ypt1=yback 1213 xpt2=xback+abs(r2x) 1214 ypt2=yback 1215 brash 1.21 xpt3=xback 1216 ypt3=yback-abs(r1y) 1217 xpt4=xback 1218 ypt4=yback+abs(r2y) 1219 1220 if(xpt1.lt.-1.0chamber_xsize)conetest=0 1221 if(xpt2.gt.chamber_xsize)conetest=0 1222 if(ypt3.lt.-1.0chamber_ysize)conetest=0 1223 if(ypt4.gt.chamber_ysize)conetest=0 1224 1225 c write(,)'(',xpt1,',',ypt1,')' 1226 c write(,)'(',xpt2,',',ypt2,')' 1227 c write(,)'(',xpt3,',',ypt3,')' 1228 c write(,)'(',xpt4,',',ypt4,')' 1229 c 1230 c write(,)chamber_xsize,chamber_ysize 1231 c 1232 c write(,*)conetest
1233 brash 1.19 1234 return 1235 end
1236 brash 1.7
1237 brash 1.22 subroutine cone_test(xin1,yin1,zin1,xin2,yin2,zin2, 1238 & theta,phi,zclose,zplane)
1239 brash 1.7
1240 brash 1.22 implicit none 1241 include 'fpp_local.h' 1242 include 'geant_local.h'
1243 brash 1.7
1244 brash 1.22 real8 xin1,yin1,zin1,xin2,yin2,zin2 1245 real8 theta,phi,zclose,zplane,h,r 1246 real8 xproj,yproj,xclose,yclose 1247 real8 xlong,xshort 1248 iconef=1 1249 c 1250 c Calculate the incident particle projection on the 1251 c back of the chamber using two points to define 1252 c the incident track and the z pos of the plane 1253 c 1254 xproj=((zplane-zin1)(xin2-xin1))/(zin2-zin1)+xin1 1255 yproj=((zplane-zin1)(yin2-yin1))/(zin2-zin1)+yin1 1256 write(,)zplane,xin1,yin1,zin1,xin2,yin2,zin2 1257 write(,)'xproj: ',xproj 1258 write(,)'yproj: ',yproj 1259 c 1260 c Calculate the x and y of the scattering point given 1261 c the z of the scattering point 1262 c 1263 xclose=((zclose-zin1)(xin2-xin1))/(zin2-zin1)+xin1 1264 yclose=((zclose-zin1)(yin2-yin1))/(zin2-zin1)+yin1 1265 brash 1.22 write(,)'xclose: ',xclose 1266 write(,)'yclose: ',yclose 1267 write(,)'zclose: ',zclose 1268 write(,)'zplane: ',zplane 1269 write(,)'theta: ',theta 1270 c 1271 c Calculate the height of the cone 1272 h=sqrt((xproj-xclose)2+(yproj-yclose)2+(zplane-zclose)*2) 1273 c 1274 c Calculate the radius of the cone 1275 r=htan(theta0.0174532925) 1276 c 1277 write(,)'height of cone= ',h 1278 write(,)'radius of cone= ',r 1279 c 1280 c Check to see if any part of the distribution is 1281 c outside of the chamber dimensions 1282 c 1283 if((xproj+r)>(41.5-5.08))then 1284 write(,)'over x' 1285 iconef=0 1286 brash 1.22 endif 1287 if((xproj-r)<(-41.5-5.08))then 1288 write(,)'under x' 1289 iconef=0 1290 endif 1291 if((yproj+r)>33.5)then 1292 write(,)'over y ' 1293 iconef=0 1294 endif 1295 if((yproj-r)<-33.5)then 1296 write(,*)'under y' 1297 iconef=0 1298 endif 1299 c 1300 return 1301 end

Analyzer/Replay: Mark Jones, Documents: Stephen Wood