/[PAMELA software]/gpamela/gphys/gfluct.F
ViewVC logotype

Annotation of /gpamela/gphys/gfluct.F

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3.1.1.1 - (hide annotations) (download) (vendor branch)
Thu Jul 11 16:02:14 2002 UTC (22 years, 4 months ago) by cafagna
Branch: v3r0
CVS Tags: v4r0, v4r1, v4r2, v4r3, firstrelease, v3r3, v3r1, v3r2
Changes since 3.1: +0 -0 lines
First GPAMELA release on CVS

1 cafagna 3.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

  ViewVC Help
Powered by ViewVC 1.1.23