cintePhoton.f

Go to the documentation of this file.

      subroutine cintephoton
      implicit none
#include  "Ztrack.h"
#include  "Ztrackv.h"
      integer i
      character*70 msg

!///////////////////
!              to see possibility of photon projectile 
!              with Primary is 'gamma' ;    see last
!              part of this file.
!      call cpredpmjet( MovedTrack.p,  Pwork, Nproduced)
!//////////

      if(intinfarray(processno)%process .eq. 'pair') then
         call cpair
      elseif(intinfarray(processno)%process .eq. 'compt') then
         call ccompt
      elseif(intinfarray(processno)%process .eq. 'photoe') then
         call cphotoee  
      elseif(intinfarray(processno)%process .eq. 'photop') then
         call cphotop
      elseif(intinfarray(processno)%process .eq. 'cohs') then
         call ccohs
      elseif(intinfarray(processno)%process .eq. 'mpair') then
         call cmpair
      else
         write(msg, *) ' process for photon; Proc#=',
     *        processno,
     *        intinfarray(processno)%process, ' undef'
         call cerrormsg(msg,  1)
         write(0, '("energy =",g12.4, "where=",i4, " w=",g12.3)') 
     *    movedtrack%p%fm%p(4), movedtrack%where, movedtrack%wgt
         write(0, * ) " coszenith=", movedtrack%vec%coszenith
         write(0, *) ' MoveStat=',movestat, 'No Of inte=',
     *       numberofinte
         do i = 1, numberofinte
            write(0,*) intinfarray(i)%process, ' dt=',
     *         intinfarray(i)%thickness,intinfarray(i)%length
         enddo
         write(0,*) 
         write(0, * ) " (dep,h)B//A==", trackbefmove%pos%depth,
     *            trackbefmove%pos%height,
     *            movedtrack%pos%depth,  movedtrack%pos%height
         stop
      endif
      end
!     *******************
      subroutine cpair
!     *******************
      implicit none

#include  "Zcode.h"
#include  "Zmass.h"
#include  "Ztrackp.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"

!
      real*8 e1, e2, u, eg, cs, sn
      real*8 teta, teta1, teta2, cos1, sin1, cos2, sin2
      integer ica
      real*8 den, cvh2den
      type(track)::aTrack
      type(coord)::dc
      type(coord)::dce
      real*8 temp
!     
      eg = movedtrack%p%fm%p(4)
      if(lpmeffect .and. eg .gt. lpmpairemin) then
!         den = cthick2den(TrackBefMove.pos.depth)
         den = cvh2den(trackbefmove%pos%height)  ! better
         call cpairerglpm(eg, den, e1)
      else
         call cpairenergy(eg, e1)
      endif

      e2=eg - e1
      if( e1 .gt. e2) then
!          store higher energy ptcl later
         temp = e1
         e1= e2
         e2 = temp
      endif
!         now e1 < e2
      
!            assign charge for e1
      call rndc(u)
      if(u .gt. .5) then
         ica=1
      else
         ica=-1
      endif
!     
      atrack = movedtrack
!      if(eg .lt. 100.e-3) then
!          take  pair angle if < 100 MeV
      if(eg .lt. 10) then
!          take  pair angle if < 10GeV
!     call cPairAng(1, eg, e1, teta1) ! teta1 < pi/2
         call cpairang(e1, masele, teta1)
          if(teta1 .lt. 0.03d0) then
             cos1 = 1. - teta1**2/2
             sin1 = teta1
          else
             cos1 = cos(teta1)
             sin1 = sin(teta1)
          endif
!
          sin2 = sin1 * sqrt(  (e1**2-masele**2)/(e2**2-masele**2) )
          if(sin2 .lt. 0.03d0) then
             cos2 = 1.- sin2**2/2
          else
             cos2 = sqrt(1.d0 - sin2**2)
          endif
          call kcossn(cs, sn)

