c***************************************************************** subroutine cosgen c generates cosmics track c c c last modified on: 21/03/2005 by R.Holzmann c***************************************************************** implicit none #include "geant321/gclist.inc" #include "geant321/gckine.inc" #include "geant321/gcflag.inc" #include "hcons.inc" #include "kinetups.inc" #include "geatarget.inc" #include "user.inc" #include "pmc.inc" real p_lab(3), buf(4), E_cine, xi, yi, fi, p real theta, phi, ptot, vertex(3), restm real rad_setup, sint, cost, sinp, cosp real e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z real x1, x2, x3, v1, v2, v3 integer jpart, itrack, nvtx, i character*20 cdum1 real dum2, dum3, dum4, dum5, dum6 real rndm c******************************************************************* c c pmc0(1) = radius of setup in mm * c c******************************************************************* rad_setup = pmc0(1) if(rad_setup.gt.10000.) rad_setup = 2000.0 ! in mm E_cine = 0.753/(1.075-rndm(1)) ! sample energy according to 1/E**2 law ! (in GeV) p = rndm(2) ! sample zenith angle according to ! cos(theta)**2 law (in rad) xi = 1.570796*(0.9+p)*p ! initial guess do i=1,5 ! find root iteratively fi = 1.0 + cos(xi) yi = xi + sin(xi) - 3.141593*p xi = (fi*xi-yi)/fi end do theta = 0.5*xi phi = 6.283185*rndm(3) ! sample azimuthal angle uniformly p = rndm(4) if(p.lt.0.539) then jpart = 5 ! make mu+ else if(p.gt.0.569) then jpart = 6 ! make mu- else jpart = 14 ! make proton end if call gfpart(jpart,cdum1,dum2,restm,dum3,dum4,dum5,dum6) ! get mass sint = sin(theta) cost = cos(theta) sinp = sin(phi) cosp = cos(phi) e1x = sint*cosp ! incidence of particle e1y = cost ! with y-axis vertical e1z = sint*sinp e2x = cost*cosp ! e2, e3 define a plane perpendicular to e1 e2y = -sint e2z = cost*sinp e3x = sinp e3y = 0.0 e3z = -cosp x1 = 1.0 i = 4 10 i = i+1 x2 = 2.0*rndm(i)-1.0 ! sample point in (e1,e2) plane x3 = 2.0*rndm(i+1)-1.0 ! into disc of radius 1 if(x2*x2+x3*x3 .ge. 1.0) goto 10 v1 = x1*e1x + x2*e2x + x3*e3x ! displace incident particle v2 = x1*e1y + x2*e2y + x3*e3y ! laterally on this disc v3 = x1*e1z + x2*e2z + x3*e3z vertex(1) = v1*rad_setup ! scale up to setup radius vertex(2) = v2*rad_setup vertex(3) = v3*rad_setup vertex(1) = vertex(1) + tgtpar(1) ! shift vertex on demand vertex(2) = vertex(2) + tgtpar(2) ! (=useful to concentrate vertex(3) = vertex(3) + tgtpar(3) ! cosmics on part of device) ptot = sqrt((restm+E_cine)**2 - restm**2) ! total momentum (GeV/c) p_lab(1) = -ptot*e1x p_lab(2) = -ptot*e1y p_lab(3) = -ptot*e1z buf(1) = 300.0 ! generator type buf(2) = 0.0 buf(3) = -1. buf(4) = 1.0 ! weight do i=1,3 vertex(i) = vertex(i)*0.10 ! go from mm to cm for GEANT enddo call gsvert(vertex,0,0,buf,0,nvtx) ! pass vertex call gskine(p_lab,jpart,nvtx,buf,4,itrack) ! pass particle call copykine(itrack,jpart,p_lab,vertex,0,0,1.,-1., + buf(1),buf(2),buf(3)) ebeam = E_cine*1000. ! cosmic particle energy bpar = 57.295*theta ! zenith angle phi_event = 57.295*phi ! azimuth return end