c********************************************************************** subroutine varangle(nvtx,xvert,deltaphi) c Generates 1 track of given particle type (or a close pair) per event. c Polar (theta) and azimuthal (phi) angles and momentum (p) can be chosen c in an interval. For close pairs, the opening angle and the orientation c can be chosen in an interval. c c Needs IKINE>0 to be set, e.g. with GEANT command kine 1 ... c c last modified on: 13/09/2007 by R.Holzmann c********************************************************************** implicit none integer nvtx real xvert(3), deltaphi #include "geant321/gclist.inc" #include "geant321/gckine.inc" #include "geant321/gcflag.inc" #include "hcons.inc" #include "kinetups.inc" #include "geatarget.inc" #include "user.inc" REAL P(3), P2(3), BUF(4), VECL6(6), PABS, PHI, THE REAL CTHEMIN, CTHEMAX, CTH, STH, CPH, SPH, CAL, SAL, CBE, SBE REAL ALPHA, ALFMIN, ALFMAX, BETA, NORM INTEGER IPHI, ITHE, IMOM, ID1, ID2 INTEGER ITRACK, I REAL PMIN, PMAX, PHIMIN, PHIMAX, THEMIN, THEMAX c****************************************************** c PKINE(1) = THEMIN , PKINE(2) = THEMAX * c PKINE(3) = PHIMIN , PKINE(4) = PHIMAX * c PKINE(5) = MOM min , PKINE(6) = MOM max * c PKINE(7) = PARTICLE TYPE * c PKINE(8) = ALPHAMIN , PKINE(9) = ALPHAMAX * c PKINE(10) = BETA * c****************************************************** THEMIN = PKINE(1) THEMAX = PKINE(2) PHIMIN = PKINE(3) PHIMAX = PKINE(4) PMIN = PKINE(5) PMAX = PKINE(6) IF (PKINE(7) .LT. 100.) THEN ! 1 id ID1 = INT(PKINE(7)+0.001) ID2 = ID1 ELSE ! decode 2 ids ID2 = PKINE(7)/1000. ID1 = PKINE(7) - ID2*1000. ENDIF ALFMIN = PKINE(8) ALFMAX = PKINE(9) BETA = PKINE(10) call grndm(vecl6,6) cthemin = cos(themin*dtor) cthemax = cos(themax*dtor) cth = cthemin + vecl6(1)*(cthemax-cthemin) sth = sqrt(1.0-cth**2) the = rtod*acos(cth) phi = phimin + vecl6(2)*(phimax-phimin) phi = phi + deltaphi ! switch sector pabs = pmin + vecl6(3)*(pmax-pmin) cph = cos(phi*dtor) sph = sin(phi*dtor) c convert now from MeV/c to GeV/c for GEANT p(1) = 0.001*pabs*sth*cph p(2) = 0.001*pabs*sth*sph p(3) = 0.001*pabs*cth buf(1) = 200.0 ! generator type buf(2) = 0.0 buf(3) = -1. buf(4) = 1.0 ! weight if(lprin(1).eq.1) + print *,'theta, phi, p in varangle.F: ', the, phi, pabs call gskine(p,id1,nvtx,buf,4,itrack) if(imacc.eq.1) then ! store for acceptance calculation call hf2(20,pabs,the,1.0) end if call copykine(itrack,id1,p,xvert,0,0,1.,-1., + buf(1),buf(2),buf(3)) ebeam = pabs ! beam energy in MeV bpar = the ! theta in degrees phi_event = phi ! use this to save phi of the one-shot event if(alfmin.lt.0.0 .or. alfmin.gt.180.0) return if(alfmax.le.0.0 .or. alfmax.gt.180.0) return ! else do partner particle (of close pair) ! with opening angle alpha and azimuth beta ! (see logbook p. 182 for details) if(alfmax.ge.alfmin) then ! uniform opening angle alpha = alfmin + vecl6(4)*(alfmax-alfmin) else ! this is like true gamma->e+e- conversion alpha = alfmin -0.456*log(vecl6(4)) endif if(beta.lt.0. .or. beta.gt.360.) beta = 360.*vecl6(5) ! make random beta cal = cos(alpha*dtor) sal = sin(alpha*dtor) cbe = cos(beta*dtor) sbe = sin(beta*dtor) norm = sqrt(sth*sth*cph*cph+cth*cth) c if() pabs = pmin + vecl6(6)*(pmax-pmin) ! |p2| .ne. |p1| p2(1) = 0.001*pabs*((sal*cbe*cth-sal*sbe*sth*sth*cph*sph)/norm + + cal*sth*cph) p2(2) = 0.001*pabs*(sal*sbe*(sth*sth*cph*cph+cth*cth)/norm + + cal*sth*sph) p2(3) = 0.001*pabs*((-sal*cbe*sth*cph-sal*sbe*cth*sth*sph)/norm + + cal*cth) call gskine(p2,id2,nvtx,buf,4,itrack) bpar = alpha ! opening angle in degrees phi_event = beta ! orientation in degrees call copykine(itrack,id2,p2,xvert,0,0,1.,-1., + buf(1),buf(2),buf(3)) end