| 1 |
* |
| 2 |
* $Id$ |
| 3 |
* |
| 4 |
* $Log$ |
| 5 |
* |
| 6 |
#if !defined(GPAMELA_NOGFLUCT) |
| 7 |
*CMZ : 3.00/00 22/12/2000 19.00.22 by Marialuigia Ambriola |
| 8 |
*CMZ : 2.01/00 06/03/2000 13.07.03 by Francesco Cafagna |
| 9 |
*CMZ : 2.00/00 03/03/2000 15.39.06 by Francesco Cafagna |
| 10 |
*CMZ : 1.02/00 18/03/97 18.25.28 by Francesco Cafagna |
| 11 |
*CMZ : 3.21/02 29/03/94 15.41.21 by S.Giani |
| 12 |
*-- Author : |
| 13 |
SUBROUTINE GFLUCT(DEMEAN,DE) |
| 14 |
C. |
| 15 |
C. ****************************************************************** |
| 16 |
C. * * |
| 17 |
C. * Subroutine to decide which method is used to simulate * |
| 18 |
C. * the straggling around the mean energy loss. * |
| 19 |
C. * * |
| 20 |
C. * * |
| 21 |
C. * DNMIN: <---------->1<-------->30<--------->50<---------> * |
| 22 |
C. * * |
| 23 |
C. * LOSS=2 : * |
| 24 |
C. * STRA=0 <----------GLANDZ-------------------><--GLANDO--> * |
| 25 |
C. * LOSS=1,3: * |
| 26 |
C. * STRA=0 <---------------------GLANDZ--------------------> * |
| 27 |
C. * * |
| 28 |
C. * STRA=1 <-----------PAI---------------------><--GLANDZ--> * |
| 29 |
#if defined(GPAMELA_ASHO) |
| 30 |
C. * STRA=2 <---PAI----><---ASHO---><----PAI----><--GLANDZ--> * |
| 31 |
C. * * |
| 32 |
#endif |
| 33 |
C. * SSTR=1: Silicon straggling, gaussian straggling * |
| 34 |
C. * * |
| 35 |
C. * DNMIN : an estimation of the number of collisions * |
| 36 |
C. * with energy close to the ionization energy * |
| 37 |
C. * (see PHYS333) * |
| 38 |
C. * * |
| 39 |
C. * Input : DEMEAN (mean energy loss) * |
| 40 |
C. * Output : DE (energy loss in the current step) * |
| 41 |
C. * * |
| 42 |
C. * ==>Called by : GTELEC,GTMUON,GTHADR * |
| 43 |
C. * * |
| 44 |
C. ****************************************************************** |
| 45 |
C. |
| 46 |
#include "gcbank.inc" |
| 47 |
#include "gcjloc.inc" |
| 48 |
#include "gconsp.inc" |
| 49 |
#include "gcmate.inc" |
| 50 |
#include "gccuts.inc" |
| 51 |
#include "gckine.inc" |
| 52 |
#include "gcmulo.inc" |
| 53 |
#include "gcphys.inc" |
| 54 |
#include "gctrak.inc" |
| 55 |
* |
| 56 |
* |
| 57 |
* *** VARIABLE FOR SILICON STRAGGLING |
| 58 |
#include "gpkey.inc" |
| 59 |
#include "gppmat.inc" |
| 60 |
INTEGER NMAT,NWNMAT |
| 61 |
CHARACTER*20 CHNMAT, |
| 62 |
REAL ANMAT,ZNMAT,DNMAT,RNMAT,ABNMAT,UBNMAT(NWMMAT),GRNDM0,GPGAUS |
| 63 |
* |
| 64 |
PARAMETER (EULER=0.577215,GAM1=EULER-1) |
| 65 |
PARAMETER (P1=.60715,P2=.67794,P3=.52382E-1,P4=.94753, |
| 66 |
+ P5=.74442,P6=1.1934) |
| 67 |
PARAMETER (DGEV=0.153536 E-3, DNLIM=50) |
| 68 |
#if defined(GPAMELA_ASHO) |
| 69 |
PARAMETER (ASHMIN=1,ASHMAX=30) |
| 70 |
#endif |
| 71 |
DIMENSION RNDM(2) |
| 72 |
FLAND(X) = P1+P6*X+(P2+P3*X)*EXP(P4*X+P5) |
| 73 |
* |
| 74 |
IF(STEP.LE.0) THEN |
| 75 |
DE=DEMEAN |
| 76 |
ELSE |
| 77 |
DEDX = DEMEAN/STEP |
| 78 |
POTI=Q(JPROB+9) |
| 79 |
IF(ISTRA.EQ.0.AND.ILOSS.NE.2) THEN |
| 80 |
CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,Q(JPROB+10)) |
| 81 |
ELSE |
| 82 |
* |
| 83 |
* *** mean ionization potential (GeV) |
| 84 |
* POTI=16E-9*Z**0.9 |
| 85 |
* |
| 86 |
GAMMA = GETOT/AMASS |
| 87 |
BETA = VECT(7)/GETOT |
| 88 |
BET2 = BETA**2 |
| 89 |
* |
| 90 |
* *** low energy transfer |
| 91 |
XI = DGEV*CHARGE**2*STEP*DENS*Z/(A*BET2) |
| 92 |
* |
| 93 |
* *** regime |
| 94 |
* *** ISTRA = 1 --> PAI + URBAN |
| 95 |
#if defined(GPAMELA_ASHO) |
| 96 |
* *** ISTRA = 2 --> PAI + URBAN + ASHO |
| 97 |
#endif |
| 98 |
DNMIN = MIN(XI,DEMEAN)/POTI |
| 99 |
* |
| 100 |
IF (ISTRA.EQ.0) THEN |
| 101 |
IF(DNMIN.GE.DNLIM) THEN |
| 102 |
* |
| 103 |
* Energy straggling using Gaussian, Landau & Vavilov theories. |
| 104 |
* |
| 105 |
* STEP = current step-length (cm) |
| 106 |
* |
| 107 |
* DELAND = DE/DX - <DE/DX> (GeV) |
| 108 |
* |
| 109 |
* Author : G.N. Patrick |
| 110 |
* |
| 111 |
IF(STEP.LT.1.E-7)THEN |
| 112 |
DELAND=0. |
| 113 |
ELSE |
| 114 |
* |
| 115 |
* Maximum energy transfer to atomic electron (GeV). |
| 116 |
ETA = BETA*GAMMA |
| 117 |
RATIO = EMASS/AMASS |
| 118 |
* |
| 119 |
* Calculate Kappa significance ratio. |
| 120 |
* EMAX=(2*EMASS*ETA**2)/(1+2*RATIO*GAMMA+RATIO**2) |
| 121 |
* CAPPA = XI/EMAX |
| 122 |
CAPPA = XI*(1+2*RATIO*GAMMA+RATIO**2)/(2*EMASS* |
| 123 |
+ ETA**2) |
| 124 |
IF (CAPPA.GE.10.) THEN |
| 125 |
* |
| 126 |
* +-----------------------------------+ |
| 127 |
* I Sample from Gaussian distribution I |
| 128 |
* +-----------------------------------+ |
| 129 |
SIGMA = XI*SQRT((1.-0.5*BET2)/CAPPA) |
| 130 |
CALL GRNDM(RNDM,2) |
| 131 |
F1 = -2.*LOG(RNDM(1)) |
| 132 |
DELAND = SIGMA*SQRT(F1)*COS(TWOPI*RNDM(2)) |
| 133 |
ELSE |
| 134 |
XMEAN = -BET2-LOG(CAPPA)+GAM1 |
| 135 |
IF (CAPPA.LT.0.01) THEN |
| 136 |
XLAMX = FLAND(XMEAN) |
| 137 |
* +---------------------------------------------------------------+ |
| 138 |
* I Sample lambda variable from Kolbig/Schorr Landau distribution I |
| 139 |
* +---------------------------------------------------------------+ |
| 140 |
* 10 CALL GRNDM(RNDM,1) |
| 141 |
* IF( RNDM(1) .GT. 0.980 ) GO TO 10 |
| 142 |
* XLAMB = GLANDR(RNDM(1)) |
| 143 |
* +---------------------------------------------------------------+ |
| 144 |
* I Sample lambda variable from James/Hancock Landau distribution I |
| 145 |
* +---------------------------------------------------------------+ |
| 146 |
10 CALL GLANDG(XLAMB) |
| 147 |
IF(XLAMB.GT.XLAMX) GO TO 10 |
| 148 |
ELSE |
| 149 |
* +---------------------------------------------------------+ |
| 150 |
* I Sample lambda variable (Landau not Vavilov) from I |
| 151 |
* I Rotondi&Montagna&Kolbig Vavilov distribution I |
| 152 |
* +---------------------------------------------------------+ |
| 153 |
CALL GRNDM(RNDM,1) |
| 154 |
XLAMB = GVAVIV(CAPPA,BET2,RNDM(1)) |
| 155 |
ENDIF |
| 156 |
* |
| 157 |
* Calculate DE/DX - <DE/DX> |
| 158 |
DELAND = XI*(XLAMB-XMEAN) |
| 159 |
ENDIF |
| 160 |
ENDIF |
| 161 |
DE = DEMEAN + DELAND |
| 162 |
ELSE |
| 163 |
CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI, |
| 164 |
+ Q(JPROB+ 10)) |
| 165 |
ENDIF |
| 166 |
ELSE IF (ISTRA.LE.2) THEN |
| 167 |
IF(DNMIN.GE.DNLIM) THEN |
| 168 |
CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI, |
| 169 |
+ Q(JPROB+ 10)) |
| 170 |
ELSE |
| 171 |
NMEC = NMEC+1 |
| 172 |
LMEC(NMEC)=109 |
| 173 |
#if defined(GPAMELA_ASHO) |
| 174 |
IF (DNMIN.GE.ASHMIN.AND.DNMIN.LT.ASHMAX .AND.ISTRA.EQ |
| 175 |
+ .2) THEN |
| 176 |
CALL GASHO(VECT(7),AMASS,STEP,DE) |
| 177 |
ELSE |
| 178 |
DE = GSTREN(GAMMA,DCUTE,STEP) |
| 179 |
ENDIF |
| 180 |
#endif |
| 181 |
#if !defined(GPAMELA_ASHO) |
| 182 |
DE = GSTREN(GAMMA,DCUTE,STEP) |
| 183 |
#endif |
| 184 |
* |
| 185 |
* *** Add brem losses to ionisation |
| 186 |
IF(ITRTYP.EQ.2) THEN |
| 187 |
JBASE = LQ(JMA-1)+2*NEK1+IEKBIN |
| 188 |
DE = DE +(1.-GEKRAT)*Q(JBASE)+GEKRAT*Q(JBASE+1) |
| 189 |
ELSEIF(ITRTYP.EQ.5) THEN |
| 190 |
JBASE = LQ(JMA-2)+NEK1+IEKBIN |
| 191 |
DE = DE +(1.-GEKRAT)*Q(JBASE)+GEKRAT*Q(JBASE+1) |
| 192 |
ENDIF |
| 193 |
ENDIF |
| 194 |
ENDIF |
| 195 |
ENDIF |
| 196 |
* |
| 197 |
* *** Add straggling specifically for Silicon By P.Papini and F.S.Cafagna |
| 198 |
* |
| 199 |
* change in de/dx with gaussian convolution |
| 200 |
* see Hall paper: NIM 220 (1984) 356. |
| 201 |
* |
| 202 |
IF(ISSTR.EQ.1) THEN |
| 203 |
* |
| 204 |
* *** First of all get a look to the DELTA2 and SIGMA2 parameter |
| 205 |
|
| 206 |
CALL GFMATE(NMAT,CHNMAT,ANMAT,ZNMAT,DNMAT,RNMAT,ABNMAT, |
| 207 |
+ UBNMAT,NWNMAT) |
| 208 |
* |
| 209 |
* *** Save user infos |
| 210 |
DELTA2 = UBNMAT(1) |
| 211 |
SIGMA2 = DELTA2*STEP |
| 212 |
GRNDM0 = GPGAUS(0.) |
| 213 |
* Gaussian energy spread... |
| 214 |
* Delta is in GeV^2/cm |
| 215 |
DE = DE + GRNDM0 * SQRT(SIGMA2) |
| 216 |
ENDIF |
| 217 |
ENDIF |
| 218 |
* |
| 219 |
END |
| 220 |
#endif |
Parent Directory
| 
Revision Log