(file) Return to gukine.f CVS log (file) (dir) Up to [HallC] / geant_gep / src

Diff for /geant_gep/src/gukine.f between version 1.2 and 1.14

version 1.2, 2005/09/22 02:13:57 version 1.14, 2006/01/26 21:52:59
Line 25 
Line 25 
       integer i,j,k,ichoice,ichoice2,iremain       integer i,j,k,ichoice,ichoice2,iremain
       real*8 rotmat,xyz(3),xyznew(3),termang       real*8 rotmat,xyz(3),xyznew(3),termang
       integer ic,ii,jj       integer ic,ii,jj
       integer*4 junk1,ikinsetting        integer*4 junk1,ikinsetting,isectrack
       real*8 xfp,yfp,tthfp,tphfp,pfp,junk2,junk3,thfp,phfp       real*8 xfp,yfp,tthfp,tphfp,pfp,junk2,junk3,thfp,phfp
       real*8 e0(10),eang(10),hang(10),targ_thick(10)       real*8 e0(10),eang(10),hang(10),targ_thick(10)
       real*8 xtgt,ytgt,thtgt,phtgt,ptgt,dptgt       real*8 xtgt,ytgt,thtgt,phtgt,ptgt,dptgt
       real yf,ff,thetaf        real yf,ff,thetaf,dptgt2
       character*80 junkline       character*80 junkline
 c c
 c The following variables are the SIMC variables c The following variables are the SIMC variables
Line 57 
Line 57 
 c      include 'option.h' c      include 'option.h'
 c c
       common/kincom/rotmat(3,3)       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 )       call grndm ( rndm , 3 )
 c c
  111  format(a80)  111  format(a80)
Line 69 
Line 83 
 c         write(*,*)'Getting kinematics setting ...' c         write(*,*)'Getting kinematics setting ...'
          open(unit=1,file='geant_kinematics.dat',type='UNKNOWN')          open(unit=1,file='geant_kinematics.dat',type='UNKNOWN')
          read(1,*)ikinsetting          read(1,*)ikinsetting
            read(1,*)isectrack
          close(unit=1)          close(unit=1)
            if(isectrack.eq.1) then
                   sectrack=.true.
            else
                   sectrack=.false.
            endif
          open(unit=1,file='hdr_gep.dat',status='old')          open(unit=1,file='hdr_gep.dat',status='old')
          do i=1,10          do i=1,10
             read(1,*)e0(i),eang(i),hang(i),targ_thick(i)             read(1,*)e0(i),eang(i),hang(i),targ_thick(i)
