Merge branch 'positron-update' into 'master'

Modify to add flag doing_positron

See merge request jpsi007/simc-file-input!3

simc/dbase.f

@@ -190,7 +190,9 @@ C DJG:
 	   doing_deutrho = (nint(targ%A).eq.2)
 	   doing_herho = (nint(targ%A).eq.3)
 	   doing_eep=.false.
-
+        else if (doing_positron) then
+	   Mh = Me
+	   doing_eep=.false.
 	else		!doing_eep if nothing else set.
 	  Mh=Mp
 	  doing_eep = .true.
@@ -302,6 +304,9 @@ C DJG:
 	      sign_hadron=-1.0
 	   endif
 
+	else if(doing_positron) then
+	   targ%Mtar_struck = Mp
+	   targ%Mrec_struck = Mp
 
 ! ... for normal production, Strike p (n), recoil partile is n(p). 
 ! ... for bound final state, use targ.Mrec if it appears to be OK (same A
@@ -411,7 +416,7 @@ C DJG:
 	    targ%angle = 0.0
 	    write(6,*) 'Forcing target angle to zero for cryotarget.'
 	  endif
-	  if (targ%can.ne.1 .and. targ%can.ne.2) stop 'bad targ.can value'
+	  if (.not.(targ%can .ge. 1 .and. targ%can.le.3)) stop 'bad targ.can value'
 	endif
 	if(sin(targ%angle) .gt. 0.85) then
 	  write(6,*) 'BAD targ.angle (0 is perp. to beam, +ve is rotated towards SOS)'
@@ -500,7 +505,7 @@ C DJG:
 	if (nint(targ%Z).eq.1) using_Coulomb=.false.  !no coulomb corr. for Z=1
 
 ! ... target
-
+        write(*,*) " Mh2 = " , Mh2
 	call target_init(using_Eloss, using_Coulomb, spec%e%theta, spec%p%theta,
      >		spec%p%P, Mh2, Ebeam, spec%e%P)
 
@@ -769,10 +774,12 @@ C DJG:
 	   else
 	      stop 'I don''t have ANY idea what (e,e''rho) we''re doing!!!'
 	   endif
-
+	else if (doing_positron) then
+	  write(6,*) ' ****--------  Doing positron  --------****'
 	else if (doing_phsp) then
 	  write(6,*) ' ****--------  PHASE SPACE - NO physics, NO radiation  --------****'
 	else
+	   write(*,*) ' doing_positron ',doing_positron
 	  stop 'I don''t have ANY idea what we''re doing!!!'
 	endif
 
@@ -886,6 +893,7 @@ C DJG:
 	ierr = regparmint('doing_hplus', doing_hplus,1)
 	ierr = regparmint('doing_2pi',doing_2pi,0)
 	ierr = regparmint('doing_rho',doing_rho,0)
+	ierr = regparmint('doing_positron',doing_positron,0)
 	ierr = regparmint('doing_decay',doing_decay,0)
 	ierr = regparmdouble('ctau',ctau,0)
 

simc/simc.f

@@ -793,6 +793,8 @@ c	  write(7,*) 'BP thingie in/out     ',shmsSTOP_BP_in,shmsSTOP_BP_out
 	  endif
 	else if (doing_phsp) then
 	  write(iun,*) '              ****--- PHASE SPACE - NO physics, NO radiation (may not work)---****'
+	else if (doing_positron) then
+	  write(iun,*) '              ****--- doing positron****'
 	else
 	  stop 'I don''t have ANY idea what we''re doing!!!'
 	endif

simc/simulate.inc

@@ -68,6 +68,7 @@
 	logical doing_hyd_elast, doing_deuterium, doing_heavy
 	logical doing_eep, doing_pion, doing_delta, doing_2pi, doing_kaon, doing_rho
 	logical doing_semi, doing_semipi, doing_semika, doing_hplus
+	logical doing_positron
 	integer*4 which_kaon, which_pion
 	logical using_cit_generation, using_Coulomb, using_Eloss
 	logical correct_Eloss, correct_raster,do_fermi
@@ -90,7 +91,7 @@
      >    spect_mode, phsp_mode, base,
      >    extra_dbase_file,tgt_field_file,using_E_arm_montecarlo,using_P_arm_montecarlo,
      >    doing_phsp, using_rad, hard_cuts, doing_hyd_elast, doing_deuterium, doing_heavy,
-     >    doing_eep, doing_pion, doing_delta, doing_2pi, doing_kaon, doing_rho, doing_semi,
+     >    doing_eep, doing_pion, doing_delta, doing_2pi, doing_kaon, doing_rho, doing_semi,doing_positron,
      >    doing_semipi, doing_semika, doing_hplus, which_kaon,
      >    which_pion, using_cit_generation, using_Coulomb, using_Eloss,
      >    correct_Eloss, correct_raster, do_fermi,using_tgt_field,

simc/target.f

@@ -13,6 +13,7 @@
 	real*8 Eloss_target, Eloss_Al,Eloss_air		! energy losses
 	real*8 Eloss_kevlar,Eloss_mylar			! (temporary)
 	real*8 z_can,t,atmp,btmp,ctmp,costmp,th_can	!for the pudding-can target.
+	real*8 ecir,ecor,entec,twall,tcm
 	logical liquid
 
 	s_Al = 0.0
@@ -36,6 +37,8 @@
 	    s_Al = s_Al + 0.0028*inch_cm
 	  else if (targ%can .eq. 2) then	!pudding can (5 mil Al, for now)
 	    s_Al = s_Al + 0.0050*inch_cm
+	  else if (targ%can .eq. 3) then	!cryo2017 10cm
+	     s_Al = s_Al + 0.005*inch_cm
 	  endif
 	endif
 
@@ -63,12 +66,12 @@
 ! .....		X0=28.6cm for Kapton, X0=28.7cm for Mylar) 
 ! ...... ASSUMES liquid targets are 2.65" wide, and have 5.0 mil Al side walls.
 
