cinteraction.f

Go to the documentation of this file.

#include "ZsubstRec.h"
!            treat interaction of MovedTrack
!
      subroutine cinteraction
      use modxsecmedia
      implicit none

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


      integer i
      real*8 Ein1, Ein2, Eout, dEabs1, deabs2, dErel1, dErel2
      real*8  dErel, dEabs, Ein

!
!          used to judge if user hook should be called 
!          after MovedTrack interacted.
!
      integer neverNEP/0/,  neverE/0/,  neverG/0/
      save  nevernep, nevere, neverg
      integer never
      integer stackpos
      save never
      integer icon
      integer,parameter::loopmax=100
      integer:: loopc



!c      record /ptcl/ fragA(maxHeavyMassN),  nonIntNucA(maxHeavyMassN)
!c      integer noOfFrag, noOfNonIntN
!
!      record /track/aTrack
!
!   **  ptcl stacking is done in each subroutine; should be changed
!           (except for hadronic interactions)
!       logic of employing events satisfying some conditions only 
!
!        generate 1 event ; data is in PWork(:)
!        call chookNEPInt(neverNEP)  ! the user may give 5 to neverNEP
!                                  ! to discard the current event
!        never = neverNEP
!        
!           
!

      loopc = 0
      do while (loopc < loopmax )
!          try until desired event is generated
!          (mainly for multiple production)
         loopc = loopc + 1
         nproduced = 0

!
!c    assume MovedTrack is not changed in interaction routine
!     but may be added asflag by chookEing if it is elec.
!      (in the case of Job ='newskel' )
!     
!
         if(movedtrack%p%code .eq. kelec) then
            call cinteelec
         elseif(movedtrack%p%code .eq. kphoton) then
            call cintephoton
         else
            call cintenep
            if(intinfarray(processno)%process .eq. 'coll') then
               movedtrack%pos%colheight = movedtrack%pos%height
            endif
         endif

         if(movedtrack%p%code .eq. kelec) then
            if(nevere .ne. 1) then
               call chookeint(nevere)
               never = nevere
            endif
         elseif(movedtrack%p%code .eq. kphoton) then
            if(neverg  .ne. 1)  then
               call chookgint(neverg)
               never = neverg
            endif
         else
            if(nevernep .ne. 1) then
               call chooknepint(nevernep)
               never = nevernep
            endif
         endif
         if( never /= 5) exit
         never = 0
      enddo   ! end of while;   
      if(btest(eabsorb(1), biteconsv-1) ) then
         if(intinfarray(processno)%process .eq. 'coll' .or.
     *        intinfarray(processno)%process .eq. 'photop' .or.
     *        intinfarray(processno)%process .eq. 'munuci' ) then

!                                    last one not used yet
            call chookeabsorbc( movedtrack, nproduced,  pwork, 0)
!               what is being done below is almost the same as
!               done above call.
            ein1 = movedtrack%p%fm%p(4)  
     *           + masn*(targetnucleonno-targetprotonno) +
     *           masp*targetprotonno
            ein2 = movedtrack%p%fm%p(4) + masp
            eout = 0.
            do i = 1, nproduced
               eout = eout + pwork(i)%fm%p(4)
            enddo

            deabs1 = eout- ein1
            derel1 = eout/ein1 -1.0
            deabs2 = eout- ein2
            derel2 = eout/ein2 -1.0
            if( abs(derel1) .lt. abs(derel2)) then
               derel=derel1
               deabs =deabs1
               ein = ein1
            else
               derel=derel2
               deabs =deabs2
               ein = ein2
!                 no mass case in Eout
               ein1 = movedtrack%p%fm%p(4)
               deabs1 = eout- ein1
               derel1 = eout/ein1 -1.0
               if(abs(deabs1) .lt. abs(deabs) ) then
                  derel = derel1
                  ein = ein1
                  deabs = deabs1
               endif
            endif
            if( abs(derel) .gt.  0.1 .or.
     *           abs(deabs) .gt.  1.e5 ) then
               write(0,*) " code=",movedtrack%p%code
               write(0,*) " chg=",movedtrack%p%charge
               write(0,*) ' Moved E=',movedtrack%p%fm%p(4)
               write(0,*) ' Ein =', ein, ' Eout=',eout
               write(0,*) ' Rerr =', eout/ein -1.0
               write(0,*) ' dEabscol= ',deabs
               write(0,*) 'ActiveModel=', activemdl
            endif
         else
