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

*
* $Id$
*
* $Log$
* Revision 1.2  1997/11/13 09:25:40  gunter
* Correction By Laszlo Urban; protect against divide by 0. if AL becomes 0.
* around GEKIN 2.93E-5.
* Reported by harald@psiclu.cern.ch (signed by keller@biomed.ee.ethz.ch)
*
* Revision 1.1.1.1  1995/10/24 10:21:22  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/03 02/08/94  19.37.36  by  S.Ravndal
*-- Author :
      SUBROUTINE GBREME
C.
C.    ******************************************************************
C.    *                                                                *
C.    *  Simulates discrete hard BREMSSTRAHLUNG by electrons.          *
C.    *                                                                *
C.    *  The secondary photon energy is sampled from a                 *
C.    *  parametrization of the bremsstrahlung calculation of          *
C.    *  Seltzer and Berger (NIM B12,p.95(1985)) for electron          *
C.    *  energies 1 keV - 10 GeV . For higher energies the             *
C.    *  parametrization agrees with the screened Bethe-Heitler        *
C.    *  bremsstrahlung spectrum.Migdal corrections are applied        *
C.    *      by default. The routine works ( together                  *
C.    *  with the routines GBRELE and GBRSGE ) without the Migdal      *
C.    *  corrections using the Patchy switch +USE,BETHE.               *
C.    *                                                                *
C.    *  NOTE :                                                        *
C.    *      BCUTE is the cut-off energy above which the photon energy *
C.    *      spectrum is generated.                                    *
C.    *                                                                *
C.    *    ==>Called by : GTELEC                                       *
C.    *       Authors    L.Urban    29/10/93 *********                 *
C.    *                                                                *
C.    *  13/11/97 bug corrected by L.Urban                             *
C.    ******************************************************************
C.
#include "geant321/gcbank.inc"
#include "geant321/gcphys.inc"
#include "geant321/gcjloc.inc"
#include "geant321/gconsp.inc"
#include "geant321/gctrak.inc"
#include "geant321/gcking.inc"
#include "geant321/gcmate.inc"
#include "geant321/gccuts.inc"
      DIMENSION RNDM(2)
      PARAMETER (CORFAC=0.805485E-10)
      PARAMETER (TLIM=0.000999)
      PARAMETER  (
     +AL00=-0.205398E+01,AL01= 0.238815E-01,AL02= 0.525483E-03,
     +AL10=-0.769748E-01,AL11=-0.691499E-01,AL12= 0.222453E-02,
     +AL20= 0.406463E-01,AL21=-0.101281E-01,AL22= 0.340919E-03,
     +A10 = 0.467733E+01,A11 =-0.619012E+00,A12 = 0.202225E-01,
     +A20 =-0.734101E+01,A21 = 0.100462E+01,A22 =-0.320985E-01,
     +A30 = 0.293119E+01,A31 =-0.403761E+00,A32 = 0.125153E-01,
     +BL00= 0.104133E+01,BL01=-0.943291E-02,BL02=-0.454758E-03,
     +BL10= 0.119253E+00,BL11= 0.407467E-01,BL12=-0.130718E-02,
     +BL20=-0.159391E-01,BL21= 0.727752E-02,BL22=-0.194405E-03,
     +B10 = 0.423071E+01,B11 =-0.610995E+00,B12 = 0.195531E-01,
     +B20 =-0.712527E+01,B21 = 0.969160E+00,B22 =-0.274255E-01,
     +B30 = 0.269925E+01,B31 =-0.363283E+00,B32 = 0.955316E-02)
C.    ------------------------------------------------------------------
C.
C             Ensure cut-off avoids infra-red catastrophe.
C
      IF (GEKIN.LE.BCUTE) GO TO 30
      KCASE=NAMEC(9)
C
*******************************>  Z3=Q(JPROB+2)
      Z3=(Z*(Z+1.))**0.3333333
      IF(Z3.LE.0.)GO TO 30
      Z32=Z3**2
      EEL1   = GETOT
      XC=BCUTE/GEKIN
      ALXC=LOG(XC)
      U=LOG(GEKIN/EMASS)
      U2=U**2
      V=LOG(Z)
