Return to gukine.f CVS log

Up to [HallC] / geant_gep / src

More updates for investigating driftmap issues ...

      subroutine gukine
c
c
c This is a subroutine to generate (x,y,theta,phi,delta-p) at
c the focal plane for proton from the H(e,e'p) reaction. 
c It makes use of a simple kinematics routine, some knowledge
c of the extended target acceptances, and polynomials for the
c transport through the spectrometer from John Lerose.
c
      implicit none
      integer       ikine,itra,istak,ivert,ipart,itrtyp,napart,ipaold
      real          pkine,amass,charge,tlife,vert,pvert
      common/gckine/ikine,pkine(10),itra,istak,ivert,ipart,itrtyp
     +      ,napart(5),amass,charge,tlife,vert(3),pvert(4),ipaold
      integer nubuf
      parameter (nubuf=10)
      real ubuf(nubuf),rndm(3)
      integer ntbeam,nttarg,nvtx
      save nvtx
      real plab(3),ptot,etot
      real*8 ang,angr,phi,phir,psir
      real mp,ea,pa,remain
      real r2,r3
      integer nt,ierri,i0,i1,i2,i3
      integer i,j,k,ichoice,ichoice2,iremain
      real*8 rotmat,xyz(3),xyznew(3),termang
      integer ic,ii,jj
      integer*4 junk1,ikinsetting,isectrack
      real*8 xfp,yfp,tthfp,tphfp,pfp,junk2,junk3,thfp,phfp
      real*8 e0(10),eang(10),hang(10),targ_thick(10)
      real*8 xtgt,ytgt,thtgt,phtgt,ptgt,dptgt
      real yf,ff,thetaf,dptgt2
      character*80 junkline
c
c The following variables are the SIMC variables
c
      real*8 x,y,z                            !(cm)
      real*8 dpp                              !delta p/p (%)
      real*8 dxdz,dydz                        !X,Y slope in spectrometer
      real*8 x_fp,y_fp,dx_fp,dy_fp            !Focal plane values to return
      real*8 p_spec,th_spec                   !spectrometer setting
      real*8 fry                              !vertical position@tgt (+y=down)  
      real*8 pathlen
      logical ms_flag                         !mult. scattering flag
      logical wcs_flag                        !wire chamber smearing flag
      logical decay_flag                      !check for particle decay
      logical ok_spec                         !true if particle makes it
      real*8 resmult,m2
c
      include 'fpp_local.h'
      include 'geant_local.h'
c
c      include 'parameter.h'
c      include 'espace_type.h'
c      include 'detector.h'
c      include 'transport.h'
c      include 'option.h'
c
      common/kincom/rotmat(3,3)
c
      nhu1=0
      nhx1=0
      nhv1=0
      nhu2=0
      nhx2=0
      nhv2=0
      nhu3=0
      nhx3=0
      nhv3=0
      nhu4=0
      nhx4=0
      nhv4=0
c
      call grndm ( rndm , 3 )
c
 111  format(a80)
      iremain=nevent/1000
      remain=nevent/1000.0
c      if(iremain.eq.remain) write(*,*)'nevent =',nevent
      write(*,*)'nevent =',nevent

      if(nevent.eq.0) then
c         write(*,*)'Getting kinematics setting ...'
	 open(unit=1,file='geant_kinematics.dat',type='UNKNOWN')
         read(1,*)ikinsetting
	 read(1,*)isectrack
	 close(unit=1)
	 if(isectrack.eq.1) then
		sectrack=.true.
	 else
		sectrack=.false.
	 endif
         open(unit=1,file='hdr_gep.dat',status='old')
         do i=1,10
            read(1,*)e0(i),eang(i),hang(i),targ_thick(i)
         enddo
         einc=e0(ikinsetting)
         hrse_ang=eang(ikinsetting)
         hrsh_ang=-1.0*hang(ikinsetting)
         trg_thk=targ_thick(ikinsetting)
         close(unit=1)
      endif

c
c Now we have the incident electron energy and the angles of 
c the two spectrometers, as well as the target thickness.
c We need to use a) kinematics routines and b) some knowledge
c of the extended target acceptances to choose the x,y,theta,phi,
c and dpmom at the target.
c
      call kincalc(einc,hrse_ang,hrsh_ang,trg_thk,xtgt,ytgt,
     $        thtgt,phtgt,ptgt,dptgt)

c
      p_spec=ptgt
c
      if(nevent.le.2)write(*,*)'******** mom = ',p_spec,' ******'