-! ... For SHMS, use SOS windows for now (just a space filler for now).
+
 
 20	continue
 	if (electron_arm.eq.1) then		!electron is in HMS
-	  s_Al = 0.016*inch_cm
-	  s_air = 15
+	  s_Al = 0.020*inch_cm
+	  s_air = 24.61
  	  s_kevlar = 0.015*inch_cm
 	  s_mylar = 0.005*inch_cm
 	  forward_path = (targ%length/2.-zpos) / abs(cos(theta+targ%angle))
@@ -91,10 +94,16 @@
 	  s_mylar = 0.010*inch_cm
 	  forward_path = (targ%length/2.-zpos) / abs(cos(theta-targ%angle))
 	else if (electron_arm.eq.5 .or. electron_arm.eq.6) then	!SHMS
-	  s_Al = 0.008*inch_cm
-	  s_air = 15
- 	  s_kevlar = 0.005*inch_cm
-	  s_mylar = 0.003*inch_cm
+C Scattering before magnets:  Approximate all scattering as occuring AT TARGET.
+C SHMS
+C  20 mil Al scattering chamber window (X0=8.89cm)
+C  57.27 cm air (X0=30420cm)
+C spectrometer entrance window
+C  10 mil Al s (X0=8.89cm)
+	  s_Al = 0.020*inch_cm + 0.010*inch_cm
+	  s_air = 57.27
+ 	  s_kevlar = 0.0
+	  s_mylar = 0.0
 	  forward_path = (targ%length/2.-zpos) / abs(cos(theta-targ%angle))
 	endif
 	s_target = forward_path
@@ -133,6 +142,26 @@ c	      stop 'z_can > can radius in target.f !!!'
 c	       stop
 	    endif
 	    s_Al = s_Al + 0.0050*inch_cm/abs(sin(target_pi/2 - (theta - th_can)))
+	  else if (targ%can .eq. 3) then	!cry02017 10cm
+            ecir=1.315*2.54           ! endcap inner radius (cm)
+            ecor=1.320*2.54           ! endcap outer radius (cm)
+            entec=targ%length-ecir        ! entrance to end cap (cm)
+            twall = (ecor-ecir)    !  Al wall
+	    
+	    tcm = (targ%length/2. + zpos)
+            if((tcm+ecir/tan(targ%angle)).lt.entec) then  ! e goes through sidewall
+               s_target=ecir/sin(targ%angle)   ! liquid target
+               s_Al=s_Al+twall/sin(targ%angle)            ! wall material
+            else
+               s_target=                              ! e goes throught end cap
+     >     (sqrt(ecir**2-((targ%length-ecir-tcm)*sin(targ%angle))**2)
+     >    +(targ%length-ecir-tcm)*cos(targ%angle)) ! liquid target
+
+              s_Al=   s_Al+                        ! wall
+     >    +(sqrt(ecor**2-((targ%length-ecir-tcm)*sin(targ%angle))**2)
+     >    -sqrt(ecir**2-((targ%length-ecir-tcm)*sin(targ%angle))**2))
+     >    *twall/(ecor-ecir)                   ! & end cap
+        endif
 	  endif
 
 	endif		
