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                 common/gctrak/vect(7),getot,gekin,vout(7),nmec,lmec(30),namec(30)
 11                +     ,nstep ,maxnst,destep,destel,safety,sleng 
 12                +     ,step  ,snext ,sfield
 13                +     ,tofg  ,gekrat,upwght,ignext,inwvol,istop ,idecad,iekbin
 14                +     , ilosl, imull,ingoto,nldown,nlevin,nlvsav,istory
 15           c     
 16                 integer nmec,lmec,namec,nstep ,maxnst,ignext,inwvol,istop
 17                +     ,idecad,iekbin,ilosl, imull,ingoto,nldown,nlevin
 18                +     ,nlvsav,istory
 19                 real  vect,getot,gekin,vout,destep,destel,safety,sleng ,step
 20                +     ,snext,sfield,tofg  ,gekrat,upwght
 21           c     end gctrak
 22 jones 1.1 *     keep,gcvolu.
 23           *--   author :
 24                 common/gcvolu/nlevel,names(15),number(15),
 25                +     lvolum(15),lindex(15),infrom,nlevmx,nldev(15),linmx(15),
 26                +     gtran(3,15),grmat(10,15),gonly(15),glx(3)
 27           c     
 28                 integer nlevel,number,lvolum,lindex,infrom,nlevmx,
 29                +     nldev,linmx
 30                 character*4 names
 31                 real gtran,grmat,gonly,glx
 32           c     end gcvolu
 33           c     
 34           *     keep,gcbank.
 35           *--   author :
 36                 integer iq,lq,nzebra,ixstor,ixdiv,ixcons,lmain,lr1
 37                 integer kwbank,kwwork,iws
 38                 real gversn,zversn,fendq,ws,q
 39           c     
 40                 parameter (kwbank=69000,kwwork=5200)
 41                 common/gcbank/nzebra,gversn,zversn,ixstor,ixdiv,ixcons,fendq(16)
 42                +     ,lmain,lr1,ws(kwbank)
 43 jones 1.1       dimension iq(2),q(2),lq(8000),iws(2)
 44                 equivalence (q(1),iq(1),lq(9)),(lq(1),lmain),(iws(1),ws(1))
 45                 common/gclink/jdigi ,jdraw ,jhead ,jhits ,jkine ,jmate ,jpart
 46                +     ,jrotm ,jrung ,jset  ,jstak ,jgstat,jtmed ,jtrack,jvertx
 47                +     ,jvolum,jxyz  ,jgpar ,jgpar2,jsklt
 48           c     
 49                 integer       jdigi ,jdraw ,jhead ,jhits ,jkine ,jmate ,jpart
 50                +     ,jrotm ,jrung ,jset  ,jstak ,jgstat,jtmed ,jtrack,jvertx
 51                +     ,jvolum,jxyz  ,jgpar,jgpar2 ,jsklt
 52           c     
 53           *     keep,gcking.
 54           *--   author :
 55                 common/gcking/kcase,ngkine,gkin(5,100),tofd(100),iflgk(100)
 56                 integer       kcase,ngkine ,iflgk
 57                 real          gkin,tofd
 58           c     end gcking
 59           c     
 60           *     keep,gckine.
 61           *--   author :
 62           *--   author :
 63                 integer       ikine,itra,istak,ivert,ipart,itrtyp,napart,ipaold
 64 jones 1.1       real          pkine,amass,charge,tlife,vert,pvert
 65                 common/gckine/ikine,pkine(10),itra,istak,ivert,ipart,itrtyp
 66                +     ,napart(5),amass,charge,tlife,vert(3),pvert(4),ipaold
 67           c     end gckine
 68           c     
 69                 integer ihset,ihdet,iset,idet,idtype,nvname,numbv
 70                 common/gcsets/ihset,ihdet,iset,idet,idtype,nvname,numbv(20)
 71                 real x1,y1,z1,lpar,v1,v2,v3,newdist,x1new,y1new,z1new
 72                 real xstr,ystr,zstr,xstrnew,ystrnew,zstrnew
 73                 real rotmat2,rotmat3,rotmat4,rotmat1

 74 brash 1.6 c

 75 jones 1.1       common/geomstep/rotmat1(3,3),rotmat2(3,3),rotmat3(3,3),
 76                &     rotmat4(3,3)
 77           c     
 78                 include 'fpp_local.h'
 79                 include 'geant_local.h'
 80           c      include 'parameter.h'
 81           c      include 'espace_type.h'
 82           c      include 'detector.h'
 83           c      include 'transport.h'
 84           c      include 'option.h'
 85           c     
 86           c     
 87                 integer i,make_hist,ioff,ihit,ieffcheck
 88                 real pt_pi,ppar_pi,arg1,arg2,rapid_pi,pchmb,hits(6)
 89                 real a,b,c,beta,z,y,straw,ypath,rndm(3),ycompare

 90 brash 1.2 c      write(6,*)'entering gustep'
 91           c      write(6,*)'inwvol =',inwvol
 92           c      write(6,*)'position =',vect(1),vect(2),vect(3)
 93           c      write(6,*)'names =',names(nlevel)

 94 jones 1.1 c     
 95           c
 96           c
 97                 if ( ngkine.gt.0. ) then
 98                    do i=1,nmec
 99                       if ( lmec(i).eq.12 ) then
