[u/mrichter/AliRoot.git] / GEANT321 / gphys / gfluct.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1999/05/18 15:55:20  fca
* AliRoot sources
*
* Revision 1.1.1.1  1995/10/24 10:21:24  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.21  by  S.Giani
*-- Author :
      SUBROUTINE GFLUCT(DEMEAN,DE)
C.
C.    ******************************************************************
C.    *                                                                *
C.    *   Subroutine to decide which method is used to simulate        *
C.    *   the straggling around the mean energy loss.                  *
C.    *                                                                *
C.    *                                                                *
C.    *   DNMIN:  <---------->1<-------->30<--------->50<--------->    *
C.    *                                                                *
C.    *   LOSS=2  :                                                    *
C.    *   STRA=0  <----------GLANDZ-------------------><--GLANDO-->    *
C.    *   LOSS=1,3:                                                    *
C.    *   STRA=0  <---------------------GLANDZ-------------------->    *
C.    *                                                                *
C.    *   STRA=1  <-----------PAI---------------------><--GLANDZ-->    *
#if defined(CERNLIB_ASHO)
C.    *   STRA=2  <---PAI----><---ASHO---><----PAI----><--GLANDZ-->    *
C.    *                                                                *
#endif
C.    *                                                                *
C.    *   DNMIN :  an estimation of the number of collisions           *
C.    *            with energy close to the ionization energy          *
C.    *            (see PHYS333)                                       *
C.    *                                                                *
C.    *   Input  : DEMEAN (mean energy loss)                           *
C.    *   Output : DE   (energy loss in the current step)              *
C.    *                                                                *
C.    *    ==>Called by : GTELEC,GTMUON,GTHADR                         *
C.    *                                                                *
C.    ******************************************************************
C.
#include "geant321/gcbank.inc"
#include "geant321/gcjloc.inc"
#include "geant321/gconsp.inc"
#include "geant321/gcmate.inc"
#include "geant321/gccuts.inc"
#include "geant321/gckine.inc"
#include "geant321/gcmulo.inc"
#include "geant321/gcphys.inc"
#include "geant321/gctrak.inc"
*
      PARAMETER (EULER=0.577215,GAM1=EULER-1)
      PARAMETER (P1=.60715,P2=.67794,P3=.52382E-1,P4=.94753,
     +           P5=.74442,P6=1.1934)
      PARAMETER (DGEV=0.153536 E-3, DNLIM=50)
#if defined(CERNLIB_ASHO)
      PARAMETER (ASHMIN=1,ASHMAX=30)
#endif
* These parameters are needed by M.Kowalski's fluctuation algorithm
      PARAMETER (FPOT=20.77E-9, EEND=10E-6, EEXPO=2.2)
      PARAMETER (XEXPO=-EEXPO+1, YEXPO=1/XEXPO)
* These parameters are needed by M.Kowalski's fluctuation algorithm
      DIMENSION RNDM(2)
      DE2(DPOT,RAN)=(DPOT**XEXPO*(1-RAN)+EEND**XEXPO*RAN)**YEXPO
      FLAND(X) = P1+P6*X+(P2+P3*X)*EXP(P4*X+P5)
*
      IF(STEP.LE.0) THEN
         DE=DEMEAN
      ELSE
         DEDX = DEMEAN/STEP
         POTI=Q(JPROB+9)
         IF(ISTRA.EQ.0.AND.(ILOSS.EQ.1.OR.ILOSS.EQ.3)) THEN
            CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,Q(JPROB+10))
         ELSEIF (ILOSS.EQ.5) THEN
* This is Marek Kowalski's fluctuation algorithm, it works only when
* the step size has been limited to one ionisation on average
            CALL GRNDM(RNDM,1)
            DE=DE2(FPOT,RNDM(1))
*
         ELSE
*
* *** mean ionization potential (GeV)
*        POTI=16E-9*Z**0.9
*
            GAMMA = GETOT/AMASS
            BETA = VECT(7)/GETOT
            BET2 = BETA**2
*
* ***    low energy transfer
            XI = DGEV*CHARGE**2*STEP*DENS*Z/(A*BET2)
*
* ***    regime
* ***    ISTRA = 1 --> PAI + URBAN
#if defined(CERNLIB_ASHO)
* ***    ISTRA = 2 --> PAI + URBAN + ASHO
#endif
            DNMIN = MIN(XI,DEMEAN)/POTI
*
            IF (ISTRA.EQ.0) THEN
               IF(DNMIN.GE.DNLIM) THEN
