| 1 | SUBROUTINE PIGEN2
|
|---|
| 2 | C
|
|---|
| 3 | C*********************************************************************
|
|---|
| 4 | C DESIGN : D. HECK IK3 FZK KARLSRUHE
|
|---|
| 5 | C DATE : JUL 31, 1989
|
|---|
| 6 | C*********************************************************************
|
|---|
| 7 | C THIS SUBROUTINE DESCRIBES THE PHOTONUCLEAR REACTION
|
|---|
| 8 | C GAMMA + NUCLEON -----> PION + PION + NUCLEON
|
|---|
| 9 | C*********************************************************************
|
|---|
| 10 | DOUBLE PRECISION BETA,DUMMY,ECM,ENUCL,GAMMA,PEIG,PTRANS
|
|---|
| 11 | *KEEP,ELABCT.
|
|---|
| 12 | COMMON /ELABCT/ ELCUT
|
|---|
| 13 | DOUBLE PRECISION ELCUT(4)
|
|---|
| 14 | *KEEP,PAM.
|
|---|
| 15 | COMMON /PAM/ PAMA,SIGNUM
|
|---|
| 16 | DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
|
|---|
| 17 | *KEEP,PARPAR.
|
|---|
| 18 | COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
|
|---|
| 19 | * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
|
|---|
| 20 | DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
|
|---|
| 21 | * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
|
|---|
| 22 | INTEGER ITYPE,LEVL
|
|---|
| 23 | *KEND.
|
|---|
| 24 | DOUBLE PRECISION PI0MSQ
|
|---|
| 25 | COMMON/PION/PI0MSQ,PITHR,PICMAS,PI0MAS,AMASK0,AMASKC,AMASPR,AMASNT
|
|---|
| 26 | *
|
|---|
| 27 | *KEEP,RANDPA.
|
|---|
| 28 | COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
|
|---|
| 29 | DOUBLE PRECISION FAC,U1,U2
|
|---|
| 30 | REAL RD(3000)
|
|---|
| 31 | INTEGER ISEED(103,10),NSEQ
|
|---|
| 32 | LOGICAL KNOR
|
|---|
| 33 | *KEEP,RUNPAR.
|
|---|
| 34 | COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
|
|---|
| 35 | * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
|
|---|
| 36 | * MONIOU,MDEBUG,NUCNUC,
|
|---|
| 37 | * CETAPE,
|
|---|
| 38 | * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
|
|---|
| 39 | * N1STTR,MDBASE,
|
|---|
| 40 | * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
|
|---|
| 41 | * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
|
|---|
| 42 | * ,GHEISH,GHESIG
|
|---|
| 43 | COMMON /RUNPAC/ DSN,HOST,USER
|
|---|
| 44 | DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
|
|---|
| 45 | REAL STEPFC
|
|---|
| 46 | INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
|
|---|
| 47 | * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
|
|---|
| 48 | * N1STTR,MDBASE
|
|---|
| 49 | INTEGER CETAPE
|
|---|
| 50 | CHARACTER*79 DSN
|
|---|
| 51 | CHARACTER*20 HOST,USER
|
|---|
| 52 |
|
|---|
| 53 | LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
|
|---|
| 54 | * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
|
|---|
| 55 | * ,GHEISH,GHESIG
|
|---|
| 56 | *KEEP,STACKE.
|
|---|
| 57 | COMMON/STACKE/ E,TIME,X,Y,Z,U,V,W,DNEAR,IQ,IGEN,IR,IOBS,LPCTE,NP
|
|---|
| 58 | DOUBLE PRECISION E(60),TIME(60)
|
|---|
| 59 | REAL X(60),Y(60),Z(60),U(60),V(60),W(60),DNEAR(60)
|
|---|
| 60 | INTEGER IQ(60),IGEN(60),IR(60),IOBS(60),LPCTE(60),NP
|
|---|
| 61 | *KEND.
|
|---|
| 62 | COMMON/UPHIIN/SINC0,SINC1,SIN0(20002),SIN1(20002)
|
|---|
| 63 | COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
|
|---|
| 64 | COMMON/ACLOCK/NCLOCK,JCLOCK
|
|---|
| 65 | C_____IF (NCLOCK.GT.JCLOCK) THEN
|
|---|
| 66 | C______WRITE(MDEBUG,* )' PIGEN2:NP=',NP,' IR=',IR(NP),' IOBS=',IOBS(NP)
|
|---|
| 67 | C______CALL AUSGB2
|
|---|
| 68 | C_____END IF
|
|---|
| 69 | IF(DEBUG)WRITE(MDEBUG,*)'PIGEN2: E=',E(NP)
|
|---|
| 70 | PEIG=E(NP)
|
|---|
| 71 | C*** NUMBERS AT THE VARIABLES MEAN :
|
|---|
| 72 | C*** 1 INCOMING GAMMA RAY
|
|---|
| 73 | C*** 2 HIT NUCLEON
|
|---|
| 74 | C*** 3 FIRST PRODUCED PION
|
|---|
| 75 | C*** 4 SECOND PRODUCED PION
|
|---|
| 76 | C*** 5 RECOILING NUCLEON
|
|---|
| 77 | CALL RMMAR(RD,2,2)
|
|---|
| 78 | RNNO81=RD(1)
|
|---|
| 79 | RNNO82=RD(2)
|
|---|
| 80 | C*** LOOK WHICH TYPE OF REACTION
|
|---|
| 81 | C*** 0.49923 IS THE FRACTION OF PROTONS IN AIR
|
|---|
| 82 | IF (RNNO81.LE.0.49923) THEN
|
|---|
| 83 | C *** HIT NUCLEON IS PROTON
|
|---|
| 84 | AMASS2=AMASPR
|
|---|
| 85 | C *** BRANCHING FOR COLLISION WITH PROTON
|
|---|
| 86 | IF (RNNO82.LE.0.3) THEN
|
|---|
| 87 | C *** PI(0) + PI(0) + PROTON
|
|---|
| 88 | IQ(NP)= 7
|
|---|
| 89 | IQ(NP+1)= 7
|
|---|
| 90 | IQ(NP+2)= 14
|
|---|
| 91 | ELSE IF(RNNO82.LE.0.6) THEN
|
|---|
| 92 | C *** PI(+) + PI(-) + PROTON
|
|---|
| 93 | IQ(NP)= 8
|
|---|
| 94 | IQ(NP+1)= 9
|
|---|
| 95 | IQ(NP+2)= 14
|
|---|
| 96 | ELSE
|
|---|
| 97 | C *** PI(+) + PI(0) + NEUTRON
|
|---|
| 98 | IQ(NP)= 8
|
|---|
| 99 | IQ(NP+1)= 7
|
|---|
| 100 | IQ(NP+2)= 13
|
|---|
| 101 | END IF
|
|---|
| 102 | ELSE
|
|---|
| 103 | C *** HIT NUCLEON IS NEUTRON
|
|---|
| 104 | C *** BRANCHING FOR COLLISION WITH NEUTRON
|
|---|
| 105 | AMASS2=AMASNT
|
|---|
| 106 | IF (RNNO82.LE.0.3) THEN
|
|---|
| 107 | C *** PI(0) + PI(0) + NEUTRON
|
|---|
| 108 | IQ(NP)= 7
|
|---|
| 109 | IQ(NP+1)= 7
|
|---|
| 110 | IQ(NP+2)= 13
|
|---|
| 111 | ELSE IF(RNNO82.LE.0.6) THEN
|
|---|
| 112 | C *** PI(+) + PI(-) + NEUTRON
|
|---|
| 113 | IQ(NP)= 8
|
|---|
| 114 | IQ(NP+1)= 9
|
|---|
| 115 | IQ(NP+2)= 13
|
|---|
| 116 | ELSE
|
|---|
| 117 | C *** PI(-) + PI(0) + PROTON
|
|---|
| 118 | IQ(NP)= 9
|
|---|
| 119 | IQ(NP+1)= 7
|
|---|
| 120 | IQ(NP+2)= 14
|
|---|
| 121 | END IF
|
|---|
| 122 | END IF
|
|---|
| 123 | C*** CALCULATE AUXILIARY PARAMETERS
|
|---|
| 124 | ECM=SQRT(AMASS2*(AMASS2+2.D0*PEIG))
|
|---|
| 125 | C*** NOTE: THE ENERGIES IN EGS ARE IN MEV, IN CORSIKA IN GEV
|
|---|
| 126 | AMASS3=PAMA(IQ(NP))*1.D3
|
|---|
| 127 | AMASS4=PAMA(IQ(NP+1))*1.D3
|
|---|
| 128 | AMASS5=PAMA(IQ(NP+2))*1.D3
|
|---|
| 129 | AUX1=(AMASS3+AMASS4)**2
|
|---|
| 130 | AUX2A=(ECM - AMASS5)**2
|
|---|
| 131 | AUX2=AUX2A-AUX1
|
|---|
| 132 | AUX3=(AMASS3+AMASS5)**2
|
|---|
| 133 | AUX4A=(ECM - AMASS4)**2
|
|---|
| 134 | AUX4=AUX4A-AUX3
|
|---|
| 135 | AUX5=AMASS3**2-AMASS4**2
|
|---|
| 136 | AUX6=ECM**2-AMASS5**2
|
|---|
| 137 | AUX7=0.5/ECM
|
|---|
| 138 | AUX8=(ECM - AMASS3)**2
|
|---|
| 139 | BETA=PEIG/(AMASS2+PEIG)
|
|---|
| 140 | GAMMA=2.*(PEIG+AMASS2)*AUX7
|
|---|
| 141 | 230 CONTINUE
|
|---|
| 142 | CALL RMMAR(RD,2,2)
|
|---|
| 143 | RNNO84=RD(1)
|
|---|
| 144 | RNNO85=RD(2)
|
|---|
| 145 | C*** ARE INVARIANT MASS SQUARES INSIDE BOUNDARY OF DALITZ PLOT?
|
|---|
| 146 | AM34SQ=AUX2*RNNO84+AUX1
|
|---|
| 147 | AM35SQ=AUX4*RNNO85+AUX3
|
|---|
| 148 | AM34I=0.5/SQRT(AM34SQ)
|
|---|
| 149 | E3STAR=(AUX5+AM34SQ)*AM34I
|
|---|
| 150 | E5STAR=(AUX6-AM34SQ)*AM34I
|
|---|
| 151 | ROOT1=SQRT(E3STAR**2-AMASS3**2)
|
|---|
| 152 | ROOT2=SQRT(E5STAR**2-AMASS5**2)
|
|---|
| 153 | C*** REJECT RANDOM NUMBERS, IF NOT INSIDE KINEMATIC BOUNDARY
|
|---|
| 154 | DISCR=AM35SQ-(E3STAR+E5STAR)**2
|
|---|
| 155 | IF((DISCR.GT.-(ROOT1-ROOT2)**2))GOTO 230
|
|---|
| 156 | IF((DISCR.LT.-(ROOT1+ROOT2)**2))GOTO 230
|
|---|
| 157 | C*** E3CM,E4CM,E5CM ARE ENERGIES IN C.M. SYSTEM
|
|---|
| 158 | E4CM=(ECM**2+AMASS4**2-AM35SQ)*AUX7
|
|---|
| 159 | E5CM=(ECM**2+AMASS5**2-AM34SQ)*AUX7
|
|---|
| 160 | C*** NOW TAKE PION WITH HIGHEST ENERGY AS PARTICLE 3
|
|---|
| 161 | E3CM=ECM-E4CM-E5CM
|
|---|
| 162 | IF (E4CM.GT.E3CM) THEN
|
|---|
| 163 | C *** INTERCHANGE PARTICLE 3 AND 4
|
|---|
| 164 | HELP=E3CM
|
|---|
| 165 | E3CM=E4CM
|
|---|
| 166 | E4CM=HELP
|
|---|
| 167 | HELP=AMASS3
|
|---|
| 168 | AMASS3=AMASS4
|
|---|
| 169 | AMASS4=HELP
|
|---|
| 170 | IHELP=IQ(NP)
|
|---|
| 171 | IQ(NP)=IQ(NP+1)
|
|---|
| 172 | IQ(NP+1)=IHELP
|
|---|
| 173 | END IF
|
|---|
| 174 | C*** P3CM,P4CM,P5CM ARE MOMENTA IN C.M. SYSTEM
|
|---|
| 175 | C*** P3SQ,P4SQ,P5SQ ARE SQUARED MOMENTA IN C.M. SYSTEM
|
|---|
| 176 | P3SQ=E3CM**2-AMASS3**2
|
|---|
| 177 | P3CM=SQRT(MAX(0.,P3SQ))
|
|---|
| 178 | P4SQ=E4CM**2-AMASS4**2
|
|---|
| 179 | P4CM=SQRT(MAX(0.,P4SQ))
|
|---|
| 180 | P5SQ=E5CM**2-AMASS5**2
|
|---|
| 181 | P5CM=SQRT(MAX(0.,P5SQ))
|
|---|
| 182 | COSA=(P5SQ-P3SQ-P4SQ)/(2.*P3CM*P4CM)
|
|---|
| 183 | SINA=-SQRT(MAX(0.,1.-COSA**2))
|
|---|
| 184 | COSB=(P4SQ-P3SQ-P5SQ)/(2.*P3CM*P5CM)
|
|---|
| 185 | SINB= SQRT(MAX(0.,1.-COSB**2))
|
|---|
| 186 | C*** NOW SELECT THE THREE INDEPENDENT ANGLES IN C.M. SYSTEM
|
|---|
| 187 | PT3=1.D3*PTRANS(DUMMY)
|
|---|
| 188 | SIN3CM=MIN(1.,PT3/P3CM)
|
|---|
| 189 | COS3CM=SQRT(1.-SIN3CM**2)
|
|---|
| 190 | CALL RMMAR(RNNO86,1,2)
|
|---|
| 191 | PSI=TWOPI*RNNO86
|
|---|
| 192 | LPSI=SINC1*PSI+SINC0
|
|---|
| 193 | SINPSI=SIN1(LPSI)*PSI+SIN0(LPSI)
|
|---|
| 194 | CPSI=PI5D2-PSI
|
|---|
| 195 | LCPSI=SINC1*CPSI+SINC0
|
|---|
| 196 | COSPSI=SIN1(LCPSI)*CPSI+SIN0(LCPSI)
|
|---|
| 197 | C*** THIRD INDEPENDENT ANGLE PHI IS CHOOSEN LATER IN SUBROUTINE UPHI
|
|---|
| 198 | C*** NOW MAKE LORENTZ-TRANSFORMATION FOR PARTICLE 3 (PION)
|
|---|
| 199 | E(NP)=GAMMA*(E3CM+BETA*P3CM*COS3CM)
|
|---|
| 200 | C*** COSTHE AND SINTHE ARE ANGLES IN LAB SYSTEM FOR PARTICLE 3 (PION)
|
|---|
| 201 | COSTHE= MIN((BETA*E3CM+P3CM*COS3CM)*GAMMA/ SQRT(MAX(0.D0,E(NP)**2
|
|---|
| 202 | *-AMASS3**2)),1.D0)
|
|---|
| 203 | SINTHE=SQRT(MAX(0.0,1.-COSTHE**2))
|
|---|
| 204 | C*** SINPHI AND COSPHI ARE NOW SET IN SUBROUTINE UPHI
|
|---|
| 205 | CALL UPHI(2,1)
|
|---|
| 206 | SINFI3=SINPHI
|
|---|
| 207 | COSFI3=COSPHI
|
|---|
| 208 | C*** NOW MAKE LORENTZ-TRANSFORMATION FOR PARTICLE 4 = PION
|
|---|
| 209 | COS4CM=COS3CM*COSA-SIN3CM*COSPSI*SINA
|
|---|
| 210 | NP=NP+1
|
|---|
| 211 | E(NP)=GAMMA*(E4CM+BETA*P4CM*COS4CM)
|
|---|
| 212 | SINT4=SQRT(MAX(0.,1.-COS4CM**2))
|
|---|
| 213 | IF (SINT4.NE.0.) THEN
|
|---|
| 214 | SINT4I =1./SINT4
|
|---|
| 215 | AUXA=COS3CM*COSPSI*SINA+SIN3CM*COSA
|
|---|
| 216 | C *** COSPHI AND SINPHI ARE IN LAB SYSTEM FOR PARTICLE 4 (PION)
|
|---|
| 217 | COSPHI=(COSFI3*AUXA-SINFI3*SINPSI*SINA)*SINT4I
|
|---|
| 218 | SINPHI=(SINFI3*AUXA+COSFI3*SINPSI*SINA)*SINT4I
|
|---|
| 219 | ELSE
|
|---|
| 220 | COSPHI=0.
|
|---|
| 221 | SINPHI=1.
|
|---|
| 222 | END IF
|
|---|
| 223 | C*** COSTHE AND SINTHE ARE IN LAB SYSTEM FOR PARTICLE 4 (PION)
|
|---|
| 224 | COSTHE= MIN((BETA*E4CM+P4CM*COS4CM)*GAMMA/ SQRT(MAX(0.D0,E(NP)**2
|
|---|
| 225 | *-AMASS4**2)),1.D0)
|
|---|
| 226 | SINTHE=SQRT(MAX(0.0,1.-COSTHE**2))
|
|---|
| 227 | CALL UPHI(3,2)
|
|---|
| 228 | C*** NOW MAKE LORENTZ-TRANSFORMATION FOR PARTICLE 5 = RECOIL NUCLEON
|
|---|
| 229 | COS5CM=COS3CM*COSB-SIN3CM*COSPSI*SINB
|
|---|
| 230 | ENUCL=GAMMA*(E5CM+BETA*P5CM*COS5CM)
|
|---|
| 231 | IF ((ENUCL-AMASS5).GT.ELCUT(1)*1000.D0) THEN
|
|---|
| 232 | C *** RECOIL NUCLEON IS ABOVE THRESHOLD AND MUST BE TREATED
|
|---|
| 233 | NP=NP+1
|
|---|
| 234 | E(NP)=ENUCL
|
|---|
| 235 | SINT5=SQRT(MAX(0.,1.-COS5CM**2))
|
|---|
| 236 | IF (SINT5.NE.0.) THEN
|
|---|
| 237 | SINT5I =1./SINT5
|
|---|
| 238 | AUXB=COS3CM*COSPSI*SINB+SIN3CM*COSB
|
|---|
| 239 | C *** COSPHI AND SINPHI ARE IN LAB SYSTEM FOR PART. 5 (NUCLEON)
|
|---|
| 240 | COSPHI=(COSFI3*AUXB-SINFI3*SINPSI*SINB)*SINT5I
|
|---|
| 241 | SINPHI=(SINFI3*AUXB+COSFI3*SINPSI*SINB)*SINT5I
|
|---|
| 242 | ELSE
|
|---|
| 243 | COSPHI=0.
|
|---|
| 244 | SINPHI=1.
|
|---|
| 245 | END IF
|
|---|
| 246 | C *** COSTHE AND SINTHE ARE IN LAB SYSTEM FOR PARTICLE 5 (NUCLEON)
|
|---|
| 247 | COSTHE=MIN((BETA*E5CM+P5CM*COS5CM)*GAMMA/SQRT(ENUCL**2-AMASS5**2)
|
|---|
| 248 | * , 1.D0)
|
|---|
| 249 | SINTHE=SQRT(MAX(0.0,1.-COSTHE**2))
|
|---|
| 250 | CALL UPHI(3,2)
|
|---|
| 251 | IF (W(NP).GT.C(29)) THEN
|
|---|
| 252 | C *** ANGLE WITH RESPECT TO X AXIS
|
|---|
| 253 | IF (U(NP)**2+V(NP)**2.GT.3.E-38) THEN
|
|---|
| 254 | ANGLEX = -ATAN2(V(NP),U(NP))
|
|---|
| 255 | ELSE
|
|---|
| 256 | ANGLEX = 0.
|
|---|
| 257 | END IF
|
|---|
| 258 | C *** ADD NUCLEON TO CORSIKA STACK
|
|---|
| 259 | SECPAR(1)=IQ(NP)
|
|---|
| 260 | SECPAR(2)=E(NP)/AMASS5
|
|---|
| 261 | SECPAR(3)=W(NP)
|
|---|
| 262 | SECPAR(4)=ANGLEX
|
|---|
| 263 | SECPAR(5)=-Z(NP)
|
|---|
| 264 | SECPAR(6)=TIME(NP)
|
|---|
| 265 | SECPAR(7)=X(NP)
|
|---|
| 266 | SECPAR(8)=-Y(NP)
|
|---|
| 267 | SECPAR(11)=1.D0
|
|---|
| 268 | SECPAR(12)=0.D0
|
|---|
| 269 | CALL TSTOUT
|
|---|
| 270 | END IF
|
|---|
| 271 | C *** ELIMINATE NUCLEON FROM EGS-STACK
|
|---|
| 272 | NP=NP-1
|
|---|
| 273 | END IF
|
|---|
| 274 | C*** END OF RECOIL NUCLEON TREATEMENT CASE
|
|---|
| 275 | RETURN
|
|---|
| 276 | END
|
|---|