100           c               write ( 6,* ) ' gustep: hadronic interaction'
101           c               write ( 6,* ) ' nevent=',nevent
102                       end if
103                    end do
104                    mylast = min(100,ngkine)
105                    do i=1,mylast
106                       iflgk(i) = 1
107                       if ( gkin(5,i).eq.4 ) iflgk(i) = 0
108                       if ( gkin(4,i).gt.0.001 ) iflgk(i)=0
109                    end do
110           c         n_2nd = n_2nd + ngkine
111                 endif
112           c     
113           c  the following call makes sure all of the secondary particles
114           c  get tracked too (provided the flag iflgk(i) for that particle
115 jones 1.1 c  was set in the loop above -- this point is not correctly or
116           c  clearly documented in the version of the geant manual i have).
117           c
118           c      if(sectrack) then
119           c         call gsking ( 0 )
120           c      endif
121           c     

122 brash 1.6 c      write(*,*)'In GUSTEP: ... names(nlevel) = ',names(nlevel)
123           

124 jones 1.1       make_hist = 0
125                 if ( inwvol.eq.1 .and. names(nlevel).eq."hall" ) then
126                    make_hist=1
127           c     
128           c     if we get here, the tracking is done for this track
129           c     (istop=1) but make a bunch of histograms before exitting
130           c     note that ipart=8 means a pi+ and 9 is a pi-
131           c     
132           c     istop = 1
133                 else if ( istop.ne.0.and. names(nlevel).eq."aira" ) then
134                    make_hist=0
135                 else if ( istop.ne.0.and. names(nlevel).eq."airb" ) then
136                    make_hist=0
137                 else if ( istop.ne.0.and. names(nlevel).eq."airc" ) then
138                    make_hist=0
139                 else if ( istop.ne.0.and. names(nlevel).eq."aird" ) then
140                    make_hist=0

141 brash 1.2       else if ( istop.ne.0.and. names(nlevel).eq."aire" ) then
142                    make_hist=0
143                 else if ( istop.ne.0.and. names(nlevel).eq."airf" ) then
144                    make_hist=0

145 brash 1.6       else if ( names(nlevel).eq."airg" ) then
146           c         write(*,*)'In airg ... inwvol = ',inwvol
147           c         write(*,*)'Z-value = ',vect(3)
148                    if (istop.ne.0) then
149                       make_hist=0
150                    endif
151                 else if ( names(nlevel).eq."HALL" ) then
152           c         write(*,*)'In HALL ... inwvol = ',inwvol
153           c         write(*,*)'Z-value = ',vect(3)
154                    if (istop.ne.0) then
155                       make_hist=0
156                    endif

157 brash 1.2       else if ( istop.ne.0.and. names(nlevel).eq."airh" ) then
158                    make_hist=0
159                 else if ( istop.ne.0.and. names(nlevel).eq."hch1" ) then
160                    make_hist=0
161                 else if ( istop.ne.0.and. names(nlevel).eq."hch2" ) then
162                    make_hist=0

163 brash 1.4       else if ( names(nlevel).eq."fch1" ) then

164 brash 1.6 c         write(*,*)'In fch1 ... inwvol = ',inwvol

165 brash 1.4          if(inwvol.eq.1) then

166 brash 1.5 c           write(6,*)'Coordinates at fch1'
167           c           write(6,*)'x=',vect(1),' y=',vect(2),' z=',vect(3)

168 brash 1.6            x1a=vect(1)
169                      y1a=vect(2)
170                      z1a=vect(3)

171 brash 1.4          endif

