1 | SUBROUTINE FIRSTPOINT_HADRON(Lmax) |
---|
2 | IMPLICIT NONE |
---|
3 | #include "event.inc" |
---|
4 | #include "detector.inc" |
---|
5 | #include "erandom.inc" |
---|
6 | #include "constants.inc" |
---|
7 | #include "book_histos.inc" |
---|
8 | #include "cwntvars.inc" |
---|
9 | |
---|
10 | REAL zmin,z_rand,r_rand,x_in_cone,y_in_cone,z_in_cone |
---|
11 | REAL H_R,THETA_R,LENGTH12,SCALAR11,SCALAR12,SCALAR22 |
---|
12 | REAL dir_unit(3),R2(3) |
---|
13 | REAL THETAMAX,Lmax |
---|
14 | REAL Xinteraction,HINT, ATM |
---|
15 | REAL*8 RNDM,aran,bran,cran,r2ran,rinv |
---|
16 | LOGICAL bFirstPoint |
---|
17 | * |
---|
18 | zmin = ISS(3)*SIN(FOV)**2 + |
---|
19 | + COS(FOV)*SQRT(R_EARTH**2-ISS(3)**2*SIN(FOV)**2) |
---|
20 | bFirstPoint = .FALSE. |
---|
21 | do while (.NOT.bFirstPoint) |
---|
22 | * generate point on the Earth surface |
---|
23 | aran=0. |
---|
24 | bran=0. |
---|
25 | cran=0. |
---|
26 | r2ran=1. |
---|
27 | do while(r2ran.gt.0.25) |
---|
28 | aran = ESAFRNDM() -0.5 |
---|
29 | bran = ESAFRNDM() -0.5 |
---|
30 | cran = ESAFRNDM() -0.5 |
---|
31 | r2ran = aran*aran+bran*bran+cran*cran |
---|
32 | enddo |
---|
33 | rinv = R_EARTH/SQRT(r2ran) |
---|
34 | x_in_cone = aran*rinv |
---|
35 | y_in_cone = bran*rinv |
---|
36 | z_in_cone = cran*rinv |
---|
37 | if (z_in_cone.GT.Zmin) bFirstPoint = .TRUE. |
---|
38 | enddo |
---|
39 | z_rand = z_in_cone |
---|
40 | r_rand = x_in_cone**2 + y_in_cone**2 |
---|
41 | |
---|
42 | PI = ACOS(-1.) |
---|
43 | RNDM = ESAFRNDM() |
---|
44 | PHI = PHIR(1) + (PHIR(2)-PHIR(1))*RNDM |
---|
45 | H_R = SQRT(r_rand**2+z_rand**2) - R_EARTH |
---|
46 | * ... SIN(2THETA) RANDOM THETA ... |
---|
47 | RNDM = ESAFRNDM() |
---|
48 | THETA_R = 0.5*ACOS(COS(2*THETAR(1)/180*PI) - |
---|
49 | + RNDM*(COS(2*THETAR(1)/180*PI)-COS(2*THETAR(2)/180*PI))) |
---|
50 | THETA = THETA_R/PI*180 |
---|
51 | * ... Find far first point |
---|
52 | dir_unit(1) = SIN(THETA_R)*COS(PHI/180.*PI) |
---|
53 | dir_unit(2) = SIN(THETA_R)*SIN(PHI/180.*PI) |
---|
54 | dir_unit(3) = COS(THETA_R) |
---|
55 | R2(1) = x_in_cone |
---|
56 | R2(2) = y_in_cone |
---|
57 | R2(3) = z_rand |
---|
58 | |
---|
59 | CALL VECTORDOT(3,dir_unit,R2,SCALAR12) |
---|
60 | CALL VECTORDOT(3,R2,R2,SCALAR22) |
---|
61 | SCALAR22 = SQRT(SCALAR22) |
---|
62 | LENGTH12 = SQRT(80*(SCALAR22+80.) + SCALAR12**2) - SCALAR12 |
---|
63 | RINT(1) = R2(1) + LENGTH12*SIN(THETA/180*PI)*COS(PHI/180.*PI) |
---|
64 | RINT(2) = R2(2) + LENGTH12*SIN(THETA/180*PI)*SIN(PHI/180.*PI) |
---|
65 | RINT(3) = R2(3) + LENGTH12*COS(THETA/180*PI) |
---|
66 | |
---|
67 | Xinteraction = 0. |
---|
68 | * ... RANDOM FIRST INTERACTION DEPTH ... |
---|
69 | CALL First_Depth_Hadron |
---|
70 | do while (Xinteraction.LT.X1) |
---|
71 | RINT(1) = RINT(1) - 0.1*SIN(THETA/180*PI)*COS(PHI/180.*PI) |
---|
72 | RINT(2) = RINT(2) - 0.1*SIN(THETA/180*PI)*SIN(PHI/180.*PI) |
---|
73 | RINT(3) = RINT(3) - 0.1*COS(THETA/180*PI) |
---|
74 | CALL VECTORDOT(3,RINT,RINT,SCALAR22) |
---|
75 | HINT = SQRT(SCALAR22) - R_EARTH |
---|
76 | Xinteraction = Xinteraction + 0.1*ATM(HINT) |
---|
77 | enddo |
---|
78 | * ... THETA,PHI ARE RETURNED IN DEGREES ... |
---|
79 | CALL ACCEPTANCE(Lmax,ACCEPTANCE_TYPE) |
---|
80 | RINT(3) = RINT(3) - R_EARTH |
---|
81 | * |
---|
82 | END |
---|
83 | |
---|