c
      call grndm(rndm,3)
      dptgt2=rndm(1)*0.200-0.100
      xfp=rndm(2)*120.00-60.00
      phfp=rndm(3)*80.0-40.0
c      phfp=3.0
c      phfp=rndm(3)*0.00001-0.000005
c
      call grndm(rndm,3)
      yfp=rndm(1)*60.00-30.00
      thfp=rndm(2)*80.0-40.0
c      thfp=3.0
c      thfp=rndm(3)*0.00001-0.000005
c
c      write(6,*)'x,phi,y,theta,ok_spec ='
c      write(6,*)xfp,phfp,yfp,thfp,ok_spec
c
      do ii=1,20
	ntuple_array(ii)=0.0
      enddo
      ntuple_array(11)=real(xfp)
      ntuple_array(12)=real(yfp)
      ntuple_array(13)=real(thfp)
      ntuple_array(14)=real(phfp)
c
      pcentral=ptgt
      pfp=dptgt2 
c    
c Finally, we convert this information over to a format that GEANT likes.
c

 1135 dpmom=pfp
      pfp=pcentral*(1.0+pfp)
      pmom=pfp
      ea=sqrt(pfp**2+938.2796**2)
      tinit=ea-938.2796
      vert(1)=xfp
      vert(2)=yfp
      vert(3)=160.0
      ntbeam  = 0.0
      nttarg  = 0.0
      ubuf(1) = 0.0
      
 1000 continue
c     
      ipart      = 14           ! geant pid  (8=pi+,9,pi-, 14 =p)
      ptot=pfp/1000.
      call grndm ( rndm , 3 )
      angr=thfp*3.14159265/180.0
      phir=phfp*3.14159265/180.0
      psir=datan(dtan(phir)*dcos(angr))

      ang=thfp
      phi=phfp
c
c these are the actual parameters for the track.  
c
      xinit=vert(1)
      yinit=vert(2)
      thinit=angr
      phiinit=phir
      
c      sptransport.l.particle.fp_h.ph=dtan(angr)
c      sptransport.l.particle.fp_h.th=dtan(phir)
c      
c      sptransport.l.particle.fp_h.x=vert(1)/100.0
c      
c      sptransport.l.particle.fp_h.y=vert(2)/100.0
c
c
c next we misalign the track to simulate misalignment of the
c entire space frame with respect to the vdc's
c
c just include the translational offsets to start
c
c	write(*,*)vert(1),vert(2),angr,phir,psir

	xofff=0.0
	yofff=0.0
	thofff=0.0
	phofff=0.0
	psofff=0.0
c
c Define the inverse Euler rotation for (thofff,phiofff,psiofff)
c
c
      rotmat(1,1)=dcos(psofff)*dcos(thofff)+dsin(psofff)*dsin(thofff)
     $     *dsin(phofff)
      rotmat(1,2)=-dcos(phofff)*dsin(thofff)
      rotmat(1,3)=-dsin(psofff)*dcos(thofff)+dcos(psofff)*dsin(thofff)
     $     *dcos(phofff)
      rotmat(2,1)=dcos(psofff)*dsin(thofff)-dsin(psofff)*dcos(thofff)
     $     *dsin(phofff)
      rotmat(2,2)=dcos(phofff)*dcos(thofff)
      rotmat(2,3)=-dsin(psofff)*dsin(thofff)-dcos(psofff)*dcos(thofff)
     $     *dsin(phofff)
      rotmat(3,1)=dsin(psofff)*dcos(phofff)
      rotmat(3,2)=dsin(phofff)
      rotmat(3,3)=dcos(psofff)*dcos(phofff)
c
c
      vert(1)=vert(1)-xofff
      vert(2)=vert(2)-yofff
c     
      xyz(1)=dsin(psir)
      xyz(2)=dcos(psir)*dsin(angr)
      xyz(3)=dcos(psir)*dcos(angr)
c
      do ic=1,3
         xyznew(ic)=xyz(1)*rotmat(ic,1)+
     &        xyz(2)*rotmat(ic,2)+xyz(3)*rotmat(ic,3)
      enddo

      angr=datan(xyznew(2)/xyznew(3))
      termang=xyznew(3)/dcos(angr)
      if(termang.gt.1.0) termang=1.0
      if(termang.lt.-1.0) termang=-1.0
      psir=dacos(termang)*abs(xyznew(1))/xyznew(1)
      phir=datan(dtan(psir)/dcos(angr))
	