172 brash 1.6          if(inwvol.eq.2) then
173                       x1b=vect(1)
174                       y1b=vect(2)
175                       z1b=vect(3)
176           
177                       call get_wire_numbers(x1a,y1a,z1a,x1b,y1b,z1b,n1u,n1x,n1v)
178           

179 brash 1.7             call get_drift_distance(x1a,y1a,z1a,x1b,y1b,z1b,
180                &           n1u,n1x,n1v,d1u,d1x,d1v)
181           

182 brash 1.6             write(*,*)'Hit in first chamber ...'
183                       write(*,*)'Wire Numbers: ',n1u,n1x,n1v
184           
185                    endif
186           
187                    if ( istop.ne.0 ) then
188                       make_hist=0
189                    endif
190                 else if ( names(nlevel).eq."fch2" ) then
191           c         write(*,*)'In fch2 ... inwvol = ',inwvol
192           c         write(*,*)'Z-value = ',vect(3)
193                    if(inwvol.eq.1) then
194                      x2a=vect(1)
195                      y2a=vect(2)
196                      z2a=vect(3)
197                    endif
198                    if(inwvol.eq.2) then
199                       x2b=vect(1)
200                       y2b=vect(2)
201                       z2b=vect(3)
202           
203 brash 1.6             call get_wire_numbers(x2a,y2a,z2a,x2b,y2b,z2b,n2u,n2x,n2v)
204           
205                       write(*,*)'Hit in second chamber ...'
206                       write(*,*)'Wire Numbers: ',n2u,n2x,n2v
207           
208                    endif
209           c
210                    if ( istop.ne.0 ) then
211                       make_hist=0
212                    endif
213                 else if ( names(nlevel).eq."fch3" ) then
214                    if(inwvol.eq.1) then
215                      x3a=vect(1)
216                      y3a=vect(2)
217                      z3a=vect(3)
218                    endif
219                    if(inwvol.eq.2) then
220                       x3b=vect(1)
221                       y3b=vect(2)
222                       z3b=vect(3)
223           
224 brash 1.6             call get_wire_numbers(x3a,y3a,z3a,x3b,y3b,z3b,n3u,n3x,n3v)
225           
226                       write(*,*)'Hit in third chamber ...'
227                       write(*,*)'Wire Numbers: ',n3u,n3x,n3v
228           
229                    endif
230           c
231                    if ( istop.ne.0 ) then
232                       make_hist=0
233                    endif     
234                 else if ( names(nlevel).eq."fch4" ) then
235                    if(inwvol.eq.1) then
236                      x4a=vect(1)
237                      y4a=vect(2)
238                      z4a=vect(3)
239                    endif
240                    if(inwvol.eq.2) then
241                       x4b=vect(1)
242                       y4b=vect(2)
243                       z4b=vect(3)
244           
245 brash 1.6             call get_wire_numbers(x4a,y4a,z4a,x4b,y4b,z4b,n4u,n4x,n4v)
246           
247                       write(*,*)'Hit in fourth chamber ...'
248                       write(*,*)'Wire Numbers: ',n4u,n4x,n4v
249           
250                    endif
251           c

252 brash 1.4          if ( istop.ne.0 ) then
253                       make_hist=0
254                    endif

255 brash 1.2       else if ( istop.ne.0.and. names(nlevel).eq."sci1" ) then
256                    make_hist=0
257                 else if ( names(nlevel).eq."anl1" ) then

258 jones 1.1          if(inwvol.eq.1) then

259 brash 1.5 c           write(6,*)'We have a hit in the fist analyzer at'
260           c           write(6,*)'x=',vect(1),' y=',vect(2),' z=',vect(3)

261 brash 1.4            xdet=vect(1)
262                      ydet=vect(2)
263                      zdet=vect(3)

264 jones 1.1          endif
265                    if ( istop.ne.0 ) then
266                       make_hist=0
267                    endif
268                 end if
269           c     
270           c     
271           c     store current track parameters (including position ) in jxyz structure.
272           c     
273                 call gsxyz
274           c     
275           c     moved histograming stuff to gulast
276           c     
277           c     write(6,*)'done in gustep'
278            9999 return

279 brash 1.6       end
280           
281                 subroutine get_wire_numbers(xa,ya,za,xb,yb,zb,nu,nx,nv)
282           

283 brash 1.7       implicit none
284           

285 brash 1.6       real*8 xa,ya,za,xb,yb,zb

286 brash 1.7       integer*4 nu,nx,nv

