source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/externals/slast/src/kernel/firstpoint_neutrino.F @ 114

Last change on this file since 114 was 114, checked in by moretto, 11 years ago

actual version of ESAF at CCin2p3
File size: 2.7 KB
Line 
1      SUBROUTINE FIRSTPOINT_NEUTRINO(Lmax)
2      IMPLICIT NONE
3#include "detector.inc"
4#include "event.inc"
5#include "erandom.inc"
6      INTEGER ID_TEST
7      LOGICAL TEST_FOV
8      DATA TEST_FOV/.FALSE./
9      SAVE BETA
10      REAL Zmin,R,BETA,BETAR,h,PI,h_r
11      REAL ATM,THETAMAX,Lmax,COSBETAMIN
12      REAL depth
13      REAL*8 RNDM
14*
15      PI = ACOS(-1.)
16      Zmin = ISS(3)*SIN(FOV)**2+
17     +       COS(FOV)*SQRT(R_EARTH**2-ISS(3)**2*SIN(FOV)**2)
18      R = SQRT(R_EARTH**2-Zmin**2)
19      BETA = R/R_EARTH        ! SIN(BETA)
20      COSBETAMIN = SQRT(1.-BETA**2)
21        IF(TEST_FOV) THEN
22          ID_TEST = 900
23          CALL HBOOK1(ID_TEST+1,'height',100,0.,50.,0.)
24          CALL HBOOK1(ID_TEST+2,'x',100,-250.,250.,0.)
25          CALL HBOOK1(ID_TEST+3,'y',100,-250.,250.,0.)
26          CALL HBOOK1(ID_TEST+4,'z',100,R_EARTH,R_EARTH+50.
27     &         ,0.)
28          CALL HBOOK1(ID_TEST+5,'BETA',100,0.,0.05,0.)
29          CALL HBOOK1(ID_TEST+6,'PHI LOCAL',100,0.,360.,0.)
30          CALL HBOOK1(ID_TEST+7,'THETA',100,0.,110.,0.)
31          CALL HBOOK1(ID_TEST+8,'PHI ',100,0.,360.,0.)
32          CALL HBOOK1(ID_TEST+9,'acceptance',10,0.,10.,0.)
33        ENDIF
34*     ... GENERATE RANDOM LOCAL HEIGHT ACCORDING TO THE AIR DENSITY
35      H_R = -7.8*LOG(ESAFRNDM())
36*     ... GENERATE RANDOM SIN(BETA) - BETA IS THE ANGLE FROM EARTH
37*     ... CENTER, AND PHI - ASIMUTAL ANGLE
38      RNDM = ESAFRNDM()
39      PHI = PHIR(1) + (PHIR(2)-PHIR(1))*RNDM
40      IF(TEST_FOV) CALL HF1(ID_TEST+6,PHI,1.)
41      RNDM = ESAFRNDM()
42      BETAR = ACOS(COSBETAMIN + (1.-COSBETAMIN)*RNDM)
43*     ... THE FIRST INTERACTION POINT ...
44      RINT(1) = (R_EARTH+H_R)*BETAR*COS(PHI/180*PI)
45      RINT(2) = (R_EARTH+H_R)*BETAR*SIN(PHI/180*PI)
46      RINT(3) = (R_EARTH+H_R)*SQRT(1.-BETAR**2)
47*     ... GENERATE RANDOM DIRECTION ...
48      RNDM = ESAFRNDM() 
49      PHI = PHIR(1) + (PHIR(2)-PHIR(1))*RNDM
50      RNDM = ESAFRNDM() 
51      THETA = 0.5*ACOS(COS(2*THETAR(1)/180*PI) -
52     +    RNDM*(COS(2*THETAR(1)/180*PI)-COS(2*THETAR(2)/180*PI)))
53*     ... COMPUTE X1 ...
54      PHI   = PHI/180*PI
55      X1 = DEPTH(RINT)
56      THETA = THETA/PI*180
57      PHI   = PHI/PI*180
58*     ... SAVE INFO FOR LATER CHECKING ...
59      IF(TEST_FOV) CALL HF1(ID_TEST+1,H_R,1.)
60      IF(TEST_FOV) CALL HF1(ID_TEST+2,RINT(1),1.)
61      IF(TEST_FOV) CALL HF1(ID_TEST+3,RINT(2),1.)
62      IF(TEST_FOV) CALL HF1(ID_TEST+4,RINT(3),1.)
63      IF(TEST_FOV) CALL HF1(ID_TEST+5,BETAR,1.)
64
65      IF(TEST_FOV) CALL HF1(ID_TEST+7,THETA,1.)
66      IF(TEST_FOV) CALL HF1(ID_TEST+8,PHI,1.)
67*     ...  THETA,PHI ARE RETURNED IN DEGREES ...
68      CALL ACCEPTANCE(Lmax,ACCEPTANCE_TYPE)
69      IF(TEST_FOV) CALL HF1(ID_TEST+9,ACCEPTANCE_TYPE,1.)
70      RINT(3) = RINT(3) - R_EARTH
71*
72      END
Note: See TracBrowser for help on using the repository browser.