COSMOS v7.655  COSMOSv7655
(AirShowerMC)
cifXObsSite.f
Go to the documentation of this file.
1  subroutine cifxobssite(nextwhere)
2 ! see if MovedTrack crosses an observation depth.
3 ! or go out of boundry.
4 ! if to be observed, reset MovedTrack information
5  implicit none
6 #include "Ztrack.h"
7 #include "Ztrackp.h"
8 #include "Ztrackv.h"
9 #include "Zincidentv.h"
10 #include "Zobs.h"
11 #include "Zobsp.h"
12 #include "Zobsv.h"
13 #include "Zearth.h"
14 #include "Zcode.h"
15 !
16  integer nextwhere ! output next where. to be used after ptcl
17  ! crosses an observation level.
18 ! for perpendicular case
19 ! cross border
20 ! 1 O BorderL< Moved.Height <BorderH
21 ! 2 X //
22 ! 3 O M.H > B.H
23 ! 4 O M.H < B.L
24 ! 5 X M.H > B.H
25 ! 6 X M.H < B.L
26 !
27  integer loc
28  real*8 leng, dedt, cosfromaxis
29  real*8 clen2thick
30 
31  logical cross, seeupper
32 
33  type(coord)::xyz1
34  type(coord)::xyz2
35  type(coord)::dircos
36  type(coord)::diffvec
37  type(coord)::unitv
38  type(coord)::unitp
39 
40  type(coord)::temp
41 
42  real*8 dummylen, r0, rp, obsrl, cosp, dt2, m2
43  real*8 coscr2, r0r
44  real*8 dedtF
45 !
46  cross = .false.
47 
48  loc = trackbefmove%where
49 
50  if(intinfarray(processno)%length .gt. 0. ) then
51 ! later use
52  call cdiffvec(trackbefmove%pos%xyz,
53  * movedtrack%pos%xyz, diffvec)
54 ! make it to the unit vector vector: Bef 1--->2 Moved
55  call c3dv2ddcos(diffvec, unitv, dummylen)
56  if( obsplane .ne. notused ) then
57  if( abs(obsplane) .eq. perpendicular ) then
58  call cifxperpen(unitv, cross, leng, seeupper)
59  else
60  call cifxhorizon(diffvec, unitv, cross, leng, seeupper)
61  endif
62  endif
63  endif
64 
65  if(cross) then
66  intinfarray(processno)%length = leng
67  intinfarray(processno)%thickness = clen2thick(
68  * trackbefmove%pos%radiallen - eradius,
69  * trackbefmove%vec%coszenith, leng)
70 
71  call cmovestreight(leng, unitv)
72  movedtrack%vec%w = unitv
73 
74 ! bef. replace, w for low energy electron
75 ! may be largely scattered and even
76 ! backward at the path end.
77 ! At Xing point with obs. depth, this
78 ! is not good.
79 ! if(MovedTrack.pos.height .le. BorderHeightL) then
80 ! write(0,*) MovedTrack.pos.height, BorderHeightL
81 ! MoveStat = BorderL
82 ! MovedTrack.where = NoOfSites+1
83 
84 
85  if( movestat .ne. borderh .and. movestat .ne. borderl) then
86  movestat = tobeobserved
87  if(seeupper) then
88  nextwhere = loc - 1
89  movedtrack%where = nextwhere
90  else
91  movedtrack%where = loc
92  nextwhere = loc + 1
93  endif
94  endif
95 
96  if(movedtrack%p%charge /= 0) then
97 ! sinc path has been cut, recalculate Energy loss
98  ! get de/dx computed previously
99  call cqelossrate(dedt, dedtf)
100  energyloss = dedt * intinfarray(processno)%thickness
101  movedtrack%p%fm%p(4) = trackbefmove%p%fm%p(4)-energyloss
102  call cresetmombyedir(movedtrack)
103 ! reset Energyloss since E in MovedTrack might have been
104 ! reset to mass if energyloss was too big.
105  energyloss = trackbefmove%p%fm%p(4)- movedtrack%p%fm%p(4)
106  else
107 ! charge = 0
108  ! nothing to do. even path is cut, momenutum etc
109 ! dose not change.
110  endif
111  else
112  nextwhere = loc
113  movedtrack%where = loc
114  endif
115  end
116 
117  subroutine cifxhorizon(diffvec, unitv, cross, leng, seeupper)
118  implicit none
119 #include "Ztrack.h"
120 #include "Ztrackp.h"
121 #include "Ztrackv.h"
122 #include "Zincidentv.h"
123 #include "Zobs.h"
124 #include "Zobsp.h"
125 #include "Zobsv.h"
126 #include "Zearth.h"
127 #include "Zcode.h"
128 !
129 
130  integer loc
131  real*8 leng
132  real*8 clen2thick
133 
134  logical cross, seeupper
135 
136  type(coord)::diffvec
137  type(coord)::unitv
138  type(coord)::unitp
139 
140 
141  real*8 dummylen, r0, rp, obsrl, cosp, dt2, m2
142  real*8 coscr2, r0r
143 !
144 
145  loc = trackbefmove%where
146 
147  r0 = obssites(loc)%pos%radiallen
148  rp = trackbefmove%pos%radiallen
149  dt2 = rp**2 - r0**2
150  if(dt2 .lt. 0.d0) then
151  if(dt2 .lt. -1.d0) then
152  write(0,*) ' loc =',loc, ' r0=',r0, ' rp=',rp,
153  * ' dt2=',dt2, ' code=',trackbefmove%p%code,
154  * ' rp%h=', trackbefmove%pos%height,
155  * obssites(loc)%pos%height
156  call cerrormsg('position info invalid',0)
157  endif
158  dt2 = abs(dt2)
159  endif
160 ! this is cos^2 of zenith for a ptcl which touches the inner sphere (loc-th sphere)
161  coscr2 = dt2 /rp**2
162 
163  call c3dv2ddcos(trackbefmove%pos%xyz, unitp, dummylen)
164  call cscalerprod(unitv, unitp, cosp)
165  cosp = -cosp
166  seeupper = cosp .le. 0
167  if(cosp .gt. 0) then
168 ! down going but not cross the lower sphere
169  seeupper = cosp**2 .le. coscr2
170  endif
171 
172 
173  if(.not. seeupper ) then
174  if(movedtrack%pos%radiallen .le. r0) then
175 ! cross ; goes down
176  cross = .true.
177  else
178 ! get moved length ^2
179  call cscalerprod(diffvec, diffvec, m2)
180  if(m2 .gt. dt2) then
181 ! double cross; first one is the nearest cross
182  cross =.true.
183  endif
184  endif
185  else
186 ! may see upper
187  if(movedtrack%pos%radiallen .ge.
188  * obssites(loc-1)%pos%radiallen) then
189 ! if loc=1, this means going outside of the borderH
190  cross =.true.
191  endif
192  endif
193 
194  if(cross) then
195 ! get length, leng, from the starting point to
196 ! the observation depth
197  call cscalerprod(unitv, trackbefmove%pos%xyz, r0r)
198  if(seeupper) then
199  obsrl = obssites(loc-1 )%pos%radiallen
200  if(loc .eq. 1) movestat=borderh
201  else
202  obsrl = obssites(loc)%pos%radiallen
203  if(loc .eq. noofsites+1 )then
204 ! this seems not needed.
205  movestat = borderl
206  endif
207  endif
208 ! get distance to the crossing point
209  if(rp .gt. obsrl ) then
210  leng = -r0r - sqrt(r0r**2 -
211  * (rp**2-obsrl**2))
212  else
213  leng = -r0r + sqrt(r0r**2 -
214  * (rp**2 - obsrl**2))
215  endif
216  leng = leng + 0.001d0
217  endif
218  end
219  subroutine cifxperpen(unitv, cross, leng, seeupper)
220  implicit none
221 #include "Ztrack.h"
222 #include "Ztrackp.h"
223 #include "Ztrackv.h"
224 #include "Zincidentv.h"
225 #include "Zobs.h"
226 #include "Zobsp.h"
227 #include "Zobsv.h"
228 #include "Zearth.h"
229 #include "Zcode.h"
230 !
231 
232 
233  real*8 leng
234  logical cross, seeupper
235 
236  real*8 clen2thick, cosfromaxis, leng2
237 
238  integer loc, i, icon
239  type(coord)::xyz1
240  type(coord)::xyz2
241  type(coord)::dircos
242  type(coord)::unitv
243 
244  loc = trackbefmove%where
245 !///////////
246 ! write(0,*) ' loc =', loc
247 !//////////
248 !
249 ! observation plane is perpendicular to primary
250 ! convert coord into 1ry system
251  call cxyz2prim(obssites(noofsites)%pos%xyz,
252  * movedtrack%pos%xyz, xyz2)
253  call cxyz2prim(obssites(noofsites)%pos%xyz,
254  * trackbefmove%pos%xyz, xyz1)
255 ! call cscalerProd(TrackBefMove.vec.w, DcAtObsXyz,
256 ! * cosfromaxis)
257  call cscalerprod(unitv, dcatobsxyz,
258  * cosfromaxis)
259 
260  seeupper = cosfromaxis .le. 0.
261  if(.not. seeupper) then
262 ! current ptcl is down going
263  if( loc .gt. noofsites ) then
264  if( movedtrack%pos%height .le. borderheightl ) then
265 ! 6
266  movestat = borderl
267  movedtrack%where = noofsites+1
268  cross = .true.
269  elseif( movedtrack%pos%height .ge. borderheighth ) then
270 ! 5
271  movestat = borderh
272  movedtrack%where = 0
273  cross = .true.
274  endif
275  else
276  if(xyz1%r(3) .gt. obssites(loc)%zpl .and.
277  * xyz2%r(3) .le. obssites(loc)%zpl) then
278  cross = .true.
279  elseif( movedtrack%pos%height .le. borderheightl ) then
280 ! escape from side 6
281  cross = .true.
282  movestat = borderl
283  movedtrack%where = noofsites + 1
284  elseif( movedtrack%pos%height .ge. borderheighth ) then
285 ! escape from side 5
286  cross = .true.
287  movestat = borderh
288  movedtrack%where = 0
289  endif
290  endif
291  else
292 ! upgooing
293  if(loc .eq. 1) then
294  if( movedtrack%pos%height .ge. borderheighth ) then
295 ! 5
296  movestat = borderh
297  movedtrack%where = 0
298  cross = .true.
299  elseif(movedtrack%pos%height .le. borderheightl ) then
300 ! 6
301  movestat = borderl
302  movedtrack%where = noofsites+1
303  cross = .true.
304  endif
305  else
306  if(xyz1%r(3) .lt. obssites(loc-1)%zpl .and.
307  * xyz2%r(3) .ge. obssites(loc-1)%zpl ) then
308  cross = .true.
309  elseif( movedtrack%pos%height .ge. borderheighth ) then
310 ! escaep from side 5
311  cross = .true.
312  movestat = borderh
313  movedtrack%where = 0
314  elseif( movedtrack%pos%height .le. borderheightl ) then
315 ! escape from side 6
316  cross = .true.
317  movestat = borderl
318  movedtrack%where = noofsites + 1
319  endif
320  endif
321  endif
322 !/////////////
323 ! write(0,*) ' see upper=',seeupper, ' cross=',cross
324 ! write(0,*) ' where=', MovedTrack.where
325 ! write(0,*) ' MoveStat=',MoveStat
326 ! write(0,*) ' BorderL=',BorderL, 'BorderH=',BorderH
327 !/////////////
328  if(.not. cross ) then
329 ! * .and. MovedTrack.pos.height .gt. BorderHeightL
330 ! * .and. MovedTrack.pos.height .lt. BorderHeightH ) then
331 ! most common case; remain in the same layer; nothing to do
332  elseif( cross .and. movestat .eq. borderl ) then
333 ! get length to the x-ing point
334  call cxplsph(
335  * trackbefmove%pos%xyz%r(1), trackbefmove%pos%xyz%r(2),
336  * trackbefmove%pos%xyz%r(3),
337  * unitv%r(1), unitv%r(2), unitv%r(3),
338  * borderheightl+eradius, leng, icon)
339 !/////////////
340 ! write(0,*) ' 1 leng=',leng, ' icon=',icon
341 !//////////////
342  leng = leng + 0.01d0
343 ! MovedTrack.where = NoOfSites + 1
344  elseif( cross .and. movestat .eq. borderh ) then
345 ! get length to the x-ing point
346  call cxplsph(
347  * trackbefmove%pos%xyz%r(1), trackbefmove%pos%xyz%r(2),
348  * trackbefmove%pos%xyz%r(3),
349  * unitv%r(1), unitv%r(2), unitv%r(3),
350  * borderheighth+eradius, leng, icon)
351 !//////////
352 ! write(0,*) ' borderH leng=',leng, 'icon=',icon
353 !//////////////
354  leng = leng + 0.01d0
355  else
356 !
357 ! get length to the x-ing point
358  call cxyz2primd(unitv, dircos)
359  if(seeupper) then
360  leng = (obssites(loc-1)%zpl - xyz1%r(3))/dircos%r(3)
361  else
362  leng = (obssites(loc)%zpl - xyz1%r(3))/dircos%r(3)
363  endif
364 !/////////
365 ! write(0,*) ' 3 leng=',leng
366 !//////////
367  leng = abs(leng)
368 ! leng = leng + abs(0.01d0/dircos.r(3))
369  leng = leng + 0.01d0
370 
371  if( movedtrack%pos%height .gt. borderheightl .and.
372  * movedtrack%pos%height .lt. borderheighth ) then
373 ! 1
374 ! most common cross case; if the track length is too long
375 ! this logic is not complete, but we neglect such case
376  elseif( movedtrack%pos%height .ge. borderheighth ) then
377 ! 3
378 ! must find shorter x-ing point
379  call cxplsph(
380  * trackbefmove%pos%xyz%r(1), trackbefmove%pos%xyz%r(2),
381  * trackbefmove%pos%xyz%r(3),
382  * unitv%r(1), unitv%r(2), unitv%r(3),
383  * borderheighth+eradius, leng2, icon)
384 !///////////
385 ! write(0,*) ' 3 leng2=',leng2, ' icon=',icon
386 !/////////////
387  if( leng2 .lt. leng ) then
388  leng = leng2
389  movestat = borderh
390  endif
391  elseif( movedtrack%pos%height .le. borderheightl ) then
392 ! 4
393 ! must find shorter x-ng point
394  call cxplsph(
395  * trackbefmove%pos%xyz%r(1), trackbefmove%pos%xyz%r(2),
396  * trackbefmove%pos%xyz%r(3),
397  * unitv%r(1), unitv%r(2), unitv%r(3),
398  * borderheightl+eradius, leng2, icon)
399 !///////////
400 ! write(0,*) ' 4 leng2=',leng2, ' icon=',icon
401 !/////////////
402  if( leng2 .lt. leng ) then
403  leng = leng2
404  movestat = borderl
405  endif
406  else
407 ! should not come
408  write(0, *) 'logic error'
409  stop
410  endif
411  endif
412  end
413  subroutine cdiffvec(r1, r2, diff)
414 ! get diff= r2-r1 as 3 D vector
415 #include "Zcoord.h"
416  type(coord)::r1
417  type(coord)::r2 ! input
418  type(coord)::diff ! output.
419 
420  integer i
421 
422  do i = 1, 3
423  diff%r(i) = r2%r(i) - r1%r(i)
424  enddo
425  end
426  subroutine cxplsph(x0, y0, z0, l, m, n, r, el, icon)
427 ! this is the same as kxplsph in Epics
428  implicit none
429  real*8 x0, y0, z0 ! input. the line passes this point
430  real*8 l, m, n ! input. direc cos. of the line
431  real*8 r ! input. radius of the sphere
432  real*8 el ! output. el>=0 distance to the
433  ! sphere from x0,y0,z0
434  integer icon ! output. icon =0. x-point exists
435  ! x0,.. is inside
436  ! icon = 1 x-point exists
437  ! x0.. is outside
438  ! =-1. no x-point
439 
440  real*8 rsqr, r0l, d
441  integer icon1, icon2
442 
443  rsqr = x0**2 + y0**2 + z0**2 -r**2
444  if(rsqr .le. 0.) then
445 ! inside
446  icon2 = 0
447  else
448  icon2 = 1
449  endif
450  r0l = x0*l + y0*m + z0*n
451  d = r0l**2 - rsqr
452  if(d .ge. 0.) then
453  d = sqrt(d)
454  el = -r0l - d
455  if(el .ge. 0.) then
456  icon1 = 0
457  else
458  el = -r0l + d
459  if(el .ge. 0.) then
460  icon1 = 0
461  else
462  icon1 = 1
463  endif
464  endif
465  else
466  icon1 = 1
467  endif
468 !
469  if(icon2 .eq. 0) then
470  icon = 0
471  elseif(icon1 .eq. 0) then
472  icon = 1
473  else
474  icon = -1
475  endif
476  end
cerrormsg
subroutine cerrormsg(msg, needrtn)
Definition: cerrorMsg.f:4
cxplsph
subroutine cxplsph(x0, y0, z0, l, m, n, r, el, icon)
Definition: cifXObsSite.f:427
notused
const int notused
Definition: Zobs.h:16
cifxhorizon
subroutine cifxhorizon(diffvec, unitv, cross, leng, seeupper)
Definition: cifXObsSite.f:118
true
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon true
Definition: cblkElemag.h:7
cxyz2prim
subroutine cxyz2prim(base, a, b)
Definition: cxyz2prim.f:5
tobeobserved
const int tobeobserved
Definition: Ztrackv.h:94
borderl
const int borderl
Definition: Ztrackv.h:97
borderh
const int borderh
Definition: Ztrackv.h:98
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
cifxperpen
subroutine cifxperpen(unitv, cross, leng, seeupper)
Definition: cifXObsSite.f:220
cifxobssite
subroutine cifxobssite(nextwhere)
Definition: cifXObsSite.f:2
cresetmombyedir
subroutine cresetmombyedir(aTrack)
Definition: cresetMom.f:19
c3dv2ddcos
subroutine c3dv2ddcos(a, b, len)
Definition: c3DV2DDCos.f:3
cxyz2primd
subroutine cxyz2primd(a, b)
Definition: cxyz2prim.f:37
false
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon ! knockon is considered Obsolete *PhotoProd false
Definition: cblkElemag.h:7
cmovestreight
subroutine cmovestreight(leng, dir)
Definition: ccompPathEnd.f:40
cdiffvec
subroutine cdiffvec(r1, r2, diff)
Definition: cifXObsSite.f:414
coord
Definition: Zcoord.h:43
perpendicular
! Zobs h header file for observation sites definition ! integer * perpendicular
Definition: Zobs.h:4
cscalerprod
subroutine cscalerprod(a, b, c)
Definition: cscalerProd.f:4