287 brash 1.6 
288                 real*8 zc,zu,zx,zv,zt,xp,yp,xu,xx,xv,yu,yx,yv,uw,xw,vw

289 brash 1.7       real*8 anu,anx,anv

290 brash 1.6 
291           c We have the (x,y,z) coordinates of the entrance (a) and exit (b) points
292           c of the track.  We can use this information to calculate the wire numbers that
293           c were hit in each plane.
294           c
295                 zc=(zb-za)/2.0+za
296                 zu=zc-1.60
297                 zx=zc
298                 zv=zc+1.60
299                 zt=(zb-za)
300                 xp=(xb-xa)/zt
301                 yp=(yb-ya)/zt

302 brash 1.7       xu=xa+xp*(zu-za)
303                 yu=ya+yp*(zu-za)
304                 xx=xa+xp*(zx-za)
305                 yx=ya+yp*(zx-za)
306                 xv=xa+xp*(zv-za)
307                 yv=ya+yp*(zv-za)

308 brash 1.6 c
309                 uw=(xu+yu)/sqrt(2.0)
310                 xw=xx
311                 vw=(-xv+yv)/sqrt(2.0)

312 brash 1.7 c
313                 write(*,*)'********************'
314                 write(*,*)xa,ya,za
315                 write(*,*)xb,yb,zb
316                 write(*,*)xu,yu,zu
317                 write(*,*)xx,yx,zx
318                 write(*,*)xv,yv,zv
319                 write(*,*)uw,xw,vw
320                 write(*,*)'********************'
321           c     
322                 anu=(-uw-3.592+104.0)/2.0
323                 anx=(-xw-5.080+84.0)/2.0
324                 anv=(vw-3.592+104.0)/2.0
325           c
326                 write(*,*)anu,anx,anv
327                 nu=anu
328                 nv=anv
329                 nx=anx
330                 if((anu-nu).ge.0.500)nu=nu+1
331                 if((anx-nx).ge.0.500)nx=nx+1
332                 if((anv-nv).ge.0.500)nv=nv+1
333 brash 1.7 c
334                 return
335                 end
336           
337                 subroutine get_drift_distance(xa,ya,za,xb,yb,zb,
338                &           nu,nx,nv,du,dx,dv)
339           
340                 implicit none
341           
342                 real*8 xa,ya,za,xb,yb,zb,du,dx,dv
343                 integer*4 nu,nx,nv
344                 real*8 tphi,ttheta
345                 real*8 uw,xw,vw
346                 real*8 c11,c12,c21,c22,d1,d2,a,zt,xt,yt
347                 real*8 xu,yu,zu
348           
349                 tphi=(xb-xa)/(zb-za)
350                 ttheta=(yb-ya)/(zb-za)
351           
352                 uw=-2.0*nu-3.592+104.0
353                 xw=-2.0*nx-5.080+84.0
354 brash 1.7       vw=2.0*nv+3.592-104.0
355                 zu=(zb+za)/2.0-1.60
356           
357                 write(*,*)uw,xw,vw
358           
359                 c11=tphi-ttheta
360                 c12=sqrt(2.0)
361                 d1=-xa+ya
362                 c21=tphi*tphi+ttheta*ttheta+1.0
363                 c22=-ttheta/sqrt(2.0)+tphi/sqrt(2.0)
364                 d2=uw*(ttheta+tphi)/sqrt(2.0)-xa*tphi-ya*ttheta+zu-za
365           
366                 a=(d1*c21-d2*c11)/(c21*c12-c22*c11)
367                 zt=(d1-c12*a)/c11
368                 write(*,*)zt,a
369           
370                 xt=xa+zt*tphi
371                 yt=ya+zt*ttheta
372                 xu=(uw-a)/sqrt(2.0)
373                 yu=(uw+a)/sqrt(2.0)
374           
375 brash 1.7       zt=zt+za
376           
377                 write(*,*)'Calculating Drift Distance ...'
378                 write(*,*)'A: ',xa,ya,za
379                 write(*,*)'B: ',xb,yb,zb
380                 write(*,*)'T: ',xt,yt,zt
381                 write(*,*)'U: ',xu,yu,zu
382           
383                 du=sqrt((xt-xu)**2+(yt-yu)**2+(zt-zu)**2)
384                 dx=0.0
385                 dv=0.0
386           
387                 write(*,*)'Drift distance = ',du

388 brash 1.6 
389                 return

390 jones 1.1       end

391 brash 1.7 
392                 
393           
394                 
395           
396                 
397                 
398                 

399 jones 1.1