*
      IF(GEKIN.LE.TLIM) THEN
         AL0=AL00+AL01*Z3+AL02*Z32
         AL1=AL10+AL11*Z3+AL12*Z32
         AL2=AL20+AL21*Z3+AL22*Z32
         AL=AL0+AL1*U+AL2*U2
         BL0=BL00+BL01*Z3+BL02*Z32
         BL1=BL10+BL11*Z3+BL12*Z32
         BL2=BL20+BL21*Z3+BL22*Z32
         BL=BL0+BL1*U+BL2*U2
         GMAX=1.+AL*XC+BL*XC**2
         IF(GEKIN.LT.0.0001) THEN
            G1=1.+AL+BL
            IF(G1.GT.GMAX) GMAX=G1
*
            IF(ABS(AL).GT.1.e-6) THEN
              X0=-BL/(2.*AL)
              IF((XC.LT.X0).AND.(X0.LT.1.)) THEN
                 G0=1.+AL*X0+BL*X0**2
                 IF(G0.GT.GMAX) GMAX=G0
              ENDIF
            ENDIF
*
         ENDIF
      ELSE
         U3=U2*U
         A1=A10+A11*Z3+A12*Z32
         A2=A20+A21*Z3+A22*Z32
         A3=A30+A31*Z3+A32*Z32
         AH=1.+A1/U+A2/U2+A3/U3
         B1=B10+B11*Z3+B12*Z32
         B2=B20+B21*Z3+B22*Z32
         B3=B30+B31*Z3+B32*Z32
         BH=0.75+B1/U+B2/U2+B3/U3
*
         F=4*V-0.55*V**2
         DEL0=136.*EMASS/(Z3*EEL1)
         EPC=XC*GEKIN/EEL1
         DC=DEL0*EPC/(1.-EPC)
         CC=42.392-F
*
         IF(DC.LE.1.) THEN
            DC2=DC**2
            F1=(42.392-7.796*DC+1.961*DC2-F)/CC
            F2=(41.734-6.484*DC+1.250*DC2-F)/CC
         ELSE
            F1=(42.24-8.368*LOG(DC+0.952)-F)/CC
            IF(F1.LT.0.) F1=0.
            F2=F1
         ENDIF
*
         GMAX=(1.-AH*EPC)*F1+BH*EPC**2*F2
      ENDIF
*
      CORR0=CORFAC*DENS*Z/A
      EPM=GEKIN/EEL1
      SC0=1.+CORR0/(EPM*EPM)
*
*     sample photon energy  ( according to 1/Ephoton)
*
   10 CALL GRNDM(RNDM,2)
*
      X=EXP(RNDM(1)*ALXC)
      EP=X*GEKIN/EEL1
*
*     Migdal correction for Ephoton->0. or no correction (Bethe)
*
#if !defined(CERNLIB_BETHE)
      CORR=SC0/(1.+CORR0/(EP*EP))
#endif
#if defined(CERNLIB_BETHE)
      CORR=1.
#endif
*
*     calculate rejection function g(x)
*
      IF(GEKIN.LE.TLIM) THEN
         G=1.+AL*X+BL*X**2
      ELSE
         D=DEL0*EP/(1.-EP)
         IF(D.LE.1.) THEN
            D2=D**2
            F1=(42.392-7.796*D+1.961*D2-F)/CC
            F2=(41.734-6.484*D+1.250*D2-F)/CC
         ELSE
            F1=(42.24-8.368*LOG(D+0.952)-F)/CC
            IF(F1.LT.0.) F1=0.
            F2=F1
         ENDIF
         G=(1.-AH*EP)*F1+BH*EP**2*F2
      ENDIF
      G=G*CORR/GMAX
      IF(RNDM(2).GT.G) GOTO 10
*
*     photon energy is sampled according to the Seltzer-Berger spectrum
*
      EGAMMA=EEL1*EP