c	write(*,*)vert(1),vert(2),angr,phir,psir

      call gsvert ( vert,ntbeam,nttarg,ubuf,0,nvtx )
      etot       = ea - 938.2796
      plab(1)    = ptot*dsin(psir)
      plab(2)    = ptot*dsin(angr)*dcos(psir)
      plab(3)    = ptot*dcos(angr)*dcos(psir)

      call hfill ( 100, etot, 0., 1. )

      call gskine ( plab,ipart,nvtx,ubuf,0,nt )

      if ( nt.le.0 ) then
          write ( 6,* ) ' gukine: error defining track'
          write ( 6,* ) '         i=',i,' nt=',nt
          stop
      end if

      nevent=nevent+1

      return

      end

      subroutine kincalc(e0,eang,hang,trg,xtgt,ytgt,
     $        thtgt,phtgt,ptgt,dptgt)

      implicit none

      real*8 e0,eang,hang,trg,xtgt,ytgt,thtgt,phtgt,ptgt,dptgt
      real*8 fg,gf
      integer i,j,k
      real*8 mt,mtg,mr,mpi,mn,mp,me,pi,mhe,alpha
      real*8 escat,pscat,pcentral,thetae,phie,thetap,phip
      real rndm(3)

      fg=3.14159265/180.0
      gf=1.0/fg

      me  = 0.511
      mpi = 139.57
      mp  = 938.2796
      mn  = 939.5731
      mhe = 2808.41
      mt  = mp
      mtg = mt/1.e3
      
      alpha=1./137.


 1432 call grndm ( rndm , 3 )
      escat=mp/(1.0+mp/e0-cos(fg*eang))
      pscat=sqrt(e0**2+escat**2-
     $     2.0*e0*escat*cos(fg*eang))
      pcentral=pscat

c     First, we assume that whichever are is the more backward will
c     determine the acceptance.  We then randomly choose the theta
c     and phi for the arm that determines the acceptance within the
c     usual full HRS acceptance and then determine from kinematics 
c     what the corresponding theta and phi are for the other arm.

      if (abs(eang).gt.abs(hang)) then
c         phie=-.025+rndm(1)*.050
c         thetae=-.009+rndm(2)*.018
         if(pcentral.ge.5000) then
c           phie=-.00+rndm(1)*.00
c           thetae=-.00+rndm(2)*.00
           phie=-.130+rndm(1)*.260
           thetae=-.065+rndm(2)*.130
	 else if(pcentral.ge.3000.and.pcentral.lt.5000) then
           phie=-.067+rndm(1)*.135
           thetae=-.034+rndm(2)*.067
         else
           phie=-.087+rndm(1)*.174
           thetae=-.044+rndm(2)*.087
	 endif
c         phie=-.00+rndm(1)*.00
c         thetae=-.00+rndm(2)*.00
         escat=mp/(1.0+mp/e0-cos(fg*eang+thetae)*cos(phie))
         pscat=sqrt(e0**2+escat**2-
     $        2.0*e0*escat*cos(fg*eang+thetae)*cos(phie))
         phip=dasin(escat*dsin(phie)/pscat)
         thetap=dasin((escat*sin(fg*eang+thetae)*cos(phie))/
     $        (pscat*cos(phip)))-fg*hang
      else
c Hadron arm defining acceptance
 1221    call grndm ( rndm, 3)
         phie=-.080+rndm(1)*.160
         thetae=-.030+rndm(2)*.060
c         phie=-.00+rndm(1)*.00
c         thetae=-.00+rndm(2)*.00
         escat=mp/(1.0+mp/e0-cos(fg*eang+thetae)*cos(phie))
         pscat=sqrt(e0**2+escat**2-
     $        2.0*e0*escat*cos(fg*eang+thetae)*cos(phie))
         phip=-1.0*dasin(escat*sin(phie)/pscat)
         thetap=dasin((escat*sin(fg*eang+thetae)*cos(phie))/
     $        (pscat*cos(phip)))-fg*hang
         if(abs(thetap).gt.0.030.or.abs(phip).gt.0.065) goto 1221
      endif

      ptgt=pscat
c	write(*,*)escat,eang,pscat,hang
      dptgt=(pcentral-pscat)/pcentral
c
c Following statement to fill just the high and low dp bins
c which have very low statistics.  This is normally commented
c out.
c
c      if(abs(dptgt).le.0.030) goto 1432
c
      thtgt=thetap
      phtgt=phip
      xtgt=0.0
      ytgt=0.0

      return 
      end