Added some user files

michel-decay/validation/user.f95

+
+
+
+                 !!!!!!!!!!!!!!!!!!!!!
+                     MODULE  USER
+                 !!!!!!!!!!!!!!!!!!!!!
+
+  use functions
+
+  implicit none
+
+!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!
+
+  integer, parameter :: nr_q = 2
+  integer, parameter :: nr_bins = 40
+  real, parameter ::                                                  &
+     min_val(nr_q) = (/ 0., 0. /)
+  real, parameter ::                                                &
+     max_val(nr_q) = (/ 1., 1. /)
+  integer mc
+
+!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!
+
+!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!
+
+    !! ============================================== !!
+    !! DO NOT EVEN THINK ABOUT CHANGING ANYTHING HERE !!
+    !! ============================================== !!
+
+  integer ::  set_zero(nr_q)
+  character (len = 6), dimension(nr_q) :: names
+  character(len=10) :: filenamesuffix
+
+
+!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!
+
+            !! ----------------------------------------- !!
+            !!     There are two versions of binning     !!
+            !!     One for computing   d \sigma/ d Q     !!
+            !!     One for computing  Q d \sigma/ d Q    !!
+            !!  choose by setting the variable bin_kind  !!
+            !! ----------------------------------------- !!
+  real (kind=prec) :: bl_div = 1._prec
+  integer :: bin_flag = 0       !!  0 for standard;  +1 for combined; 
+  integer :: bin_kind = 0       !!  0 for d \sig/dQ; +1 for Q d \sig/dQ;
+
+
+  contains
+
+
+
+
+  SUBROUTINE FIX_MU
+
+  !! ==== Specify the scale mu AND musq==mu**2 ==== !!
+
+  mu = Mm
+
+  musq = mu**2
+
+  END SUBROUTINE FIX_MU
+
+
+
+  SUBROUTINE USERINIT
+  ! This is called without arguments once as soon as McMule
+  ! starts and has read all other configuration, i.e. you can
+  ! access which_piece and flavour. Use this to read any
+  ! further information from the user (like cut configuration
+  ! etc). You do not have to print the hashes - this is
+  ! already taken care of - but you are very much invited to
+  ! include information of what it is you are doing
+  !
+  ! If you are using the cut channel of the menu, you may need to set
+  ! the filenamesuffix variable which is appended to the name of the
+  ! VEGAS file.
+
+  ! Example for reading a cut:
+  ! integer cut
+  ! read*,cut
+  ! write(filenamesuffix,'(I2)') cut
+
+  print*, "This calculates the muon decay in the configuration"
+  print*, "(presumably) used by [hep-ph/0505069] and definitely"
+  print*, "used by our paper."
+
+  read*,mc
+  write(filenamesuffix,'(I1)') mc
+
+  END SUBROUTINE
+
+
+
+
+  FUNCTION QUANT(q1,q2,q3,q4,q5,q6,q7)
+
+  real (kind=prec), intent(in) :: q1(4),q2(4),q3(4),q4(4)
+  real (kind=prec), optional :: q5(4),q6(4),q7(4)
+  real (kind=prec) :: qq5(4),qq6(4),qq7(4)
+  real (kind=prec) :: peL(4),peS(4),peE(4)
+  real (kind=prec) :: quant(nr_q)
+  real (kind=prec) :: ez(4),elec(4),etot,mtot
+  real (kind=prec) :: cthL, cthS, cthE, e_gg
+  integer:: survive5, survive6
+
+  qq5 = 0._prec
+  qq6 = 0._prec
+  qq7 = 0._prec
+
+  if(present(q5)) qq5=q5
+  if(present(q6)) qq6=q6
+  if(present(q7)) qq7=q7
+
+!  pol1 = (/ 0._prec, 0._prec, 0.85_prec, 0._prec /)
+  pol1 = (/ 0._prec, 0._prec, 0._prec, 0._prec /)
+
+  !! ==== keep the line below in any case ==== !!
+
+  set_zero = 1
+  call fix_mu
+
+  if (which_piece .eq. "m2ennee0") then
+    names(1) = 'xh'
+    names(2) = 'xs'
+    if (qq5(4).gt.q2(4)) then
+      quant(1) = 2*qq5(4) / Mm
+      quant(2) = 2* q2(4) / Mm
+    else
+      quant(2) = 2*qq5(4) / Mm
+      quant(1) = 2* q2(4) / Mm
+    endif
+  else
+    names(1) = 'xe'
+    quant(1) = 2* q2(4) / Mm
+     
+    if (mc .eq. 1) then
+        e_gg = 0._prec
+        if(cos_th(q2,qq5)>0.8_prec) e_gg = e_gg + qq5(4)
+        if(cos_th(q2,qq6)>0.8_prec) e_gg = e_gg + qq6(4)
+        if(e_gg>10._prec) set_zero=0
+    endif
+  endif
+
+    ! move variable out of bounds if cut has been applied
+  quant = quant + 2*max_val*(1-set_zero)
+
+  END FUNCTION QUANT
+
+
+
+
+
+
+                 !!!!!!!!!!!!!!!!!!!!!!!
+                     END MODULE  USER
+                 !!!!!!!!!!!!!!!!!!!!!!!
+
+
+