*
*  Energy straggling using Gaussian, Landau & Vavilov theories.
*
*  STEP   =  current step-length (cm)
*
*  DELAND =  DE/DX - <DE/DX>     (GeV)
*
*  Author      : G.N. Patrick
*
                  IF(STEP.LT.1.E-7)THEN
                     DELAND=0.
                  ELSE
*
*     Maximum energy transfer to atomic electron (GeV).
                     ETA = BETA*GAMMA
                     RATIO = EMASS/AMASS
*
*     Calculate Kappa significance ratio.
*                 EMAX=(2*EMASS*ETA**2)/(1+2*RATIO*GAMMA+RATIO**2)
*                 CAPPA = XI/EMAX
                     CAPPA = XI*(1+2*RATIO*GAMMA+RATIO**2)/(2*EMASS*
     +               ETA**2)
                     IF (CAPPA.GE.10.) THEN
*
*     +-----------------------------------+
*     I Sample from Gaussian distribution I
*     +-----------------------------------+
                        SIGMA = XI*SQRT((1.-0.5*BET2)/CAPPA)
                        CALL GRNDM(RNDM,2)
                        F1 = -2.*LOG(RNDM(1))
                        DELAND = SIGMA*SQRT(F1)*COS(TWOPI*RNDM(2))
                     ELSE
                        XMEAN = -BET2-LOG(CAPPA)+GAM1
                        IF (CAPPA.LT.0.01) THEN
                           XLAMX = FLAND(XMEAN)
*     +---------------------------------------------------------------+
*     I Sample lambda variable from Kolbig/Schorr Landau distribution I
*     +---------------------------------------------------------------+
*  10                   CALL GRNDM(RNDM,1)
*                       IF( RNDM(1) .GT. 0.980 ) GO TO 10
*                       XLAMB = GLANDR(RNDM(1))
*     +---------------------------------------------------------------+
*     I Sample lambda variable from James/Hancock Landau distribution I
*     +---------------------------------------------------------------+
   10                      CALL GLANDG(XLAMB)
                           IF(XLAMB.GT.XLAMX) GO TO 10
                        ELSE
*            +---------------------------------------------------------+
*            I Sample lambda variable (Landau not Vavilov) from        I
*            I Rotondi&Montagna&Kolbig Vavilov distribution            I
*            +---------------------------------------------------------+
                           CALL GRNDM(RNDM,1)
                           XLAMB = GVAVIV(CAPPA,BET2,RNDM(1))
                        ENDIF
*
*     Calculate DE/DX - <DE/DX>
                        DELAND = XI*(XLAMB-XMEAN)
                     ENDIF
                  ENDIF
                  DE = DEMEAN + DELAND
               ELSE
                  CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,
     +            Q(JPROB+ 10))
               ENDIF
            ELSE IF (ISTRA.LE.2) THEN
               IF(DNMIN.GE.DNLIM) THEN
                  CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,
     +            Q(JPROB+ 10))
               ELSE
                  NMEC = NMEC+1
                  LMEC(NMEC)=109
#if defined(CERNLIB_ASHO)
                  IF (DNMIN.GE.ASHMIN.AND.DNMIN.LT.ASHMAX .AND.ISTRA.EQ
     +            .2) THEN
                     CALL GASHO(VECT(7),AMASS,STEP,DE)
                  ELSE
                     DE = GSTREN(GAMMA,DCUTE,STEP)
                  ENDIF
#endif
#if !defined(CERNLIB_ASHO)
                  DE = GSTREN(GAMMA,DCUTE,STEP)
#endif
*
* ***   Add brem losses to ionisation
                  IF(ITRTYP.EQ.2) THEN
                     JBASE = LQ(JMA-1)+2*NEK1+IEKBIN
                     DE = DE +(1.-GEKRAT)*Q(JBASE)+GEKRAT*Q(JBASE+1)
                  ELSEIF(ITRTYP.EQ.5) THEN
                     JBASE = LQ(JMA-2)+NEK1+IEKBIN
                     DE = DE +(1.-GEKRAT)*Q(JBASE)+GEKRAT*Q(JBASE+1)
                  ENDIF
               ENDIF
            ENDIF
         ENDIF
      ENDIF
*
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
24b41ee6	5	* Revision 1.1.1.1 1999/05/18 15:55:20 fca
	6	* AliRoot sources
	7	*