C
C        CUT ON ENERGY THRESHOLD ?
C
      IF((IBREM.NE.1).OR.(EGAMMA.LE.CUTGAM)) THEN
         DESTEP = DESTEP + EGAMMA
         GO TO 20
      ENDIF
C
C             Generate emitted photon angles with respect to a Z-axis
C             defined along parent track. PHI is generated isotropically
C             and THETA is assigned a universal angular distribution
C
      EMASS1 = EMASS
      THETA  = GBTETH(EEL1, EMASS1, EP)*EMASS/EEL1
      SINTH  = SIN(THETA)
      COSTH  = COS(THETA)
      CALL GRNDM(RNDM,1)
      PHI    = TWOPI*RNDM(1)
      COSPHI = COS(PHI)
      SINPHI = SIN(PHI)
C
C             Polar co-ordinates to momentum components.
C
      NGKINE = NGKINE + 1
      GKIN(1,1)=EGAMMA*SINTH*COSPHI
      GKIN(2,1)=EGAMMA*SINTH*SINPHI
      GKIN(3,1)=EGAMMA*COSTH
      GKIN(4,1)=EGAMMA
      GKIN(5,1)=1.
      TOFD(1)  =0.
      GPOS(1,1) = VECT(1)
      GPOS(2,1) = VECT(2)
      GPOS(3,1) = VECT(3)
C
C             Rotate photon into GEANT system
C
      CALL GVROT(VECT(4),GKIN)
C
C             Correct track for lost energy
C
   20 CONTINUE
      GEKIN = GEKIN - EGAMMA
      GETOT = GEKIN + EMASS
      VECT(7)=SQRT (ABS((GETOT+EMASS)*GEKIN))
      CALL GEKBIN
C
C             Update probabilities
C
   30 CALL GRNDM(RNDM,1)
      ZINTBR=-LOG(RNDM(1))