!              possible process; compton, mscat, bscat,
!              anihi, decay, photoe, brems, pair cohs
            ein = movedtrack%p%fm%p(4)
            if(intinfarray(processno)%process .ne. 'decay' .and.
     *           intinfarray(processno)%process .ne. 'brems' .and.
     *           intinfarray(processno)%process .ne. 'pair'  .and.
     *           intinfarray(processno)%process .ne. 'cohs' )  then
               ein = ein + masele
            endif
            if(ein .gt. movedtrack%p%mass) then
               eout = 0.
               do i = 1, nproduced
                  eout = eout + pwork(i)%fm%p(4)
               enddo
               deabs = eout- ein
               derel = eout/ein -1.0
               if( abs(derel) .gt. 0.2 )  then
!            if( abs(dEabs) .gt. 1.e5) then
                  if( abs(deabs) .gt. 1.e5 ) then
                     write(0,*) '****************************'
                  else
                     write(0,*) '----------------------------'
                  endif
                  write(0,*) 'proc=', intinfarray(processno)%process
                  write(0,*) 'code=',movedtrack%p%code, ' charge=',
     *                 movedtrack%p%charge, ' E=',ein
                  write(0,*) 'dEabs= ', deabs, derel, nproduced
                  do i = 1, nproduced
                     write(0,*) i, pwork(i)%code, pwork(i)%fm%p(4)
                  enddo
               endif
            endif
         endif
      endif
!///////////////////////      

      if(onedim .eq. 0) then
!            3 dimensional
!                 stack the leading ptcl  first (to save stack area)
         call cmoveptcl3(movedtrack, pwork, nproduced, nstacked)
      else
         movedtrack%vec = incidentcopy%vec
         call cmoveptcl1(movedtrack, pwork, nproduced, nstacked)
      endif


      if(never .eq. 0 .or. never .eq. 1 ) then
!          user may set never=3 
      elseif(never .eq. 3) then
!          don't follow this and  child; reset stackpos
         call cgetcurrentstackpos(stackpos)
         stackpos=stackpos-nstacked
         call cresetstackpos(stackpos)
      elseif(never .eq. 4) then
!             discard this event generated by the current primary
!             clear stack
         call cinitstack
      else
         call cerrormsg('return value from chookE,G,NEPInt wrong', 1)
         write(0,*)  ' never=', never
         stop
      endif

      end
!          following is remnant of never=2 and ad-hoc model
!          when chookNEPInt is called before push is called.
!      The reason that we put chookNEPInt interface after push
!      is to have easy interface for skeleton making
!        
!      elseif(never .eq. 2) then
!              save only fragments and non interacting nucleons
!         if(MovedTrack.p.code .eq. kgnuc ) then
!c               get fragment and non interacting nuc.
!            call cqHvyIntF(fragA, noOfFrag)
!            call cqHvyIntNIN(nonIntNucA, noOfNonIntN)
!c
!            if(OneDim .eq. 0) then
!                call cmovePtcl3(MovedTrack, fragA, noOfFrag)
!                call cmovePtcl3(MovedTrack, nonIntNucA, noOfNonIntN)
!             else
!                MovedTrack.vec = IncidentCopy.vec
!                call cmovePtcl1(MovedTrack, fragA, noOfFrag)
!                call cmovePtcl1(MovedTrack, nonIntNucA, noOfNonIntN)
!             endif
!          endif
!
!     ************************************
!          move partcles in a given array to stack
!          3 dimensional case.

      subroutine  cmoveptcl3(iTrack, pw, n, npush)
      implicit none