fe4da5cc	8	* Revision 1.1.1.1 1995/10/24 10:21:24 cernlib
	9	* Geant
	10	*
	11	*
	12	#include "geant321/pilot.h"
	13	*CMZ : 3.21/02 29/03/94 15.41.21 by S.Giani
	14	*-- Author :
	15	SUBROUTINE GFLUCT(DEMEAN,DE)
	16	C.
	17	C. ******************************************************************
	18	C. * *
	19	C. * Subroutine to decide which method is used to simulate *
	20	C. * the straggling around the mean energy loss. *
	21	C. * *
	22	C. * *
	23	C. * DNMIN: <---------->1<-------->30<--------->50<---------> *
	24	C. * *
	25	C. * LOSS=2 : *
	26	C. * STRA=0 <----------GLANDZ-------------------><--GLANDO--> *
	27	C. * LOSS=1,3: *
	28	C. * STRA=0 <---------------------GLANDZ--------------------> *
	29	C. * *
	30	C. * STRA=1 <-----------PAI---------------------><--GLANDZ--> *
	31	#if defined(CERNLIB_ASHO)
	32	C. * STRA=2 <---PAI----><---ASHO---><----PAI----><--GLANDZ--> *
	33	C. * *
	34	#endif
	35	C. * *
	36	C. * DNMIN : an estimation of the number of collisions *
	37	C. * with energy close to the ionization energy *
	38	C. * (see PHYS333) *
	39	C. * *
	40	C. * Input : DEMEAN (mean energy loss) *
	41	C. * Output : DE (energy loss in the current step) *
	42	C. * *
	43	C. * ==>Called by : GTELEC,GTMUON,GTHADR *
	44	C. * *
	45	C. ******************************************************************
	46	C.
	47	#include "geant321/gcbank.inc"
	48	#include "geant321/gcjloc.inc"
	49	#include "geant321/gconsp.inc"
	50	#include "geant321/gcmate.inc"
	51	#include "geant321/gccuts.inc"
	52	#include "geant321/gckine.inc"
	53	#include "geant321/gcmulo.inc"
	54	#include "geant321/gcphys.inc"
	55	#include "geant321/gctrak.inc"
	56	*
	57	PARAMETER (EULER=0.577215,GAM1=EULER-1)
	58	PARAMETER (P1=.60715,P2=.67794,P3=.52382E-1,P4=.94753,
	59	+ P5=.74442,P6=1.1934)
	60	PARAMETER (DGEV=0.153536 E-3, DNLIM=50)
	61	#if defined(CERNLIB_ASHO)
	62	PARAMETER (ASHMIN=1,ASHMAX=30)
	63	#endif
24b41ee6	64	* These parameters are needed by M.Kowalski's fluctuation algorithm
	65	PARAMETER (FPOT=20.77E-9, EEND=10E-6, EEXPO=2.2)
	66	PARAMETER (XEXPO=-EEXPO+1, YEXPO=1/XEXPO)
	67	* These parameters are needed by M.Kowalski's fluctuation algorithm
fe4da5cc	68	DIMENSION RNDM(2)
24b41ee6	69	DE2(DPOT,RAN)=(DPOT*XEXPO(1-RAN)+EEND*XEXPORAN)**YEXPO
fe4da5cc	70	FLAND(X) = P1+P6X+(P2+P3X)EXP(P4X+P5)
	71	*
	72	IF(STEP.LE.0) THEN
	73	DE=DEMEAN
	74	ELSE
	75	DEDX = DEMEAN/STEP
	76	POTI=Q(JPROB+9)
24b41ee6	77	IF(ISTRA.EQ.0.AND.(ILOSS.EQ.1.OR.ILOSS.EQ.3)) THEN
fe4da5cc	78	CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,Q(JPROB+10))
24b41ee6	79	ELSEIF (ILOSS.EQ.5) THEN
	80	* This is Marek Kowalski's fluctuation algorithm, it works only when
	81	* the step size has been limited to one ionisation on average
	82	CALL GRNDM(RNDM,1)
	83	DE=DE2(FPOT,RNDM(1))
	84	*