C
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.2 1997/11/13 09:25:40 gunter
	6	* Correction By Laszlo Urban; protect against divide by 0. if AL becomes 0.
	7	* around GEKIN 2.93E-5.
	8	* Reported by harald@psiclu.cern.ch (signed by keller@biomed.ee.ethz.ch)
	9	*
	10	* Revision 1.1.1.1 1995/10/24 10:21:22 cernlib
	11	* Geant
	12	*
	13	*
	14	#include "geant321/pilot.h"
	15	*CMZ : 3.21/03 02/08/94 19.37.36 by S.Ravndal
	16	*-- Author :
	17	SUBROUTINE GBREME
	18	C.
	19	C. ******************************************************************
	20	C. * *
	21	C. * Simulates discrete hard BREMSSTRAHLUNG by electrons. *
	22	C. * *
	23	C. * The secondary photon energy is sampled from a *
	24	C. * parametrization of the bremsstrahlung calculation of *
	25	C. * Seltzer and Berger (NIM B12,p.95(1985)) for electron *
	26	C. * energies 1 keV - 10 GeV . For higher energies the *
	27	C. * parametrization agrees with the screened Bethe-Heitler *
	28	C. * bremsstrahlung spectrum.Migdal corrections are applied *
	29	C. * by default. The routine works ( together *
	30	C. * with the routines GBRELE and GBRSGE ) without the Migdal *
	31	C. * corrections using the Patchy switch +USE,BETHE. *
	32	C. * *
	33	C. * NOTE : *
	34	C. * BCUTE is the cut-off energy above which the photon energy *
	35	C. * spectrum is generated. *
	36	C. * *
	37	C. * ==>Called by : GTELEC *
	38	C. * Authors L.Urban 29/10/93 ********* *
	39	C. * *
	40	C. * 13/11/97 bug corrected by L.Urban *
	41	C. ******************************************************************
	42	C.
	43	#include "geant321/gcbank.inc"
	44	#include "geant321/gcphys.inc"
	45	#include "geant321/gcjloc.inc"
	46	#include "geant321/gconsp.inc"
	47	#include "geant321/gctrak.inc"
	48	#include "geant321/gcking.inc"
	49	#include "geant321/gcmate.inc"
	50	#include "geant321/gccuts.inc"
	51	DIMENSION RNDM(2)
	52	PARAMETER (CORFAC=0.805485E-10)
	53	PARAMETER (TLIM=0.000999)
	54	PARAMETER (
	55	+AL00=-0.205398E+01,AL01= 0.238815E-01,AL02= 0.525483E-03,
	56	+AL10=-0.769748E-01,AL11=-0.691499E-01,AL12= 0.222453E-02,
	57	+AL20= 0.406463E-01,AL21=-0.101281E-01,AL22= 0.340919E-03,
	58	+A10 = 0.467733E+01,A11 =-0.619012E+00,A12 = 0.202225E-01,
	59	+A20 =-0.734101E+01,A21 = 0.100462E+01,A22 =-0.320985E-01,
	60	+A30 = 0.293119E+01,A31 =-0.403761E+00,A32 = 0.125153E-01,
	61	+BL00= 0.104133E+01,BL01=-0.943291E-02,BL02=-0.454758E-03,
	62	+BL10= 0.119253E+00,BL11= 0.407467E-01,BL12=-0.130718E-02,
	63	+BL20=-0.159391E-01,BL21= 0.727752E-02,BL22=-0.194405E-03,
	64	+B10 = 0.423071E+01,B11 =-0.610995E+00,B12 = 0.195531E-01,
65	+B20 =-0.712527E+01,B21 = 0.969160E+00,B22 =-0.274255E-01,
66	+B30 = 0.269925E+01,B31 =-0.363283E+00,B32 = 0.955316E-02)
67	C. ------------------------------------------------------------------
68	C.
69	C Ensure cut-off avoids infra-red catastrophe.
70	C
71	IF (GEKIN.LE.BCUTE) GO TO 30
72	KCASE=NAMEC(9)
73	C
74	*******************************> Z3=Q(JPROB+2)
75	Z3=(Z(Z+1.))*0.3333333
76	IF(Z3.LE.0.)GO TO 30
77	Z32=Z3**2
78	EEL1 = GETOT
79	XC=BCUTE/GEKIN
80	ALXC=LOG(XC)
81	U=LOG(GEKIN/EMASS)
82	U2=U**2
83	V=LOG(Z)
84	*
85	IF(GEKIN.LE.TLIM) THEN
86	AL0=AL00+AL01Z3+AL02Z32
87	AL1=AL10+AL11Z3+AL12Z32
88	AL2=AL20+AL21Z3+AL22Z32
89	AL=AL0+AL1U+AL2U2
90	BL0=BL00+BL01Z3+BL02Z32
91	BL1=BL10+BL11Z3+BL12Z32
92	BL2=BL20+BL21Z3+BL22Z32
93	BL=BL0+BL1U+BL2U2
94	GMAX=1.+ALXC+BLXC**2
95	IF(GEKIN.LT.0.0001) THEN
96	G1=1.+AL+BL
97	IF(G1.GT.GMAX) GMAX=G1
98	*
99	IF(ABS(AL).GT.1.e-6) THEN
100	X0=-BL/(2.*AL)
101	IF((XC.LT.X0).AND.(X0.LT.1.)) THEN
102	G0=1.+ALX0+BLX0**2
103	IF(G0.GT.GMAX) GMAX=G0
104	ENDIF
105	ENDIF
106	*
107	ENDIF
108	ELSE
109	U3=U2*U
110	A1=A10+A11Z3+A12Z32
111	A2=A20+A21Z3+A22Z32
112	A3=A30+A31Z3+A32Z32
113	AH=1.+A1/U+A2/U2+A3/U3
114	B1=B10+B11Z3+B12Z32
115	B2=B20+B21Z3+B22Z32
116	B3=B30+B31Z3+B32Z32
117	BH=0.75+B1/U+B2/U2+B3/U3
118	*
119	F=4V-0.55V**2
120	DEL0=136.EMASS/(Z3EEL1)
121	EPC=XC*GEKIN/EEL1
122	DC=DEL0*EPC/(1.-EPC)
123	CC=42.392-F
124	*
125	IF(DC.LE.1.) THEN
126	DC2=DC**2
127	F1=(42.392-7.796DC+1.961DC2-F)/CC
128	F2=(41.734-6.484DC+1.250DC2-F)/CC
129	ELSE
130	F1=(42.24-8.368*LOG(DC+0.952)-F)/CC
131	IF(F1.LT.0.) F1=0.
132	F2=F1
133	ENDIF
134	*
135	GMAX=(1.-AHEPC)F1+BHEPC2F2
136	ENDIF
137	*
138	CORR0=CORFACDENSZ/A
139	EPM=GEKIN/EEL1
140	SC0=1.+CORR0/(EPM*EPM)
141	*
142	* sample photon energy ( according to 1/Ephoton)
143	*
144	10 CALL GRNDM(RNDM,2)
145	*
146	X=EXP(RNDM(1)*ALXC)
147	EP=X*GEKIN/EEL1
148	*
149	* Migdal correction for Ephoton->0. or no correction (Bethe)
150	*
151	#if !defined(CERNLIB_BETHE)
152	CORR=SC0/(1.+CORR0/(EP*EP))
153	#endif
154	#if defined(CERNLIB_BETHE)
155	CORR=1.
156	#endif
157	*
158	* calculate rejection function g(x)
159	*
160	IF(GEKIN.LE.TLIM) THEN
161	G=1.+ALX+BLX**2
162	ELSE
163	D=DEL0*EP/(1.-EP)
164	IF(D.LE.1.) THEN
165	D2=D**2
166	F1=(42.392-7.796D+1.961D2-F)/CC
167	F2=(41.734-6.484D+1.250D2-F)/CC
168	ELSE
169	F1=(42.24-8.368*LOG(D+0.952)-F)/CC
170	IF(F1.LT.0.) F1=0.
171	F2=F1
172	ENDIF
173	G=(1.-AHEP)F1+BHEP2F2
174	ENDIF
175	G=G*CORR/GMAX
176	IF(RNDM(2).GT.G) GOTO 10
177	*
178	* photon energy is sampled according to the Seltzer-Berger spectrum
179	*
180	EGAMMA=EEL1*EP
181	C
182	C CUT ON ENERGY THRESHOLD ?
183	C
184	IF((IBREM.NE.1).OR.(EGAMMA.LE.CUTGAM)) THEN
185	DESTEP = DESTEP + EGAMMA
186	GO TO 20
187	ENDIF
188	C
189	C Generate emitted photon angles with respect to a Z-axis
190	C defined along parent track. PHI is generated isotropically
191	C and THETA is assigned a universal angular distribution
192	C
193	EMASS1 = EMASS
194	THETA = GBTETH(EEL1, EMASS1, EP)*EMASS/EEL1
195	SINTH = SIN(THETA)
196	COSTH = COS(THETA)
197	CALL GRNDM(RNDM,1)
198	PHI = TWOPI*RNDM(1)
199	COSPHI = COS(PHI)
200	SINPHI = SIN(PHI)
201	C
202	C Polar co-ordinates to momentum components.
203	C
204	NGKINE = NGKINE + 1
205	GKIN(1,1)=EGAMMASINTHCOSPHI
206	GKIN(2,1)=EGAMMASINTHSINPHI
207	GKIN(3,1)=EGAMMA*COSTH
208	GKIN(4,1)=EGAMMA
209	GKIN(5,1)=1.
210	TOFD(1) =0.
211	GPOS(1,1) = VECT(1)
212	GPOS(2,1) = VECT(2)
213	GPOS(3,1) = VECT(3)
214	C
215	C Rotate photon into GEANT system
216	C
217	CALL GVROT(VECT(4),GKIN)
218	C
219	C Correct track for lost energy
220	C
221	20 CONTINUE
222	GEKIN = GEKIN - EGAMMA
223	GETOT = GEKIN + EMASS
224	VECT(7)=SQRT (ABS((GETOT+EMASS)*GEKIN))
225	CALL GEKBIN
226	C
227	C Update probabilities
228	C
229	30 CALL GRNDM(RNDM,1)
230	ZINTBR=-LOG(RNDM(1))
231	C
232	END