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
286 brash 1.8
287 brash 1.12 c write(,)'Hit in first chamber ...' 288 c write(,)'Wire Numbers: ',n2ua,n2xa,n2va
289 brash 1.6 290 endif 291 c 292 if ( istop.ne.0 ) then 293 make_hist=0 294 endif 295 else if ( names(nlevel).eq."fch3" ) then 296 if(inwvol.eq.1) then 297 x3a=vect(1) 298 y3a=vect(2) 299 z3a=vect(3) 300 endif 301 if(inwvol.eq.2) then 302 x3b=vect(1) 303 y3b=vect(2) 304 z3b=vect(3) 305
306 brash 1.13 call get_wire_numbers(x3a,y3a,z3a,x3b,y3b,z3b,nhu,nhx,nhv,
307 brash 1.12 & n3ua,n3xa,n3va,n3ub,n3xb,n3vb)
308 brash 1.11
309 brash 1.13 nhu3=nhu 310 nhx3=nhx 311 nhv3=nhv 312 313 if(nhu.eq.1.and.nhv.eq.1.and.nhx.eq.1) then
314 brash 1.12 call get_drift_distance_ejb(x3a,y3a,z3a,x3b,y3b,z3b,
315 brash 1.14 & n3ua,n3xa,n3va,d3ue,d3xe,d3ve,idflag)
316 brash 1.13 call get_drift_distance(x3a,y3a,z3a,x3b,y3b,z3b, 317 & n3ua,n3xa,n3va,d3u,d3x,d3v)
318 brash 1.12 else 319 call get_drift_distance_ejb(x3a,y3a,z3a,x3b,y3b,z3b,
320 brash 1.14 & n3ua,n3xa,n3va,d3ue,d3xe,d3ve,idflag) 321 c write(,)'Drift Distance 3a: ',d3ue,d3xe,d3ve
322 brash 1.12 call get_drift_distance_ejb(x3a,y3a,z3a,x3b,y3b,z3b,
323 brash 1.14 & n3ub,n3xb,n3vb,d3ue,d3xe,d3ve,idflag) 324 c write(,)'Drift Distance 3b: ',d3ue,d3xe,d3ve
325 brash 1.13 call get_drift_distance(x3a,y3a,z3a,x3b,y3b,z3b, 326 & n3ua,n3xa,n3va,d3u,d3x,d3v)
327 brash 1.14 c write(,)'Drift Distance 3c: ',d3u,d3x,d3v
328 brash 1.12 endif 329
330 brash 1.8
331 brash 1.12 c write(,)'Hit in first chamber ...' 332 c write(,)'Wire Numbers: ',n3ua,n3xa,n3va
333 brash 1.6 334 endif 335 c 336 if ( istop.ne.0 ) then 337 make_hist=0 338 endif 339 else if ( names(nlevel).eq."fch4" ) then 340 if(inwvol.eq.1) then 341 x4a=vect(1) 342 y4a=vect(2) 343 z4a=vect(3) 344 endif 345 if(inwvol.eq.2) then 346 x4b=vect(1) 347 y4b=vect(2) 348 z4b=vect(3) 349
350 brash 1.13 call get_wire_numbers(x4a,y4a,z4a,x4b,y4b,z4b,nhu,nhx,nhv,
351 brash 1.12 & n4ua,n4xa,n4va,n4ub,n4xb,n4vb)
352 brash 1.11
353 brash 1.13 nhu4=nhu 354 nhx4=nhx 355 nhv4=nhv 356 357 if(nhu.eq.1.and.nhv.eq.1.and.nhx.eq.1) then
358 brash 1.12 call get_drift_distance_ejb(x4a,y4a,z4a,x4b,y4b,z4b,
359 brash 1.14 & n4ua,n4xa,n4va,d4ue,d4xe,d4ve,idflag)
360 brash 1.13 call get_drift_distance(x4a,y4a,z4a,x4b,y4b,z4b, 361 & n4ua,n4xa,n4va,d4u,d4x,d4v)
362 brash 1.12 else 363 call get_drift_distance_ejb(x4a,y4a,z4a,x4b,y4b,z4b,
364 brash 1.14 & n4ua,n4xa,n4va,d4ue,d4xe,d4ve,idflag) 365 c write(,)'Drift Distance 4a: ',d4ue,d4xe,d4ve
366 brash 1.12 call get_drift_distance_ejb(x4a,y4a,z4a,x4b,y4b,z4b,
367 brash 1.14 & n4ub,n4xb,n4vb,d4ue,d4xe,d4ve,idflag) 368 c write(,)'Drift Distance 4b: ',d4ue,d4xe,d4ve
369 brash 1.13 call get_drift_distance(x4a,y4a,z4a,x4b,y4b,z4b, 370 & n4ua,n4xa,n4va,d4u,d4x,d4v)
371 brash 1.14 c write(,)'Drift Distance 4c: ',d4u,d4x,d4v
372 brash 1.12 endif 373
374 brash 1.8
375 brash 1.12 c write(,)'Hit in first chamber ...' 376 c write(,)'Wire Numbers: ',n3ua,n3xa,n3va
377 brash 1.6 378 endif 379 c
380 brash 1.4 if ( istop.ne.0 ) then 381 make_hist=0 382 endif
383 brash 1.2 else if ( istop.ne.0.and. names(nlevel).eq."sci1" ) then 384 make_hist=0 385 else if ( names(nlevel).eq."anl1" ) then
386 jones 1.1 if(inwvol.eq.1) then
387 brash 1.5 c write(6,)'We have a hit in the fist analyzer at' 388 c write(6,)'x=',vect(1),' y=',vect(2),' z=',vect(3)
389 brash 1.4 xdet=vect(1) 390 ydet=vect(2) 391 zdet=vect(3)
392 jones 1.1 endif 393 if ( istop.ne.0 ) then 394 make_hist=0 395 endif 396 end if 397 c 398 c 399 c store current track parameters (including position ) in jxyz structure. 400 c 401 call gsxyz 402 c 403 c moved histograming stuff to gulast 404 c 405 c write(6,*)'done in gustep' 406 9999 return
407 brash 1.6 end 408
409 brash 1.13 subroutine get_wire_numbers(xa,ya,za,xb,yb,zb,nhu,nhx,nhv,
410 brash 1.12 % nu1,nx1,nv1,nu2,nx2,nv2)
411 brash 1.6
412 brash 1.7 implicit none 413
414 brash 1.6 real*8 xa,ya,za,xb,yb,zb
415 brash 1.12 integer4 nu1,nx1,nv1 416 integer4 nu2,nx2,nv2
417 brash 1.13 integer*4 nhu,nhx,nhv
418 brash 1.6 419 real*8 zc,zu,zx,zv,zt,xp,yp,xu,xx,xv,yu,yx,yv,uw,xw,vw
420 brash 1.7 real*8 anu,anx,anv
421 brash 1.6
422 brash 1.12 nu2=0 423 nx2=0 424 nv2=0 425 c
426 brash 1.6 c We have the (x,y,z) coordinates of the entrance (a) and exit (b) points 427 c of the track. We can use this information to calculate the wire numbers that 428 c were hit in each plane. 429 c 430 zc=(zb-za)/2.0+za 431 zu=zc-1.60 432 zx=zc 433 zv=zc+1.60 434 zt=(zb-za) 435 xp=(xb-xa)/zt 436 yp=(yb-ya)/zt
437 brash 1.12 c
438 brash 1.14 c Project to the FRONT of the "cell" associated with each plane. 439 c 440 xu=xa+xp(zu-za-0.8) 441 yu=ya+yp(zu-za-0.8) 442 xx=xa+xp(zx-za-0.8) 443 yx=ya+yp(zx-za-0.8) 444 xv=xa+xp(zv-za-0.8) 445 yv=ya+yp(zv-za-0.8) 446 c 447 c xu=xa 448 c yu=ya 449 c xx=xa 450 c yx=ya 451 c xv=xa 452 c yv=ya
453 brash 1.6 c 454 uw=(xu+yu)/sqrt(2.0) 455 xw=xx 456 vw=(-xv+yv)/sqrt(2.0)
457 brash 1.7 c
458 brash 1.8 c write(,)'********************'
459 brash 1.10 c write(,)'A: ',xa,ya,za 460 c write(,)'B: ',xb,yb,zb 461 c write(,)'U: ',xu,yu,zu 462 c write(,)'X: ',xx,yx,zx 463 c write(,)'V: ',xv,yv,zv 464 c write(,)'W: ',uw,xw,vw
465 brash 1.14 c write(,)'********************'
466 brash 1.7 c 467 anu=(-uw-3.592+104.0)/2.0 468 anx=(-xw-5.080+84.0)/2.0 469 anv=(vw-3.592+104.0)/2.0 470 c
471 brash 1.10 c write(,)'Wires: ',anu,anx,anv 472 c write(,)'********************'
473 brash 1.12 nu1=anu 474 nv1=anv 475 nx1=anx 476 if((anu-nu1).ge.0.500)nu1=nu1+1 477 if((anx-nx1).ge.0.500)nx1=nx1+1 478 if((anv-nv1).ge.0.500)nv1=nv1+1 479 c write(,)'using front of chamber: ',nu1,nx1,nv1 480 c
481 brash 1.14 c Now project to the BACK of the "cell" associated with each plane.
482 brash 1.12 c
483 brash 1.14 xu=xa+xp(zu-za+0.8) 484 yu=ya+yp(zu-za+0.8) 485 xx=xa+xp(zx-za+0.8) 486 yx=ya+yp(zx-za+0.8) 487 xv=xa+xp(zv-za+0.8) 488 yv=ya+yp(zv-za+0.8) 489 c 490 c xu=xb 491 c yu=yb 492 c xx=xb 493 c yx=yb 494 c xv=xb 495 c yv=yb
496 brash 1.12 c 497 uw=(xu+yu)/sqrt(2.0) 498 xw=xx 499 vw=(-xv+yv)/sqrt(2.0)
500 brash 1.7 c
501 brash 1.12 c write(,)'*******************' 502 c write(,)'A: ',xa,ya,za 503 c write(,)'B: ',xb,yb,zb 504 c write(,)'U: ',xu,yu,zu 505 c write(,)'X: ',xx,yx,zx 506 c write(,*)'V: ',xv,yv,zv
507 brash 1.10 c write(,)'W: ',uw,xw,vw 508 c write(,)'********************'
509 brash 1.12 c 510 anu=(-uw-3.592+104.0)/2.0 511 anx=(-xw-5.080+84.0)/2.0 512 anv=(vw-3.592+104.0)/2.0 513 c 514 c write(,)'Wires: ',anu,anx,anv 515 c write(,)'********************' 516 nu2=anu 517 nv2=anv 518 nx2=anx 519 if((anu-nu2).ge.0.500)nu2=nu2+1 520 if((anx-nx2).ge.0.500)nx2=nx2+1 521 if((anv-nv2).ge.0.500)nv2=nv2+1
522 brash 1.13 nhu=1 523 nhx=1 524 nhv=1 525 if(nu1.ne.nu2)nhu=2 526 if(nx1.ne.nx2)nhx=2 527 if(nv1.ne.nv2)nhv=2
528 brash 1.7
529 brash 1.12 return 530 end
531 brash 1.8
532 brash 1.12 c subroutine get_drift_distance_ejb(xa,ya,za,xb,yb,zb, 533 c & nu,nx,nv,du,dx,dv) 534 c 535 c implicit none 536 c 537 c real8 xa,ya,za,xb,yb,zb,du,dx,dv 538 c integer4 nu,nx,nv 539 c real8 tphi,ttheta 540 c real8 uw,xw,vw 541 c real8 c11,c12,c21,c22,d1,d2,a,zt,xt,yt 542 c real8 xu,yu,zu 543 c real8 xv,yv,zv 544 c real8 xx,yx,zx 545 c 546 c tphi=(xb-xa)/(zb-za) 547 c ttheta=(yb-ya)/(zb-za) 548 c 549 c uw=-2.0nu-3.592+104.0 550 c xw=-2.0nx-5.080+84.0 551 c vw=2.0nv+3.592-104.0 552 c zu=(zb+za)/2.0-1.60 553 brash 1.12 c zv=(zb+za)/2.0+1.60 554 c zx=(zb+za)/2.0 555 c 556 cc write(,)'W: ',uw,xw,vw 557 cc write(,)'******************' 558 c 559 c c11=tphi-ttheta 560 c c12=sqrt(2.0) 561 c d1=-xa+ya 562 c c21=tphitphi+tthetattheta+1.0 563 c c22=-ttheta/sqrt(2.0)+tphi/sqrt(2.0) 564 c d2=uw(ttheta+tphi)/sqrt(2.0)-xatphi-yattheta+zu-za 565 c 566 c a=(d1c21-d2c11)/(c21c12-c22c11) 567 c zt=(d1-c12a)/c11 568 cc write(,)'U Plane zt,a: ',zt,a 569 c 570 c xt=xa+zttphi 571 c yt=ya+ztttheta 572 c xu=(uw-a)/sqrt(2.0) 573 c yu=(uw+a)/sqrt(2.0) 574 brash 1.12 c 575 c zt=zt+za 576 c 577 c du=sqrt((xt-xu)2+(yt-yu)2+(zt-zu)2) 578 c 579 c c11=tphi+ttheta 580 c c12=-sqrt(2.0) 581 c d1=-xa-ya 582 c c21=tphitphi+tthetattheta+1.0 583 c c22=-ttheta/sqrt(2.0)-tphi/sqrt(2.0) 584 c d2=vw(ttheta+tphi)/sqrt(2.0)-xatphi-yattheta+zv-za 585 c 586 c a=(d1c21-d2c11)/(c21c12-c22c11) 587 c zt=(d1-c12a)/c11 588 cc write(,)'V Plane zt,a: ',zt,a 589 c 590 c xt=xa+zttphi 591 c yt=ya+ztttheta 592 c xv=-(vw-a)/sqrt(2.0) 593 c yv=(vw+a)/sqrt(2.0) 594 c 595 brash 1.12 c zt=zt+za 596 c 597 c dv=sqrt((xt-xv)2+(yt-yv)2+(zt-zv)2) 598 c 599 c c11=tphi 600 c c12=-1.0 601 c d1=-ya 602 c c21=tphitphi+tthetattheta+1.0 603 c c22=-ttheta 604 c d2=xwtphi-xatphi-yattheta+zx-za 605 c 606 c a=(d1c21-d2c11)/(c21c12-c22c11) 607 c zt=(d1-c12a)/c11 608 cc write(,)'X Plane zt,a: ',zt,a 609 c 610 c xt=xa+zttphi 611 c yt=ya+ztttheta 612 c xx=xw 613 c yx=a 614 c 615 c zt=zt+za 616 brash 1.12 c 617 c dx=sqrt((xt-xx)2+(yt-yx)2+(zt-zx)*2) 618 c 619 cc if((du.gt.1.00)) then
620 brash 1.11 c write(,)'Calculating Drift Distance ...' 621 c write(,)'A: ',xa,ya,za 622 c write(,)'B: ',xb,yb,zb
623 brash 1.12 cc write(,)'U: ',xu,yu,zu 624 cc write(,)'V: ',xv,yv,zv 625 cc write(,)'X: ',xx,yx,zx
626 brash 1.11 c write(,)'U Drift distance = ',du 627 c write(,)'V Drift distance = ',dv 628 c write(,)'X Drift distance = ',dx
629 brash 1.12 cc endif 630 c 631 c return 632 c end
633 brash 1.11 634 635 subroutine get_drift_distance(xa,ya,za,xb,yb,zb, 636 & nu,nx,nv,distu,distx,distv) 637 c 638 c This subroutine uses the entrance (a) and exit (b) 639 c points of a chamber to define the line of the track. 640 c It uses the wire number and wire direction to define 641 c the wire line. It then uses these to build parallel 642 c planes by computing a normal vector to both of the lines. 643 c Finally, it calculates the distance between these two 644 c planes which is the distance of shortest approach. 645 c 646 implicit none 647 c 648 real8 xa,ya,za,xb,yb,zb 649 real8 distu,distx,distv,xp,yp,zc 650 integer4 nu,nx,nv 651 real8 uw,xw,vw 652 real8 zt,xt,yt 653 real8 xu,yu,zu 654 brash 1.11 real8 xv,yv,zv 655 real8 xx,yx,zx 656 c 657 c the direction of each of the lines 658 c vect1=track vectu=u wire
659 brash 1.14 real8 vect1(1:3), vectu(1:3) 660 real8 vectx(1:3), vectv(1:3)
661 brash 1.11 c 662 c the normal vector to both lines
663 brash 1.14 real8 normu(1:3) 664 real8 normx(1:3) 665 real*8 normv(1:3)
666 brash 1.11 c 667 c the coefficients of the plane
668 brash 1.14 real8 pvectu(1:4) 669 real8 pvectx(1:4) 670 real*8 pvectv(1:4)
671 brash 1.11 c 672 vect1(1)=xb-xa 673 vect1(2)=yb-ya 674 vect1(3)=zb-za 675 c 676 zu=(zb+za)/2.0-1.60 677 zv=(zb+za)/2.0+1.60 678 zx=(zb+za)/2.0 679 c 680 c use line number to calculate distance relative to 681 c wire plane, then convert to x and y 682 c write(,)'nu nx nv',nu,nx,nv 683 uw=-2.0nu-3.592+104.0 684 xw=-2.0nx-5.080+84.0 685 vw=2.0*nv+3.592-104.0 686 xu=uw/sqrt(2.0) 687 yu=uw/sqrt(2.0) 688 xx=xw 689 yx=0 690 xv=-vw/sqrt(2.0) 691 yv=vw/sqrt(2.0)
692 brash 1.12 c write(,)'uw xu yu zu',uw,xu,yu,zu
693 brash 1.11 c write(,)'xw xx yx zx',xw,xx,yx,zx 694 c write(,)'vw xv yv zv',vw,xv,yv,zv 695 c 696 c define direction vector for wires, will be the same 697 c for each wire in a given plane, and is known 698 c for each plane 699 vectu(1)=1.0/sqrt(2.0) 700 vectu(2)=1.0/sqrt(2.0) 701 vectu(3)=0.0 702 vectx(1)=0.0 703 vectx(2)=1.0 704 vectx(3)=0.0 705 vectv(1)=-1.0/sqrt(2.0) 706 vectv(2)=1.0/sqrt(2.0) 707 vectv(3)=0.0 708 c
709 brash 1.12 c write(,)'distance calculations .....' 710 c write(,)xa,ya,za 711 c write(,)vect1(1),vect1(2),vect1(3) 712 c write(,)xu,yu,zu 713 c write(,)vectu(1),vectu(2),vectu(3) 714 c write(,)xx,yx,zx 715 c write(,)vectx(1),vectx(2),vectx(3) 716 c write(,)xv,yv,zv 717 c write(,)vectv(1),vectv(2),vectv(3) 718 c write(,)'distance calculations .....' 719 c
720 brash 1.11 c cross product 721 normu(1)=vect1(2)vectu(3)-vect1(3)vectu(2) 722 normu(2)=vect1(1)vectu(3)-vect1(3)vectu(1) 723 normu(3)=vect1(1)vectu(2)-vect1(2)vectu(1) 724 normx(1)=vect1(2)vectx(3)-vect1(3)vectx(2) 725 normx(2)=vect1(1)vectx(3)-vect1(3)vectx(1) 726 normx(3)=vect1(1)vectx(2)-vect1(2)vectx(1) 727 normv(1)=vect1(2)vectv(3)-vect1(3)vectv(2) 728 normv(2)=vect1(1)vectv(3)-vect1(3)vectv(1) 729 normv(3)=vect1(1)vectv(2)-vect1(2)vectv(1) 730 c 731 pvectu(1)=normu(1) 732 pvectu(2)=normu(2) 733 pvectu(3)=normu(3) 734 pvectu(4)=normu(1)(-xu)+normu(2)(-yu)+normu(3)(-zu) 735 pvectx(1)=normx(1) 736 pvectx(2)=normx(2) 737 pvectx(3)=normx(3) 738 pvectx(4)=normx(1)(-xx)+normx(2)(-yx)+normx(3)(-zx) 739 pvectv(1)=normv(1) 740 pvectv(2)=normv(2) 741 brash 1.11 pvectv(3)=normv(3) 742 pvectv(4)=normv(1)(-xv)+normv(2)(-yv)+normv(3)(-zv) 743 c 744 c distance formula 745 distu=(pvectu(1)xa+pvectu(2)ya+pvectu(3)za+pvectu(4)) 746 & /sqrt(normu(1)2+normu(2)2+normu(3)*2) 747 distx=(pvectx(1)xa+pvectx(2)ya+pvectx(3)za+pvectx(4)) 748 & /sqrt(normx(1)2+normx(2)2+normx(3)*2) 749 distv=(pvectv(1)xa+pvectv(2)ya+pvectv(3)za+pvectv(4)) 750 & /sqrt(normv(1)2+normv(2)2+normv(3)*2) 751 c write(,*)'Drift distance: ',distu, distx, distv
752 brash 1.14 c 753 c write(,)'Brads routine ....' 754 c write(,)xa,ya,za 755 c write(,)xb,yb,zb 756 c write(,)nu,nx,nv 757 c write(,)uw,xw,vw 758 c write(,)'Drift distance: ',distu, distx, distv 759
760 brash 1.11 return 761 end 762
763 brash 1.12 subroutine get_drift_distance_ejb(xa,ya,za,xb,yb,zb,
764 brash 1.14 & nu,nx,nv,distu,distx,distv,idflag)
765 brash 1.12 c 766 c Author: Ed Brash - December 15th, 2005 767 c Yet another attempt at a full drift distance calculation 768 c 769 implicit none 770 c 771 real8 xa,ya,za,xb,yb,zb 772 real8 distu,distx,distv,xp,yp,zc 773 integer4 nu,nx,nv 774 real8 uw,xw,vw 775 real8 zt,xt,yt 776 real8 xu,yu,zu 777 real8 xv,yv,zv 778 real8 xx,yx,zx
779 brash 1.14 logical idflag
780 brash 1.12 c 781 c the direction of each of the lines 782 c vect1=track vectu=u wire 783 real8 vect1(1:3), vectu(1:3) 784 real8 vectx(1:3), vectv(1:3) 785 786 c the difference vector between the defining points 787 real8 du(1:3),dx(1:3),dv(1:3) 788 c 789 c the normal vector to both lines 790 real8 normu(1:3) 791 real8 normx(1:3) 792 real8 normv(1:3) 793 real8 normumag,normxmag,normvmag 794 c 795 c the coefficients of the distance vector 796 real8 dvectu(1:4) 797 real8 dvectx(1:4) 798 real8 dvectv(1:4) 799 c
800 brash 1.14 idflag=.false.
801 brash 1.12 vect1(1)=xb-xa 802 vect1(2)=yb-ya 803 vect1(3)=zb-za 804 c 805 zu=(zb+za)/2.0-1.60 806 zv=(zb+za)/2.0+1.60 807 zx=(zb+za)/2.0 808 c 809 c use line number to calculate distance relative to 810 c wire plane, then convert to x and y 811 c write(,)'nu nx nv',nu,nx,nv 812 uw=-2.0nu-3.592+104.0 813 xw=-2.0nx-5.080+84.0 814 vw=2.0nv+3.592-104.0 815 xu=uw/sqrt(2.0) 816 yu=uw/sqrt(2.0) 817 xx=xw 818 yx=0 819 xv=-vw/sqrt(2.0) 820 yv=vw/sqrt(2.0) 821 c write(,)'uw xu yu zu',uw,xu,yu,zu 822 brash 1.12 c write(,)'xw xx yx zx',xw,xx,yx,zx 823 c write(,)'vw xv yv zv',vw,xv,yv,zv 824 c 825 c define direction vector for wires, will be the same 826 c for each wire in a given plane, and is known 827 c for each plane 828 vectu(1)=1.0/sqrt(2.0) 829 vectu(2)=-1.0/sqrt(2.0) 830 vectu(3)=0.0 831 vectx(1)=0.0 832 vectx(2)=1.0 833 vectx(3)=0.0 834 vectv(1)=1.0/sqrt(2.0) 835 vectv(2)=1.0/sqrt(2.0) 836 vectv(3)=0.0 837 c 838 c write(,)'distance calculations .....' 839 c write(,)xa,ya,za 840 c write(,)vect1(1),vect1(2),vect1(3) 841 c write(,)xu,yu,zu 842 c write(,)vectu(1),vectu(2),vectu(3) 843 brash 1.12 c write(,)xx,yx,zx 844 c write(,)vectx(1),vectx(2),vectx(3) 845 c write(,)xv,yv,zv 846 c write(,)vectv(1),vectv(2),vectv(3) 847 c write(,)'distance calculations .....' 848 c 849 c cross product 850 normu(1)=vect1(2)vectu(3)-vect1(3)vectu(2) 851 normu(2)=vect1(3)vectu(1)-vect1(1)vectu(3) 852 normu(3)=vect1(1)vectu(2)-vect1(2)vectu(1) 853 normx(1)=vect1(2)vectx(3)-vect1(3)vectx(2) 854 normx(2)=vect1(3)vectx(1)-vect1(1)vectx(3) 855 normx(3)=vect1(1)vectx(2)-vect1(2)vectx(1) 856 normv(1)=vect1(2)vectv(3)-vect1(3)vectv(2) 857 normv(2)=vect1(3)vectv(1)-vect1(1)vectv(3) 858 normv(3)=vect1(1)vectv(2)-vect1(2)vectv(1) 859 c write(,)normu(1),normu(2),normu(3) 860 c 861 normumag=sqrt(normu(1)2+normu(2)2+normu(3)2) 862 normxmag=sqrt(normx(1)2+normx(2)2+normx(3)2) 863 normvmag=sqrt(normv(1)2+normv(2)2+normv(3)*2) 864 brash 1.12 normu(1)=normu(1)/normumag 865 normu(2)=normu(2)/normumag 866 normu(3)=normu(3)/normumag
867 brash 1.14 normx(1)=normx(1)/normxmag 868 normx(2)=normx(2)/normxmag 869 normx(3)=normx(3)/normxmag 870 normv(1)=normv(1)/normvmag 871 normv(2)=normv(2)/normvmag 872 normv(3)=normv(3)/normvmag
873 brash 1.12 c write(,)normumag 874 c 875 du(1)=xa-xu 876 du(2)=ya-yu 877 du(3)=za-zu 878 dx(1)=xa-xx 879 dx(2)=ya-yx 880 dx(3)=za-zx 881 dv(1)=xa-xv 882 dv(2)=ya-yv 883 dv(3)=za-zv 884 c 885 c 886 c distance formula 887 distu=du(1)normu(1)+du(2)normu(2)+du(3)normu(3) 888 distx=dx(1)normx(1)+dx(2)normx(2)+dx(3)normx(3) 889 distv=dv(1)normv(1)+dv(2)normv(2)+dv(3)*normv(3) 890 c
891 brash 1.14 if(distu.gt.1.0.or.distx.gt.1.0.or.distv.gt.1.0) 892 & idflag=.true. 893 if(distu.gt.1.28.or.distx.gt.1.28.or.distv.gt.1.28)then 894 write(,)'Problem Child !!!' 895 write(,)'distance calculations .....' 896 write(,)xa,ya,za 897 c write(,)vect1(1),vect1(2),vect1(3) 898 write(,)xu,yu,zu 899 c write(,)vectu(1),vectu(2),vectu(3) 900 write(,)xx,yx,zx 901 c write(,)vectx(1),vectx(2),vectx(3) 902 write(,)xv,yv,zv 903 c write(,)vectv(1),vectv(2),vectv(3) 904 write(,)'normalization factors' 905 write(,)normu(1),normu(2),normu(3) 906 write(,)normumag 907 write(,)normx(1),normx(2),normx(3) 908 write(,)normxmag 909 write(,)normv(1),normv(2),normv(3) 910 write(,)normvmag 911 write(,)'Drift distance: ',distu, distx, distv 912 brash 1.14 endif
913 brash 1.12 c 914 return 915 end 916
917 brash 1.11 918 subroutine calc_theta_phi(xin1,yin1,zin1,xin2,yin2,zin2, 919 & xsc1,ysc1,zsc1,xsc2,ysc2,zsc2,theta,phi) 920 c 921 implicit none
922 brash 1.12 include 'fpp_local.h' 923 include 'geant_local.h'
924 brash 1.11 c 925 real8 xin1,yin1,zin1,xin2,yin2,zin2 926 real8 xsc1,ysc1,zsc1,xsc2,ysc2,zsc2,theta,phi 927 real*8 ftheta, fphi, fpsi
928 brash 1.12 real*8 lin,lout,theta_ejb,phi_ejb
929 brash 1.11 c 930 real invect(1:3) 931 real scvect(1:3) 932 real scvect2(1:3)
933 brash 1.12 real in(1:3) 934 real out(1:3) 935 real scat(1:3)
936 brash 1.11 c 937 invect(1)=xin2-xin1 938 invect(2)=yin2-yin1 939 invect(3)=zin2-zin1 940 scvect(1)=xsc2-xsc1 941 scvect(2)=ysc2-ysc1 942 scvect(3)=zsc2-zsc1
943 brash 1.12 c write(,)'INCOMING: ',invect(1),invect(2),invect(3) 944 c write(,)'SCATTERED: ',scvect(1),scvect(2),scvect(3) 945 c 946 c EJB calculation of theta and phi 947 c 948 in(1)=invect(1)/invect(3) 949 in(2)=invect(2)/invect(3) 950 in(3)=invect(3)/invect(3) 951 out(1)=scvect(1)/scvect(3) 952 out(2)=scvect(2)/scvect(3) 953 out(3)=scvect(3)/scvect(3) 954 lin=sqrt(in(1)2+in(2)2+in(3)2) 955 lout=sqrt(out(1)2+out(2)2+out(3)2) 956 scat(1)=out(1)-in(1) 957 scat(2)=out(2)-in(2) 958 scat(3)=out(3) 959 x_ejb=scat(1) 960 y_ejb=scat(2) 961 z_ejb=scat(3) 962 if(scat(1).ge.0.0.and.scat(2).gt.0.0)then 963 phi_ejb=atan(scat(1)/scat(2))57.2957795 964 brash 1.12 else if(scat(1).ge.0.0.and.scat(2).lt.0.0)then 965 phi_ejb=atan(scat(1)/scat(2))57.2957795+180.00 966 else if(scat(1).le.0.0.and.scat(2).lt.0.0)then 967 phi_ejb=atan(scat(1)/scat(2))57.2957795+180.00 968 else if(scat(1).le.0.0.and.scat(2).gt.0.0)then 969 phi_ejb=atan(scat(1)/scat(2))57.2957795+360.00 970 endif 971 c 972 theta_ejb=acos((in(1)out(1)+in(2)out(2)+in(3)out(3))/ 973 & (linlout))57.2957795 974 c write(,)'EJB Incoming Vector = ',in(1),in(2),in(3) 975 c write(,)'EJB Outgoing Vector = ',out(1),out(2),out(3) 976 c write(,)'EJB Scattered Vector = ',scat(1),scat(2),scat(3) 977 c write(,*)'EJB Thetas = ',theta_ejb,phi_ejb
978 brash 1.11 c
979 brash 1.12 c end EJB calculation 980 c 981 982
983 brash 1.11 ftheta=acos(invect(3)/sqrt(invect(1)2+invect(3)2)) 984 fphi=acos(invect(3)/sqrt(invect(2)2+invect(3)2)) 985 fpsi=acos(sqrt(invect(2)2+invect(3)2)/sqrt(invect(1)2 986 & +invect(2)2+invect(3)**2)) 987 c
988 brash 1.12 c write(,)'ftheta, fphi, fpsi', 989 c & ftheta57.296,fphi57.296,fpsi*57.296
990 brash 1.11 scvect2(1)=scvect(1)cos(fpsi)-sin(fpsi)(scvect(2)sin(fphi) 991 & +scvect(3)cos(fphi)) 992 scvect2(2)=scvect(2)cos(fphi)-scvect(3)sin(fphi) 993 scvect2(3)=scvect(1)sin(fpsi)+cos(fpsi)(scvect(2)sin(fphi) 994 & +scvect(3)cos(fphi)) 995 c
996 brash 1.12 c write(,)'SCATTERED 2: ',scvect2(1),scvect2(2),scvect2(3)
997 brash 1.11 theta=atan(sqrt(scvect2(1)2+scvect2(2)2)/scvect2(3))57.2957795 998 phi=atan(scvect2(1)/scvect2(2))57.2957795 999 if (scvect2(1).lt.0.0.and.scvect2(2).gt.0.0) 1000 & phi=phi+360.00 1001 if (scvect2(1).lt.0.0.and.scvect2(2).lt.0.0) 1002 & phi=phi+180.00 1003 if (scvect2(1).gt.0.0.and.scvect2(2).lt.0.0) 1004 & phi=phi+180.00
1005 brash 1.12 1006 c write(,)'Theta,phi =',theta,phi
1007 brash 1.11 1008 c 1009 return 1010 end 1011
1012 brash 1.7 1013 1014 1015 1016 1017 1018 1019
1020 jones 1.1

Analyzer/Replay: Mark Jones, Documents: Stephen Wood