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