Skip to content
Snippets Groups Projects
Commit b39343c0 authored by Sylvester Joosten's avatar Sylvester Joosten
Browse files

disable multiple scattering but not wire chamber smearing for muons

parent 9c6eca3d
No related branches found
No related tags found
No related merge requests found
Pipeline #55894 failed
...@@ -1321,6 +1321,7 @@ c enddo ...@@ -1321,6 +1321,7 @@ c enddo
real*8 eloss_E_arm, eloss_P_arm, r, beta, dangles(2), dang_in(2) real*8 eloss_E_arm, eloss_P_arm, r, beta, dangles(2), dang_in(2)
logical success logical success
logical ok_E_arm, ok_P_arm logical ok_E_arm, ok_P_arm
logical ms_flag, wcs_flag !Mult. scat. and wire chamber smearing
type(event):: orig, recon type(event):: orig, recon
type (event_main):: main type (event_main):: main
...@@ -1352,6 +1353,11 @@ c enddo ...@@ -1352,6 +1353,11 @@ c enddo
frx = 0.0 frx = 0.0
endif endif
! Disable multiple scattering for muons, but not the wire chamber
! smearing
ms_flag = mc_smear .and. (.not. doing_muons)
wcs_flag = mc_smear
!BEAM MULTIPLE SCATTERING: NOT YET IMPLEMENTED, JUST GETTING IT READY! !BEAM MULTIPLE SCATTERING: NOT YET IMPLEMENTED, JUST GETTING IT READY!
! ... multiple scattering of the beam. Generate angles for beam deflection, ! ... multiple scattering of the beam. Generate angles for beam deflection,
...@@ -1359,7 +1365,7 @@ c enddo ...@@ -1359,7 +1365,7 @@ c enddo
! ... The yptar offset goes directly to yptar of both arms (small angle approx). ! ... The yptar offset goes directly to yptar of both arms (small angle approx).
! ... xptar is multiplied by cos(theta) to get the xptar offsets. ! ... xptar is multiplied by cos(theta) to get the xptar offsets.
if (mc_smear) then if (ms_flag) then
call target_musc(orig%Ein,beta_electron,main%target%teff(1),dang_in) call target_musc(orig%Ein,beta_electron,main%target%teff(1),dang_in)
else else
dang_in(1)=0.0 !yp offset, goes directly to yp of both arms dang_in(1)=0.0 !yp offset, goes directly to yp of both arms
...@@ -1386,7 +1392,7 @@ c enddo ...@@ -1386,7 +1392,7 @@ c enddo
! ... multiple scattering ! ... multiple scattering
if (mc_smear) then if (ms_flag) then
beta = orig%p%p/orig%p%E beta = orig%p%p/orig%p%E
call target_musc(orig%p%p, beta, main%target%teff(3), dangles) call target_musc(orig%p%p, beta, main%target%teff(3), dangles)
else else
...@@ -1468,27 +1474,27 @@ C DJG moved this to the last part of generate!!! ...@@ -1468,27 +1474,27 @@ C DJG moved this to the last part of generate!!!
if (hadron_arm.eq.1) then if (hadron_arm.eq.1) then
call mc_hms(spec%p%P, spec%p%theta, delta_P_arm, x_P_arm, call mc_hms(spec%p%P, spec%p%theta, delta_P_arm, x_P_arm,
> y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp, > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp,
> m2, mc_smear, mc_smear, doing_decay, > m2, ms_flag, wcs_flag, doing_decay,
> ntup%resfac, fry, ok_P_arm, pathlen) > ntup%resfac, fry, ok_P_arm, pathlen)
else if (hadron_arm.eq.2) then else if (hadron_arm.eq.2) then
call mc_sos(spec%p%P, spec%p%theta, delta_P_arm, x_P_arm, call mc_sos(spec%p%P, spec%p%theta, delta_P_arm, x_P_arm,
> y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp, > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp,
> m2, mc_smear, mc_smear, doing_decay, > m2, ms_flag, wcs_flag, doing_decay,
> ntup%resfac, fry, ok_P_arm, pathlen) > ntup%resfac, fry, ok_P_arm, pathlen)
else if (hadron_arm.eq.3) then else if (hadron_arm.eq.3) then
call mc_hrsr(spec%p%P, spec%p%theta, delta_P_arm, x_P_arm, call mc_hrsr(spec%p%P, spec%p%theta, delta_P_arm, x_P_arm,
> y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp, > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp,
> m2, mc_smear, mc_smear, doing_decay, > m2, ms_flag, wcs_flag, doing_decay,
> ntup%resfac, fry, ok_P_arm, pathlen) > ntup%resfac, fry, ok_P_arm, pathlen)
else if (hadron_arm.eq.4) then else if (hadron_arm.eq.4) then
call mc_hrsl(spec%p%P, spec%p%theta, delta_P_arm, x_P_arm, call mc_hrsl(spec%p%P, spec%p%theta, delta_P_arm, x_P_arm,
> y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp, > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp,
> m2, mc_smear, mc_smear, doing_decay, > m2, ms_flag, wcs_flag, doing_decay,
> ntup%resfac, fry, ok_P_arm, pathlen) > ntup%resfac, fry, ok_P_arm, pathlen)
else if (hadron_arm.eq.5 .or. hadron_arm.eq.6) then else if (hadron_arm.eq.5 .or. hadron_arm.eq.6) then
call mc_shms(spec%p%P, spec%p%theta, delta_P_arm, x_P_arm, call mc_shms(spec%p%P, spec%p%theta, delta_P_arm, x_P_arm,
> y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp, > y_P_arm, z_P_arm, dx_P_arm, dy_P_arm, xfp, dxfp, yfp, dyfp,
> m2, mc_smear, mc_smear, doing_decay, > m2, ms_flag, wcs_flag, doing_decay,
> ntup%resfac, fry, ok_P_arm, pathlen, hadron_arm, use_first_cer) > ntup%resfac, fry, ok_P_arm, pathlen, hadron_arm, use_first_cer)
endif endif
...@@ -1580,7 +1586,7 @@ C recon%p%delta = (recon%p%P-spec%p%P)/spec%p%P*100. ...@@ -1580,7 +1586,7 @@ C recon%p%delta = (recon%p%P-spec%p%P)/spec%p%P*100.
! ... multiple scattering ! ... multiple scattering
if (mc_smear) then if (mc_smear .and. (.not. doing_muons)) then
call target_musc(orig%e%p, beta_electron, main%target%teff(2), dangles) call target_musc(orig%e%p, beta_electron, main%target%teff(2), dangles)
else else
dangles(1)=0.0 dangles(1)=0.0
...@@ -1647,31 +1653,34 @@ C recon%p%delta = (recon%p%P-spec%p%P)/spec%p%P*100. ...@@ -1647,31 +1653,34 @@ C recon%p%delta = (recon%p%P-spec%p%P)/spec%p%P*100.
! ok flag ! ok flag
m2 = me2 m2 = me2
if (doing_muons) then
m2 = Mmu2
endif
pathlen = 0.0 pathlen = 0.0
if (electron_arm.eq.1) then if (electron_arm.eq.1) then
call mc_hms(spec%e%P, spec%e%theta, delta_E_arm, x_E_arm, call mc_hms(spec%e%P, spec%e%theta, delta_E_arm, x_E_arm,
> y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp, > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp,
> me2, mc_smear, mc_smear, .false., > m2, ms_flag, wcs_flag, .false.,
> tmpfact, fry, ok_E_arm, pathlen) > tmpfact, fry, ok_E_arm, pathlen)
else if (electron_arm.eq.2) then else if (electron_arm.eq.2) then
call mc_sos(spec%e%P, spec%e%theta, delta_E_arm, x_E_arm, call mc_sos(spec%e%P, spec%e%theta, delta_E_arm, x_E_arm,
> y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp, > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp,
> me2, mc_smear, mc_smear, .false., > m2, ms_flag, wcs_flag, .false.,
> tmpfact, fry, ok_E_arm, pathlen) > tmpfact, fry, ok_E_arm, pathlen)
else if (electron_arm.eq.3) then else if (electron_arm.eq.3) then
call mc_hrsr(spec%e%P, spec%e%theta, delta_E_arm, x_E_arm, call mc_hrsr(spec%e%P, spec%e%theta, delta_E_arm, x_E_arm,
> y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp, > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp,
> me2, mc_smear, mc_smear, .false., > m2, ms_flag, wcs_flag, .false.,
> tmpfact, fry, ok_E_arm, pathlen) > tmpfact, fry, ok_E_arm, pathlen)
else if (electron_arm.eq.4) then else if (electron_arm.eq.4) then
call mc_hrsl(spec%e%P, spec%e%theta, delta_E_arm, x_E_arm, call mc_hrsl(spec%e%P, spec%e%theta, delta_E_arm, x_E_arm,
> y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp, > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp,
> me2, mc_smear, mc_smear, .false., > m2, ms_flag, wcs_flag, .false.,
> tmpfact, fry, ok_E_arm, pathlen) > tmpfact, fry, ok_E_arm, pathlen)
else if (electron_arm.eq.5 .or. electron_arm.eq.6) then else if (electron_arm.eq.5 .or. electron_arm.eq.6) then
call mc_shms(spec%e%P, spec%e%theta, delta_E_arm, x_E_arm, call mc_shms(spec%e%P, spec%e%theta, delta_E_arm, x_E_arm,
> y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp, > y_E_arm, z_E_arm, dx_E_arm, dy_E_arm, xfp, dxfp, yfp, dyfp,
> me2, mc_smear, mc_smear, .false., > m2, ms_flag, wcs_flag, .false.,
> tmpfact, fry, ok_E_arm, pathlen, electron_arm, use_first_cer) > tmpfact, fry, ok_E_arm, pathlen, electron_arm, use_first_cer)
else if (electron_arm.eq.7 .or. electron_arm .eq. 8) then else if (electron_arm.eq.7 .or. electron_arm .eq. 8) then
if (abs(spec%p%phi-pi/2) .eq. 10.) then if (abs(spec%p%phi-pi/2) .eq. 10.) then
...@@ -1681,7 +1690,7 @@ C recon%p%delta = (recon%p%P-spec%p%P)/spec%p%P*100. ...@@ -1681,7 +1690,7 @@ C recon%p%delta = (recon%p%P-spec%p%P)/spec%p%P*100.
endif endif
call mc_calo(spec%e%p, spec%e%theta, delta_e_arm, x_e_arm, call mc_calo(spec%e%p, spec%e%theta, delta_e_arm, x_e_arm,
> y_e_arm, z_e_arm, dx_e_arm, dy_e_arm, xfp, dxfp, yfp, dyfp, > y_e_arm, z_e_arm, dx_e_arm, dy_e_arm, xfp, dxfp, yfp, dyfp,
> m2, mc_smear, mc_smear, doing_decay, > m2, ms_flag, wcs_flag, doing_decay,
> ntup%resfac, frx, fry, ok_e_arm, pathlen, using_tgt_field, > ntup%resfac, frx, fry, ok_e_arm, pathlen, using_tgt_field,
> zhadron,electron_arm,drift_to_cal) > zhadron,electron_arm,drift_to_cal)
endif endif
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment