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

Analyzer/Replay: Mark Jones, Documents: Stephen Wood