cEfield0.f

Go to the documentation of this file.

      module modefield
        implicit none
!         next one is only for printing
        namelist /efparam/ howefield, defofr, myef

        type  efield
          sequence
           real(8)::ef(3)=0.
           real(8)::h1=0., h2=0.
           real(8)::r1=0., r2=0.
           real(8)::t1=0., t2=0.
        end type efield
      
        integer,parameter:: nmaxefield=5
!           HowEfield and myEf are namelist inputable
!           parameters
        integer,save:: howefield=0  
               ! 0-->no Efield assumed
               ! 1-->simple Efield defined by myEf
               ! 2-->User supplies Efield subroutine cmyEfield.f

        character(len=1),save:: defofr ='h'  ! define distance to the shower axis
                      !  h--> horizontal distance
                      !  p--> perpendicular to the shower axis.
        type(efield),save:: myef(nmaxefield)
        integer,save:: nefield = 0  ! acutal # of E flied 

        logical,save:: checkh(nmaxefield) = .false.
        logical,save:: checkr(nmaxefield) = .false.
        logical,save:: checkt(nmaxefield) = .false.
        logical,save:: useh=.false.
        logical,save:: user=.false.
        logical,save:: uset=.false.



        real(8),save::efxyz(3,nmaxefield) !copy of
                  ! myEf(:)%Ef(3) but in E-xyz coordinate

      contains
!       Nest two sub are used only at test time
      subroutine cefieldread
      implicit none
      read(*, efparam) 
      end subroutine cefieldread

      subroutine cefieldwrite
      implicit none
        write(0, efparam) 
      end      subroutine cefieldwrite

      subroutine cefieldcheck
      implicit none
      integer:: i

      if(howefield /= 1 ) return  !!!!!
      do i =  nmaxefield, 1, -1   ! last non zero Ef 
         if(dot_product(myef(i)%Ef(:),myef(i)%Ef(:))  > 0 ) then
            nefield = i
            exit
         endif
      enddo
      checkh(1:nefield) =
     *     myef(1:nefield)%H1**2+myef(1:nefield)%H2**2 > 0.
      checkr(1:nefield) =
     *     myef(1:nefield)%R1**2+myef(1:nefield)%R2**2 > 0.
      checkt(1:nefield) =
     *     myef(1:nefield)%T1**2+myef(1:nefield)%T2**2 > 0.
      do i = 1, nefield
         useh = useh .or. checkh(i)
         user = user .or. checkr(i)
         uset = uset .or. checkt(i)
      enddo
      useh = useh .and. .not. uset  ! useT has priority
      
      if( nefield >= 1 ) then
              ! check if height is from low to heigh
         do i = 1, nefield
            if( myef(i)%H1 > myef(i)%H2) then
               write(0,*)  i,'-th Efiled height error'
               write(0,*)
     *              'H1=', myef(i)%H1,' >', myef(i)%H2,'=H2'
               stop
            endif
            if( myef(i)%R1 > myef(i)%R2) then
               write(0,*)  i,'-th Efiled Radius error'
               write(0,*)
     *              'R1=', myef(i)%R1, ' >',  myef(i)%R2,'=R2'
               stop
            endif
            if( myef(i)%T1 > myef(i)%T2) then
               write(0,*)  i,'-th Efiled Time error'
               write(0,*)
     *              'T1=', myef(i)%T1, ' >',  myef(i)%T2,'=T2'
               stop
            endif
            if( nefield > i ) then
               if( myef(i+1)%H1 < myef(i)%H2 ) then
                  write(0,*)
     *                 i,'-th Ef height region<=',myef(i)%H2
                  write(0,*)
     *                 ' but is not < next region=',myef(i+1)%H1
                  stop
               endif
            endif                  
         enddo
      endif
      write(0,*) ' # of Efield(s)=',nefield
      end subroutine cefieldcheck


      subroutine cefielddet2xyz
      ! convert Ef in det system into E-xyz system
      ! should be called once before starting tracking
      implicit none
      integer::i
      do i = 1, nefield
         call cdet2xyzd(myef(i)%Ef(:), efxyz(:,i))
      enddo
      end  subroutine cefielddet2xyz

      subroutine cefield0(aTrack, Efout)
!      Electric field is given depending on particle position
!      and/or time.  
!      H          T                      R
!    given   not given --> H is used   R at that H is used. if R is out 
!                                      of range, E=0. else E at H
!                                      if R is not given  E at H
!  not given    given  --> T is used   same as above, use T instead of H.
!    given      given  --> T is used    //
!  not gvien   not given               use first R,.if it exists
!                                      if in range, use E(1) else E=0
!                                      if R absent, use E(1)
      implicit none