mu-e-scattering/muone/user.f95

+
+
+
+                 !!!!!!!!!!!!!!!!!!!!!
+                     MODULE  USER
+                 !!!!!!!!!!!!!!!!!!!!!
+
+  use functions
+
+  implicit none
+
+!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!
+
+  integer, parameter :: nr_q = 11
+  integer, parameter :: nr_bins = 200
+  integer :: hardcut
+  real, parameter :: &
+     min_val(nr_q) = (/ 0.           ,0.               ,0.510998928    ,0.          ,-142894.  , &
+                        -142894.     ,0.510998928**2   ,0.             ,0.          ,0.        , &
+                        0.           /)
+  real, parameter :: &
+     max_val(nr_q) = (/ pi*5./180    ,pi*5./7/180      ,139819.        ,150.*10**3  ,-1021.48  , &
+                        -1021.48     ,250.**2           ,200.     ,     200.   ,     200.      , &
+                        0.001/)
+
+!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!
+
+!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!
+
+    !! ============================================== !!
+    !! DO NOT EVEN THINK ABOUT CHANGING ANYTHING HERE !!
+    !! ============================================== !!
+
+  integer ::  set_zero(nr_q)
+  character (len = 6), dimension(nr_q) :: names
+  character(len=10) :: filenamesuffix
+
+
+!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!==!
+
+            !! ----------------------------------------- !!
+            !!     There are two versions of binning     !!
+            !!     One for computing   d \sigma/ d Q     !!
+            !!     One for computing  Q d \sigma/ d Q    !!
+            !!  choose by setting the variable bin_kind  !!
+            !! ----------------------------------------- !!
+  real (kind=prec) :: bl_div = 1._prec
+  integer :: bin_flag = 0       !!  0 for standard;  +1 for combined;
+  integer :: bin_kind = 0       !!  0 for d \sig/dQ; +1 for Q d \sig/dQ;
+
+
+  contains
+
+
+
+
+  SUBROUTINE FIX_MU
+
+  !! ==== Specify the scale mu AND musq==mu**2 ==== !!
+
+  mu = Mm
+
+  musq = mu**2
+
+  END SUBROUTINE FIX_MU
+
+
+
+  SUBROUTINE USERINIT
+  ! This is called without arguments once as soon as McMule
+  ! starts and has read all other configuration, i.e. you can
+  ! access which_piece and flavour. Use this to read any
+  ! further information from the user (like cut configuration
+  ! etc). You do not have to print the hashes - this is
+  ! already taken care of - but you are very much invited to
+  ! include information of what it is you are doing
+  !
+  ! If you are using the cut channel of the menu, you may need to set
+  ! the filenamesuffix variable which is appended to the name of the
+  ! VEGAS file.
+
+  print*, "with acocut"
+  print*, "hard cut for r, rr and rv is 10^-(filenamesuffix+1)"
+  print*, "coll cut for r, rr and rv is 10^-(filenamesuffix+3)"
+
+  ! Example for reading a cut:
+  ! integer cut
+   read*, hardcut
+   write(filenamesuffix,'(I1,A6)') hardcut,"acocut"
+
+  END SUBROUTINE
+
+
+
+
+  FUNCTION QUANT(q1,q2,q3,q4,q5,q6,q7)
+
+  real (kind=prec), intent(in) :: q1(4),q2(4),q3(4),q4(4)
+  real (kind=prec), optional :: q5(4),q6(4),q7(4)
+  real (kind=prec) :: qq5(4),qq6(4),q1Rest(4),q2Rest(4),q3Rest(4),q4Rest(4)
+  real (kind=prec) :: quant(nr_q)
+  real (kind=prec) :: mgg2,theta_e,theta_m,tee,tmm,Ee,Em,pTe,pTm,pTem,phi_e,phi_m
+
+  qq5 = 0._prec
+  qq6 = 0._prec
+
+  if(present(q5)) qq5=q5
+  if(present(q6)) qq6=q6
+
+!  pol1 = (/ 0._prec, 0._prec, 0.85_prec, 0._prec /)
+  pol1 = (/ 0._prec, 0._prec, 0._prec, 0._prec /)
+
+  !! ==== keep the line below in any case ==== !!
+  call fix_mu
+
+  q1Rest = boost_rf(q1,q1)
+  q2Rest = boost_rf(q1,q2)
+  q3Rest = boost_rf(q1,q3)
+  q4Rest = boost_rf(q1,q4)
+  theta_e = acos(cos_th(q2Rest,q3Rest))
+  theta_m = acos(cos_th(q2Rest,q4Rest))
+  phi_e = phi(q3Rest)
+  phi_m = phi(q4Rest)
+  tee = sq(q1-q3)
+  tmm = sq(q2-q4)
+  Ee = q3Rest(4)
+  Em = q4Rest(4)
+  mgg2 = sq(q3+qq5+qq6)
+  pTe = PT(q3Rest)
+  pTm = PT(q4Rest)
+  pTem = PT(q3Rest+q4Rest)
+
+  set_zero = 1
+  !energy&angular cut
+  if(Ee<1000..OR.theta_e>100.E-3.OR.theta_m>100.E-3) set_zero = 0
+  !acoplanarity cut
+  if(abs(pi-abs(phi_e-phi_m))>3.5E-3) set_zero = 0
+  names(1) = "thetae"
+  names(2) = "thetam"
+  names(3) = "Ee"
+  names(4) = "Em"
+  names(5) = "tee"
+  names(6) = "tmm"
+  names(7) = "mgg2"
+  names(8) = "pTe"
+  names(9) = "pTm"
+  names(10) = "pTem"
+  names(11) = "thetam_ref"
+
+  quant(1) = theta_e
+  quant(2) = theta_m
+  quant(3) = Ee
+  quant(4) = Em
+  quant(5) = tee
+  quant(6) = tmm
+  quant(7) = mgg2
+  quant(8) = pTe
+  quant(9) = pTm
+  quant(10) = pTem
+  quant(11) = theta_m
+
+  ! move variable out of bounds if cut has been applied
+  quant = quant + 2*max_val*(1-set_zero)
+
+  END FUNCTION QUANT
+
+
+
+
+
+
+                 !!!!!!!!!!!!!!!!!!!!!!!
+                     END MODULE  USER
+                 !!!!!!!!!!!!!!!!!!!!!!!
+
+
+