cinteraction.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	cinteraction
subroutine	cmoveptcl3 (iTrack, pw, n, npush)
subroutine	cthinstack (stackloc, n, iTrack, nout)
subroutine	ctrickycount (iTrack, aTrack, pw, i)
subroutine	cmoveptcl1 (iTrack, pw, n, npush)
subroutine	cqinteptcl (ptclA, num)

Function/Subroutine Documentation

cinteraction()

subroutine cinteraction

(

)

Definition at line 5 of file cinteraction.f.

References cerrormsg(), cgetcurrentstackpos(), chookeabsorbc(), chookeint(), chookgint(), chooknepint(), cinitstack(), cinteelec(), cintenep(), cintephoton(), cmoveptcl1(), cmoveptcl3(), cresetstackpos(), kelec, kphoton, masele, masn, and masp.

Referenced by ctracking().

      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

Here is the call graph for this function:

Here is the caller graph for this function:

cmoveptcl1()

subroutine cmoveptcl1	(	type(track)	iTrack,
		type(ptcl), dimension(n)	pw,
		integer	n,
		integer	npush
	)

Definition at line 400 of file cinteraction.f.

References cgetcurrentstackpos(), cgetzenith(), cpush(), cpxyzp(), cscalerprod(), cthinstack(), and ctrickycount().

Referenced by cinteraction().

      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

Here is the call graph for this function:

Here is the caller graph for this function:

cmoveptcl3()

subroutine cmoveptcl3	(	type(track)	iTrack,
		type(ptcl), dimension(n)	pw,
		integer	n,
		integer	npush
	)

Definition at line 238 of file cinteraction.f.

References cgetcurrentstackpos(), cpush(), cresetdirec(), cthinstack(), and ctrickycount().

Referenced by cinteraction().

      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

Here is the call graph for this function:

Here is the caller graph for this function:

cqinteptcl()

subroutine cqinteptcl	(	type(ptcl), dimension(*)	ptclA,
		integer	num
	)

Definition at line 482 of file cinteraction.f.

References kgnuc.

      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

cthinstack()

subroutine cthinstack	(	integer	stackloc,
		integer	n,
		type(track)	iTrack,
		integer	nout
	)

Definition at line 309 of file cinteraction.f.

References cresetstackpos(), and cthinning().

Referenced by cmoveptcl1(), and cmoveptcl3().

#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)

Here is the call graph for this function:

Here is the caller graph for this function:

ctrickycount()

subroutine ctrickycount	(	type(track)	iTrack,
		type(track)	aTrack,
		type(ptcl), dimension(*)	pw,
		integer	i
	)

Definition at line 335 of file cinteraction.f.

References d3, false, kelec, and true.

Referenced by cmoveptcl1(), and cmoveptcl3().

      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

Here is the caller graph for this function: