subroutine gukine c c kinematics for the primary interaction and tracks c called by GTRIG c implicit none include 'event.cmn' include 'guser.cmn' include 'gcflag.cmn' include 'gconst.cmn' include 'gcunit.cmn' include 'gmail.cmn' include 'ukine.cmn' include 'montb.cmn' include 'geninfo.cmn' include 'uevent.cmn' include 'vertices.cmn' logical first/.true./ integer ipart, nvtx, nt real vrtx(3), ppar(3), pmom, etot real phigh/100./ real rmass real rndm(3), ubuf(1) real pelec(4), mupol(3), imu real costh, phi integer i,j,k, ierr integer irct integer modder/1/ c c event number message c if(mod(ievent,modder).eq.0) then write(6,*) ' Event:',ievent if(ievent.eq.modder*10) modder = modder*10 end if c c initialize number of vertices and tracks c nvtx = 0 nt = 0 c c initialize a few user parameter c iuneu = 0 iumod = 0 iprnt = 0 call vzero(iupar,nparmax) nphsum = 0 nphpmt = 0 call vzero(npmt,nset) call vzero(uvert,3) call vzero(upini,3*nparmax) plepton = 0. call vzero(ustop,3) call vzero(upfin,3) ssleng = 0. c c c write random number at beginning of event if wanted c if(iswit(1).gt.0) then write(lout,*) ' Random number at the beginning of event' write(lout,*) nrndm(1),nrndm(2) end if c c branch according to event type flag c c if ipar > 0 this is a single particle with fixed momentum PABS c if ipar = 0 use the file of generated events c if ipar < 0 cosmic ray particle with momentum between pabs, phigh c if ipar < -1000 cosmic ray particle entering from side only c if (ipar.gt.0) then ! single particle; fixed mom; random dir nupar = 1 iupar(1) = ipar ubuf(1) = ipar ubuf(1) = 1 plepton = pabs c call nu_vert(vrtx) c fix position and direction for the moment vrtx(1) = 0. vrtx(2) = 0. vrtx(3) = 0. c call randir(ppar) ppar(1) = 1. ppar(2) = 0. ppar(3) = 0. if (ipar.eq.200) then c c generate michel electron with prper spectrum c do i=1,3 mupol(i) = verdir(i,ievent) enddo call mudecay(pelec,mupol,1) vrtx(1)=verpos(1,ievent) vrtx(2)=verpos(2,ievent) vrtx(3)=verpos(3,ievent) c ppar(1) = verdir(1,ievent) c ppar(2) = verdir(2,ievent) c ppar(3) = verdir(3,ievent) pabs = sqrt(pelec(1)**2+pelec(2)**2+pelec(3)**2) ppar(1) = pelec(1)/pabs ppar(2) = pelec(2)/pabs ppar(3) = pelec(3)/pabs endif pupar(1) = pabs do i=1,3 uvert(i) = vrtx(i) upini(i,1) = ppar(i) ppar(i) = pabs*ppar(i) end do call gsvert(vrtx,0,0,ubuf,1,nvtx) call gskine(ppar,ipar,nvtx,ubuf,1,nt) else if (ipar.lt.0) then ! cosmic ray particle 672 continue nupar = 1 iupar(1) = iabs(ipar) call cosgeo(vrtx,ppar,pmom) pupar(1) = pmom if ((pmom.lt.pabs).or.(pmom.gt.phigh)) goto 672 ! momentum cut if (ipar.lt.-1000) then if(abs(vrtx(3)).gt.725.) goto 672 ! demand vertex on side end if ipart = iabs(ipar) if(ipart.gt.1000) ipart = ipart-1000 ubuf(1) = ipart ubuf(1) = 1 plepton = pmom do i=1,3 uvert(i) = vrtx(i) upini(i,1) = ppar(i) ppar(i) = pmom*ppar(i) end do call gsvert(vrtx,0,0,ubuf,1,nvtx) call gskine(ppar,ipart,nvtx,ubuf,1,nt) else if(ipar.eq.0) then ! read mc generated events from file c if(first) then c first = .false. c open(unit=62,file='mc_gen.dat',form='UNFORMATTED', c + status='OLD', err=998) c end if 10 continue ! generator event loop call readev(ierr) if(ierr.ne.0) goto 888 c c generate event location within the fiducial volume c or use the one from the file c if(ifidu.eq.1) then call nu_vert(vrtx) ubuf(1) = 4. ! neutrino call gsvert(vrtx,0,0,ubuf,1,nvtx) else ubuf(1) = 4. vrtx(1) = ppos(1,1) vrtx(2) = ppos(1,2) vrtx(3) = ppos(1,3) call gsvert(vrtx,0,0,ubuf,1,nvtx) endif c c get reaction type c ireact = iparm(1) pneu(1) = float(iparm(3))/1000. pneu(2) = float(iparm(4))/1000. pneu(3) = float(iparm(5))/1000. pneu(4) = float(iparm(6))/1000. orgneu(1) = float(iparm(7))/1000. orgneu(2) = float(iparm(8))/1000. orgneu(3) = float(iparm(9))/1000. rctq2 = float(iparm(2))/1000. ifreepr = iparm(22) c c now depending on the reaction type store the most muon like particle c actually it doesn't matter which reaction. By convention the most c muon like particle is the first one c c c now loop through the generated particles and put them on the c stack c do j = 1, npart c c first get the energy and the momentum c etot = pkin(j) ! in MeV etot = etot/1000. ! convert to GeV rmass = qmass(j)/1000. ! convert to GeV pmom = sqrt(etot**2 - rmass**2) ppar(1) = pdir(j,3)*pmom ! z -> x ppar(2) = pdir(j,1)*pmom ppar(3) = pdir(j,2)*pmom c c get the particle type c ipart = 0 ubuf(1) = 0 if ((qmass(j).eq.0.).and.(ipchg(j).eq.0)) then ipart = 1 ! photon if(ireact.eq.11 .or. ireact.eq.13) then ! flag pi0 gammas from NC if(j.le.2) ubuf(1) = j+50 else if( ireact.eq.4) then ! flag pi0 gammas from CC if(j.le.3) ubuf(1) = j+50 end if end if if ((qmass(j).eq.0.511).and.(ipchg(j).eq.-1)) then ipart = 3 ! electron end if if ((qmass(j).eq.0.511).and.(ipchg(j).eq.1)) then ipart = 2 ! positron end if if ((qmass(j).eq.105.659).and.(ipchg(j).eq.-1)) then ipart = 6 ! muon- end if if ((qmass(j).eq.105.659).and.(ipchg(j).eq.1)) then ipart = 5 ! muon+ end if if ((qmass(j).eq.134.963).and.(ipchg(j).eq.0)) then ipart = 7 ! pi0 end if if ((qmass(j).eq.139.567).and.(ipchg(j).eq.1)) then ipart = 8 ! pi+ end if if ((qmass(j).eq.139.567).and.(ipchg(j).eq.-1)) then ipart = 9 ! pi- end if if ((qmass(j).eq.938.280).and.(ipchg(j).eq.1)) then ipart = 14 ! proton end if if ((qmass(j).eq.939.573).and.(ipchg(j).eq.0)) then ipart = 13 ! neutron end if c write(6,*) ipart, qmass(j),ipchg(j), etot,pmom, ppar if(ipart.eq.0) then goto 20 ! unknown particle type skip end if c c stcik it in c if(ubuf(1).eq.0) then ubuf(1) = ipart end if if(j.eq.1) then plepton = pmom do i=1,3 uvert(i) = vrtx(i) upini(i,j) = ppar(i)/pmom ! remember the z->x end do end if iupar(j) = ipart pupar(j) = pmom call gskine(ppar, ipart, nvtx, ubuf,1,nt) 20 continue end do ! loop over generatoed particles if(nt.eq.0) then write(lout,*) ' No tracks on the stack in this event' end if end if ! ipar.eq.0 return 888 continue write(lout,*) ' No more events in generator file' ievent = ievent -1 idevt = idevt -1 ieorun = 1 ieotri = 1 return 998 continue write(lout,*)' Open file mc_gen.dat error ' stop end