@@ -157,8 +186,8 @@ c	       stop
 
 30	continue
 	if (hadron_arm.eq.1) then		!proton in HMS
-	  s_Al = 0.016*inch_cm
-	  s_air = 15
+	  s_Al = 0.020*inch_cm
+	  s_air = 24.61
  	  s_kevlar = 0.015*inch_cm
 	  s_mylar = 0.005*inch_cm
 	  forward_path = (targ%length/2.-zpos) / abs(cos(theta+targ%angle))
@@ -185,10 +214,16 @@ c	       stop
 	  s_mylar = 0.010*inch_cm
 	  forward_path = (targ%length/2.-zpos) / abs(cos(theta-targ%angle))
 	else if (hadron_arm.eq.5 .or. hadron_arm.eq.6) then	!SHMS
-	  s_Al = 0.008*inch_cm
-	  s_air = 15
- 	  s_kevlar = 0.005*inch_cm
-	  s_mylar = 0.003*inch_cm
+C Scattering before magnets:  Approximate all scattering as occuring AT TARGET.
+C SHMS
+C  20 mil Al scattering chamber window (X0=8.89cm)
+C  57.27 cm air (X0=30420cm)
+C spectrometer entrance window
+C  10 mil Al s (X0=8.89cm)
+	  s_Al = 0.020*inch_cm + 0.010*inch_cm
+	  s_air = 57.27
+ 	  s_kevlar = 0.0
+	  s_mylar = 0.0
 	  forward_path = (targ%length/2.-zpos) / abs(cos(theta-targ%angle))
 	endif
 
@@ -202,6 +237,7 @@ c	       stop
 	      s_target = side_path
 	      s_Al = s_Al + 0.005*inch_cm / abs(sin(theta))
 	    endif
+	    s_Al = s_Al + 0.0050*inch_cm/abs(sin(target_pi/2 - (theta - th_can)))
 	  else if (targ%can .eq. 2) then	!pudding can (5 mil Al, for now)
 
 ! this is ugly.  Solve for z position where particle intersects can.  The
@@ -227,6 +263,26 @@ c	      stop 'z_can > can radius in target.f !!!'
 c	      stop
 	    endif
 	    s_Al = s_Al + 0.0050*inch_cm/abs(sin(target_pi/2 - (theta - th_can)))
+	  else if (targ%can .eq. 3) then	!cry02017 10cm
+            ecir=1.315*2.54           ! endcap inner radius (cm)
+            ecor=1.320*2.54           ! endcap outer radius (cm)
+            entec=targ%length-ecir        ! entrance to end cap (cm)
+            twall = (ecor-ecir)    !  Al wall
+	    
+	    tcm = (targ%length/2. + zpos)
+            if((tcm+ecir/tan(targ%angle)).lt.entec) then  ! e goes through sidewall
+               s_target=ecir/sin(targ%angle)   ! liquid target
+               s_Al=s_Al+twall/sin(targ%angle)            ! wall material
+            else
+               s_target=                              ! e goes throught end cap
+     >     (sqrt(ecir**2-((targ%length-ecir-tcm)*sin(targ%angle))**2)
+     >    +(targ%length-ecir-tcm)*cos(targ%angle)) ! liquid target
+
+              s_Al=   s_Al+                        ! wall
+     >    +(sqrt(ecor**2-((targ%length-ecir-tcm)*sin(targ%angle))**2)
+     >    -sqrt(ecir**2-((targ%length-ecir-tcm)*sin(targ%angle))**2))
+     >    *twall/(ecor-ecir)                   ! & end cap
+        endif
 	  endif
 
 	endif
@@ -273,7 +329,7 @@ c	      stop
 	logical	liquid
 
 	real*8 zero
-	parameter (zero=0.0e0)	!real*8 zero for subroutine calls
+	parameter (zero=0.0e0)	!double precision zero for subroutine calls
 
 !Given limiting values for the electron/proton angles, the z-position in the
 !target, and beta for the proton, determine min and max losses in target (and