!         teta = masele/eg
!         teta1=teta* e2/eg
!         teta2=teta* e1/eg
!         cos1=1. - teta1**2/2
!         cos2=1. - teta2**2/2
!               sample direction cos. of 1st
!         sin1=teta1
!         sin2=teta2
!
         dc%r(1) = cs * sin1
         dc%r(2) = sn * sin1
         dc%r(3)=  cos1
         call ctransvectz(movedtrack%vec%w, dc, dce)
         call cmkptc(kelec, 0, ica, atrack%p)
         atrack%p%fm%p(4) = e1
         call csetdircos(dce, atrack)
         call ce2p(atrack)
         nproduced = nproduced + 1
         pwork(nproduced) = atrack%p
!            another electron (higher E)

         dc%r(1) = -cs*sin2
         dc%r(2) = -sn*sin2
         dc%r(3) = cos2
         call ctransvectz(movedtrack%vec%w, dc, dce)
         atrack%p%fm%p(4) = e2
         call cmkptc(kelec, 0, -ica,  atrack%p)
         call csetdircos(dce, atrack)
         call ce2p(atrack)
         nproduced = nproduced + 1
         pwork(nproduced) = atrack%p
      else
!          neglect pair angle
         atrack%p%fm%p(4) = e1
         call ce2p(atrack)
         call cmkptc(kelec, 0, ica, atrack%p)
         nproduced = nproduced + 1
         pwork(nproduced) = atrack%p
!     
         atrack%p%fm%p(4) = e2
         call ce2p(atrack)
         call cmkptc(kelec, 0, -ica, atrack%p)
         nproduced = nproduced + 1
         pwork(nproduced) = atrack%p
      endif
      end
!     ***********
      subroutine ccompt
!     ***********
      implicit none
!----      include '../Particle/Zcode.h'
#include  "Zcode.h"
!----      include 'Ztrack.h'
#include  "Ztrack.h"
!----      include 'Ztrackv.h'
#include  "Ztrackv.h"
!
      type(track)::aTrack
      real*8 eg, e1, cs, sn, cosg, cose
      real*8 sine, tmp, sing
      type(coord)::dc
      type(coord)::dce
      type(coord)::dcg
!  
      call ccomptea(movedtrack%p%fm%p(4), eg, e1, cosg, cose)

      call kcossn(cs,sn)
!           electron direction
      tmp=max(1.d0-cose*cose, 0.d0)
      sine=sqrt(tmp)
      dc%r(1)=cs*sine
      dc%r(2)=sn*sine
      dc%r(3)=cose
      call ctransvectz(movedtrack%vec%w, dc, dce)
      atrack = movedtrack
      call cmkptc(kelec, 0, -1, atrack%p)
      atrack%p%fm%p(4) = e1
      call csetdircos(dce, atrack)
      call ce2p(atrack)
      nproduced = nproduced + 1
      pwork(nproduced) = atrack%p
!            gamma dicrection
      tmp=max(1.d0-cosg*cosg, 0.d0)
      sing=sqrt(tmp)
      dc%r(1) = -cs*sing
      dc%r(2) = -sn*sing
      dc%r(3) = cosg
      call ctransvectz(movedtrack%vec%w, dc, dcg)
      atrack%p%fm%p(4) = eg
      call cmkptc(kphoton, kcasg, 0, atrack%p)
      call csetdircos(dcg, atrack)
      call ce2p(atrack)
      nproduced = nproduced + 1
      pwork(nproduced) = atrack%p


      end
!     ***********
      subroutine cphotop
!     ***********
      use modxsecmedia
      implicit none
#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Ztrackp.h"
!
!
      integer ngen
      character(8):: whichcode 
      nproduced = 0
      if( howphotop == 0 ) then  ! will not happen
         ngen = 1
         pwork(nproduced+1) =  movedtrack%p
      elseif( howphotop == 1 ) then
         whichcode = "sofia"
      elseif( howphotop == 2 ) then
         if( movedtrack%p%fm%p(4) < 2.5 )  then
            whichcode ="current"
         else
            whichcode ="sofia"
         endif
      elseif( howphotop == 3 ) then
         if( movedtrack%p%fm%p(4) < 2.5 )  then
            whichcode ="sofia"
         else
            whichcode ="current"
         endif
      elseif( howphotop == 4 ) then
         whichcode = "current"
      else
         write(0,*) 'HowPhotoP =', howphotop, ' invalid '
         stop
      endif

      if( whichcode == "sofia" ) then
         call csofia( targetnucleonno, targetprotonno, 
     *      movedtrack%p, pwork(nproduced+1), ngen)
      elseif( whichcode == "current" ) then
         call cgphad(targetnucleonno, targetprotonno, 
     *        movedtrack%p, pwork(nproduced+1), ngen)
      else
         write(0,*) ' setting mistake of  whichcode=',
     *    whichcode, ' in cphotop'
         stop
      endif
      nproduced = nproduced + ngen
         
      end
