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