Line 92 
Line 112 
      $        thtgt,phtgt,ptgt,dptgt)      $        thtgt,phtgt,ptgt,dptgt)
  
 c c
 c At this point, after using polynomials from Lerose, we should  
 c have x,y,theta,phi, and dpmom at the focal plane.  
 c  
 c      write(*,*)'Back from kincalc ...'  
 c      write(*,*)xtgt,ytgt,thtgt,phtgt,dptgt  
 c      x(1)=tan(phtgt)  
 c      x(2)=ytgt  
 c      x(3)=tan(thtgt)  
 c      x(4)=dptgt  
       y_target=ytgt*100.0  
 c  
 c Here is where we call the SIMC routines to calculate the focal plane  
 c quantities from the target quantities.  
 c  
 c      xfp=xfinal(x,4)*100.0  
 c      phfp=atan(thetaf(x,4))/3.14159265*180.0  
 c      yfp=yf(x,4)*100.0  
 c      thfp=atan(ff(x,4))/3.14159265*180.0  
       p_spec=ptgt       p_spec=ptgt
       th_spec=hrsh_ang  
       dpp=dptgt*100.00  
       x=xtgt  
       y=ytgt  
       z=0.0  
       dxdz=tan(phtgt*3.14159265/180.0)  
       dydz=tan(thtgt*3.14159265/180.0)  
       write(*,*)p_spec,th_spec,dpp  
       write(*,*)x,y,z  
       write(*,*)dxdz,dydz  
 c  
 c Call to SIMC routine  
 c  
       call mc_hms (p_spec, th_spec, dpp, x, y, z, dxdz, dydz,  
      >          x_fp, dx_fp, y_fp, dy_fp, m2,  
      >          ms_flag, wcs_flag, decay_flag, resmult, fry,  
      > ok_spec, pathlen)  
       write(*,*)p_spec,th_spec,dpp  
       write(*,*)x,y,z  
       write(*,*)dxdz,dydz  
       write(6,*)x_fp,dx_fp,y_fp,dy_fp  
       write(*,*)m2,ms_flag,wcs_flag,decay_flag  
       write(*,*)resmult,fry,ok_spec,pathlen  
   
       xfp=x_fp  
       phfp=atan(dx_fp)/3.14159265*180.0  
       yfp=y_fp  
       thfp=atan(dy_fp)/3.14159265*180.0  
   
 c c
       pcentral=ptgt        if(nevent.le.2)write(*,*)'******** mom = ',p_spec,' ******'
       pfp=dptgt  
 c 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 c
 c       write(*,*)'Focal plane quantities ...'        call grndm(rndm,3)
 c       write(*,*)xfp,yfp,thfp,phfp,pcentral,pfp        yfp=rndm(1)*60.00-30.00
         thfp=rndm(2)*80.0-40.0
 c      xfp=xtgt      ! in cm  c      thfp=3.0
 c      yfp=ytgt      ! in cm  c      thfp=rndm(3)*0.00001-0.000005
 c      thfp=thtgt/3.14159265*180.0   ! in degrees  c
 c      phfp=phtgt/3.14159265*180.0   ! in degrees  c      write(6,*)'x,phi,y,theta,ok_spec ='
 c      pcentral=ptgt ! in MeV/c  c      write(6,*)xfp,phfp,yfp,thfp,ok_spec
 c      pfp=dptgt     ! fraction of central momentum  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
 c Finally, we convert this information over to a format that GEANT likes. c Finally, we convert this information over to a format that GEANT likes.
 c c
Line 165 
Line 153 
       tinit=ea-938.2796       tinit=ea-938.2796
       vert(1)=xfp       vert(1)=xfp
       vert(2)=yfp       vert(2)=yfp
       vert(3)=0.0        vert(3)=160.0
       ntbeam  = 0.0       ntbeam  = 0.0
       nttarg  = 0.0       nttarg  = 0.0
       ubuf(1) = 0.0       ubuf(1) = 0.0
Line 310 
Line 298 
 c     what the corresponding theta and phi are for the other arm. c     what the corresponding theta and phi are for the other arm.
  
       if (abs(eang).gt.abs(hang)) then       if (abs(eang).gt.abs(hang)) then
          phie=-.065+rndm(1)*.130  c         phie=-.025+rndm(1)*.050
          thetae=-.030+rndm(2)*.060  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))          escat=mp/(1.0+mp/e0-cos(fg*eang+thetae)*cos(phie))
          pscat=sqrt(e0**2+escat**2-          pscat=sqrt(e0**2+escat**2-
      $        2.0*e0*escat*cos(fg*eang+thetae)*cos(phie))      $        2.0*e0*escat*cos(fg*eang+thetae)*cos(phie))
Line 321 
Line 323 
       else       else
 c Hadron arm defining acceptance c Hadron arm defining acceptance
  1221    call grndm ( rndm, 3)  1221    call grndm ( rndm, 3)
          phie=-.065+rndm(1)*.130           phie=-.080+rndm(1)*.160
          thetae=-.030+rndm(2)*.060          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))          escat=mp/(1.0+mp/e0-cos(fg*eang+thetae)*cos(phie))
          pscat=sqrt(e0**2+escat**2-          pscat=sqrt(e0**2+escat**2-
      $        2.0*e0*escat*cos(fg*eang+thetae)*cos(phie))      $        2.0*e0*escat*cos(fg*eang+thetae)*cos(phie))
Line 333 
Line 337 
       endif       endif
  
       ptgt=pscat       ptgt=pscat
   c       write(*,*)escat,eang,pscat,hang
       dptgt=(pcentral-pscat)/pcentral       dptgt=(pcentral-pscat)/pcentral
 c c
 c Following statement to fill just the high and low dp bins c Following statement to fill just the high and low dp bins

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.14
Analyzer/Replay: Mark Jones, Documents: Stephen Wood
 Powered by
ViewCVS 0.9.2-cvsgraph-1.4.0