!     ****************
      subroutine ccohs
!     ****************
!        coherent scattering
      implicit none
#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
!     ******************
!        since coherent scattering is effective at
!        low energies where angular distribution can be
!        approximated by (1+cos^2) dcos, we simply use
   
      type(track)::aTrack
      type(coord)::w
      real*8 cosg, tmp, cs, sn, sing
!             sample scattering angle from (1+cos^2)dcos
      call ksamprsa(cosg)
      tmp=1.d0-cosg*cosg
      sing = sqrt(tmp)
      call kcossn(cs,sn)
      
      w%r(1) = cs*sing
      w%r(2) = sn*sing
      w%r(2) = cosg
      call ctransvectz(movedtrack%vec%w,  w, w)
      atrack = movedtrack
!        energy unchaged;
      call csetdircos(w, atrack)
      call ce2p(atrack)
      nproduced = nproduced + 1
      pwork(nproduced) = atrack%p
      end
!     ****************
      subroutine cphotoee
!     ****************
!        photo electric effect
      implicit none
#include  "Zcode.h"
#include  "Zmass.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"

!     ******************
      type(track)::aTrack
      type(coord)::w
      real*8  cs, sn, sing

      real*8 eout, eg, rEg, rEe, cost, a, tmp
!          essentially no shell down to 1keV.
      eg = movedtrack%p%fm%p(4)
      eout = masele + eg
      
      atrack = movedtrack
      reg=eg/masele
      ree=(eout-masele) /masele
      if(ree .le. 0.) then
         cost=1.
      else
         a = ( targetatomicn/137.0/137.0 + 2.0*ree + reg**2 )/
     *    (2.0*reg*sqrt(2.0*ree) )
         call ksamppeang(a, cost)
      endif
      
      tmp=1.d0-cost*cost
      sing = sqrt(tmp)
      call kcossn(cs,sn)
      
      w%r(1) = cs*sing
      w%r(2) = sn*sing
      w%r(2) = cost
      call ctransvectz(movedtrack%vec%w,  w, w)
      atrack%p%fm%p(4) = eout
      call cmkptc(kelec, 0, -1, atrack%p)
      call csetdircos(w, atrack)
      call ce2p(atrack)
      nproduced = nproduced + 1
      pwork(nproduced) = atrack%p
      end
!     *******************
      subroutine cmpair
!         magnetic pair creation
!     *******************
      implicit none
#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
!
      real*8 e1, e2, u
      integer ica, nc
      type(track)::aTrack
!       
      call cmpaire(xai, e2, nc)
!           e2 is higher energy fraction; change to real energy
!         store later in  working array, then higher one is
!        stored first in the stack to save the memory.
      e2 = movedtrack%p%fm%p(4) * e2
      e1=movedtrack%p%fm%p(4) - e2
!            assign charge for e1
      call rndc(u)
      if(u .gt. .5) then
         ica=1
      else
         ica=-1
      endif
!     
      atrack = movedtrack
      atrack%p%fm%p(4) = e1
      call ce2p(atrack)
      call cmkptc(kelec, 0, ica, atrack%p)
      nproduced = nproduced + 1
      pwork(nproduced) = atrack%p

      atrack%p%fm%p(4) = e2
      call ce2p(atrack)
      call cmkptc(kelec, 0, -ica, atrack%p)
      nproduced = nproduced + 1
      pwork(nproduced) = atrack%p
      end
      
!/////////////
      subroutine cpredpmjet(pj,  a, ntp)
#include  "Zair.h"
#include "Zptcl.h"
      type(ptcl)::pj
      integer ntp
      type(ptcl)::a(*)
      targetnucleonno=14
      targetprotonno= 7 
      call cdpmjet( pj, targetnucleonno, targetprotonno,
     *   a, ntp)
      write(0,*) ' ntp =',ntp
      stop
      end
!/////////////