fe4da5cc	85	ELSE
	86	*
	87	* *** mean ionization potential (GeV)
	88	* POTI=16E-9Z*0.9
	89	*
	90	GAMMA = GETOT/AMASS
	91	BETA = VECT(7)/GETOT
	92	BET2 = BETA**2
	93	*
	94	* *** low energy transfer
	95	XI = DGEVCHARGE2STEPDENSZ/(A*BET2)
	96	*
	97	* *** regime
	98	* *** ISTRA = 1 --> PAI + URBAN
	99	#if defined(CERNLIB_ASHO)
	100	* *** ISTRA = 2 --> PAI + URBAN + ASHO
	101	#endif
	102	DNMIN = MIN(XI,DEMEAN)/POTI
	103	*
	104	IF (ISTRA.EQ.0) THEN
	105	IF(DNMIN.GE.DNLIM) THEN
	106	*
	107	* Energy straggling using Gaussian, Landau & Vavilov theories.
	108	*
	109	* STEP = current step-length (cm)
	110	*
	111	* DELAND = DE/DX - <DE/DX> (GeV)
	112	*
	113	* Author : G.N. Patrick
	114	*
	115	IF(STEP.LT.1.E-7)THEN
	116	DELAND=0.
	117	ELSE
	118	*
	119	* Maximum energy transfer to atomic electron (GeV).
	120	ETA = BETA*GAMMA
	121	RATIO = EMASS/AMASS
	122	*
	123	* Calculate Kappa significance ratio.
	124	* EMAX=(2EMASSETA*2)/(1+2RATIOGAMMA+RATIO*2)
	125	* CAPPA = XI/EMAX
	126	CAPPA = XI(1+2RATIOGAMMA+RATIO2)/(2EMASS*
	127	+ ETA**2)
	128	IF (CAPPA.GE.10.) THEN
	129	*
	130	* +-----------------------------------+
	131	* I Sample from Gaussian distribution I
	132	* +-----------------------------------+
	133	SIGMA = XISQRT((1.-0.5BET2)/CAPPA)
	134	CALL GRNDM(RNDM,2)
	135	F1 = -2.*LOG(RNDM(1))
	136	DELAND = SIGMASQRT(F1)COS(TWOPI*RNDM(2))
	137	ELSE
	138	XMEAN = -BET2-LOG(CAPPA)+GAM1
	139	IF (CAPPA.LT.0.01) THEN
	140	XLAMX = FLAND(XMEAN)
	141	* +---------------------------------------------------------------+
	142	* I Sample lambda variable from Kolbig/Schorr Landau distribution I
	143	* +---------------------------------------------------------------+
	144	* 10 CALL GRNDM(RNDM,1)
	145	* IF( RNDM(1) .GT. 0.980 ) GO TO 10
	146	* XLAMB = GLANDR(RNDM(1))
	147	* +---------------------------------------------------------------+
	148	* I Sample lambda variable from James/Hancock Landau distribution I
149	* +---------------------------------------------------------------+
150	10 CALL GLANDG(XLAMB)
151	IF(XLAMB.GT.XLAMX) GO TO 10
152	ELSE
153	* +---------------------------------------------------------+
154	* I Sample lambda variable (Landau not Vavilov) from I
155	* I Rotondi&Montagna&Kolbig Vavilov distribution I
156	* +---------------------------------------------------------+
157	CALL GRNDM(RNDM,1)
158	XLAMB = GVAVIV(CAPPA,BET2,RNDM(1))
159	ENDIF
160	*
161	* Calculate DE/DX - <DE/DX>
162	DELAND = XI*(XLAMB-XMEAN)
163	ENDIF
164	ENDIF
165	DE = DEMEAN + DELAND
166	ELSE
167	CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,
168	+ Q(JPROB+ 10))
169	ENDIF
170	ELSE IF (ISTRA.LE.2) THEN
171	IF(DNMIN.GE.DNLIM) THEN
172	CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,
173	+ Q(JPROB+ 10))
174	ELSE
175	NMEC = NMEC+1
176	LMEC(NMEC)=109
177	#if defined(CERNLIB_ASHO)
178	IF (DNMIN.GE.ASHMIN.AND.DNMIN.LT.ASHMAX .AND.ISTRA.EQ
179	+ .2) THEN
180	CALL GASHO(VECT(7),AMASS,STEP,DE)
181	ELSE
182	DE = GSTREN(GAMMA,DCUTE,STEP)
183	ENDIF
184	#endif
185	#if !defined(CERNLIB_ASHO)
186	DE = GSTREN(GAMMA,DCUTE,STEP)
187	#endif
188	*
189	* *** Add brem losses to ionisation
190	IF(ITRTYP.EQ.2) THEN
191	JBASE = LQ(JMA-1)+2*NEK1+IEKBIN
192	DE = DE +(1.-GEKRAT)Q(JBASE)+GEKRATQ(JBASE+1)
193	ELSEIF(ITRTYP.EQ.5) THEN
194	JBASE = LQ(JMA-2)+NEK1+IEKBIN
195	DE = DE +(1.-GEKRAT)Q(JBASE)+GEKRATQ(JBASE+1)
196	ENDIF
197	ENDIF
198	ENDIF
199	ENDIF
200	ENDIF
201	*
202	END