!                put n ptcls in pw into stack.
!          if ThinSampling, 
#include  "Zcode.h"
#include  "Ztrackp.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
!
      integer n
      type(track)::iTrack  ! input. incident ptcl.
      type(ptcl)::pw(n)
      integer npush   ! output. actual number of ptcls put in stack.
                   ! in case of ThinSampling, this may be <= Nproduced. 

      type(track)::aTrack  
      integer i
      integer loc1
      integer nact


      atrack = itrack
      npush = 0
      call cgetcurrentstackpos(loc1)  ! upto loc1 is already filled 

      do i =  n, 1, -1          ! move leading ptcl first
#ifdef SUBSTREC
         atrack%p = pw(i)
#else
         atrack%p%fm%p = pw(i)%fm%p
         atrack%p%mass = pw(i)%mass
         atrack%p%code = pw(i)%code
         atrack%p%subcode = pw(i)%subcode
         atrack%p%charge = pw(i)%charge
#endif

!               reset direction cos and related stuffs
         call cresetdirec(atrack)


#if LABELING == 1
!
!              whennever an interaction occur,  update labelcounter
!             if info>0, clear the timer and infor counters.
!
         if(atrack%info .gt. 0) then
!cc               aTrack.info = 0
!cc               aTrack.t = 0.
         endif
         labelcounter = labelcounter + 1
         atrack%label = labelcounter
#elif LABELING == 2
!            the above simple counter may be replaced by the
!            next sophisticated one. 
!            the same one is in the 1dim mode move routine below.
!
         call ctrickycount(itrack, atrack, pw, i)
!
#endif
         npush = npush + 1
         call cpush(atrack)
      enddo 
      if(thinsampling .and. npush .ge. 2 ) then
!         if(ThinSampling .and. npush .ge. 2  .and.
!     *   IntInfArray(ProcessNo).process .ne. 'photop') then
         call cthinstack(loc1+1, npush, itrack, nact)
!               among npush from loc1+1, nact is accepted
!               adjust npush
         npush = nact
      endif
      end
      subroutine cthinstack(stackloc, n, iTrack, nout)
#include  "Zcode.h"
#include  "Ztrackp.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
#include  "Zstackv.h"
! 
      integer stackloc  ! first loc of stack where tracks of current int.
      integer n
      type(track)::iTrack  ! input. incident ptcl. of the interaction
      integer nout

      call cthinning(stack(stackloc), n, itrack, nout)
      call cresetstackpos(stackloc-1+nout)
      end




#if LABELING == 2
!       ************************************
!             This routine updates label counters.
!          but for the survival particle form brems and knock-on
!          the label counter is not updated. 
!          For those ptcl with info > 0, timer and info counters
!          are cleared. 
      subroutine ctrickycount(iTrack, aTrack, pw, i)
      implicit none
#include  "Zcode.h"
#include  "Ztrackp.h"
#include  "Ztrack.h"
#include  "Ztrackv.h"
!
      integer n
      type(track)::iTrack  ! input. incident ptcl.
      type(track)::aTrack  ! input/output. i-th  secondary track 
      type(ptcl)::pw(*)    ! input.  secondary pool
      integer i              ! input.  i-th secondary (index)

      logical reset

      reset = .true.

      if( itrack%pos%height .gt. 40.d3 ) then
         if(intinfarray(processno)%process .eq. 'brems' ) then
!             for brem, electron has the same label as the incident
            if( atrack%p%fm%p(4) .lt. itrack%p%fm%p(4)*0.8) then
               atrack%label = itrack%label
               reset = .false.
            endif
         elseif(intinfarray(processno)%process .eq. 'knock' .or.
     *        intinfarray(processno)%process .eq. 'mscat' .or.
     *         intinfarray(processno)%process .eq. 'bscat' ) then
!                for the knock-on, survival particle has the same
!                label.
            if(itrack%p%code .ne. kelec  .or.
     *           itrack%p%charge .ne. -1) then
!                     knockon by p,mu,pi..e+ (not by e-)
               if(atrack%p%code .ne. kelec) then
!                       survival one has the same label
                  atrack%label = itrack%label
                  reset = .false.
               endif
            else
