| fe4da5cc |
1 | * |
| 2 | * $Id$ |
| 3 | * |
| 4 | * $Log$ |
| 5 | * Revision 1.1.1.1 1995/10/24 10:21:25 cernlib |
| 6 | * Geant |
| 7 | * |
| 8 | * |
| 9 | #include "geant321/pilot.h" |
| 10 | *CMZ : 3.21/04 22/02/95 10.38.36 by S.Giani |
| 11 | *-- Author : |
| 12 | SUBROUTINE GLANDZ(Z,STEP,P,E,DEDX,DE,POTI,POTIL) |
| 13 | C. |
| 14 | C. ****************************************************************** |
| 15 | C. * * |
| 16 | C. * Energy straggling using a Monte Carlo model. * |
| 17 | C. * It can be used with or without delta ray generation. * |
| 18 | C. * * |
| 19 | C. * It is a NEW VERSION of the model , which reproduces * |
| 20 | C. * the experimental data rather well. * |
| 21 | C. * * |
| 22 | C. * Input : STEP = current step-length (cm) * |
| 23 | C. * * |
| 24 | C. * Output: DE = actual energy loss (Gev) * |
| 25 | C. * ( NOT the fluctuation DE/DX-<DE/DX> !) * |
| 26 | C. * * |
| 27 | C. * ==> Called by : GTELEC,GTHADR,GTMUON * |
| 28 | C. * * |
| 29 | C. * Author : L.Urban * |
| 30 | C. * Date : 28.04.1988 Last update : 1.02.90 * |
| 31 | C. * * |
| 32 | C. ****************************************************************** |
| 33 | C. |
| 34 | #include "geant321/gconsp.inc" |
| 35 | #include "geant321/gcflag.inc" |
| 36 | #include "geant321/gckine.inc" |
| 37 | #include "geant321/gccuts.inc" |
| 38 | PARAMETER (MAXRND=100) |
| 39 | DIMENSION RNDM(MAXRND),APOIS(3),NPOIS(3),FPOIS(3) |
| 40 | #if !defined(CERNLIB_SINGLE) |
| 41 | DOUBLE PRECISION E1,E2,P2,B2,BG2,TM,DEC,W,WW,WWW,RR |
| 42 | DOUBLE PRECISION X,AK,ALFA,EA,SA,AKL,REALK,FPOIS |
| 43 | DOUBLE PRECISION ONE,HALF,ZERO |
| 44 | PARAMETER (HMXINT=2**30) |
| 45 | #endif |
| 46 | #if defined(CERNLIB_SINGLE) |
| 47 | PARAMETER (HMXINT=2**45) |
| 48 | #endif |
| 49 | PARAMETER (ONE=1,HALF=ONE/2,ZERO=0) |
| 50 | PARAMETER (RCD=0.40, RCD1=1.-RCD, PROBLM=0.01) |
| 51 | PARAMETER( C1=4, C2=16) |
| 52 | * |
| 53 | * ------------------------------------------------------------------ |
| 54 | * |
| 55 | * POTI=1.6E-8*Z**0.9 |
| 56 | * POTIL=LOG(POTI) |
| 57 | * |
| 58 | IF(Z.GT.2.) THEN |
| 59 | F2=2./Z |
| 60 | ELSE |
| 61 | F2=0. |
| 62 | ENDIF |
| 63 | F1=1.-F2 |
| 64 | * |
| 65 | E2 = 1.E-8*Z*Z |
| 66 | E2L= LOG(E2) |
| 67 | E1L= (POTIL-F2*E2L)/F1 |
| 68 | E1 = EXP(E1L) |
| 69 | * |
| 70 | * |
| 71 | P2=P*P |
| 72 | B2=P2/(E*E) |
| 73 | BG2=P2/(AMASS*AMASS) |
| 74 | IF(ITRTYP.EQ.2) THEN |
| 75 | TM=P2/(E+EMASS) |
| 76 | IF(CHARGE.LT.0.) TM=TM/2. |
| 77 | TM=TM-POTI |
| 78 | IF(TM.GT.DCUTE) TM=DCUTE |
| 79 | ELSE |
| 80 | TM=EMASS*P2/(0.5*AMASS*AMASS+EMASS*E) |
| 81 | TM=TM-POTI |
| 82 | IF(TM.GT.DCUTM) TM=DCUTM |
| 83 | ENDIF |
| 84 | * |
| 85 | * |
| 86 | * *** Protection against negative TM --------------------- |
| 87 | * TM can be negative only for heavy particles with a very |
| 88 | * low kinetic energy (e.g. for proton with T 100-300 kev) |
| 89 | TM=MAX(TM,ZERO) |
| 90 | * |
| 91 | W = TM+POTI |
| 92 | WW = W/POTI |
| 93 | WWW= 2.*EMASS*BG2 |
| 94 | WL = LOG(WWW) |
| 95 | CSB=STEP*RCD1*DEDX/(WL-POTIL-B2) |
| 96 | APOIS(1)=CSB*F1*(WL-E1L-B2)/E1 |
| 97 | APOIS(2)=CSB*F2*(WL-E2L-B2)/E2 |
| 98 | * |
| 99 | IF(TM.GT.0.) THEN |
| 100 | APOIS(3)=RCD*DEDX*STEP*TM/(POTI*W*LOG(WW)) |
| 101 | ELSE |
| 102 | APOIS(1)=APOIS(1)/RCD1 |
| 103 | APOIS(2)=APOIS(2)/RCD1 |
| 104 | APOIS(3)=0. |
| 105 | ENDIF |
| 106 | * |
| 107 | * calculate the probability of the zero energy loss |
| 108 | * |
| 109 | APSUM=APOIS(1)+APOIS(2)+APOIS(3) |
| 110 | IF(APSUM.LT.0.) THEN |
| 111 | PROB=1. |
| 112 | ELSEIF(APSUM.LT.50.) THEN |
| 113 | PROB=EXP(-APSUM) |
| 114 | ELSE |
| 115 | PROB=0. |
| 116 | ENDIF |
| 117 | * |
| 118 | * |
| 119 | * do it differently if prob > problm <==================== |
| 120 | IF(PROB.GT.PROBLM) THEN |
| 121 | E0=1.E-8 |
| 122 | EMEAN=DEDX*STEP |
| 123 | IF(TM.LE.0.) THEN |
| 124 | * excitation only .... |
| 125 | APOIS(1)=EMEAN/E0 |
| 126 | * |
| 127 | CALL GPOISS(APOIS,NPOIS,1) |
| 128 | FPOIS(1)=NPOIS(1) |
| 129 | DE=FPOIS(1)*E0 |
| 130 | * |
| 131 | ELSE |
| 132 | * ionization only .... |
| 133 | EM=TM+E0 |
| 134 | APOIS(1)=EMEAN*(EM-E0)/(EM*E0*LOG(EM/E0)) |
| 135 | CALL GPOISS(APOIS,NPOIS,1) |
| 136 | NN=NPOIS(1) |
| 137 | DE=0. |
| 138 | * |
| 139 | IF(NN.GT.0) THEN |
| 140 | RCORR=1. |
| 141 | IF(NN.GT.MAXRND) THEN |
| 142 | RCORR=FLOAT(NN)/MAXRND |
| 143 | NN=MAXRND |
| 144 | * |
| 145 | ENDIF |
| 146 | W=(EM-E0)/EM |
| 147 | CALL GRNDM(RNDM,NN) |
| 148 | DO 10 I=1,NN |
| 149 | DE=DE+E0/(1.-W*RNDM(I)) |
| 150 | 10 CONTINUE |
| 151 | DE=RCORR*DE |
| 152 | * |
| 153 | ENDIF |
| 154 | ENDIF |
| 155 | GOTO 999 |
| 156 | ENDIF |
| 157 | * |
| 158 | IF(MAX(APOIS(1),APOIS(2),APOIS(3)).LT.HMXINT) THEN |
| 159 | CALL GPOISS(APOIS,NPOIS,3) |
| 160 | FPOIS(1)=NPOIS(1) |
| 161 | FPOIS(2)=NPOIS(2) |
| 162 | FPOIS(3)=NPOIS(3) |
| 163 | ELSE |
| 164 | DO 20 JPOIS=1, 3 |
| 165 | IF(APOIS(JPOIS).LT.HMXINT) THEN |
| 166 | CALL GPOISS(APOIS(JPOIS),NPOIS(JPOIS),1) |
| 167 | FPOIS(JPOIS)=NPOIS(JPOIS) |
| 168 | ELSE |
| 169 | CALL GRNDM(RNDM,2) |
| 170 | FPOIS(JPOIS)=ABS(SQRT(-2.*LOG(RNDM(1)*ONE) |
| 171 | + *APOIS(JPOIS))*SIN(TWOPI*RNDM(2)*ONE)+APOIS(JPOIS)) |
| 172 | ENDIF |
| 173 | 20 CONTINUE |
| 174 | ENDIF |
| 175 | |
| 176 | * |
| 177 | * Now we have all three numbers in REAL/DOUBLE |
| 178 | * variables. REALK is actually an INTEGER that now may |
| 179 | * exceed the machine representation limit for integers. |
| 180 | * |
| 181 | DE=FPOIS(1)*E1+FPOIS(2)*E2 |
| 182 | * |
| 183 | * smear to avoid peaks in the energy loss (note: E1<<E2) |
| 184 | * |
| 185 | IF(DE.GT.0.) THEN |
| 186 | CALL GRNDM(RNDM,1) |
| 187 | DE=DE+E1*(1.-2.*RNDM(1)) |
| 188 | ENDIF |
| 189 | * |
| 190 | ALFA=1. |
| 191 | REALK=0 |
| 192 | DEC=0. |
| 193 | * |
| 194 | ANC=FPOIS(3) |
| 195 | IF(ANC.GE.C2) THEN |
| 196 | R=ANC/(C2+ANC) |
| 197 | AN=ANC*R |
| 198 | SN=C1*R |
| 199 | CALL GRNDM(RNDM,2) |
| 200 | RR=SQRT(-2.*LOG(RNDM(1))) |
| 201 | PHI=TWOPI*RNDM(2) |
| 202 | X=RR*COS(PHI) |
| 203 | AK=AN+SN*X |
| 204 | ALFA=WW*(C2+ANC)/(C2*WW+ANC) |
| 205 | EA=AK*POTI*ALFA*LOG(ALFA)/(ALFA-1.) |
| 206 | SA=SQRT(ABS(AK*ALFA*POTI*POTI-EA*EA/AK)) |
| 207 | AKL=(EA-C1*SA)/POTI |
| 208 | IF(AK.LE.AKL) THEN |
| 209 | X=RR*SIN(PHI) |
| 210 | DEC=EA+SA*X |
| 211 | REALK=AK+HALF-MOD(AK+HALF,ONE) |
| 212 | ELSE |
| 213 | ALFA=1. |
| 214 | ENDIF |
| 215 | ENDIF |
| 216 | NN=NINT(FPOIS(3)-REALK) |
| 217 | IF(NN.GT.MAXRND) THEN |
| 218 | W=1.-ALFA/WW |
| 219 | WW=POTI*ALFA |
| 220 | * |
| 221 | * Here we take a gaussian distribution to avoid loosing |
| 222 | * too much time in computing |
| 223 | * |
| 224 | AVERAG=-LOG(1.-W)/W |
| 225 | SIGMA =SQRT(NN*(1./(1.-W)-AVERAG**2)) |
| 226 | CALL GRNDM(RNDM,2) |
| 227 | DEC = DEC+WW*(NN*AVERAG+SIGMA*SQRT(-2.*LOG(RNDM(1)))* |
| 228 | + SIN(TWOPI*RNDM(2))) |
| 229 | DE=DE+DEC |
| 230 | ELSEIF(NN.GT.0) THEN |
| 231 | W=1.-ALFA/WW |
| 232 | WW=POTI*ALFA |
| 233 | CALL GRNDM(RNDM,NN) |
| 234 | DO 30 I=1,NN |
| 235 | DEC=DEC+WW/(1.-W*RNDM(I)) |
| 236 | 30 CONTINUE |
| 237 | DE=DE+DEC |
| 238 | ENDIF |
| 239 | * |
| 240 | 999 END |
git://git.uio.no / u / mrichter / AliRoot.git / blame