* * $Id: gfluct.F,v 3.1.1.1 2002/07/11 16:02:14 cafagna Exp $ * * $Log: gfluct.F,v $ * Revision 3.1.1.1 2002/07/11 16:02:14 cafagna * First GPAMELA release on CVS * * #if !defined(GPAMELA_NOGFLUCT) *CMZ : 3.00/00 22/12/2000 19.00.22 by Marialuigia Ambriola *CMZ : 2.01/00 06/03/2000 13.07.03 by Francesco Cafagna *CMZ : 2.00/00 03/03/2000 15.39.06 by Francesco Cafagna *CMZ : 1.02/00 18/03/97 18.25.28 by Francesco Cafagna *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(GPAMELA_ASHO) C. * STRA=2 <---PAI----><---ASHO---><----PAI----><--GLANDZ--> * C. * * #endif C. * SSTR=1: Silicon straggling, gaussian straggling * 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 "gcbank.inc" #include "gcjloc.inc" #include "gconsp.inc" #include "gcmate.inc" #include "gccuts.inc" #include "gckine.inc" #include "gcmulo.inc" #include "gcphys.inc" #include "gctrak.inc" * * * *** VARIABLE FOR SILICON STRAGGLING #include "gpkey.inc" #include "gppmat.inc" INTEGER NMAT,NWNMAT CHARACTER*20 CHNMAT, REAL ANMAT,ZNMAT,DNMAT,RNMAT,ABNMAT,UBNMAT(NWMMAT),GRNDM0,GPGAUS * 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(GPAMELA_ASHO) PARAMETER (ASHMIN=1,ASHMAX=30) #endif DIMENSION RNDM(2) 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.NE.2) THEN CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,Q(JPROB+10)) 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(GPAMELA_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 - (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 - 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(GPAMELA_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(GPAMELA_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 * * *** Add straggling specifically for Silicon By P.Papini and F.S.Cafagna * * change in de/dx with gaussian convolution * see Hall paper: NIM 220 (1984) 356. * IF(ISSTR.EQ.1) THEN * * *** First of all get a look to the DELTA2 and SIGMA2 parameter CALL GFMATE(NMAT,CHNMAT,ANMAT,ZNMAT,DNMAT,RNMAT,ABNMAT, + UBNMAT,NWNMAT) * * *** Save user infos DELTA2 = UBNMAT(1) SIGMA2 = DELTA2*STEP GRNDM0 = GPGAUS(0.) * Gaussian energy spread... * Delta is in GeV^2/cm DE = DE + GRNDM0 * SQRT(SIGMA2) ENDIF ENDIF * END