!                  electron; make the higher one has the same label
               if(pw(1)%fm%p(4) .gt. pw(2)%fm%p(4)) then
                  if(i .eq. 1) then
                     atrack%label = itrack%label
                     reset = .false.
                  endif
               else
                  if(i .eq. 2) then
                     atrack%label = itrack%label
                     reset = .false.
                  endif
               endif
            endif                     
         endif
      endif
      if(reset) then
         labelcounter = labelcounter +1
         atrack%label = labelcounter
!c         aTrack.t = 0.          ! timer reset
!c         aTrack.info = 0        ! cross counter reset
      endif
      end
#endif
!     ************************************
!          move partcles in a given array to stack
!          1 dimensional case.

      subroutine  cmoveptcl1(iTrack, pw, n, npush)
      implicit none

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


      integer n
      type(track)::iTrack   ! input. incident track
      type(ptcl)::pw(n)
      integer  npush   !  output.  actuall number of ptcls put in stack.
      real*8 temp, p

      integer i
      
      type(track)::aTrack 
      integer loc1
      integer nact

      atrack = itrack

      call cgetcurrentstackpos(loc1)

      do i =  n, 1, -1          ! move leading ptcl last
#ifdef SUBSTREC
         atrack%p = pw(i)
#else
         atrack%p%fm%p = pw(i)%fm%p
         atrack%p%mass = pw(i)%mass
         atrack%p%code = pw(i)%code
         atrack%p%subcode = pw(i)%subcode
         atrack%p%charge = pw(i)%charge
#endif
!            see if angle of particle is larger than a lmit
         call cscalerprod(atrack%p%fm%p, dcatobsxyz, temp)
         call cpxyzp(atrack%p%fm, p)
         if(p .gt. 0.) then
            temp = temp/p
         else
            temp = 1.
         endif
         if(temp .gt. backanglimit)  then
!              only take some limitted angle particles
!            call cresetMom(aTrack)  which is
            atrack%p%fm%p(1) = p * atrack%vec%w%r(1)
            atrack%p%fm%p(2) = p * atrack%vec%w%r(2)
            atrack%p%fm%p(3) = p * atrack%vec%w%r(3)

            call cgetzenith(atrack, atrack%vec%coszenith)


#if LABELING == 1
!                  whennever secondary particles are generated,
!               each of them get an updated label cocunter
!               if the particle has crossed the highest level
!               (info > 0),  timer and info counter is cleared
!
            if( atrack%info .gt. 0) then
!cc                  aTrack.info = 0
!c                  aTrack.t = 0.
            endif
            labelcounter = labelcounter + 1
            atrack%label = labelcounter
#elif LABELING == 2
!                   this may be used if a tricky count is needed
!               in stead of above  counting
            call ctrickycount(itrack, atrack, pw, i)
#endif
            npush = npush + 1  
            call cpush(atrack)
         endif
      enddo
      if(thinsampling .and. npush .gt. 0 ) then
         call cthinstack(loc1+1, npush, itrack, nact)
         npush = nact
      endif
      end

!     ****************************************************************
      subroutine cqinteptcl(ptclA, num)
      implicit none
!          inquire the particle information that made interactions
!        to produce secondary particles.
!        If "MovedTrack" is a heavy,  ptclA will get interacting nucleons
!        otherwise, ptclA will have MovedTrack.p itself.
!
!           
#include  "Zcode.h"
#include  "Ztrack.h"
#include  "Ztrackp.h"
#include  "Ztrackv.h"

      type(ptcl)::ptclA(*)   ! output. interacted particles. max size
                               ! should be maxHeavyMassN (= 56 =Fe)
      integer num              ! output. number of ptcls in ptclA
!
!
!c      if(MovedTrack.p.code .ge. kdeut .and. 
!c     *    MovedTrack.p.code .le. khvymax ) then
      if( movedtrack%p%code .eq. kgnuc ) then
         call cqhvyintin(ptcla, num)
      else
         num = 1
#ifdef SUBSTREC
         ptcla(1) = movedtrack%p
#else
         ptcla(1)%fm%p = movedtrack%p%fm%p
         ptcla(1)%mass = movedtrack%p%mass
         ptcla(1)%code = movedtrack%p%code
         ptcla(1)%subcode = movedtrack%p%subcode
         ptcla(1)%charge = movedtrack%p%charge
#endif
      endif
      end