#include  "Zglobalc.h"
#include "Ztrack.h"
#include "Zincidentv.h"
      type(track):: atrack  ! input. Track info in E-xyz.
         ! position and "where" info are used.
      real(8),intent(out):: Efout(3) ! Elecric field
                       ! in E-xyz system.  V/m.
      integer:: i, hindex, tindex, out
      logical,save:: first=.true.
      real(8)::time
      
      real(8)::r
      real(8)::Zero(3)=(/0.,0.,0./)

      if( first ) then   ! convert Ef in det into E-xyz
         call cefieldcheck
         call cefielddet2xyz
         first = .false.
!!!!!!!
!         call cEfieldWrite
!!!!!!!!!
      endif

      hindex = 0
      if(useh) then
         do i = 1, nefield
            if(checkh(i)) then
               if( atrack%pos%height >= myef(i)%H1 .and.
     *              atrack%pos%height <= myef(i)%H2 ) then
                  hindex = i
                  exit
               endif
            endif
         enddo
         if( hindex == 0) then  ! H is out of range
            efout(:) = zero(:)
         else
              !  see if R is in the range
            call cefieldcheckr(atrack, hindex, out)
            if(out == 0 ) then
             !  Efout(:) = myEf(hindex)%Ef(:)  ! bug
               efout(:) = efxyz(:,hindex)
            else
               efout(:) = zero(:)
            endif
         endif
      elseif(uset) then
         tindex = 0
         do i = 1, nefield
            if(checkt(i)) then
               time = atrack%t/c*tonsec
               if( time >= myef(i)%T1 .and.
     *             time <= myef(i)%T2 ) then
                  tindex = i
                  exit
               endif
            endif
         enddo
         if( tindex == 0 ) then  ! T is out of range
            efout(:) = zero(:)
         else
              !  see if R is in the range
            call cefieldcheckr(atrack, tindex, out)
            if(out == 0 ) then
!               Efout(:) = myEf(tindex)%Ef(:)  bug
               efout(:) = efxyz(:,tindex)
            else
               efout(:) = zero(:)
            endif
         endif
      else  ! H,T is not used. 
         i = 1   ! use first data if possible
         if( checkr(i) ) then   ! check first R if any
            call cefieldcheckr(atrack, i, out)
            if(out == 0 ) then
!               Efout(:) = myEf(i)%Ef(:)
               efout(:) = efxyz(:,i)
            else
               efout(:) = zero(:)
            endif
         else  ! no check of R.  assume first Ef is applied
!            Efout(:) = myEf(i)%Ef(:)
            efout(:) = efxyz(:,i)
         endif
      endif
      end   subroutine cefield0
      
      subroutine cefieldcheckr(aTrack, i, out)
      implicit none
#include "Ztrack.h"
      type(track):: atrack  ! input
      integer,intent(in):: i ! R or H index
      integer,intent(out)::out ! 0 aTrack is in the range or no check is needed
                   !  else out of range
      real(8):: r

      if( checkr(i) ) then
                ! core distance
         if(defofr == 'p' ) then ! perpedicular distance
            call cgetpcoredist(atrack, r)
         elseif(defofr == 'h') then ! horizontal distance
            call cgethcoredist(atrack, r)
         else
            write(0,*) ' DefofR =',defofr, ' invalid'
            stop
         endif
         if( r >= myef(i)%R1 .and.
     *        r <= myef(i)%R2 ) then
            out = 0
         else                   ! R is out of range
            out = 1
         endif
      else
         out = 0
      endif
      end subroutine cefieldcheckr


      end module modefield

      subroutine cgetpcoredist(aTrack, r)
      implicit none
#include "Ztrack.h"
#include "Zincidentv.h"
#include  "Zobs.h"
#include "Zobsv.h"
      type(track):: atrack  ! input. in E-xyz system
      real(8),intent(out):: r ! Perpendicular core distance

      real(8):: rdet(3)  !   ptcl pos in det-sys

      call cxyz2det(obssites(atrack%where)%pos%xyz,
     *              atrack%pos%xyz%r(:), rdet(:))
      r = abs(dot_product(rdet(:), angleatobscopy%r(:)))  ! Angle..  is dwn.g
      end      subroutine cgetpcoredist

      subroutine cgethcoredist(aTrack, r)
      implicit none
#include "Ztrack.h"
#include "Zincidentv.h"
#include "Zobs.h"
#include "Zobsv.h"
      type(track):: atrack  ! input. in E-xyz system
      real(8),intent(out):: r ! Perpendicular core distance

      real(8):: rdet(3)  !   ptcl pos in det-sys
      real(8):: Raxis(3) ! shower core at the same height
             !         as the ptcl


      call cxyz2det(obssites(atrack%where)%pos%xyz,
     *               atrack%pos%xyz%r(:), rdet(:))
      raxis(:) = 
     *  (atrack%pos%height - obssites(atrack%where)%pos%height) /
     *  angleatobscopy%r(3)     ! length from origin( <0 )
     *  * angleatobscopy%r(:)   ! vector (now upgoing) since
                       ! AngleAtObsCopy is downogin.
      r = sqrt(sum( ( rdet(:)-raxis(:))**2 ))
      end      subroutine cgethcoredist