* * * Revision 1.1.1.1 1996/09/24 12:14:40 hneumann * GEANT base and includes imported * * Revision 1.1.1.1 1995/10/24 10:21:23 cernlib * Geant * * #include "geant321/pilot.h" *CMZ : 3.21/02 29/03/94 15.41.21 by S.Giani *-- Author : SUBROUTINE GDALIT(XM0,XM1,XM2,XM3,PCM) C. C. ****************************************************************** C. * * C. * This routine is based on GDECA3 which * C. * simulates three-body decays in center-of-mass. * C. * GDECA3 was Written by Jurgen Schukraft * C. * and Vincent Hedberg for R807. * C. * Adapted for Geant3 by Sverker Johansson (Aug-1987). * C. * (Quite extensive modifications; comments etc. are not always * C. * consistent with the code anymore. ) * C. * * C. * In GDECA3 the decay is done according to pure phase-space. * C. * This would result in wrong opening angles of * C. * the dielectrons from Pi0-Dalitz decay in LAB =45 deg * C. * Ebeam=1AGeV. * C. * but it has to be =15 deg * ************************************************************************ C. C. * Input: XM0: Mass of decaying Particle: XM0=Pi_0 mass C. * XM1: Gamma mass=0 C. * XM2, XM3: Electron mass me C. * Output: PCM(4,3) Four Momentum vectors of Gamma, C. * Electron and Positron C. * C. ****************************************************************** C. * This routine simulates PI0-DALITZ decay Pi0 -> Gamma e+e- * C. * in the center-of-mass system according to the * C. * Kroll-Wada Decay rate * C. * (N.Kroll & W.Wada, Phys.Rev.98 (1955)1355): * C. * * C. * d Gamma/dxdy ~ (1-xinv)^3*{RC/xinv+yinv^2+1}/2 * C. * * C. * with: * C. * RC: =4*(me/mpi0)^2 (=4*epsilon) * C. * xinv: The squared invariant mass of the e+e- pair, * c. * normalized to the pion mass: RC < xinv < 1 * C. * yinv: Variable describing the energy distribution between * C. * electron and positron; in contrast to xinv * C. * not lorentz invariant; only * C. * for symmetry reasons called yinv: * C. * |yinv|Gamma * C. * \e- * C. * * C. * ==>Called by : GDECAY * C. * * C. * Author H.SCHOEN (NEUMANN) ********* * C. * comments: 13/05/1996 H.S. * C. * * C. * ==>Called by : GDECAY * C. * Author S. Johansson ********* * C. * * C. ****************************************************************** C. *************************************************************** ******* original name : ******************* ******SUBROUTINE D E C 3 B (ID,IP,K1,K2,K3)******************* *************************************************************** # include "geant321/gconst.inc" C IMPLICIT NONE REAL PI2,CT1,ST1,XM3,GINV REAL XINV,YINV,BETINV,RC,Q1,Q2,Q3,G1,G2,EM,XM0,XM1,XM2 REAL GAM,E1,E2,E3,PGSQ,CTST,STST,F1,F2 REAL PCM(4,3) REAL RNDM(3) * PI2 = 2. * 3.141592654 C-------------------------------------------------------- C SIMULATION OF SECONDARY PARTICLE MOMENTA C IN THE REST FRAME OF PARENT PARTICLE C-------------------------------------------------------- ****** * lower limit for XINV **** RC = (2*XM2/XM0)**2 ***** * upper limit for the absolute valu of YINV ***** 10 CALL GRNDM(RNDM,3) ****** * random number xinv ****** XINV = RNDM(1) ****** * new random number for XINV below lower limit RC ****** IF (XINV.LE.RC) GOTO 10 BETINV = SQRT(1.-RC/XINV) ****** * random number yinv ****** YINV = BETINV*(2.*RNDM(2)-1.) GINV = RNDM(3)*2./RC GAM = XM0 * (1.-XINV)**3 * (1.+YINV**2+RC/XINV)/XINV IF (GINV.GT.GAM) GOTO 10 ****** * Energy for e+, e- (E2, E3) and Gamma (E1) ****** E2 = XM0*(XINV+1.+YINV*(1.-XINV))/4. E3 = XM0*(XINV+1.-YINV*(1.-XINV))/4. E1 = XM0-E2-E3 ****** * Absolute Momentum e+/e- (Q3,Q2) and Gamma Q1 ****** Q3 = SQRT(E3*E3-XM2*XM3) Q2 = SQRT(E2*E2-XM2*XM3) Q1 = E1 ****** * PGAMMA squared: PGSQ = Q1**2=E1**2 ****** PGSQ = (XM0*(1.-XINV)/2.)**2 ****** * Cos(theta*) * stands for virtual gamma ****** CTST = (PGSQ-Q3**2-Q2**2)/(2.*Q3*Q2) ****** * Sin(theta*) ****** STST = SQRT(ABS(1.-CTST**2)) CALL GRNDM(RNDM,3) ***** * Random polar angle t1 for one of the leptons im CMS * CT1 = cos(t1) * ST1 = sin(t1) ***** CT1 = 2.*RNDM(1)-1. ST1 = SQRT(ABS(1.-CT1*CT1)) ***** * random azimuth angle f1 for this lepton ***** F1 = PI2*RNDM(2) ***** * Fourmomentum vector of this lepton ***** PCM(1,3) = Q3*ST1*COS(F1) PCM(2,3) = Q3*ST1*SIN(F1) PCM(3,3) = Q3*CT1 PCM(4,3) = E3 ***** * random azimuthal f2 for the other lepton. ***** F2 = PI2*RNDM(3) PCM(1,2) = Q2*( STST*SIN(F2)*CT1*COS(F1) + STST*COS(F2)*SIN(F1) & + CTST*ST1*COS(F1)) PCM(2,2) = Q2*( STST*SIN(F2)*CT1*SIN(F1) - STST*COS(F2)*COS(F1) & + CTST*ST1*SIN(F1)) PCM(3,2) = Q2*(-STST*SIN(F2)*ST1 + CTST*CT1) PCM(4,2) = E2 PCM(1,1) = -PCM(1,3)-PCM(1,2) PCM(2,1) = -PCM(2,3)-PCM(2,2) PCM(3,1) = -PCM(3,3)-PCM(3,2) PCM(4,1) = E1 c Print *,'Modified gdalit: Dalitz-decay' RETURN END