* * $Id: gukine.F,v 1.14 2004/12/08 14:43:24 davidl Exp $ * * Revision 1.1.1.1 1995/10/24 10:21:52 cernlib * Geant * * c To enable beam motion spreading, define the beam box size below (cm) c#define BEAM_BOX_SIZE 5 #include "geant321/pilot.h" *CMZ : 3.21/02 29/03/94 15.41.35 by S.Giani *-- Author : SUBROUTINE GUKINE * ************************************************************************ * * * Generates Kinematics for primary tracks * * * ************************************************************************ * #include "geant321/gcunit.inc" #include "geant321/gcflag.inc" #include "geant321/gckine.inc" #include "geant321/gconsp.inc" #include "geant321/gcscan.inc" #include "geant321/gcomis.inc" * DIMENSION VERTEX(3),PLAB(3) DIMENSION RNDM(20) * include 'cobrems.inc' integer coSwitch real freqMaximum common /coherentGen/coSwitch,freqMaximum data coSwitch,freqMaximum/0,0/ save /coherentGen/ real xMinimum,Theta02,CollimPos,RadiaPos parameter (xMinimum=1e-6) parameter (Theta02=1.8) parameter (CollimPos=-2200.0,RadiaPos=0.0) integer nubuf real ubuf(10) * * ----------------------------------------------------------------- * ev = IDEVT do i=1,10 ev = ev/10. if (ev.lt.10) goto 2 enddo 2 if (int(ev).eq.ev) then write(LOUT,*) IDEVT," events simulated" endif * * Try input from MonteCarlo generator first * itry = nextInput() if (itry .eq. 0) then itry = loadInput() elseif (itry .ne. 9) then ieorun = 1 ieotri = 1 return * * Try coherent bremsstrahlung beam generation next * elseif (E.gt.0) then call GRNDM(RNDM,7) phim = RNDM(1)*TWOPI rhom = mospread*sqrt(-2*log(RNDM(2))) thxMosaic = rhom*cos(phim) thyMosaic = rhom*sin(phim) phib = RNDM(3)*TWOPI varBeam = (emit/spot)**2 rhob = sqrt(-2*varBeam*log(RNDM(4))) thxBeam = rhob*cos(phib) thyBeam = rhob*sin(phib) phis = RNDM(5)*TWOPI varMS = sigma2MS(t*RNDM(6)) rhos = sqrt(-2*varMS*log(RNDM(7))) thxMS = rhos*cos(phis) thyMS = rhos*sin(phis) cos45 = 1/sqrt(2d0) rotate(1,1) = 0 rotate(1,2) = cos45 !point (1,0,0) along beam rotate(1,3) = -cos45 !point (0,1,1) vertically rotate(2,1) = 0 rotate(2,2) = cos45 rotate(2,3) = cos45 rotate(3,1) = 1 rotate(3,2) = 0 rotate(3,3) = 0 call rotmat(rotate,thxBeam+thxMS-thx-thxMosaic,0d0,0d0) call rotmat(rotate,0d0,thyBeam+thyMS-thy-thyMosaic,0d0) do i=1,10000 call GRNDM(RNDM,4) phi = RNDM(1)*TWOPI x = xMinimum**RNDM(2) theta2 = Theta02*RNDM(3)/(1-RNDM(3)+1e-30) coSwitch = coSwitch+1 if (coSwitch/2*2.eq.coSwitch) then !try coherent generation freq = dNcdxdp(x,phi) f = freq*RNDM(3) do ip=1,q2points if (f.le.q2weight(ip)) then theta2 = q2theta2(ip) goto 5 endif enddo 5 continue ppol = polarization(x,theta2) else !try incoherent generation freq = dNidxdt2(x,theta2) freq = freq*(Theta02+theta2)**2/(TWOPI*Theta02) ppol = 0 endif freq = x*freq if (freq.gt.freqMaximum) then freqMaximum = freq*1.2 elseif (freq.ge.freqMaximum*RNDM(4)) then goto 50 endif enddo write(5,*) 'Gukine: photon beam generator failed, giving up' 50 continue theta = sqrt(theta2)*(me/E) thetaX = thxMS+theta*cos(phi) thetaY = thyMS+theta*sin(phi) PLAB(3) = E*x PLAB(1) = PLAB(3)*thetaX PLAB(2) = PLAB(3)*thetaY call GRNDM(RNDM,2) phic = RNDM(1)*TWOPI rhoc = spot*sqrt(-2*log(RNDM(2))) * VERTEX(1) = (rhoc*cos(phic) + D*thetaX)*100 * VERTEX(2) = (rhoc*sin(phic) + D*thetaY)*100 VERTEX(1) = (rhoc*cos(phic) *100) ! + D*thetaX)*100 VERTEX(2) = (rhoc*sin(phic) *100) ! + D*thetaY)*100 VERTEX(3) = RadiaPos !CollimPos ubuf(1) = ppol nubuf = 1 #if defined BEAM_BOX_SIZE call GRNDM(RNDM,2) ubuf(2) = (RNDM(1)-0.5)*BEAM_BOX_SIZE ubuf(3) = (RNDM(2)-0.5)*BEAM_BOX_SIZE VERTEX(1) = VERTEX(1) + ubuf(2) VERTEX(2) = VERTEX(2) + ubuf(3) nubuf = 3 #endif call GSVERT(VERTEX,0,0,ubuf,nubuf,NVERT) call GSKINE(PLAB,1,NVERT,0,0,NT) call storeInput(IDRUN,IDEVT,NT); c post in the post-bremsstrahlung electron as well PLAB(3) = E-PLAB(3) PLAB(1) = -PLAB(1) PLAB(2) = -PLAB(2) call GSKINE(PLAB,3,NVERT,0,0,NT) * * If all else fails, do automatic single-track generation * else VERTEX(1)=PKINE(4) VERTEX(2)=PKINE(5) VERTEX(3)=PKINE(6) if (vertex(1).eq.0.and.vertex(2).eq.0.and.vertex(3).eq.0) then call GRNDM(rndm,2) phiv = rndm(1)*TWOPI rhov = sqrt(-2*log(rndm(2))) vertex(1)=0.17*rhov*cos(phiv) vertex(2)=0.05*rhov*sin(phiv) endif #if defined BEAM_BOX_SIZE call GRNDM(RNDM,2) VERTEX(1) = VERTEX(1) + (RNDM(1)-0.5)*BEAM_BOX_SIZE VERTEX(2) = VERTEX(2) + (RNDM(2)-0.5)*BEAM_BOX_SIZE #endif IF(IKINE.GT.100)THEN IK=IKINE-100 THETA=PKINE(2)*DEGRAD PHI=PKINE(3)*DEGRAD ELSE IK=IKINE CALL GRNDM(RNDM,2) THETA=PI*RNDM(1) PHI=TWOPI*RNDM(2) ENDIF CALL GRNDM(RNDM,2) PMAG = PKINE(1) !*rndm(1) PLAB(1) = PMAG*SIN(THETA)*COS(PHI) PLAB(2) = PMAG*SIN(THETA)*SIN(PHI) PLAB(3) = PMAG*COS(THETA) CALL GSVERT(VERTEX,0,0,0.0,0,NVERT) CALL GSKINE(PLAB,IK,NVERT,0,0,NT) call storeInput(IDRUN,IDEVT,NT); endif * * Kinematic debug (controled by ISWIT(1)) * IF(IDEBUG.EQ.1.AND.ISWIT(1).EQ.1) THEN CALL GPRINT('VERT',0) CALL GPRINT('KINE',0) ENDIF * END