diff --git a/simc/dbase.f b/simc/dbase.f
index e107c7e6209322b8194d9e4b10ff79e42e1af8f3..fd841127995a60a4da10a35d857ce08067ff9a27 100644
--- a/simc/dbase.f
+++ b/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)
diff --git a/simc/simc.f b/simc/simc.f
index 8d3832bae51dd5fe6a5a3b6852d7021ff87695f3..8c8f8945cacc34e642736701361256518ed89d2c 100644
--- a/simc/simc.f
+++ b/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
diff --git a/simc/simulate.inc b/simc/simulate.inc
index 75049b0ad36b055730b9aba7a3ed8cc149de9838..87894005dc86430089d84426b16909fbfd869a40 100644
--- a/simc/simulate.inc
+++ b/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,
diff --git a/simc/target.f b/simc/target.f
index 4994dab62fc193099fd69674719ff2a16029c864..2dfb428ce2b5d4e24c386b96df99b276daa86cde 100644
--- a/simc/target.f
+++ b/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