diff --git a/gphys/gdalet.F b/gphys/gdalet.F new file mode 100644 index 0000000..d97ce07 --- /dev/null +++ b/gphys/gdalet.F @@ -0,0 +1,161 @@ +c------------------------------------------------------------------------------ + real function vdmform(x,rmas,isw) ! VDM formfactor + real*4 x, rmas, dalit, rhom ! masses in MeV + real*4 f, delta + parameter (rhom=770.0) +C + vdmform = 0.0 + if (x.ge.rmas) return + if (x.lt.2.0*0.511) return + + f = sqrt(1.-4.*(0.511/x)**2) * (1.+2.*(0.511/x)**2) + if (isw.eq.1) then ! pseudo-scalar meson -> gamma e+ e- + dalit = ((1.-(x/rmas)**2)**3/(x/rmas))*1./(1.-(x/rhom)**2)**2 + vdmform = 8.2e-4*dalit*f + else if (isw.eq.2) then ! vector meson -> pi0 e+ e- + delta = (135.0/rmas)**2 + bracket = (1.+(x/rmas)**2/(1.-delta))**2 + + - 4.*(x/rmas)**2/(1.-delta)**2 + if(bracket.le.0.0) return + ff2 = 0.64**4/((0.65**2-(x/1000.)**2)**2 + (0.042*0.65)**2) + dalit = (bracket**1.5)/(x/rmas)*ff2 + vdmform = 8.2e-4*dalit*f + end if +C + return + end +c------------------------------------------------------------------------------ + real function qedform(x,rmas,isw) ! QED formfactor + real*4 x, rmas, dalit ! masses in MeV + real*4 f, delta +C + qedform = 0.0 + if (x.ge.rmas) return + if (x.lt.2.0*0.511) return + + f = sqrt(1.-4.*(0.511/x)**2) * (1.+2.*(0.511/x)**2) + if (isw.eq.1) then ! pseudo-scalar meson + dalit = ((1.-(x/rmas)**2)**3/(x/rmas)) + qedform = 8.2e-4*dalit*f + else if (isw.eq.2) then ! vector meson + delta = (135.0/rmas)**2 + bracket = (1.+(x/rmas)**2/(1.-delta))**2 + + - 4.*(x/rmas)**2/(1.-delta)**2 + if(bracket.le.0.0) return + dalit = (bracket**1.5)/(x/rmas) + qedform = 8.2e-4*dalit*f + end if +C + return + end +c------------------------------------------------------------------------------ + +c------------------------------------------------------------------------------ +c Generation of 2 particle mass according to the VDM form factor +c for Dalitz decay of (id) 1=pi0 2=eta to (e+e-) gamma and 3=omega to pi0 e+e- +c------------------------------------------------------------------------------ + + real function dalitz_evgen(id) + real xmin(3), xmax(3), wmax(3), mesonmass(3), rndm(2) + integer isw(3) + data mesonmass/134.9766,547.75,782.59/ + data isw/1,1,2/ + data xstep/0.1/ + logical maxfound(3) + data maxfound/3*.FALSE./ + save maxfound, wmax + + if (.NOT. maxfound(id)) then + wmax(id) = 0.0 + xmas = 0.0 + 1 xmas = xmas + xstep + ff = vdmform(xmas,mesonmass(id),isw(id)) + if (ff .gt. wmax(id)) wmax(id) = ff + if (xmas .lt. mesonmass(id)) goto 1 + maxfound(id) = .TRUE. + end if + + 2 call grndm(rndm,2) + rnd1 = rndm(1)*mesonmass(id) + rnd2 = rndm(2)*wmax(id) + ff = vdmform(rnd1,mesonmass(id),isw(id)) + if (ff .lt. rnd2) goto 2 + + dalitz_evgen = rnd1/1000. + + return + end + +c------------------------------------------------------------------------------ + + SUBROUTINE GDALET(XM0,XM1,XM2,XM3,PCM,IDH) +C. +C. ****************************************************************** +C. * * +C. * Dalitz decay of pi0, eta * +C. * * +C. * * +C. * ==>Called by : GDECAY * +C. * mod. from gdec3; including eta formfactor calc. in VDM) * +C. * W.Schoen 12.8.96 * +C. * tritouille par R.H. 24.3.97 and 28.8.97 and 22.8.2000 * +C. ****************************************************************** + IMPLICIT NONE + REAL FMASS, PI2, CT1, ST1, GINV, XINV, YINV + REAL BETINV, RC, Q1, Q2, Q3, G1, G2, EM, XM0, XM1, XM2, XM3 + REAL GAM, E1, E2, E3, PGSQ, CTST, STST, F1, F2 + REAL PCM(4,3) + REAL XMETA, RNDM(4) + INTEGER IDH + REAL dalitz_evgen +* + + PI2 = 2.*3.141592654 +C-------------------------------------------------------- +C SIMULATION OF SECONDARY PARTICLE MOMENTA +C IN THE REST FRAME OF PARENT PARTICLE +C-------------------------------------------------------- + RC = (2.*XM2/(XM0-XM1))**2 +C XINV = HRNDM1(IDH)/1000. ! convert to GeV here + + CALL GRNDM(RNDM,4) + XINV = dalitz_evgen(IDH) + XINV = (XINV/(XM0-XM1))**2 + YINV = SQRT(ABS(1. - RC/XINV))*(2.*RNDM(2) - 1.) + IF(RNDM(1).LT.0.5) YINV = -YINV + E2 = (XM0-XM1)*(XINV+1.+YINV*(1.-XINV))/4. + E3 = (XM0-XM1)*(XINV+1.-YINV*(1.-XINV))/4. + E1 = XM0-E2-E3 + + Q3 = SQRT(E3*E3-XM2*XM3) + Q2 = SQRT(E2*E2-XM2*XM3) + Q1 = SQRT(E1*E1-XM1*XM1) + PGSQ = ( (XM0**2-XM1**2)*(1.-XINV)/(2.*XM0) )**2 + CTST = (PGSQ-Q3**2-Q2**2)/(2.*Q3*Q2) + STST = SQRT(ABS(1.-CTST**2)) + + CT1 = 2.*RNDM(2)-1. + ST1 = SQRT(ABS(1.-CT1*CT1)) + + F1 = PI2*RNDM(3) + + PCM(1,3) = Q3*ST1*COS(F1) + PCM(2,3) = Q3*ST1*SIN(F1) + PCM(3,3) = Q3*CT1 + PCM(4,3) = E3 + F2 = PI2*RNDM(4) + + 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 + + RETURN + END diff --git a/gphys/gdecay.F b/gphys/gdecay.F index 043b9c6..3340936 100644 --- a/gphys/gdecay.F +++ b/gphys/gdecay.F @@ -1,7 +1,10 @@ * -* $Id$ +* $Id: gdecay.F,v 1.2 2003/11/28 11:23:56 brun Exp $ * * $Log: gdecay.F,v $ +* Revision 1.2 2003/11/28 11:23:56 brun +* New version of geant321 with all geant3 routines renamed from G to G3 +* * Revision 1.1.1.1 2002/07/24 15:56:25 rdm * initial import into CVS * @@ -148,8 +151,38 @@ SUBROUTINE G3DECAY GO TO 99 41 DMASS = RMASS END IF - CALL G3DECA3(DMASS,XM(1),XM(2),XM(3),PCM) + + + IF (IPART.EQ.7) THEN ! pi0 + IF (MXX.EQ.30201) THEN ! pi0 -> gamma e+e- (=30201) + CALL GDALET(DMASS,XM(1),XM(2),XM(3),PCM,1) + ELSE ! phase space + CALL G3DECA3(DMASS,XM(1),XM(2),XM(3),PCM) + ENDIF + ENDIF + + IF (IPART.EQ.17) THEN ! eta + IF(MXX.EQ.30201) THEN ! eta -> gamma e+e- (=30201) + CALL GDALET(DMASS,XM(1),XM(2),XM(3),PCM,2) + ELSE ! phase space + CALL G3DECA3(DMASS,XM(1),XM(2),XM(3),PCM) + ENDIF + ENDIF + + IF (IPART.EQ.33) THEN ! omega + IF(MXX.EQ.30207) THEN ! omega -> pi0 e+e- + CALL GDALET(DMASS,XM(1),XM(2),XM(3),PCM,3) + ELSE ! phase space + CALL G3DECA3(DMASS,XM(1),XM(2),XM(3),PCM) + ENDIF + ENDIF + + IF (IPART.NE.7.AND.IPART.NE.17.AND.IPART.NE.33) THEN ! other particle + CALL G3DECA3(DMASS,XM(1),XM(2),XM(3),PCM) + ENDIF + ENDIF + C C LORENTZ boost into LAB system defined along parent vector C followed by rotation back into GEANT system.