source: trunk/MagicSoft/Simulation/Corsika/Mmcs/annih.f

Last change on this file was 286, checked in by harald, 25 years ago
This is the start point for further developments of the Magic Monte Carlo Simulation written by Jose Carlos Gonzales. Now it is under control of one CVS repository for the whole collaboration. Everyone should use this CVS repository for further developments.
File size: 3.1 KB
Line 
1 SUBROUTINE ANNIH
2C VERSION 4.00 -- 26 JAN 1986/1900
3C******************************************************************
4C GAMMA SPECTRUM FOR TWO GAMMA IN-FLIGHT POSITRON ANNIHILATION.
5C USING SCHEME BASED ON HEITLER'S P269-270 FORMULAE
6C THIS ROUTINE SHOULD GIVE THE CORRECT DISTRIBUTION, BUT MORE
7C THOUGHT COULD BE PUT INTO DEVISING A FASTER SCHEME. HOWEVER,
8C SINCE POSITRON ANNIHILATION IN FLIGHT IS RELATIVELY INFREQUENT
9C THIS MAY NOT BE WORTHWHILE.
10C******************************************************************
11 DOUBLE PRECISION PAVIP
12*KEEP,RANDPA.
13 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
14 DOUBLE PRECISION FAC,U1,U2
15 REAL RD(3000)
16 INTEGER ISEED(103,10),NSEQ
17 LOGICAL KNOR
18*KEEP,RUNPAR.
19 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
20 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
21 * MONIOU,MDEBUG,NUCNUC,
22 * CETAPE,
23 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
24 * N1STTR,MDBASE,
25 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
26 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
27 * ,GHEISH,GHESIG
28 COMMON /RUNPAC/ DSN,HOST,USER
29 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
30 REAL STEPFC
31 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
32 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
33 * N1STTR,MDBASE
34 INTEGER CETAPE
35 CHARACTER*79 DSN
36 CHARACTER*20 HOST,USER
37
38 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
39 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
40 * ,GHEISH,GHESIG
41*KEEP,STACKE.
42 COMMON/STACKE/ E,TIME,X,Y,Z,U,V,W,DNEAR,IQ,IGEN,IR,IOBS,LPCTE,NP
43 DOUBLE PRECISION E(60),TIME(60)
44 REAL X(60),Y(60),Z(60),U(60),V(60),W(60),DNEAR(60)
45 INTEGER IQ(60),IGEN(60),IR(60),IOBS(60),LPCTE(60),NP
46*KEND.
47 COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
48 DOUBLE PRECISION PZERO,PRM,PRMT2,RMI,VC
49 COMMON/USEFUL/PZERO,PRM,PRMT2,RMI,VC,RM,MEDIUM,MEDOLD,IBLOBE,ICALL
50 COMMON/ACLOCK/NCLOCK,JCLOCK
51C_____IF (NCLOCK.GT.JCLOCK) THEN
52C______WRITE(MDEBUG,* )' ANNIH: NP=',NP,' IR=',IR(NP),' IOBS=',IOBS(NP)
53C______CALL AUSGB2
54C_____END IF
55 PAVIP=E(NP)+PRM
56 AVIP=PAVIP
57 A=AVIP*RMI
58 AI=1.0/A
59 G=A-1.0
60 T=G-1.0
61 P=SQRT(A*T)
62 POT=P/T
63 EP0I=(A+P)
64331 CONTINUE
65 CALL RMMAR(RD,2,2)
66 RNNO01=RD(1)
67 RNNO02=RD(2)
68 EP=EXP(RNNO01*LOG(EP0I-1.0))/EP0I
69 REJF=1.0-EP+AI*AI*(2.0*G-1.0/EP)
70 IF((RNNO02.LE.REJF))GO TO332
71 GO TO 331
72332 CONTINUE
73 ESG1=AVIP*MAX(EP,1.-EP)
74 E(NP)=ESG1
75 E(NP+1)=PAVIP-E(NP)
76 ESG2=E(NP+1)
77 IQ(NP)=1
78 COSTHE=(ESG1-RM)*POT/ESG1
79 SINTHE=SQRT(MAX(1.0-COSTHE*COSTHE,0.))
80 CALL UPHI(2,1)
81 NP=NP+1
82 IQ(NP)=1
83 COSTHE=(ESG2-RM)*POT/ESG2
84 SINTHE=SQRT(MAX(1.0-COSTHE*COSTHE,0.))
85 CALL UPHI(3,2)
86 RETURN
87 END
Note: See TracBrowser for help on using the repository browser.