Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F11222276
S95tables.f90
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
75 KB
Subscribers
None
S95tables.f90
View Options
! This file is part of HiggsBounds
! -KW
!******************************************************************
module
S95tables
!******************************************************************
use
S95tables_type1
use
S95tables_type2
use
usefulbits
,
only
:
Hneut
,
Hplus
implicit none
private
public
::
nLEPtable1
,
nLEPtable2
,
nTEVtable1
,
nTEVtable2
,
&
&
S95_t1
,
S95_t2
,
&
&
setup_S95tables
,
deallocate_S95tables
,
&
&
calcfact
,
outputproc
,
&
&
check_against_bound
,
&
&
convolve_chisq_with_gaussian
,
&
&
S95_t1_or_S95_t2_idfromelementnumber
,
&
&
f_from_t2
,
f_from_t3
,
f_from_slices_t2
integer
::
nLEPtable1
,
nLEPtable2
,
nTEVtable1
,
nTEVtable2
!-----------------------------------------------------------------------------------
! delta_M*_* determine how close in mass particles should be before their masses are combined
! nb. delta*_* is only used for tables of type 1 at the moment
! for some analyses, S95_t1(x)%deltax has already been set in S95tables_type1.f90
! for analyses where S95_t1(x)%deltax is *not* has already set:
! we set
! LEP neutral Higgs tables to have deltax=delta_Mh_LEP
! Tevatron neutral Higgs tables to have deltax=delta_Mh_TEV
! LEP charged Higgs tables to have deltax=delta_Mhplus_LEP
! Tevatron charged Higgs tables to have deltax=delta_Mhplus_TEV
! other LEP tables to have deltax=delta_M_LEP_default
! other Tevatron tables to have deltax=delta_M_TEV_default
! setting delta_M*_TEV and/or delta_M*_LEP to zero turns this feature off
! DO NOT CHANGE THESE VALUES BEFORE READING THE MANUAL
! even when it is appropriate to add the cross sections, we would recommend using delta_Mh_LEP<=2.0D0 or delta_Mh_TEV<=10.0
double precision
,
parameter
::
delta_Mh_LEP
=
0.0D0
double precision
,
parameter
::
delta_Mh_TEV
=
0.0D0
double precision
,
parameter
::
delta_Mhplus_LEP
=
0.0D0
double precision
,
parameter
::
delta_Mhplus_TEV
=
0.0D0
double precision
,
parameter
::
delta_M_LEP_default
=
0.0D0
!where delta_M_LEP/TEV is not specified below,
double precision
,
parameter
::
delta_M_TEV_default
=
0.0D0
!these values will be used
!double precision, parameter :: delta_Mh_LEP=15.0D0 !crazy values - for debugging only
!double precision, parameter :: delta_Mh_TEV=15.0D0 !crazy(ish) values - for debugging only
!-----------------------------------------------------------------------------------
! eps determines how strict the Standard Model test is
double precision
,
parameter
::
eps
=
2.0D-2
!double precision, parameter :: eps=1.0D3 !crazy value - for debugging only
!-----------------------------------------------------------------------------------
!table type 1-----------------------------
type
(
table1
),
allocatable
::
S95_t1
(:)
!------------------------------------------
!table type 2------------------------------
type
(
table2
),
allocatable
::
S95_t2
(:)
!------------------------------------------
contains
!**********************************************************
subroutine
setup_S95tables
! Allocates and calls subroutines to fill S95_t1, S95_t2
! (which store the experimental data)
! Sets delta_M_TEV,delta_M_LEP which govern how close Higgs
! need to be in mass before HiggsBounds combines their cross sections
!**********************************************************
use
usefulbits
,
only
:
debug
,
&
&
TEVt1Mhmax
,
TEVt1Mhmin
,
LEPt12Mhmax
,
LEPt12Mhmin
,
&
&
TEVrange_Mhmin
,
TEVrange_Mhmax
,
LEPrange_Mhmin
,
LEPrange_Mhmax
use
theory_BRfunctions
,
only
:
BRSMt1Mhmax
,
BRSMt1Mhmin
use
theory_XS_SM_functions
,
only
:
tevXS_SM_functions_xmin
,
tevXS_SM_functions_xmax
implicit none
!-----------------------------------internal
integer
::
i
double precision
::
max_tev_delta_Mh
,
max_lep_delta_Mh
!-------------------------------------------
! these numbers have to be changed appropriately every time a table is added
! or taken away:
nLEPtable1
=
20
! table type 1 involves 1 mass
nTEVtable1
=
63
nLEPtable2
=
14
! table type 2 involves 2 masses
nTEVtable2
=
0
allocate
(
S95_t1
(
nLEPtable1
+
nTEVtable1
))
allocate
(
S95_t2
(
nLEPtable2
+
nTEVtable2
))
call
initializetables_type1_blank
(
S95_t1
)
call
initializetables_type2_blank
(
S95_t2
)
call
initializetables1
(
S95_t1
)
call
initializetables2
(
S95_t2
)
! checks that the numbers of LEP and TEV tables set in S95_t1
! corresponds to the number set here
if
(
(
count
(
S95_t1
%
expt
.
eq
.
'LEP'
).
ne
.
nLEPtable1
)
&
.
or
.(
count
(
S95_t1
%
expt
.
eq
.
'CDF'
)
+
count
(
S95_t1
%
expt
.
eq
.
' D0'
)
+
count
(
S95_t1
%
expt
.
eq
.
'TCB'
).
ne
.
nTEVtable1
))
then
stop
'error in setup_S95tables (a)'
endif
! checks that the numbers of LEP and TEV tables set in S95_t2
! corresponds to the number set here
if
(
(
count
(
S95_t2
%
expt
.
eq
.
'LEP'
).
ne
.
nLEPtable2
)
&
.
or
.(
count
(
S95_t2
%
expt
.
eq
.
'CDF'
)
+
count
(
S95_t2
%
expt
.
eq
.
' D0'
)
+
count
(
S95_t2
%
expt
.
eq
.
'TCB'
).
ne
.
nTEVtable2
))
then
stop
'error in setup_S95tables (b)'
endif
! checks that none of the id's are repeated in S95_t1
do
i
=
lbound
(
S95_t1
,
dim
=
1
),
ubound
(
S95_t1
,
dim
=
1
)
if
(
count
(
S95_t1
%
id
.
eq
.
S95_t1
(
i
)%
id
)
&
+
count
(
S95_t2
%
id
.
eq
.
S95_t1
(
i
)%
id
).
ne
.
1
)
then
stop
'error in setup_S95tables (c)'
endif
enddo
! checks that none of the id's are repeated in S95_t2
do
i
=
lbound
(
S95_t2
,
dim
=
1
),
ubound
(
S95_t2
,
dim
=
1
)
if
(
count
(
S95_t1
%
id
.
eq
.
S95_t2
(
i
)%
id
)&
+
count
(
S95_t2
%
id
.
eq
.
S95_t2
(
i
)%
id
).
ne
.
1
)
then
stop
'error in setup_S95tables (d)'
endif
enddo
! looks for the min and max values of Mh (neutral Higgs) in the LEP tables
! will be used in input.f90 to work out which SM para need to be calculated
! initial (impossible) values
LEPt12Mhmin
=
1.0D6
LEPt12Mhmax
=
0.0D0
do
i
=
lbound
(
S95_t1
,
dim
=
1
),
ubound
(
S95_t1
,
dim
=
1
)
if
((
S95_t1
(
i
)%
expt
.
eq
.
'LEP'
).
and
.(
S95_t1
(
i
)%
particle_x
.
eq
.
Hneut
))
then
if
(
S95_t1
(
i
)%
xmax
.
gt
.
LEPt12Mhmax
)
LEPt12Mhmax
=
S95_t1
(
i
)%
xmax
if
(
S95_t1
(
i
)%
xmin
.
lt
.
LEPt12Mhmin
)
LEPt12Mhmin
=
S95_t1
(
i
)%
xmin
endif
enddo
do
i
=
lbound
(
S95_t2
,
dim
=
1
),
ubound
(
S95_t2
,
dim
=
1
)
if
((
S95_t2
(
i
)%
expt
.
eq
.
'LEP'
).
and
.(
S95_t2
(
i
)%
particle_x1
.
eq
.
Hneut
)
)
then
if
(
S95_t2
(
i
)%
xmin1
.
lt
.
LEPt12Mhmin
)
LEPt12Mhmin
=
S95_t2
(
i
)%
xmin1
if
(
S95_t2
(
i
)%
xmax1
.
gt
.
LEPt12Mhmax
)
LEPt12Mhmax
=
S95_t2
(
i
)%
xmax1
endif
if
((
S95_t2
(
i
)%
expt
.
eq
.
'LEP'
).
and
.(
S95_t2
(
i
)%
particle_x2
.
eq
.
Hneut
)
)
then
if
(
S95_t2
(
i
)%
xmin2
.
lt
.
LEPt12Mhmin
)
LEPt12Mhmin
=
S95_t2
(
i
)%
xmin2
if
(
S95_t2
(
i
)%
xmax2
.
gt
.
LEPt12Mhmax
)
LEPt12Mhmax
=
S95_t2
(
i
)%
xmax2
endif
enddo
! looks for the min and max values of Mh (neutral Higgs) in the Tevatron tables
! will be used in theo_SM.f90 to work out which SM para need to be calculated
! initial (impossible) values
TEVt1Mhmin
=
1.0D6
TEVt1Mhmax
=
0.0D0
do
i
=
lbound
(
S95_t1
,
dim
=
1
),
ubound
(
S95_t1
,
dim
=
1
)
if
(
(
S95_t1
(
i
)%
expt
.
eq
.
'CDF'
)
&
.
or
.(
S95_t1
(
i
)%
expt
.
eq
.
' D0'
)
&
.
or
.(
S95_t1
(
i
)%
expt
.
eq
.
'TCB'
))
then
if
(
S95_t1
(
i
)%
particle_x
.
eq
.
Hneut
)
then
if
(
S95_t1
(
i
)%
xmax
.
gt
.
TEVt1Mhmax
)
TEVt1Mhmax
=
S95_t1
(
i
)%
xmax
if
(
S95_t1
(
i
)%
xmin
.
lt
.
TEVt1Mhmin
)
TEVt1Mhmin
=
S95_t1
(
i
)%
xmin
endif
endif
enddo
if
(
delta_M_LEP_default
.
gt
.
2.1d0
)
write
(
*
,
*
)
'WARNING: delta_M_LEP_default.gt.2.1d0'
if
(
delta_M_TEV_default
.
gt
.
1
0.1d0
)
write
(
*
,
*
)
'WARNING: delta_M_TEV_default.gt.10.1d0'
if
(
delta_Mh_LEP
.
gt
.
2.1d0
)
write
(
*
,
*
)
'WARNING: delta_Mh_LEP.gt.2.1d0'
if
(
delta_Mh_TEV
.
gt
.
1
0.1d0
)
write
(
*
,
*
)
'WARNING: delta_Mh_TEV.gt.10.1d0'
if
(
delta_Mhplus_LEP
.
gt
.
2.1d0
)
write
(
*
,
*
)
'WARNING: delta_Mhplus_LEP.gt.2.1d0'
if
(
delta_Mhplus_TEV
.
gt
.
1
0.1d0
)
write
(
*
,
*
)
'WARNING: delta_Mhplus_TEV.gt.10.1d0'
if
(
eps
.
gt
.
5.0d-1
)
write
(
*
,
*
)
'WARNING: eps.gt.5.0d-1'
! here, we set
! all LEP neutral Higgs tables to be have the same deltax (delta_Mh_LEP)
! all Tevatron neutral Higgs tables to be have the same deltax (delta_Mh_TEV)
! all LEP charged Higgs tables to be have the same deltax (delta_Mhplus_LEP)
! all Tevatron charged Higgs tables to be have the same deltax (delta_Mhplus_TEV)
! all other LEP tables to be have the same deltax (delta_M_LEP_default)
! all other Tevatron tables to be have the same deltax (delta_M_TEV_default)
! unless they have already even set to a value > -0.5D0 (in subroutine initializetables1)
do
i
=
lbound
(
S95_t1
,
dim
=
1
),
ubound
(
S95_t1
,
dim
=
1
)
if
(
S95_t1
(
i
)%
deltax
.
lt
.
-
0.5D0
)
then
!deltax has not been set yet
if
(
(
S95_t1
(
i
)%
expt
.
eq
.
'CDF'
)
&
.
or
.(
S95_t1
(
i
)%
expt
.
eq
.
' D0'
)
&
.
or
.(
S95_t1
(
i
)%
expt
.
eq
.
'TCB'
))
then
if
(
S95_t1
(
i
)%
particle_x
.
eq
.
Hneut
)
then
S95_t1
(
i
)%
deltax
=
delta_Mh_TEV
elseif
(
S95_t1
(
i
)%
particle_x
.
eq
.
Hplus
)
then
S95_t1
(
i
)%
deltax
=
delta_Mhplus_TEV
else
S95_t1
(
i
)%
deltax
=
delta_M_TEV_default
endif
else
if
(
S95_t1
(
i
)%
particle_x
.
eq
.
Hneut
)
then
S95_t1
(
i
)%
deltax
=
delta_Mh_LEP
elseif
(
S95_t1
(
i
)%
particle_x
.
eq
.
Hplus
)
then
S95_t1
(
i
)%
deltax
=
delta_Mhplus_LEP
else
S95_t1
(
i
)%
deltax
=
delta_M_LEP_default
endif
endif
endif
enddo
do
i
=
lbound
(
S95_t2
,
dim
=
1
),
ubound
(
S95_t2
,
dim
=
1
)
S95_t2
(
i
)%
deltax
=
1.0D-7
!don't change this - the relevent part of code was taken out to improve efficiency
enddo
! finds the maximum delta_Mh for the Tevatron tables
! will be used in theo_SM.f90 to work out which SM para need to be calculated
max_tev_delta_Mh
=
-
1.0D0
do
i
=
lbound
(
S95_t1
,
dim
=
1
),
ubound
(
S95_t1
,
dim
=
1
)
if
(
(
S95_t1
(
i
)%
expt
.
eq
.
'CDF'
)
&
.
or
.(
S95_t1
(
i
)%
expt
.
eq
.
' D0'
)
&
.
or
.(
S95_t1
(
i
)%
expt
.
eq
.
'TCB'
))
then
if
(
S95_t1
(
i
)%
particle_x
.
eq
.
Hneut
)
then
if
(
S95_t1
(
i
)%
deltax
.
gt
.
max_tev_delta_Mh
)
then
max_tev_delta_Mh
=
S95_t1
(
i
)%
deltax
endif
endif
endif
enddo
!do i=lbound(S95_t2,dim=1),ubound(S95_t2,dim=1) !not needed until code is adapted to allow a nonzero deltax for type 2 tables
! if( (S95_t2(i)%expt.eq.'CDF') &
! .or.(S95_t2(i)%expt.eq.' D0') &
! .or.(S95_t2(i)%expt.eq.'TCB'))then
! if(S95_t2(i)%deltax.gt.max_tev_delta_Mh)then
! max_tev_delta_Mh=S95_t2(i)%deltax
! endif
! endif
!enddo
! finds the maximum delta_Mh for the LEP tables
! will be used in theo_SM.f90 to work out which SM para need to be calculated
max_lep_delta_Mh
=
-
1.0D0
do
i
=
lbound
(
S95_t1
,
dim
=
1
),
ubound
(
S95_t1
,
dim
=
1
)
if
(
S95_t1
(
i
)%
expt
.
eq
.
'LEP'
)
then
if
(
S95_t1
(
i
)%
particle_x
.
eq
.
Hneut
)
then
if
(
S95_t1
(
i
)%
deltax
.
gt
.
max_lep_delta_Mh
)
then
max_lep_delta_Mh
=
S95_t1
(
i
)%
deltax
endif
endif
endif
enddo
!do i=lbound(S95_t2,dim=1),ubound(S95_t2,dim=1) !not needed until code is adapted to allow a nonzero deltax for type 2 tables
! if( S95_t2(i)%expt.eq.'LEP')then
! if(S95_t2(i)%deltax.gt.max_lep_delta_Mh)then
! max_lep_delta_Mh=S95_t2(i)%deltax
! endif
! endif
!enddo
if
(
debug
)
write
(
*
,
*
)
'max_tev_delta_Mh'
,
max_tev_delta_Mh
if
(
debug
)
write
(
*
,
*
)
'max_lep_delta_Mh'
,
max_lep_delta_Mh
TEVrange_Mhmin
=
TEVt1Mhmin
-
max_tev_delta_Mh
TEVrange_Mhmax
=
TEVt1Mhmax
+
max_tev_delta_Mh
LEPrange_Mhmin
=
LEPt12Mhmin
-
max_lep_delta_Mh
LEPrange_Mhmax
=
LEPt12Mhmax
+
max_lep_delta_Mh
!we need XS_SM_functions and BRfunctions to have a big enough range to cover the tables
if
(
TEVrange_Mhmax
.
gt
.
tevXS_SM_functions_xmax
)
then
write
(
*
,
*
)
TEVt1Mhmax
,
max_tev_delta_Mh
,
tevXS_SM_functions_xmax
stop
'need to extend upper range of XS_SM_functions or reduce delta_M_(LEP/TEV)'
endif
if
(
TEVrange_Mhmin
.
lt
.
tevXS_SM_functions_xmin
)
then
stop
'need to extend lower range of XS_SM_functions'
endif
if
(
(
TEVrange_Mhmax
.
gt
.
BRSMt1Mhmax
).
or
.(
LEPrange_Mhmax
.
gt
.
BRSMt1Mhmax
))
then
stop
'need to extend upper range of BRfunctions or reduce delta_M_(LEP/TEV)'
endif
if
(
TEVrange_Mhmin
.
lt
.
BRSMt1Mhmin
)
then
stop
'need to extend lower range of BRfunctions'
endif
if
(
LEPrange_Mhmin
.
lt
.
BRSMt1Mhmin
)
then
if
(
LEPt12Mhmin
.
lt
.
BRSMt1Mhmin
)
then
stop
'need to extend lower range of BRfunctions 1'
else
LEPrange_Mhmin
=
BRSMt1Mhmin
endif
endif
end subroutine
setup_S95tables
!******************************************************************
subroutine
calcfact
(
proc
,
t
,
cfact
,
Mi_av
,
Mj_av
,
nc
)
!******************************************************************
!for table type1, calls calcfact_t1
!for table type2, calls calcfact_t2
!*************************************** ***************************
use
usefulbits
,
only
:
dataset
,
listprocesses
!internal
implicit none
!--------------------------------------input
type
(
dataset
)
::
t
type
(
listprocesses
)
::
proc
!-----------------------------------output
double precision
::
cfact
,
Mj_av
,
Mi_av
integer
::
nc
!-------------------------------------------
select case
(
proc
%
ttype
)
case
(
1
)
call
calcfact_t1
(
proc
%
tlist
,
proc
%
findj
,
t
,
cfact
,
Mi_av
,
nc
)
Mj_av
=
Mi_av
case
(
2
)
call
calcfact_t2
(
proc
%
tlist
,
proc
%
findj
,
proc
%
findi
,
t
,
cfact
,
Mi_av
,
Mj_av
,
nc
)
case
default
stop
'wrong input to function calcfact in module channels'
end select
end subroutine
calcfact
!******************************************************************
subroutine
outputproc
(
proc
,
k
,
descrip
,
specific
)
!******************************************************************
!for table type1, calls outputproc_t1
!for table type2, calls outputproc_t2
!if neither table applies, writes message
!k is where output goes
!specific=1 if the specific process should be printed
!e.g. ee->h1Z->bbZ
!whereas specific==0 if the generic process should be printed
!e.g. ee->hiZ->bbZ
!******************************************************************
use
usefulbits
,
only
:
listprocesses
!input
implicit none
!--------------------------------------input
integer
,
intent
(
in
)
::
k
,
specific
integer
::
i
,
j
type
(
listprocesses
),
intent
(
in
)
::
proc
character
(
LEN
=
200
)
::
descrip
!-------------------------------------------
select case
(
specific
)
case
(
0
)
i
=
0
j
=
0
case
(
1
)
i
=
proc
%
findi
j
=
proc
%
findj
case
default
end select
select case
(
proc
%
ttype
)
case
(
0
)
descrip
=
'none of the processes apply'
case
(
1
)
call
outputproc_t1
(
proc
%
tlist
,
i
,
k
,
descrip
)
case
(
2
)
call
outputproc_t2
(
proc
%
tlist
,
i
,
j
,
k
,
descrip
)
case
default
stop
'wrong input to subroutine outputproc in module channels'
end select
end subroutine
outputproc
!******************************************************************
subroutine
check_against_bound
(
proc
,
fact
,
Mi_av
,
Mj_av
,
ratio
,
predobs
)
!******************************************************************
!for table type1, calls interpolate_tabletype1
!for table type2, calls interpolate_tabletype2
!******************************************************************
use
usefulbits
,
only
:
dataset
,
listprocesses
!internal
use
interpolate
implicit none
!--------------------------------------input
integer
::
predobs
double precision
::
fact
,
Mi_av
,
Mj_av
type
(
listprocesses
)
::
proc
!-------------------------------------output
double precision
::
ratio
!-----------------------------------internal
double precision
::
Mi
,
Mj
,
interpol
!-------------------------------------------
Mi
=
Mi_av
Mj
=
Mj_av
if
(
fact
.
ge
.
0.0D0
)
then
select case
(
proc
%
ttype
)
case
(
1
)
call
interpolate_tabletype1
(
Mi
,
S95_t1
(
proc
%
tlist
),
predobs
,
interpol
)
case
(
2
)
call
interpolate_tabletype2
(
Mi
,
Mj
,
S95_t2
(
proc
%
tlist
),
predobs
,
interpol
)
case
default
write
(
*
,
*
)
'wrong input to subroutine check_against_bound in module channels'
stop
end select
endif
if
(
interpol
.
ge
.
0
)
then
ratio
=
fact
/
interpol
else
ratio
=
-
1.0D0
endif
end subroutine
check_against_bound
!**********************************************************
subroutine
calcfact_t1
(
c
,
jj
,
t
,
cfact_t1
,
M_av
,
nc
)
!**********************************************************
! calculates fact for table type 1
! Takes in to account how Standard Model-like parameter point is
! and whether there are any Higgs with slightly higher masses which
! can be combined with his result
! note: numerator and denominator are worked out separately
!**********************************************************
use
usefulbits
,
only
:
dataset
,
np
,
div
,
TEVt1Mhmax
,
TEVt1Mhmin
,
vvsmall
use
theory_BRfunctions
use
theory_XS_SM_functions
implicit none
!--------------------------------------input
type
(
dataset
)
::
t
integer
::
c
,
jj
!-----------------------------------output
double precision
::
cfact_t1
,
M_av
integer
::
nc
!-------------------------------------------
integer
::
f
,
j
double precision
::
M_tot
double precision
::
BR_Hbb_SM_av
,
BR_HWW_SM_av
!double precision :: BR_Htautau_SM_av
!double precision :: BR_Hgaga_SM_av
double precision
::
tev_XS_HW_SM_av
,
tev_XS_HZ_SM_av
double precision
::
tev_XS_H_SM_av
!double precision :: tev_XS_H_SM_9674_av,tev_XS_H_SM_9713_av
double precision
::
tev_XS_Hb_SM_av
!double precision :: tev_XS_vbf_SM_av
double precision
::
tev_XS_ttH_SM_av
double precision
::
BR_Zll
,
BR_Znunu
,
BR_Wlnu
,
BR_Ztautau
double precision
::
BR_Whad
,
BR_Zhad
double precision
,
allocatable
::
mass
(:),
fact
(:)
logical
::
TevTable
,
HneutTevTableOutOfRange
integer
,
allocatable
::
model_like
(:)
integer
::
npart
!number of particles
!source: PDG (Yao et al. J Phys G 33 (2006))
BR_Zll
=
3.363D-2
+
3.366D-2
!BR_Zll = sum(l=e,mu), BR(Z ->l+ l-)
BR_Znunu
=
20
D
-
2
!BR_Znunu = BR(Z ->nu_l nu_l-bar) ('invisible')
BR_Wlnu
=
1
0.75D-2
+
1
0.57D-2
!BR_Wlnu = sum(l=e,mu),
!BR(W+ ->l+ nu_l) = BR(W- ->l- + nu_l-bar)
BR_Whad
=
6
7.6D-2
BR_Zhad
=
6
9.91D-2
BR_Ztautau
=
3.370D-2
npart
=
np
(
S95_t1
(
c
)%
particle_x
)
allocate
(
mass
(
npart
),
fact
(
npart
),
model_like
(
npart
))
mass
(:)
=
t
%
particle
(
S95_t1
(
c
)%
particle_x
)%
M
(:)
fact
=
0.0D0
cfact_t1
=
0.0D0
model_like
=
0
!now calculate numerator of 'fact'
do
j
=
1
,
npart
if
(
(
abs
(
mass
(
jj
)
-
mass
(
j
)).
le
.
S95_t1
(
c
)%
deltax
)
&
&
.
and
.(
mass
(
jj
).
le
.
mass
(
j
)
)
)
then
select case
(
S95_t1
(
c
)%
id
)
!these can be compacted, but not doing it yet... doing it at the
!same time as revamping the model_likeness test
case
(
142
)
fact
(
j
)
=
t
%
lep
%
XS_hjZ_ratio
(
j
)
*
t
%
BR_hjbb
(
j
)
!notice that this is not absolute XS
case
(
143
)
fact
(
j
)
=
t
%
lep
%
XS_hjZ_ratio
(
j
)
*
t
%
BR_hjtautau
(
j
)
!notice that this is not absolute XS
case
(
300
)
fact
(
j
)
=
t
%
lep
%
XS_hjZ_ratio
(
j
)
!notice that this is not absolute XS
case
(
400
,
401
,
402
,
403
)
fact
(
j
)
=
t
%
lep
%
XS_hjZ_ratio
(
j
)
*
t
%
BR_hjinvisible
(
j
)
!notice that this is not absolute XS
case
(
500
)
fact
(
j
)
=
t
%
lep
%
XS_hjZ_ratio
(
j
)
*
t
%
BR_hjgaga
(
j
)
!notice that this is not absolute XS
case
(
600
)
fact
(
j
)
=
t
%
lep
%
XS_hjZ_ratio
(
j
)
&
!notice that this is not absolute XS
&
*
(
t
%
BR_hjss
(
j
)
+
t
%
BR_hjcc
(
j
)
+
t
%
BR_hjbb
(
j
)
+
t
%
BR_hjgg
(
j
))
case
(
711
,
713
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
721
,
723
,
741
,
743
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
731
,
733
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
801
,
811
,
821
)
fact
(
j
)
=
t
%
lep
%
XS_HpjHmj_ratio
(
j
)
*
(
t
%
BR_Hpjcs
(
j
)
+
t
%
BR_Hpjcb
(
j
))
**
2.0D0
!notice that this is not absolute XS
case
(
802
)
fact
(
j
)
=
t
%
lep
%
XS_HpjHmj_ratio
(
j
)
*
(
t
%
BR_Hpjcs
(
j
)
+
t
%
BR_Hpjcb
(
j
))
*
2.0D0
*
t
%
BR_Hpjtaunu
(
j
)
!notice that this is not absolute XS
case
(
803
,
813
)
fact
(
j
)
=
t
%
lep
%
XS_HpjHmj_ratio
(
j
)
*
t
%
BR_Hpjtaunu
(
j
)
**
2.0D0
!notice that this is not absolute XS
case
(
8742
,
4493
,
9475
,
5482
,
5570
,
5876
,
1024
,
9889
,
3534
,
6089
,
10235
,
3047
,
3564
)
fact
(
j
)
=
t
%
tev
%
XS_hjZ_ratio
(
j
)
*
t
%
tev
%
XS_HZ_SM
(
j
)
*
t
%
BR_hjbb
(
j
)
case
(
8958
,
5489
,
5624
,
9236
,
3930
,
6039
,
3216
,
10433
)
fact
(
j
)
=
t
%
tev
%
XS_hj_ratio
(
j
)
*
t
%
tev
%
XS_H_SM
(
j
)
*
t
%
BR_hjWW
(
j
)
case
(
5757
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
8957
,
5472
,
9219
,
9463
,
9596
,
5828
,
1970
,
3493
,
5972
,
3155
,
5613
,
9868
,
10068
,
6092
,
10217
,
10239
,
0874
)
fact
(
j
)
=
t
%
tev
%
XS_hjW_ratio
(
j
)
*
t
%
tev
%
XS_HW_SM
(
j
)
*
t
%
BR_hjbb
(
j
)
case
(
5485
,
7307
,
5873
)
fact
(
j
)
=
t
%
tev
%
XS_hjW_ratio
(
j
)
*
t
%
tev
%
XS_HW_SM
(
j
)
*
t
%
BR_hjWW
(
j
)
case
(
9071
,
2491
,
5740
,
5980
,
1014
,
3363
)
fact
(
j
)
=
t
%
tev
%
XS_hj_ratio
(
j
)
*
t
%
tev
%
XS_H_SM
(
j
)
*
t
%
BR_hjtautau
(
j
)
case
(
8961
,
0598
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
9284
,
5503
,
5726
,
10105
)
fact
(
j
)
=
t
%
tev
%
XS_hjb_ratio
(
j
)
*
t
%
tev
%
XS_Hb_SM
(
j
)
*
t
%
BR_hjbb
(
j
)
/
0.9D0
/
2.0D0
case
(
6083
)
fact
(
j
)
=
t
%
tev
%
XS_hjb_ratio
(
j
)
*
t
%
tev
%
XS_Hb_SM
(
j
)
*
t
%
BR_hjtautau
(
j
)
/
0.1D0
/
2.0D0
case
(
1514
,
5601
,
5737
)
fact
(
j
)
=
(
t
%
tev
%
XS_hjZ_ratio
(
j
)
*
t
%
tev
%
XS_HZ_SM
(
j
)
&
&
+
t
%
tev
%
XS_hjW_ratio
(
j
)
*
t
%
tev
%
XS_HW_SM
(
j
)
&
&
+
t
%
tev
%
XS_hj_ratio
(
j
)
*
t
%
tev
%
XS_H_SM
(
j
)
&
&
+
t
%
tev
%
XS_vbf_ratio
(
j
)
*
t
%
tev
%
XS_vbf_SM
(
j
)
&
&
)
&
&
*
t
%
BR_hjgaga
(
j
)
case
(
5858
,
6177
,
1887
,
10065
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
7081
,
9166
,
9483
,
5586
,
9642
,
1266
,
0432
,
9891
,
5285
,
3935
,
6087
,
6170
,
10212
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
10010
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
9248
,
10133
,
10439
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
4800
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
5845
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
9290
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
5984
,
9714
,
6006
,
4481
,
6095
,
10432
,
6179
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
3556
,
1931
)
fact
(
j
)
=
t
%
tev
%
XS_hjb_ratio
(
j
)
*
t
%
tev
%
XS_Hb_c2_SM
(
j
)
*
t
%
BR_hjbb
(
j
)
case
(
9465
,
5871
,
9022
,
0710
,
9887
,
4162
,
10102
,
4468
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
!case(9023)
! call model_likeness(j,S95_t1(c)%id,t,model_like(j),fact(j))
case
(
9674
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
9897
,
9999
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
0024
,
5985
,
0968
,
5974
)
fact
(
j
)
=
t
%
tev
%
XS_hjb_ratio
(
j
)
*
t
%
tev
%
XS_Hb_c3_SM
(
j
)
*
t
%
BR_hjtautau
(
j
)
case
(
5739
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
fact
(
j
)
=
fact
(
j
)
*
t
%
tev
%
XS_ttH_SM
(
j
)
*
t
%
BR_Hbb_SM
(
j
)
!to get the normalisation right
case
(
0611
)
fact
(
j
)
=
t
%
tev
%
XS_hj_ratio
(
j
)
*
t
%
tev
%
XS_H_SM
(
j
)
*
t
%
BR_hjZga
(
j
)
case
(
1269
,
1270
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
1811
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
1812
,
8353
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
6008
,
9998
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
6082
,
6182
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
6096
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
(
6091
)
call
model_likeness
(
j
,
S95_t1
(
c
)%
id
,
t
,
model_like
(
j
),
fact
(
j
))
case
default
stop
'wrong input to function calcfact_t1 in module S95tables'
end select
endif
enddo
!TevTable=True if it's a tevatron table we're looking at
TevTable
=
((
S95_t1
(
c
)%
expt
.
eq
.
'CDF'
).
or
.(
S95_t1
(
c
)%
expt
.
eq
.
' D0'
))
&
&
.
or
.(
S95_t1
(
c
)%
expt
.
eq
.
'TCB'
)
if
(
fact
(
jj
).
le
.
vvsmall
)
then
!A !Higgs jj doesn't contribute - wait until another call of this subroutine before
!looking at nearby masses
M_av
=
mass
(
jj
)
nc
=
0
cfact_t1
=
0.0D0
else
!A
!find M_av (only using higgs which have non-zero fact):
f
=
0
M_tot
=
0.0D0
do
j
=
1
,
npart
if
(
fact
(
j
).
gt
.
vvsmall
)
then
f
=
f
+
1
M_tot
=
M_tot
+
mass
(
j
)
endif
enddo
nc
=
f
!f will always be > 0 because we've already made sure that fact(jj)>0.0D0
M_av
=
M_tot
/
dble
(
nc
)
! HneutTevTableOutOfRange=True if it's a tevatron table we're looking at but mass isn't in the range of
! any of the tevatron tables (means there's no need to calculate BR_*_SM_av, tev_*_SM_av, so it might be out of range of
! the subroutines that calculate them)
HneutTevTableOutOfRange
=
.
False
.
! It's now M_av that gets looked up in the expt tables
! so no point calculating cfact_t1 if M_av out of the range of the expt tables
! but this condition MUST be more strict than the one used in theo_manip to calculate t%tev%XS_*_SM(j)
if
(
S95_t1
(
c
)%
particle_x
.
eq
.
Hneut
)
then
if
(
(
M_av
.
lt
.
TEVt1Mhmin
).
or
.(
M_av
.
gt
.
TEVt1Mhmax
)
)
then
HneutTevTableOutOfRange
=
.
True
.
endif
endif
if
(.
not
.
TevTable
)
then
!B
cfact_t1
=
sum
(
fact
)
elseif
(
S95_t1
(
c
)%
particle_x
.
ne
.
Hneut
)
then
!B
cfact_t1
=
sum
(
fact
)
elseif
(
HneutTevTableOutOfRange
)
then
!B
cfact_t1
=
0.0D0
else
!B
if
(
f
.
eq
.
1
)
then
!have already calculated these in theo_manip to save time
BR_Hbb_SM_av
=
t
%
BR_Hbb_SM
(
jj
)
!BR_Htautau_SM_av = t%BR_Htautau_SM(jj)
BR_HWW_SM_av
=
t
%
BR_HWW_SM
(
jj
)
!BR_Hgaga_SM_av = t%BR_Hgaga_SM(jj)
tev_XS_HW_SM_av
=
t
%
tev
%
XS_HW_SM
(
jj
)
tev_XS_HZ_SM_av
=
t
%
tev
%
XS_HZ_SM
(
jj
)
tev_XS_H_SM_av
=
t
%
tev
%
XS_H_SM
(
jj
)
tev_XS_Hb_SM_av
=
t
%
tev
%
XS_Hb_SM
(
jj
)
!tev_XS_vbf_SM_av = t%tev%XS_vbf_SM(jj)
tev_XS_ttH_SM_av
=
t
%
tev
%
XS_ttH_SM
(
jj
)
!tev_XS_H_SM_9713_av = t%tev%XS_H_SM_9713(jj)
!tev_XS_H_SM_9674_av = t%tev%XS_H_SM_9674(jj)
else
BR_Hbb_SM_av
=
BRSM_Hbb
(
M_av
)
!BR_Htautau_SM_av = BRSM_Htautau(M_av)
BR_HWW_SM_av
=
BRSM_HWW
(
M_av
)
!BR_Hgaga_SM_av = BRSM_Hgaga(M_av)
tev_XS_HW_SM_av
=
XS_tev_qq_HW_SM
(
M_av
)
tev_XS_HZ_SM_av
=
XS_tev_qq_HZ_SM
(
M_av
)
tev_XS_H_SM_av
=
XS_tev_gg_H_SM
(
M_av
)
+
XS_tev_bb_H_SM
(
M_av
)
tev_XS_Hb_SM_av
=
XS_tev_bg_Hb_SM
(
M_av
)
!tev_XS_vbf_SM_av = XS_tev_vbf_SM(M_av)
tev_XS_ttH_SM_av
=
XS_tev_ttH_SM
(
M_av
)
!tev_XS_H_SM_9713_av = XS_tev_gg_H_SM_9713(M_av)+XS_tev_bb_H_SM(M_av)
!tev_XS_H_SM_9674_av = XS_tev_gg_H_SM_9674(M_av)+XS_tev_bb_H_SM(M_av)
endif
! now include denominator of 'fact'
select case
(
S95_t1
(
c
)%
id
)
case
(
8742
,
5482
,
5570
,
4493
,
9475
,
5876
,
1024
,
9889
,
3534
,
6089
,
10235
,
3047
,
3564
)
fact
=
fact
/
tev_XS_HZ_SM_av
/
BR_Hbb_SM_av
case
(
8958
,
5489
,
5624
,
9236
,
3930
)
do
j
=
1
,
npart
fact
(
j
)
=
div
(
fact
(
j
)
/
tev_XS_H_SM_av
,
BR_HWW_SM_av
,
0.0D0
,
0.0D0
)
enddo
case
(
8957
,
5472
,
9219
,
9463
,
9596
,
5828
,
1970
,
3493
,
5972
,
3155
,
5613
,
9868
,
10068
,
6092
,
10217
,
10239
,
0874
)
fact
=
fact
/
tev_XS_HW_SM_av
/
BR_Hbb_SM_av
case
(
5503
,
9284
,
5726
,
10105
)
fact
=
fact
/
tev_XS_Hb_SM_av
case
(
6083
)
fact
=
fact
/
tev_XS_Hb_SM_av
case
(
7307
,
5873
)
do
j
=
1
,
npart
fact
(
j
)
=
div
(
fact
(
j
)
/
tev_XS_HW_SM_av
,
BR_HWW_SM_av
,
0.0D0
,
0.0D0
)
enddo
case
(
8961
,
0598
,
10010
,
9290
,
9674
,
9897
,
9999
)
case
(
7081
,
9166
,
9483
,
5586
,
9642
,
1266
,
0432
,
9891
,
5285
,
3935
,
6087
,
6170
,
10212
)
case
(
9248
,
5845
,
4800
,
5858
,
6177
,
1887
,
10065
,
10133
,
10439
)
case
(
9465
,
5871
,
9022
,
9023
,
0710
,
9887
)
case
(
5984
,
9714
,
4162
,
10102
,
6006
,
4481
,
4468
,
5757
,
6095
,
10432
,
6179
)
case
(
5739
)
fact
=
fact
/
tev_XS_ttH_SM_av
/
BR_Hbb_SM_av
case
(
6008
,
9998
)
case
(
6082
,
6182
)
case
(
6096
)
case
(
6091
)
case
(
5485
,
9071
,
2491
,
5601
,
1514
,
3556
,
5740
,
0024
,
5980
,
1014
,
5985
,
&
&
0611
,
3363
,
6039
,
3216
,
10433
,
0968
,
5974
,
1931
)
case
default
stop
'error calculating denom. in calcfact_t1'
end select
cfact_t1
=
sum
(
fact
)
endif
!B
endif
!A
deallocate
(
mass
)
deallocate
(
fact
)
deallocate
(
model_like
)
end subroutine
calcfact_t1
!**********************************************************
subroutine
calcfact_t2
(
c
,
jj
,
ii
,
t
,
cfact_t2
,
Mi_av
,
Mj_av
,
nc
)
!**********************************************************
!calculates fact for table type 2
!**********************************************************
use
usefulbits
,
only
:
dataset
,
np
,
vsmall
implicit none
!--------------------------------------input
type
(
dataset
)
::
t
integer
::
c
,
jj
,
ii
!-----------------------------------output
double precision
::
cfact_t2
,
Mi_av
,
Mj_av
integer
::
nc
!-------------------------------------------
integer
::
f
,
i
,
j
,
np1
,
np2
double precision
::
fact
,
eps2
double precision
,
allocatable
::
massi
(:),
massj
(:)
eps2
=
0.02D0
np1
=
np
(
S95_t2
(
c
)%
particle_x1
)
np2
=
np
(
S95_t2
(
c
)%
particle_x2
)
allocate
(
massi
(
np1
))
allocate
(
massj
(
np2
))
massi
(:)
=
t
%
particle
(
S95_t2
(
c
)%
particle_x1
)%
M
(:)
massj
(:)
=
t
%
particle
(
S95_t2
(
c
)%
particle_x2
)%
M
(:)
Mi_av
=
massi
(
ii
)
Mj_av
=
massj
(
jj
)
fact
=
0.0D0
cfact_t2
=
0.0D0
f
=
0
j
=
jj
i
=
ii
! do j=1,nHneut
! if ( (abs(mass(jj)-mass(j)).lt.S95_t2(c)%deltax) &
! & .and.( mass(jj).le.mass(j) ) )then
! do i=1,nHneut
! if( ( (abs(mass(ii)-mass(i)).lt.S95_t2(c)%deltax)&
! & .and.( mass(ii).le.mass(i) ) ) &
! & .and.( mass(j).ge.mass(i) ) ) then
if
(
massj
(
j
).
ge
.
massi
(
i
)
)
then
f
=
f
+
1
select case
(
S95_t2
(
c
)%
id
)
case
(
150
)
fact
=
test_appl
(
t
%
lep
%
XS_hjZ_ratio
(
j
)
*
t
%
BR_hjhihi
(
j
,
i
)
*
t
%
BR_hjbb
(
i
)
**
2.0D0
)
!table 15 hep-ex/0602042 XS ratio
case
(
160
)
fact
=
test_appl
(
t
%
lep
%
XS_hjZ_ratio
(
j
)
*
t
%
BR_hjhihi
(
j
,
i
)
*
t
%
BR_hjtautau
(
i
)
**
2.0D0
)
!table 16 hep-ex/0602042 XS ratio
case
(
180
)
fact
=
test_appl
(
t
%
lep
%
XS_hjhi_ratio
(
j
,
i
)
*
t
%
BR_hjbb
(
j
)
*
t
%
BR_hjbb
(
i
))
!table 18 hep-ex/0602042 XS ratio
case
(
190
)
fact
=
test_appl
(
t
%
lep
%
XS_hjhi_ratio
(
j
,
i
)
*
t
%
BR_hjtautau
(
j
)
*
t
%
BR_hjtautau
(
i
))
!table 19 hep-ex/0602042 XS ratio
case
(
200
)
fact
=
test_appl
(
t
%
lep
%
XS_hjhi_ratio
(
j
,
i
)
*
t
%
BR_hjhihi
(
j
,
i
)
*
t
%
BR_hjbb
(
i
)
**
3.0D0
)
!table 20 hep-ex/0602042 XS ratio
case
(
210
)
fact
=
test_appl
(
t
%
lep
%
XS_hjhi_ratio
(
j
,
i
)
*
t
%
BR_hjhihi
(
j
,
i
)
*
t
%
BR_hjtautau
(
i
)
**
3.0D0
)
!table 21 hep-ex/0602042 XS ratio
case
(
220
)
fact
=
test_appl
(
t
%
lep
%
XS_hjZ_ratio
(
j
)
*
t
%
BR_hjhihi
(
j
,
i
)
*
t
%
BR_hjbb
(
i
)
*
t
%
BR_hjtautau
(
i
))
!table 22 hep-ex/0602042 XS ratio
case
(
230
)
fact
=
test_appl
(
t
%
lep
%
XS_hjhi_ratio
(
j
,
i
)
*
t
%
BR_hjbb
(
j
)
*
t
%
BR_hjtautau
(
i
))
!table 23 hep-ex/0602042 XS ratio
case
(
240
)
fact
=
test_appl
(
t
%
lep
%
XS_hjhi_ratio
(
j
,
i
)
*
t
%
BR_hjtautau
(
j
)
*
t
%
BR_hjbb
(
i
))
!table 24 hep-ex/0602042 XS ratio
case
(
905
)
fact
=
test_appl
(
t
%
lep
%
XS_CpjCmj
(
j
)
*
t
%
BR_CjqqNi
(
j
,
i
)
**
2.0D0
)
!fig 5 hep-ex/0401026 absolute XS in fb
case
(
906
)
fact
=
test_appl
(
t
%
lep
%
XS_CpjCmj
(
j
)
*
t
%
BR_CjqqNi
(
j
,
i
)
*
t
%
BR_CjlnuNi
(
j
,
i
))
!fig 6 hep-ex/0401026 absolute XS in fb
case
(
907
)
fact
=
test_appl
(
t
%
lep
%
XS_CpjCmj
(
j
)
*
t
%
BR_CjlnuNi
(
j
,
i
)
**
2.0D0
)
!fig 7 hep-ex/0401026 absolute XS in fb
case
(
908
)
fact
=
test_appl
(
t
%
lep
%
XS_CpjCmj
(
j
))
!fig 8 hep-ex/0401026 absolute XS in fb
case
(
909
)
fact
=
test_appl
(
t
%
lep
%
XS_NjNi
(
j
,
i
)
*
t
%
BR_NjqqNi
(
j
,
i
))
!fig 9 hep-ex/0401026 absolute XS in fb
case
(
910
)
fact
=
test_appl
(
t
%
lep
%
XS_NjNi
(
j
,
i
))
!fig 10 hep-ex/0401026 absolute XS in fb
case
default
stop
'wrong input to function calcfact_t2 in module S95tables'
end select
else
fact
=
0.0D0
endif
cfact_t2
=
cfact_t2
+
fact
! enddo
! endif
! enddo
nc
=
f
!if(f.gt.1)stop'f.ne.1 in calcfact_t2: this has not yet been thoroughly checked'
!*****
deallocate
(
massi
)
deallocate
(
massj
)
contains
!********************************************************
function
test_appl
(
x
)
!********************************************************
!if hi->hj+hj is needed, check it's kinematically possible
!********************************************************
implicit none
!--------------------------------------input
double precision
::
x
!-----------------------------------function
double precision
::
test_appl
!-------------------------------------------
select case
(
S95_t2
(
c
)%
id
)
case
(
150
,
160
,
180
,
190
,
200
,
210
,
220
,
230
,
240
)
if
(
S95_t2
(
c
)%
needs_M2_gt_2M1
.
and
.(
massj
(
j
).
lt
.
2.0D0
*
massi
(
i
)))
then
test_appl
=
0.0D0
else
test_appl
=
x
endif
case
(
905
,
906
,
907
,
909
)
if
(
abs
(
minval
(
massi
)
-
massi
(
i
)).
gt
.
vsmall
)
then
!checking that lightest neutralino in process is lightest neutralino in model
test_appl
=
0.0D0
else
test_appl
=
x
endif
case
(
908
)
if
(
abs
(
t
%
BR_CjWNi
(
j
,
i
)
-
1.0D0
)
.
gt
.
eps2
)
then
test_appl
=
0.0D0
elseif
(
abs
(
minval
(
massi
)
-
massi
(
i
)).
gt
.
vsmall
)
then
!checking that lightest neutralino in process is lightest neutralino in model
test_appl
=
0.0D0
else
test_appl
=
x
endif
case
(
910
)
if
(
abs
(
t
%
BR_NjZNi
(
j
,
i
)
-
1.0D0
)
.
gt
.
eps2
)
then
test_appl
=
0.0D0
elseif
(
abs
(
minval
(
massi
)
-
massi
(
i
)).
gt
.
vsmall
)
then
!checking that lightest neutralino in process is lightest neutralino in model
test_appl
=
0.0D0
else
test_appl
=
x
endif
case
default
stop
'error in function test_appl'
end select
end function
test_appl
end subroutine
calcfact_t2
!********************************************************
subroutine
outputproc_t1
(
tlistn
,
jj
,
k
,
descrip
)
!********************************************************
! uses information about the process to output a description
! for processes using table type 1
! note: at the moment, np(x) (and so ii and jj) needs to be 1 digit long i.e. nH<10
!********************************************************
implicit none
!--------------------------------------input
integer
::
tlistn
integer
::
jj
,
k
!-----------------------------------internal
character
(
LEN
=
1
)
::
j
character
(
LEN
=
45
)
::
label
character
(
LEN
=
200
)
::
descrip
!-------------------------------------------
if
(
jj
.
ne
.
0
)
then
write
(
j
,
'(I1)'
)
jj
else
j
=
'j'
endif
if
(
k
.
eq
.
21
)
then
label
=
''
!no need to lable each line in Key.dat
else
label
=
'('
//
trim
(
S95_t1
(
tlistn
)%
label
)
//
')'
endif
descrip
=
''
select case
(
S95_t1
(
tlistn
)%
id
)
case
(
142
)
descrip
=
' (ee)->(h'
//
j
//
')Z->(b b-bar)Z '
//
label
case
(
143
)
descrip
=
' (ee)->(h'
//
j
//
')Z->(tau tau)Z '
//
label
case
(
300
)
descrip
=
' (ee)->(h'
//
j
//
')Z->(...)Z '
//
label
case
(
400
,
401
,
402
,
403
)
descrip
=
' (ee)->(h'
//
j
//
')Z->(invisible)Z '
//
label
case
(
500
)
descrip
=
' (ee)->(h'
//
j
//
')Z->(gamma gamma)Z '
//
label
case
(
600
)
descrip
=
' (ee)->(h'
//
j
//
')Z->(2 jets)Z '
//
label
case
(
711
)
descrip
=
' (ee)->b b-bar(h'
//
j
//
')->b b-bar(b b-bar) where h'
//
j
//
' is CP even '
//
label
case
(
713
)
descrip
=
' (ee)->b b-bar(h'
//
j
//
')->b b-bar(b b-bar) where h'
//
j
//
' is CP odd '
//
label
case
(
721
,
741
)
descrip
=
' (ee)->b b-bar(h'
//
j
//
')->b b-bar(tau tau) where h'
//
j
//
' is CP even '
//
label
case
(
723
,
743
)
descrip
=
' (ee)->b b-bar(h'
//
j
//
')->b b-bar(tau tau) where h'
//
j
//
' is CP odd '
//
label
case
(
731
)
descrip
=
' (ee)->tau tau(h'
//
j
//
')->tau tau(tau tau) where h'
//
j
//
' is CP even '
//
label
case
(
733
)
descrip
=
' (ee)->tau tau(h'
//
j
//
')->tau tau(tau tau) where h'
//
j
//
' is CP odd '
//
label
case
(
801
,
811
,
821
)
descrip
=
' (ee)->(H'
//
j
//
'+)(H'
//
j
//
'-)->4 quarks '
//
label
case
(
802
)
descrip
=
' (ee)->(H'
//
j
//
'+)(H'
//
j
//
'-)->(2 quarks) tau nu '
//
label
case
(
803
,
813
)
descrip
=
' (ee)->(H'
//
j
//
'+)(H'
//
j
//
'-)->tau nu tau nu'
//
label
case
(
5482
,
5570
,
8742
,
4493
,
9475
,
5876
,
1024
,
9889
,
3534
,
6089
,
10235
,
3047
,
3564
)
descrip
=
' (p p-bar)->Z(h'
//
j
//
')->l l (b b-bar) '
//
label
case
(
9236
,
3930
,
8958
,
6039
,
3216
,
10433
)
descrip
=
' (p p-bar)->h'
//
j
//
'->W W '
//
label
case
(
9219
,
9463
,
5472
,
8957
,
9596
,
5828
,
1970
,
3493
,
5972
,
3155
,
5613
,
9868
,
10068
,
6092
,
10217
,
10239
,
0874
)
descrip
=
' (p p-bar)->W(h'
//
j
//
')->l nu (b b-bar) '
//
label
case
(
5489
)
descrip
=
' (p p-bar)->h'
//
j
//
'->W W->e mu '
//
label
case
(
5624
)
descrip
=
' (p p-bar)->h'
//
j
//
'->W W->l l '
//
label
case
(
5757
)
descrip
=
' (p p-bar)->h'
//
j
//
'/VBF->W W->l l where h'
//
j
//
' is SM-like '
//
label
case
(
5485
,
5873
)
descrip
=
' (p p-bar)->W(h'
//
j
//
')->W W W->l l nu nu '
//
label
case
(
9071
,
2491
,
5740
,
5980
,
1014
,
3363
)
descrip
=
' (p p-bar)->h'
//
j
//
'->tau tau '
//
label
case
(
8961
,
9465
,
9290
,
9713
,
9674
,
0598
,
9897
,
9998
,
9999
,
6008
,
6096
)
descrip
=
' (p p-bar)->h'
//
j
//
'+... where h'
//
j
//
' is SM-like '
//
label
case
(
9284
,
5503
,
5726
,
3556
,
10105
,
1931
)
descrip
=
' (p p-bar)->h'
//
j
//
'(b/b-bar)->(b b-bar) (b/b-bar) '
//
label
case
(
7307
)
descrip
=
' (p p-bar)->W(h'
//
j
//
')->W W W '
//
label
case
(
5601
,
5737
,
1514
)
descrip
=
' (p p-bar)->h'
//
j
//
'+...->gamma gamma+... '
//
label
case
(
5858
,
6177
,
1887
,
10065
)
descrip
=
' (p p-bar)->h'
//
j
//
'+...->gamma gamma+... where h'
//
j
//
' is SM-like '
//
label
case
(
7081
,
9166
,
9483
,
5586
,
9642
,
1266
,
0432
,
9891
,
5285
,
3935
,
6087
,
6170
,
10212
)
descrip
=
' (p p-bar)->V h'
//
j
//
'-> (b b-bar) +missing Et where h'
//
j
//
' is SM-like '
//
label
case
(
6091
)
descrip
=
' (p p-bar)->V h'
//
j
//
'-> ll + X where h'
//
j
//
' is SM-like '
//
label
case
(
10010
)
descrip
=
' (p p-bar)->V (h'
//
j
//
')/VBF-> (b b-bar) q q where h'
//
j
//
' is SM-like '
//
label
case
(
9248
,
10133
,
10439
)
descrip
=
' (p p-bar)->h'
//
j
//
'+...->tau tau +... where h'
//
j
//
' is SM-like '
//
label
case
(
4800
)
descrip
=
' (p p-bar)->h'
//
j
//
'+...->tau tau (2 jets) where h'
//
j
//
' is SM-like '
//
label
case
(
5845
)
descrip
=
' (p p-bar)->h'
//
j
//
'+...->tau tau (2 jets) where h'
//
j
//
' is SM-like '
//
label
case
(
5871
)
descrip
=
' (p p-bar)->h'
//
j
//
'+...->W W +... ->l l nu nu +... where h'
//
j
//
' is SM-like '
//
label
case
(
6082
,
6182
)
descrip
=
' (p p-bar)->h'
//
j
//
'+...->V V +... ->e mu missing Et +... where h'
//
j
//
' is SM-like '
//
label
case
(
5984
,
9714
,
6006
,
9022
,
9023
,
0710
,
9887
,
4162
,
10102
,
4481
,
4468
,
6095
,
10432
,
6179
)
descrip
=
' (p p-bar)->h'
//
j
//
'+...->W W +... where h'
//
j
//
' is SM-like '
//
label
case
(
0024
,
5985
,
0968
,
5974
,
6083
)
descrip
=
' (p p-bar)->h'
//
j
//
'(b/b-bar)->(tau tau) (b/b-bar) '
//
label
case
(
5739
)
descrip
=
' (p p-bar)->t t-bar h'
//
j
//
'->t t-bar b b-bar '
//
label
case
(
0611
)
descrip
=
' (p p-bar)->h'
//
j
//
'->Z gamma '
//
label
case
(
1811
)
descrip
=
' t->(H'
//
j
//
'+)b->(2 quarks) b '
//
label
case
(
1812
,
8353
)
descrip
=
' t->(H'
//
j
//
'+)b->tau nu b '
//
label
case
(
1269
)
descrip
=
' t->(H'
//
j
//
'+)b->(c s) b (lower mass region) '
//
label
case
(
1270
)
descrip
=
' t->(H'
//
j
//
'+)b->(c s) b (higher mass region) '
//
label
case
default
stop
'wrong input to function outputproc_t1 in module S95tables (1)'
end select
end subroutine
outputproc_t1
!********************************************************
subroutine
outputproc_t2
(
tlistn
,
ii
,
jj
,
k
,
descrip
)
!********************************************************
! uses information about the process to output a description
! for processes using table type 1
! note: at the moment, np(x) (and so ii and jj) needs to be 1 digit long i.e. np(x)<10
!********************************************************
implicit none
!--------------------------------------input
integer
::
tlistn
integer
::
ii
,
jj
,
k
!-----------------------------------internal
character
(
LEN
=
1
)
::
j
,
i
character
(
LEN
=
45
)
::
label
character
(
LEN
=
200
)
::
descrip
!-------------------------------------------
if
((
ii
.
ne
.
0
).
and
.(
jj
.
ne
.
0
))
then
write
(
i
,
'(I1)'
)
ii
write
(
j
,
'(I1)'
)
jj
else
i
=
'i'
j
=
'j'
endif
if
(
k
.
eq
.
21
)
then
label
=
''
!no need to lable each line in Key.dat
else
label
=
'('
//
trim
(
S95_t2
(
tlistn
)%
label
)
//
')'
endif
select case
(
S95_t2
(
tlistn
)%
id
)
case
(
150
)
descrip
=
' (ee)->(h'
//
j
//
'->h'
//
i
//
' h'
//
i
//
')Z->(b b b b)Z '
//
label
case
(
160
)
descrip
=
' (ee)->(h'
//
j
//
'->h'
//
i
//
' h'
//
i
//
')Z->(tau tau tau tau)Z '
//
label
case
(
180
)
descrip
=
' (ee)->(h'
//
j
//
' h'
//
i
//
')->(b b b b) '
//
label
case
(
190
)
descrip
=
' (ee)->(h'
//
j
//
' h'
//
i
//
')->(tau tau tau tau) '
//
label
case
(
200
)
descrip
=
' (ee)->(h'
//
j
//
'->h'
//
i
//
' h'
//
i
//
')h'
//
i
//
'->(b b b b)b b '
//
label
case
(
210
)
descrip
=
' (ee)->(h'
//
j
//
'->h'
//
i
//
' h'
//
i
//
')h'
//
i
//
'->(tau tau tau tau)tau tau '
//
label
case
(
220
)
descrip
=
' (ee)->(h'
//
j
//
'->h'
//
i
//
' h'
//
i
//
')Z->(b b)(tau tau)Z '
//
label
case
(
230
)
descrip
=
' (ee)->(h'
//
j
//
'->b b)(h'
//
i
//
'->tau tau) '
//
label
case
(
240
)
descrip
=
' (ee)->(h'
//
j
//
'->tau tau)(h'
//
i
//
'->b b) '
//
label
case
(
905
)
descrip
=
' (ee)->(C'
//
j
//
'+)(C'
//
j
//
'-)-> (q q N'
//
i
//
') (q q N'
//
i
//
') '
//
label
case
(
906
)
descrip
=
' (ee)->(C'
//
j
//
'+)(C'
//
j
//
'-)-> q q l nu N'
//
i
//
' N'
//
i
//
' '
//
label
case
(
907
)
descrip
=
' (ee)->(C'
//
j
//
'+)(C'
//
j
//
'-)-> (l nu N'
//
i
//
') (l nu N'
//
i
//
') '
//
label
case
(
908
)
descrip
=
' (ee)->(C'
//
j
//
'+)(C'
//
j
//
'-) with all C'
//
j
//
' decaying to W + N'
//
i
//
' '
//
label
case
(
909
)
descrip
=
' (ee)->(N'
//
j
//
') N'
//
i
//
'-> (q q N'
//
i
//
') N'
//
i
//
' '
//
label
case
(
910
)
descrip
=
' (ee)->N'
//
j
//
' N'
//
i
//
' with all N'
//
j
//
' decaying to Z + N'
//
i
//
' '
//
label
case
default
stop
'wrong input to function outputproc_t2 in module S95tables (2)'
end select
end subroutine
outputproc_t2
!******************************************************************
subroutine
model_likeness
(
j
,
id
,
t
,
model_like
,
sigmaXbr
)
!*****************************************************************
! Tests how Standard Model-like a parameter point is
! 0 means Mi.ge.MSingleLim (treat as single channel)
! 1 passes the SM-like test and Mi.lt.MSingleLim
! -1 fails the SM-like test and Mi.lt.MSingleLim
use
usefulbits
,
only
:
dataset
,
div
implicit none
!--------------------------------------input
type
(
dataset
)
::
t
integer
::
id
,
j
!-----------------------------------internal
integer
::
ns
,
nb
,
n
double precision
,
allocatable
::
XS_rat
(:),
BR_rat
(:)
integer
::
model_like
,
testSMratios
double precision
::
sigmaXbr
integer
::
is
,
ib
double precision
::
s
,
b
double precision
,
allocatable
::
dsbys
(:),
dbbyb
(:)
logical
::
correct_properties
correct_properties
=
.
True
.
n
=
t1elementnumberfromid
(
S95_t1
,
id
)
select case
(
id
)
case
(
711
,
713
,
721
,
723
,
731
,
733
,
741
,
743
,
5739
)
!the only meodel-likeness test these have is CP, so we can
!have a non-zero deltax
case
default
if
(
S95_t1
(
n
)%
deltax
.
gt
.
0.0D0
)
then
write
(
*
,
*
)
'hello id='
,
id
,
'deltax='
,
S95_t1
(
n
)%
deltax
stop
'error in subroutine model_likeness (1)'
endif
end select
select case
(
id
)
case
(
8961
,
0598
)
ns
=
3
;
allocate
(
XS_rat
(
ns
))
nb
=
2
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
tev
%
XS_hjW_ratio
(
j
)
XS_rat
(
2
)
=
t
%
tev
%
XS_hj_ratio
(
j
)
XS_rat
(
3
)
=
t
%
tev
%
XS_hjZ_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjbb
(
j
)
/
t
%
BR_Hbb_SM
(
j
)
BR_rat
(
2
)
=
div
(
t
%
BR_hjWW
(
j
)
,
t
%
BR_HWW_SM
(
j
)
,
0.0D0
,
1.0D9
)
case
(
10010
)
ns
=
3
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
tev
%
XS_hjW_ratio
(
j
)
XS_rat
(
2
)
=
t
%
tev
%
XS_vbf_ratio
(
j
)
XS_rat
(
3
)
=
t
%
tev
%
XS_hjZ_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjbb
(
j
)
/
t
%
BR_Hbb_SM
(
j
)
case
(
9290
)
ns
=
4
;
allocate
(
XS_rat
(
ns
))
nb
=
4
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
tev
%
XS_hjW_ratio
(
j
)
XS_rat
(
2
)
=
t
%
tev
%
XS_hj_ratio
(
j
)
XS_rat
(
3
)
=
t
%
tev
%
XS_hjZ_ratio
(
j
)
XS_rat
(
4
)
=
t
%
tev
%
XS_vbf_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjbb
(
j
)
/
t
%
BR_Hbb_SM
(
j
)
BR_rat
(
2
)
=
div
(
t
%
BR_hjWW
(
j
)
,
t
%
BR_HWW_SM
(
j
)
,
0.0D0
,
1.0D9
)
BR_rat
(
3
)
=
t
%
BR_hjtautau
(
j
)
/
t
%
BR_Htautau_SM
(
j
)
BR_rat
(
4
)
=
t
%
BR_hjgaga
(
j
)
/
t
%
BR_Hgaga_SM
(
j
)
case
(
9674
,
9897
,
9999
)
ns
=
4
;
allocate
(
XS_rat
(
ns
))
nb
=
3
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
tev
%
XS_hjW_ratio
(
j
)
XS_rat
(
2
)
=
t
%
tev
%
XS_hj_ratio
(
j
)
XS_rat
(
3
)
=
t
%
tev
%
XS_hjZ_ratio
(
j
)
XS_rat
(
4
)
=
t
%
tev
%
XS_vbf_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjbb
(
j
)
/
t
%
BR_Hbb_SM
(
j
)
BR_rat
(
2
)
=
div
(
t
%
BR_hjWW
(
j
)
,
t
%
BR_HWW_SM
(
j
)
,
0.0D0
,
1.0D9
)
BR_rat
(
3
)
=
t
%
BR_hjtautau
(
j
)
/
t
%
BR_Htautau_SM
(
j
)
case
(
7081
,
9166
,
9483
,
5586
,
9642
,
1266
,
0432
,
9891
,
5285
,
3935
,
6087
,
6170
,
10212
)
ns
=
2
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
tev
%
XS_hjW_ratio
(
j
)
XS_rat
(
2
)
=
t
%
tev
%
XS_hjZ_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjbb
(
j
)
/
t
%
BR_Hbb_SM
(
j
)
case
(
9248
,
10133
,
10439
)
ns
=
4
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
tev
%
XS_hjW_ratio
(
j
)
XS_rat
(
2
)
=
t
%
tev
%
XS_hj_ratio
(
j
)
XS_rat
(
3
)
=
t
%
tev
%
XS_hjZ_ratio
(
j
)
XS_rat
(
4
)
=
t
%
tev
%
XS_vbf_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjtautau
(
j
)
/
t
%
BR_Htautau_SM
(
j
)
case
(
4800
)
ns
=
4
;
allocate
(
XS_rat
(
ns
))
nb
=
3
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
tev
%
XS_hjW_ratio
(
j
)
XS_rat
(
2
)
=
t
%
tev
%
XS_hj_ratio
(
j
)
XS_rat
(
3
)
=
t
%
tev
%
XS_hjZ_ratio
(
j
)
XS_rat
(
4
)
=
t
%
tev
%
XS_vbf_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjtautau
(
j
)
/
t
%
BR_Htautau_SM
(
j
)
BR_rat
(
2
)
=
t
%
BR_hjbb
(
j
)
/
t
%
BR_Hbb_SM
(
j
)
BR_rat
(
3
)
=
(
t
%
BR_hjss
(
j
)
+
t
%
BR_hjcc
(
j
)
+
t
%
BR_hjbb
(
j
)
+
t
%
BR_hjgg
(
j
)
)
&
&
/
t
%
BR_Hjets_SM
(
j
)
case
(
5845
)
ns
=
4
;
allocate
(
XS_rat
(
ns
))
nb
=
2
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
tev
%
XS_hjW_ratio
(
j
)
XS_rat
(
2
)
=
t
%
tev
%
XS_hj_ratio
(
j
)
XS_rat
(
3
)
=
t
%
tev
%
XS_hjZ_ratio
(
j
)
XS_rat
(
4
)
=
t
%
tev
%
XS_vbf_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjtautau
(
j
)
/
t
%
BR_Htautau_SM
(
j
)
BR_rat
(
2
)
=
(
t
%
BR_hjss
(
j
)
+
t
%
BR_hjcc
(
j
)
+
t
%
BR_hjbb
(
j
)
+
t
%
BR_hjgg
(
j
)
)
&
&
/
t
%
BR_Hjets_SM
(
j
)
case
(
5858
,
6177
,
1887
,
10065
)
ns
=
4
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
tev
%
XS_hjW_ratio
(
j
)
XS_rat
(
2
)
=
t
%
tev
%
XS_hj_ratio
(
j
)
XS_rat
(
3
)
=
t
%
tev
%
XS_hjZ_ratio
(
j
)
XS_rat
(
4
)
=
t
%
tev
%
XS_vbf_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjgaga
(
j
)
/
t
%
BR_Hgaga_SM
(
j
)
case
(
9465
,
5871
,
9022
,
9023
,
0710
,
9887
,
5984
,
9714
,
4162
,
10102
,
6006
,
4481
,
4468
,
6095
,
10432
,
6179
)
ns
=
4
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
tev
%
XS_hjW_ratio
(
j
)
XS_rat
(
2
)
=
t
%
tev
%
XS_hj_ratio
(
j
)
XS_rat
(
3
)
=
t
%
tev
%
XS_hjZ_ratio
(
j
)
XS_rat
(
4
)
=
t
%
tev
%
XS_vbf_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjWW
(
j
)
/
t
%
BR_HWW_SM
(
j
)
case
(
6082
,
6182
)
ns
=
4
;
allocate
(
XS_rat
(
ns
))
nb
=
2
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
tev
%
XS_hjW_ratio
(
j
)
XS_rat
(
2
)
=
t
%
tev
%
XS_hj_ratio
(
j
)
XS_rat
(
3
)
=
t
%
tev
%
XS_hjZ_ratio
(
j
)
XS_rat
(
4
)
=
t
%
tev
%
XS_vbf_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjWW
(
j
)
/
t
%
BR_HWW_SM
(
j
)
BR_rat
(
2
)
=
t
%
BR_hjZZ
(
j
)
/
t
%
BR_HZZ_SM
(
j
)
case
(
6091
)
ns
=
2
;
allocate
(
XS_rat
(
ns
))
nb
=
2
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
tev
%
XS_hjW_ratio
(
j
)
XS_rat
(
2
)
=
t
%
tev
%
XS_hjZ_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjWW
(
j
)
/
t
%
BR_HWW_SM
(
j
)
BR_rat
(
2
)
=
t
%
BR_hjZZ
(
j
)
/
t
%
BR_HZZ_SM
(
j
)
case
(
6008
,
9998
)
ns
=
5
;
allocate
(
XS_rat
(
ns
))
nb
=
5
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
tev
%
XS_hjW_ratio
(
j
)
XS_rat
(
2
)
=
t
%
tev
%
XS_hj_ratio
(
j
)
XS_rat
(
3
)
=
t
%
tev
%
XS_hjZ_ratio
(
j
)
XS_rat
(
4
)
=
t
%
tev
%
XS_vbf_ratio
(
j
)
XS_rat
(
5
)
=
t
%
tev
%
XS_tthj_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjbb
(
j
)
/
t
%
BR_Hbb_SM
(
j
)
BR_rat
(
2
)
=
div
(
t
%
BR_hjWW
(
j
)
,
t
%
BR_HWW_SM
(
j
)
,
0.0D0
,
1.0D9
)
BR_rat
(
3
)
=
t
%
BR_hjtautau
(
j
)
/
t
%
BR_Htautau_SM
(
j
)
BR_rat
(
4
)
=
t
%
BR_hjgaga
(
j
)
/
t
%
BR_Hgaga_SM
(
j
)
BR_rat
(
5
)
=
(
t
%
BR_hjss
(
j
)
+
t
%
BR_hjcc
(
j
)
+
t
%
BR_hjbb
(
j
)
+
t
%
BR_hjgg
(
j
)
)
&
&
/
t
%
BR_Hjets_SM
(
j
)
case
(
6096
)
ns
=
5
;
allocate
(
XS_rat
(
ns
))
nb
=
6
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
tev
%
XS_hjW_ratio
(
j
)
XS_rat
(
2
)
=
t
%
tev
%
XS_hj_ratio
(
j
)
XS_rat
(
3
)
=
t
%
tev
%
XS_hjZ_ratio
(
j
)
XS_rat
(
4
)
=
t
%
tev
%
XS_vbf_ratio
(
j
)
XS_rat
(
5
)
=
t
%
tev
%
XS_tthj_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjbb
(
j
)
/
t
%
BR_Hbb_SM
(
j
)
BR_rat
(
2
)
=
div
(
t
%
BR_hjWW
(
j
)
,
t
%
BR_HWW_SM
(
j
)
,
0.0D0
,
1.0D9
)
BR_rat
(
3
)
=
t
%
BR_hjtautau
(
j
)
/
t
%
BR_Htautau_SM
(
j
)
BR_rat
(
4
)
=
t
%
BR_hjgaga
(
j
)
/
t
%
BR_Hgaga_SM
(
j
)
BR_rat
(
5
)
=
(
t
%
BR_hjss
(
j
)
+
t
%
BR_hjcc
(
j
)
+
t
%
BR_hjbb
(
j
)
+
t
%
BR_hjgg
(
j
)
)
&
&
/
t
%
BR_Hjets_SM
(
j
)
BR_rat
(
6
)
=
t
%
BR_hjZZ
(
j
)
/
t
%
BR_HZZ_SM
(
j
)
case
(
5757
)
ns
=
2
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
tev
%
XS_hj_ratio
(
j
)
XS_rat
(
2
)
=
t
%
tev
%
XS_vbf_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjWW
(
j
)
/
t
%
BR_HWW_SM
(
j
)
case
(
5739
)
ns
=
1
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
tev
%
XS_tthj_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjbb
(
j
)
/
t
%
BR_Hbb_SM
(
j
)
if
(
t
%
CP_value
(
j
).
eq
.
1
)
then
! analysis only applies if higgs is CP even
else
correct_properties
=
.
False
.
endif
case
(
711
)
ns
=
1
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
lep
%
XS_bbhj_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjbb
(
j
)
!note *not* normalised to SM
if
(
t
%
CP_value
(
j
).
eq
.
1
)
then
! analysis only applies if higgs is CP even
else
correct_properties
=
.
False
.
endif
case
(
713
)
ns
=
1
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
lep
%
XS_bbhj_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjbb
(
j
)
!note *not* normalised to SM
if
(
t
%
CP_value
(
j
).
eq
.
-
1
)
then
! analysis only applies if higgs is CP odd
else
correct_properties
=
.
False
.
endif
case
(
721
)
ns
=
1
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
lep
%
XS_bbhj_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjtautau
(
j
)
!note *not* normalised to SM
if
(
t
%
CP_value
(
j
).
eq
.
1
)
then
! analysis only applies if higgs is CP even
else
correct_properties
=
.
False
.
endif
case
(
723
)
ns
=
1
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
lep
%
XS_bbhj_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjtautau
(
j
)
!note *not* normalised to SM
if
(
t
%
CP_value
(
j
).
eq
.
-
1
)
then
! analysis only applies if higgs is CP odd
else
correct_properties
=
.
False
.
endif
case
(
731
)
ns
=
1
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
lep
%
XS_tautauhj_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjtautau
(
j
)
!note *not* normalised to SM
if
(
t
%
CP_value
(
j
).
eq
.
1
)
then
! analysis only applies if higgs is CP even
else
correct_properties
=
.
False
.
endif
case
(
733
)
ns
=
1
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
lep
%
XS_tautauhj_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjtautau
(
j
)
!note *not* normalised to SM
if
(
t
%
CP_value
(
j
).
eq
.
-
1
)
then
! analysis only applies if higgs is CP odd
else
correct_properties
=
.
False
.
endif
case
(
741
)
ns
=
1
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
lep
%
XS_bbhj_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjtautau
(
j
)
!note *not* normalised to SM
if
(
t
%
CP_value
(
j
).
eq
.
1
)
then
! analysis only applies if higgs is CP even
else
correct_properties
=
.
False
.
endif
case
(
743
)
ns
=
1
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
lep
%
XS_bbhj_ratio
(
j
)
BR_rat
(
1
)
=
t
%
BR_hjtautau
(
j
)
!note *not* normalised to SM
if
(
t
%
CP_value
(
j
).
eq
.
-
1
)
then
! analysis only applies if higgs is CP odd
else
correct_properties
=
.
False
.
endif
case
(
1811
)
ns
=
1
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
BR_tHpjb
(
j
)
BR_rat
(
1
)
=
t
%
BR_Hpjcs
(
j
)
if
(
(
t
%
BR_tHpjb
(
j
)
+
t
%
BR_tWpb
).
le
.
0.98D0
)
then
correct_properties
=
.
False
.
elseif
((
t
%
BR_Hpjcs
(
j
)
+
t
%
BR_Hpjtaunu
(
j
)).
le
.
0.98D0
)
then
correct_properties
=
.
False
.
endif
case
(
1812
,
8353
)
ns
=
1
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
BR_tHpjb
(
j
)
BR_rat
(
1
)
=
t
%
BR_Hpjtaunu
(
j
)
if
(
(
t
%
BR_tHpjb
(
j
)
+
t
%
BR_tWpb
).
le
.
0.98D0
)
then
correct_properties
=
.
False
.
elseif
((
t
%
BR_Hpjcs
(
j
)
+
t
%
BR_Hpjtaunu
(
j
)).
le
.
0.98D0
)
then
correct_properties
=
.
False
.
endif
case
(
1269
,
1270
)
ns
=
1
;
allocate
(
XS_rat
(
ns
))
nb
=
1
;
allocate
(
BR_rat
(
nb
))
XS_rat
(
1
)
=
t
%
BR_tHpjb
(
j
)
BR_rat
(
1
)
=
t
%
BR_Hpjcs
(
j
)
if
(
(
t
%
BR_tHpjb
(
j
)
+
t
%
BR_tWpb
).
le
.
0.98D0
)
then
correct_properties
=
.
False
.
elseif
(
t
%
BR_Hpjcs
(
j
)
.
le
.
0.98D0
)
then
correct_properties
=
.
False
.
endif
!case(801,802,803,811,813,821)
!ns = 1; allocate(XS_rat(ns))
! nb = 1; allocate(BR_rat(nb))
! XS_rat(1) =
!BR_rat(1) =
!if((t%BR_Hpjcb(j)+t%BR_Hpjcs(j)+t%BR_Hpjtaunu(j)).gt.0.98D0)then
!else
! correct_properties=.False.
!endif
case
default
write
(
*
,
*
)
'hello id='
,
id
stop
'error in subroutine model_likeness (2)'
end select
ns
=
ubound
(
XS_rat
,
dim
=
1
)
nb
=
ubound
(
BR_rat
,
dim
=
1
)
s
=
sum
(
XS_rat
)
/
ns
b
=
sum
(
BR_rat
)
/
nb
allocate
(
dsbys
(
ns
))
do is
=
1
,
ns
!dsbys(is)= abs(div((XS_rat(is) -s),s, 0.0D0,1.0D9))
dsbys
(
is
)
=
div
((
XS_rat
(
is
)
-
s
),
s
,
0.0D0
,
1.0D9
)
!want to allow e.g. a lower sigma
!to compensate for a higher Br
enddo
allocate
(
dbbyb
(
nb
))
do
ib
=
1
,
nb
!dbbyb(ib)= abs(div((BR_rat(ib) -b),b, 0.0D0,1.0D9))
dbbyb
(
ib
)
=
div
((
BR_rat
(
ib
)
-
b
),
b
,
0.0D0
,
1.0D9
)
!want to allow e.g. a lower sigma
!to compensate for a higher Br
enddo
testSMratios
=
1
!passes the SM-like ratios test
do is
=
1
,
ns
do
ib
=
1
,
nb
if
(
abs
(
dsbys
(
is
)
+
dbbyb
(
ib
)
+
dsbys
(
is
)
*
dbbyb
(
ib
)
).
gt
.
eps
)
then
testSMratios
=
-
1
!fails the SM-like ratios test
endif
enddo
enddo
if
(
testSMratios
.
lt
.
0
)
correct_properties
=
.
False
.
if
(
correct_properties
)
then
model_like
=
1
!passes the model-likeness test
sigmaXbr
=
s
*
b
else
model_like
=
-
1
!fails the model-likeness test
sigmaXbr
=
0.0D0
endif
deallocate
(
dsbys
)
deallocate
(
dbbyb
)
deallocate
(
XS_rat
)
deallocate
(
BR_rat
)
end subroutine
model_likeness
!***********************************************************
subroutine
fill_blank_ft1_dat
(
ft1
,
ft1_sep
,
vmasslower
,
vmasshigher
,
vmass_xmin
,
vmass_xmax
,
vmass_sep
,
valueoutsidetable
)
! don't forget to deallocate f_t1%dat
use
usefulbits
,
only
:
small
implicit none
integer
::
ilower
,
ihigher
double precision
,
intent
(
in
)
::
ft1_sep
,
vmasslower
,
vmasshigher
,
vmass_xmin
,
vmass_xmax
,
vmass_sep
,
valueoutsidetable
type
(
table1
)
::
ft1
if
(
abs
(
vmass_xmin
-
vmass_xmax
).
lt
.
small
)
stop
'problem in f_from_t3 (4)'
ft1
%
sep
=
ft1_sep
! we want f_t1%xmin to be lower than x1lower
if
((
vmasslower
-
vmass_xmin
).
ge
.
0.0D0
)
then
ilower
=
int
((
vmasslower
-
vmass_xmin
)
/
vmass_sep
)
+
1
else
!off lower edge of table
ilower
=
int
((
vmasslower
-
vmass_xmin
)
/
vmass_sep
)
+
1
-
1
!-1 since int rounds up for negative numbers
endif
ihigher
=
int
((
vmasshigher
-
vmass_xmin
)
/
vmass_sep
)
+
2
! we want f_t1%xmax to be higher than x1higher
ft1
%
xmin
=
dble
(
ilower
-
1
)
*
vmass_sep
+
vmass_xmin
ft1
%
xmax
=
dble
(
ihigher
-
1
)
*
vmass_sep
+
vmass_xmin
ft1
%
nx
=
nint
((
ft1
%
xmax
-
ft1
%
xmin
)
/
ft1
%
sep
)
+
1
allocate
(
ft1
%
dat
(
ft1
%
nx
,
1
))
ft1
%
dat
(:,
1
)
=
valueoutsidetable
end subroutine
fill_blank_ft1_dat
!***********************************************************
subroutine
f_from_t1
(
t1
,
vmasslower
,
vmasshigher
,
sepmultfactor
,
datcomp
,
&
&
f_t1
,
valueoutsidetable
)
! Fills the f_t1 array with the information from a t1 array
!
! Do not forget to deallocate f_t1%dat later on
!***********************************************************
use
interpolate
use
usefulbits
,
only
:
small
implicit none
!--------------------------------------input
type
(
table1
),
intent
(
in
)
::
t1
double precision
,
intent
(
in
)
::
vmasslower
,
vmasshigher
,
valueoutsidetable
double precision
,
intent
(
in
)
::
sepmultfactor
integer
,
intent
(
in
)
::
datcomp
!-----------------------------------output
type
(
table1
),
intent
(
out
)
::
f_t1
!-----------------------------------internal
integer
::
i
,
ilower
,
ihigher
double precision
::
interpol
double precision
::
vmass
,
vmass_xmin
,
vmass_xmax
,
vmass_sep
!-------------------------------------------
if
(
vmasslower
.
gt
.
vmasshigher
)
then
stop
'problem in f_from_t1 (1)'
endif
f_t1
%
id
=
t1
%
id
f_t1
%
deltax
=
t1
%
deltax
vmass_xmin
=
t1
%
xmin
vmass_xmax
=
t1
%
xmax
vmass_sep
=
t1
%
sep
f_t1
%
sep
=
t1
%
sep
*
sepmultfactor
call
fill_blank_ft1_dat
(
f_t1
,
f_t1
%
sep
,
vmasslower
,
vmasshigher
,
vmass_xmin
,
vmass_xmax
,
vmass_sep
,
valueoutsidetable
)
do
i
=
1
,
ubound
(
f_t1
%
dat
,
dim
=
1
)
vmass
=
dble
(
i
-
1
)
*
f_t1
%
sep
+
f_t1
%
xmin
if
(
vmass
.
lt
.
vmass_xmin
-
small
)
then
f_t1
%
dat
(
i
,
1
)
=
valueoutsidetable
elseif
(
vmass
.
gt
.
vmass_xmax
+
small
)
then
f_t1
%
dat
(
i
,
1
)
=
valueoutsidetable
else
call
interpolate_tabletype1
(
vmass
,
t1
,
datcomp
,
interpol
)
f_t1
%
dat
(
i
,
1
)
=
interpol
endif
enddo
end subroutine
f_from_t1
!***********************************************************
subroutine
f_from_t2
(
t2
,
m1_at_ref_point_1
,
m2_at_ref_point_1
,
m1_at_ref_point_2
,
m2_at_ref_point_2
,
&
&
vmassm1orm2
,
vmasslower
,
vmasshigher
,
sepmultfactor
,
datcomp
,
&
&
f_t1
,
valueoutsidetable
)
! Fills the f_t1 array with the information from a t2 array along a line
! m2 = line_grad*m1 + line_const
!
! Do not forget to deallocate f_t1%dat later on
!***********************************************************
use
interpolate
use
usefulbits
,
only
:
small
implicit none
!--------------------------------------input
type
(
table2
),
intent
(
in
)
::
t2
double precision
,
intent
(
in
)
::
m1_at_ref_point_1
,
m2_at_ref_point_1
,
m1_at_ref_point_2
,
m2_at_ref_point_2
double precision
,
intent
(
in
)
::
vmasslower
,
vmasshigher
,
valueoutsidetable
double precision
,
intent
(
in
)
::
sepmultfactor
integer
,
intent
(
in
)
::
datcomp
,
vmassm1orm2
!-----------------------------------output
type
(
table1
),
intent
(
out
)
::
f_t1
!-----------------------------------internal
type
(
table1
)
::
t1
double precision
::
line_grad
,
line_const
integer
::
i
,
ilower
,
ihigher
logical
::
const_m1
,
const_m2
integer
::
const_m1_i
,
const_m2_j
logical
::
on_m1_gridline
,
on_m2_gridline
double precision
::
interpol
,
mass1
,
mass2
double precision
::
m1bit
,
m2bit
double precision
::
vmass
,
vmass_xmin
,
vmass_xmax
,
vmass_sep
integer
::
ftype_selection
(
1
)
integer
::
n
!-------------------------------------------
if
(
vmasslower
.
gt
.
vmasshigher
)
then
stop
'problem in f_from_t2 (1)'
endif
if
(
abs
(
m1_at_ref_point_1
-
m1_at_ref_point_2
).
lt
.
small
)
then
const_m1
=
.
True
.
!line_grad is not needed
!line_const is not needed
else
const_m1
=
.
False
.
line_grad
=
(
m2_at_ref_point_1
-
m2_at_ref_point_2
)
/
(
m1_at_ref_point_1
-
m1_at_ref_point_2
)
line_const
=
(
m1_at_ref_point_1
*
m2_at_ref_point_2
-
m1_at_ref_point_2
*
m2_at_ref_point_1
)
&
&
/
(
m1_at_ref_point_1
-
m1_at_ref_point_2
)
endif
if
(
abs
(
m2_at_ref_point_1
-
m2_at_ref_point_2
).
lt
.
small
)
then
const_m2
=
.
True
.
else
const_m2
=
.
False
.
endif
f_t1
%
id
=
t2
%
id
f_t1
%
deltax
=
t2
%
deltax
select case
(
vmassm1orm2
)
case
(
1
)
if
(
const_m1
)
stop
'problem in f_from_t2 (3a)'
vmass_xmin
=
t2
%
xmin1
vmass_xmax
=
t2
%
xmax1
vmass_sep
=
t2
%
sep1
f_t1
%
sep
=
t2
%
sep1
*
sepmultfactor
case
(
2
)
if
(
const_m2
)
stop
'problem in f_from_t2 (3b)'
vmass_xmin
=
t2
%
xmin2
vmass_xmax
=
t2
%
xmax2
vmass_sep
=
t2
%
sep2
f_t1
%
sep
=
t2
%
sep2
*
sepmultfactor
case
default
stop
'problem in f_from_t2 (3)'
end select
call
fill_blank_ft1_dat
(
f_t1
,
f_t1
%
sep
,
vmasslower
,
vmasshigher
,
vmass_xmin
,
vmass_xmax
,
vmass_sep
,
valueoutsidetable
)
on_m1_gridline
=
.
False
.
if
(
const_m1
)
then
const_m1_i
=
nint
(
(
m1_at_ref_point_1
-
t2
%
xmin1
)
/
t2
%
sep1
)
+
1
m1bit
=
m1_at_ref_point_1
-
(
dble
(
const_m1_i
-
1
)
*
t2
%
sep1
+
t2
%
xmin1
)
/
t2
%
sep1
if
(
m1bit
.
lt
.
small
)
on_m1_gridline
=
.
True
.
endif
on_m2_gridline
=
.
False
.
if
(
const_m2
)
then
const_m2_j
=
nint
(
(
m2_at_ref_point_1
-
t2
%
xmin2
)
/
t2
%
sep2
)
+
1
m2bit
=
m2_at_ref_point_1
-
(
dble
(
const_m2_j
-
1
)
*
t2
%
sep2
+
t2
%
xmin2
)
/
t2
%
sep2
if
(
m2bit
.
lt
.
small
)
on_m2_gridline
=
.
True
.
endif
ftype_selection
(
1
)
=
datcomp
if
(
on_m1_gridline
)
then
call
fill_t1_from_t2
(
t2
,
2
,
const_m1_i
,
ftype_selection
,
t1
)
call
f_from_t1
(
t1
,
vmasslower
,
vmasshigher
,
sepmultfactor
,
datcomp
,
&
&
f_t1
,
valueoutsidetable
)
deallocate
(
t1
%
dat
)
elseif
(
on_m2_gridline
)
then
call
fill_t1_from_t2
(
t2
,
1
,
const_m2_j
,
ftype_selection
,
t1
)
call
f_from_t1
(
t1
,
vmasslower
,
vmasshigher
,
sepmultfactor
,
datcomp
,
&
&
f_t1
,
valueoutsidetable
)
deallocate
(
t1
%
dat
)
else
do
i
=
1
,
ubound
(
f_t1
%
dat
,
dim
=
1
)
vmass
=
dble
(
i
-
1
)
*
f_t1
%
sep
+
f_t1
%
xmin
if
(
t2
%
nx2
.
eq
.
1
)
then
mass1
=
vmass
mass2
=
t2
%
xmin2
elseif
(
vmassm1orm2
.
eq
.
1
)
then
mass1
=
vmass
mass2
=
mass1
*
line_grad
+
line_const
else
mass2
=
vmass
if
(
const_m1
)
then
mass1
=
m1_at_ref_point_1
else
mass1
=
(
mass2
-
line_const
)
/
line_grad
endif
endif
if
(
vmass
.
lt
.
vmass_xmin
-
small
)
then
f_t1
%
dat
(
i
,
1
)
=
valueoutsidetable
elseif
(
vmass
.
gt
.
vmass_xmax
+
small
)
then
f_t1
%
dat
(
i
,
1
)
=
valueoutsidetable
elseif
((
t2
%
needs_M2_gt_2M1
).
and
.(
2.0D0
*
mass1
>
mass2
))
then
f_t1
%
dat
(
i
,
1
)
=
valueoutsidetable
elseif
((.
not
.(
t2
%
needs_M2_gt_2M1
)).
and
.(
mass1
>
mass2
).
and
.(
t2
%
nx2
.
gt
.
1
))
then
f_t1
%
dat
(
i
,
1
)
=
valueoutsidetable
else
call
interpolate_tabletype2
(
mass1
,
mass2
,
t2
,
datcomp
,
interpol
)
f_t1
%
dat
(
i
,
1
)
=
interpol
endif
enddo
endif
end subroutine
f_from_t2
!******************************************************************
subroutine
f_from_slices_t2
(
slices_t2
,
m1_at_ref_point_1
,
m2_at_ref_point_1
,
m1_at_ref_point_2
,
m2_at_ref_point_2
,
z
,
&
&
vmassm1orm2
,
vmasslower
,
vmasshigher
,
sepmultfactor
,
datcomp
,
&
&
f_t1
,
valueoutsidetable
)
!******************************************************************
! fill the f_t1 array with the information from a t3 array at constant sf along a line
! m2 = line_grad*m1 + line_const
! do not forget to deallocate dat
use
S95tables_type3
use
interpolate
use
usefulbits
,
only
:
small
implicit none
type
(
table2
),
intent
(
in
)
::
slices_t2
(
2
)
type
(
table1
)
::
f_t1
double precision
,
intent
(
in
)
::
m1_at_ref_point_1
,
m2_at_ref_point_1
,
m1_at_ref_point_2
,
m2_at_ref_point_2
double precision
,
intent
(
in
)
::
z
,
vmasslower
,
vmasshigher
,
valueoutsidetable
double precision
,
intent
(
in
)
::
sepmultfactor
double precision
::
line_grad
,
line_const
integer
,
intent
(
in
)
::
datcomp
,
vmassm1orm2
integer
::
i
,
ilower
,
ihigher
,
a
logical
::
const_m1
,
const_m2
double precision
::
interpol
,
mass1
,
mass2
double precision
::
vmass
,
vmass_xmin
,
vmass_xmax
,
vmass_sep
integer
::
ilow
double precision
::
z_below
,
z_above
if
(
vmasslower
.
gt
.
vmasshigher
)
then
stop
'problem in f_from_slices_t2 (1)'
endif
if
(
abs
(
m1_at_ref_point_1
-
m1_at_ref_point_2
).
lt
.
small
)
then
const_m1
=
.
True
.
else
const_m1
=
.
False
.
endif
if
(
abs
(
m2_at_ref_point_1
-
m2_at_ref_point_2
).
lt
.
small
)
then
const_m2
=
.
True
.
else
const_m2
=
.
False
.
endif
! check if mass is within z range of table:
if
(
.
not
.
(
(
z
.
ge
.
slices_t2
(
1
)%
z
-
small
).
and
.(
z
.
le
.
slices_t2
(
2
)%
z
+
small
)
)
)
then
!#1! written in convoluted way to get the NaNs
f_t1
%
id
=
slices_t2
(
1
)%
id
f_t1
%
deltax
=
slices_t2
(
1
)%
deltax
if
((
slices_t2
(
1
)%
nx2
.
eq
.
1
).
or
.(
vmassm1orm2
.
eq
.
1
))
then
if
(
const_m1
)
stop
'problem in f_from_slices_t2 (1a)'
vmass_xmin
=
slices_t2
(
1
)%
xmin1
vmass_sep
=
slices_t2
(
1
)%
sep1
f_t1
%
sep
=
slices_t2
(
1
)%
sep1
*
sepmultfactor
else
if
(
const_m2
)
stop
'problem in f_from_slices_t2 (1b)'
vmass_xmin
=
slices_t2
(
1
)%
xmin2
vmass_sep
=
slices_t2
(
1
)%
sep2
f_t1
%
sep
=
slices_t2
(
1
)%
sep2
*
sepmultfactor
endif
call
fill_blank_ft1_dat
(
f_t1
,
f_t1
%
sep
,
vmasslower
,
vmasshigher
,
vmass_xmin
,
vmass_xmax
,
vmass_sep
,
valueoutsidetable
)
else
!#1
z_below
=
slices_t2
(
1
)%
z
z_above
=
slices_t2
(
2
)%
z
if
(
abs
(
z_below
-
z
).
lt
.
small
)
then
!z is the same as z_below !#2
call
f_from_t2
(
slices_t2
(
1
),
m1_at_ref_point_1
,
m2_at_ref_point_1
,
m1_at_ref_point_2
,
m2_at_ref_point_2
,
&
&
vmassm1orm2
,
vmasslower
,
vmasshigher
,
sepmultfactor
,
1
,
&
&
f_t1
,
valueoutsidetable
)
elseif
(
abs
(
z_above
-
z
).
lt
.
small
)
then
!z is the same as z_above !#2
call
f_from_t2
(
slices_t2
(
2
),
m1_at_ref_point_1
,
m2_at_ref_point_1
,
m1_at_ref_point_2
,
m2_at_ref_point_2
,
&
&
vmassm1orm2
,
vmasslower
,
vmasshigher
,
sepmultfactor
,
1
,
&
&
f_t1
,
valueoutsidetable
)
else
!#2
if
(
const_m1
)
then
!line_grad is not needed
!line_const is not needed
else
line_grad
=
(
m2_at_ref_point_1
-
m2_at_ref_point_2
)
/
(
m1_at_ref_point_1
-
m1_at_ref_point_2
)
line_const
=
(
m1_at_ref_point_1
*
m2_at_ref_point_2
-
m1_at_ref_point_2
*
m2_at_ref_point_1
)
&
&
/
(
m1_at_ref_point_1
-
m1_at_ref_point_2
)
endif
f_t1
%
id
=
slices_t2
(
1
)%
id
f_t1
%
deltax
=
slices_t2
(
1
)%
deltax
if
((
slices_t2
(
1
)%
nx2
.
eq
.
1
).
or
.(
vmassm1orm2
.
eq
.
1
))
then
vmass_xmin
=
slices_t2
(
1
)%
xmin1
vmass_xmax
=
slices_t2
(
1
)%
xmax1
vmass_sep
=
slices_t2
(
1
)%
sep1
f_t1
%
sep
=
slices_t2
(
1
)%
sep1
*
sepmultfactor
else
if
(
const_m2
)
stop
'problem in f_from_slices_t2 (3b)'
vmass_xmin
=
slices_t2
(
1
)%
xmin2
vmass_xmax
=
slices_t2
(
1
)%
xmax2
vmass_sep
=
slices_t2
(
1
)%
sep2
f_t1
%
sep
=
slices_t2
(
1
)%
sep2
*
sepmultfactor
endif
call
fill_blank_ft1_dat
(
f_t1
,
f_t1
%
sep
,
vmasslower
,
vmasshigher
,
vmass_xmin
,
vmass_xmax
,
vmass_sep
,
valueoutsidetable
)
do
i
=
1
,
ubound
(
f_t1
%
dat
,
dim
=
1
)
vmass
=
dble
(
i
-
1
)
*
f_t1
%
sep
+
f_t1
%
xmin
if
(
slices_t2
(
1
)%
nx2
.
eq
.
1
)
then
mass1
=
vmass
mass2
=
slices_t2
(
1
)%
xmin2
else
select case
(
vmassm1orm2
)
case
(
1
)
mass1
=
vmass
mass2
=
mass1
*
line_grad
+
line_const
case
(
2
)
mass2
=
vmass
if
(
const_m1
)
then
mass1
=
m1_at_ref_point_1
else
mass1
=
(
mass2
-
line_const
)
/
line_grad
endif
case
default
stop
'problem in f_from_slices_t2 (4b)'
end select
endif
if
(
vmass
.
lt
.
vmass_xmin
-
small
)
then
f_t1
%
dat
(
i
,
1
)
=
valueoutsidetable
elseif
(
vmass
.
gt
.
vmass_xmax
+
small
)
then
f_t1
%
dat
(
i
,
1
)
=
valueoutsidetable
elseif
((
slices_t2
(
1
)%
nx2
.
gt
.
1
).
and
.(
slices_t2
(
1
)%
needs_M2_gt_2M1
).
and
.(
2.0D0
*
mass1
>
mass2
))
then
f_t1
%
dat
(
i
,
1
)
=
valueoutsidetable
elseif
((
slices_t2
(
1
)%
nx2
.
gt
.
1
).
and
.(.
not
.(
slices_t2
(
1
)%
needs_M2_gt_2M1
)).
and
.(
mass1
>
mass2
))
then
f_t1
%
dat
(
i
,
1
)
=
valueoutsidetable
else
call
interpolate_slices_t2
(
mass1
,
mass2
,
z
,
slices_t2
,
datcomp
,
interpol
)
f_t1
%
dat
(
i
,
1
)
=
interpol
endif
enddo
endif
!#2
endif
!#1
end subroutine
f_from_slices_t2
!******************************************************************
subroutine
f_from_t3
(
t3
,
m1_at_ref_point_1
,
m2_at_ref_point_1
,
m1_at_ref_point_2
,
m2_at_ref_point_2
,
z
,
&
&
vmassm1orm2
,
vmasslower
,
vmasshigher
,
sepmultfactor
,
datcomp
,
&
&
f_t1
,
valueoutsidetable
)
!******************************************************************
! fill the f_t1 array with the information from a t3 array at constant sf along a line
! m2 = line_grad*m1 + line_const
! do not forget to deallocate dat
use
S95tables_type3
use
interpolate
use
usefulbits
,
only
:
small
implicit none
type
(
table3
),
intent
(
in
)
::
t3
type
(
table2
)
::
slices_t2
(
2
)
type
(
table1
)
::
f_t1
double precision
,
intent
(
in
)
::
m1_at_ref_point_1
,
m2_at_ref_point_1
,
m1_at_ref_point_2
,
m2_at_ref_point_2
double precision
,
intent
(
in
)
::
z
,
vmasslower
,
vmasshigher
,
valueoutsidetable
double precision
,
intent
(
in
)
::
sepmultfactor
double precision
::
line_grad
,
line_const
integer
,
intent
(
in
)
::
datcomp
,
vmassm1orm2
integer
::
i
,
ilower
,
ihigher
,
a
logical
::
const_m1
,
const_m2
double precision
::
interpol
,
mass1
,
mass2
double precision
::
vmass
,
vmass_xmin
,
vmass_xmax
,
vmass_sep
integer
::
ilow
,
c_zi
(
2
),
ftype_selection
(
1
)
double precision
::
z_below
,
z_above
if
(
vmasslower
.
gt
.
vmasshigher
)
then
stop
'problem in f_from_t3 (1)'
endif
if
(
abs
(
m1_at_ref_point_1
-
m1_at_ref_point_2
).
lt
.
small
)
then
const_m1
=
.
True
.
else
const_m1
=
.
False
.
endif
if
(
abs
(
m2_at_ref_point_1
-
m2_at_ref_point_2
).
lt
.
small
)
then
const_m2
=
.
True
.
else
const_m2
=
.
False
.
endif
! check if mass is within z range of table:
if
(
.
not
.
(
(
z
.
ge
.
t3
%
zmin
-
small
).
and
.(
z
.
le
.
t3
%
zmax
+
small
)
)
)
then
!#1! written in convoluted way to get the NaNs
f_t1
%
id
=
t3
%
id
f_t1
%
deltax
=
t3
%
deltax
if
((
t3
%
nx2
.
eq
.
1
).
or
.(
vmassm1orm2
.
eq
.
1
))
then
if
(
const_m1
)
stop
'problem in f_from_t3 (1a)'
vmass_xmin
=
t3
%
xmin1
vmass_sep
=
t3
%
sep1
f_t1
%
sep
=
t3
%
sep1
*
sepmultfactor
else
if
(
const_m2
)
stop
'problem in f_from_t3 (1b)'
vmass_xmin
=
t3
%
xmin2
vmass_sep
=
t3
%
sep2
f_t1
%
sep
=
t3
%
sep2
*
sepmultfactor
endif
call
fill_blank_ft1_dat
(
f_t1
,
f_t1
%
sep
,
vmasslower
,
vmasshigher
,
vmass_xmin
,
vmass_xmax
,
vmass_sep
,
valueoutsidetable
)
else
!#1
ilow
=
int
((
z
-
t3
%
zmin
)
/
t3
%
zsep
)
+
1
z_below
=
dble
(
ilow
-
1
)
*
t3
%
zsep
+
t3
%
zmin
z_above
=
z_below
+
t3
%
zsep
if
(
abs
(
z_below
-
z
).
lt
.
small
)
then
!z is the same as z_below !#2
c_zi
=
ilow
elseif
(
abs
(
z_above
-
z
).
lt
.
small
)
then
!z is the same as z_above !#2
c_zi
=
ilow
+
1
else
!#2
c_zi
(
1
)
=
ilow
c_zi
(
2
)
=
ilow
+
1
endif
!#2
ftype_selection
(
1
)
=
datcomp
call
fill_slices_t2_from_slices_of_t3
(
t3
,
c_zi
,
ftype_selection
,
slices_t2
)
call
f_from_slices_t2
(
slices_t2
,
m1_at_ref_point_1
,
m2_at_ref_point_1
,
m1_at_ref_point_2
,
m2_at_ref_point_2
,
z
,
&
&
vmassm1orm2
,
vmasslower
,
vmasshigher
,
sepmultfactor
,
datcomp
,
&
&
f_t1
,
valueoutsidetable
)
do
a
=
1
,
2
deallocate
(
slices_t2
(
a
)%
dat
)
enddo
endif
!#1
end subroutine
f_from_t3
!************************************************************
subroutine
convolve_chisq_with_gaussian
(
t1
,
datcomp
,
sigma
,
mass
,
result
)
!************************************************************
! intergrate exp(-t1%dat(xi,1)/2)*exp(-(massx-mass)^2/(2*sigma^2))/sqrt(2*pi*sigma^2) w.r.t. x
! between xlower and xhigher
! then do -2.0D0*log to get result
! negative data points are invalid. They are set to zero.
use
usefulbits
,
only
:
vsmall
,
vvsmall
,
pi
!internal
use
interpolate
use
S95tables_type1
implicit none
type
(
table1
),
intent
(
in
)
::
t1
integer
,
intent
(
in
)
::
datcomp
double precision
,
intent
(
in
)
::
sigma
,
mass
double precision
,
intent
(
out
)
::
result
!-----------------------------------internal
integer
::
i
,
ilow
,
ihigh
,
j
,
divisions
,
d
,
n
,
ntot
,
a
double precision
::
runningtotal
,
massx
,
datvalue
,
newsep
double precision
,
allocatable
::
newdat
(:)
double precision
::
big_number_instead_of_infinity
double precision
::
dati
,
datiplus1
!-------------------------------------------
if
((
datcomp
.
lt
.
lbound
(
t1
%
dat
,
dim
=
2
)).
or
.(
datcomp
.
gt
.
ubound
(
t1
%
dat
,
dim
=
2
)))
then
stop
'wrong datcomp inputted to subroutine convolve_with_gaussian'
elseif
(
t1
%
nx
.
le
.
1
)
then
stop
'wrong t1%nx inputted to subroutine convolve_with_gaussian (2)'
elseif
(
sigma
.
le
.
vsmall
)
then
stop
'wrong sigma inputted to subroutine convolve_with_gaussian'
elseif
(
abs
(
t1
%
sep
).
le
.
vsmall
)
then
stop
'wrong t1%sep inputted to subroutine convolve_with_gaussian'
endif
big_number_instead_of_infinity
=
1.0D5
divisions
=
5
!do i=1,t1%nx
! if(t1%dat(i,datcomp).ge.big_number_instead_of_infinity)t1%dat(i,datcomp)=1.0D20
!enddo
n
=
0
if
(
minval
(
t1
%
dat
(:,
datcomp
)).
lt
.
1.0D4
)
then
ilow
=
lbound
(
t1
%
dat
,
dim
=
1
)
ihigh
=
ubound
(
t1
%
dat
,
dim
=
1
)
if
(
ilow
.
eq
.
ihigh
)
stop
'problem in subroutine convolve_with_gaussian'
newsep
=
t1
%
sep
/
dble
(
divisions
)
ntot
=
divisions
*
(
ihigh
-
ilow
)
+
1
allocate
(
newdat
(
ntot
))
newdat
=
0.0D0
do
i
=
ilow
,
ihigh
dati
=
t1
%
dat
(
i
,
datcomp
)
if
(
dati
.
ge
.
0.0D0
)
then
n
=
n
+
1
massx
=
dble
(
i
-
1
)
*
t1
%
sep
+
t1
%
xmin
datvalue
=
dati
newdat
(
n
)
=
exp
(
-
datvalue
/
2.0D0
)
&
&
*
exp
(
-
(
massx
-
mass
)
**
2.0D0
/
(
2.0D0
*
sigma
**
2.0D0
))
/
sqrt
(
2.0D0
*
pi
*
sigma
**
2.0D0
)
if
(
i
.
lt
.
ihigh
)
then
datiplus1
=
t1
%
dat
(
i
+
1
,
datcomp
)
if
(
datiplus1
.
ge
.
0.0D0
)
then
do
j
=
2
,
divisions
-
1
n
=
n
+
1
massx
=
dble
(
i
-
1
)
*
t1
%
sep
+
t1
%
xmin
+
dble
(
j
-
1
)
*
newsep
!do a=1,ihigh
! write(*,*)a,dble(a-1)*t1%sep+t1%xmin ,t1%dat(a,datcomp)
!enddo
datvalue
=
dati
+
((
datiplus1
-
dati
)
/
t1
%
sep
)
*
dble
(
j
-
1
)
*
newsep
if
(
datvalue
.
lt
.
0.0D0
)
then
!these are invalid point or places outside range of table
datvalue
=
0.0D0
endif
newdat
(
n
)
=
exp
(
-
datvalue
/
2.0D0
)
&
&
*
exp
(
-
(
massx
-
mass
)
**
2.0D0
/
(
2.0D0
*
sigma
**
2.0D0
))
/
sqrt
(
2.0D0
*
pi
*
sigma
**
2.0D0
)
enddo
else
do
j
=
2
,
divisions
-
1
n
=
n
+
1
enddo
endif
else
!negative data points are invalid
do
j
=
2
,
divisions
-
1
n
=
n
+
1
enddo
endif
!massx=dble(i-1)*t1%sep+t1%xmin
!newdat(i)=exp(-t1%dat(i,datcomp)/2.0D0) &
! & *exp(-(massx-mass)**2.0D0/(2.0D0*sigma**2.0D0))/sqrt(2.0D0*pi*sigma**2.0D0)
else
do
j
=
1
,
divisions
-
1
n
=
n
+
1
enddo
endif
enddo
!intergrate with trapezium rule
runningtotal
=
0.5D0
*
(
newdat
(
1
)
+
newdat
(
ntot
))
if
((
ntot
).
gt
.
1
)
then
do
n
=
2
,
ntot
-
1
runningtotal
=
runningtotal
+
newdat
(
n
)
enddo
endif
deallocate
(
newdat
)
if
(
abs
(
runningtotal
).
le
.
vvsmall
)
then
result
=
big_number_instead_of_infinity
else
result
=
-
2.0D0
*
log
(
runningtotal
*
newsep
)
endif
if
(
result
.
gt
.
2
2.4D0
)
then
!corresponds to clsb=1.0D-6, which is the lowest clsb that ppchi2 can take as input
result
=
big_number_instead_of_infinity
endif
else
result
=
big_number_instead_of_infinity
endif
end subroutine
convolve_chisq_with_gaussian
!************************************************************
function
S95_t1_or_S95_t2_idfromelementnumber
(
ttype
,
tlist
)
!************************************************************
implicit none
integer
::
S95_t1_or_S95_t2_idfromelementnumber
integer
,
intent
(
in
)
::
tlist
integer
,
intent
(
in
)
::
ttype
select case
(
ttype
)
case
(
1
)
S95_t1_or_S95_t2_idfromelementnumber
=
S95_t1
(
tlist
)%
id
case
(
2
)
S95_t1_or_S95_t2_idfromelementnumber
=
S95_t2
(
tlist
)%
id
case
default
stop
'wrong input to function S95_t1_or_S95_t2_idfromelementnumber'
end select
end function
S95_t1_or_S95_t2_idfromelementnumber
!************************************************************
function
S95_t1_or_S95_t2_elementnumberfromid
(
ttype
,
id
)
!************************************************************
use
S95tables_type1
,
only
:
t1elementnumberfromid
use
S95tables_type2
,
only
:
t2elementnumberfromid
implicit none
integer
,
intent
(
in
)
::
id
integer
,
intent
(
in
)
::
ttype
integer
::
S95_t1_or_S95_t2_elementnumberfromid
select case
(
ttype
)
case
(
1
)
S95_t1_or_S95_t2_elementnumberfromid
=
t1elementnumberfromid
(
S95_t1
,
id
)
case
(
2
)
S95_t1_or_S95_t2_elementnumberfromid
=
t2elementnumberfromid
(
S95_t2
,
id
)
case
default
stop
'problem with function S95_t1_or_S95_t2_elementnumberfromid'
end select
end function
S95_t1_or_S95_t2_elementnumberfromid
!************************************************************
subroutine
deallocate_S95tables
!************************************************************
implicit none
!-----------------------------------internal
integer
x
!-------------------------------------------
do
x
=
lbound
(
S95_t1
,
dim
=
1
),
ubound
(
S95_t1
,
dim
=
1
)
deallocate
(
S95_t1
(
x
)%
dat
)
enddo
do
x
=
lbound
(
S95_t2
,
dim
=
1
),
ubound
(
S95_t2
,
dim
=
1
)
deallocate
(
S95_t2
(
x
)%
dat
)
enddo
deallocate
(
S95_t1
)
deallocate
(
S95_t2
)
end subroutine
deallocate_S95tables
!***********************************************************
end module
S95tables
!************************************************************
File Metadata
Details
Attached
Mime Type
text/plain
Expires
Wed, May 14, 11:41 AM (10 h, 49 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5111487
Default Alt Text
S95tables.f90 (75 KB)
Attached To
rHIGGSBOUNDSSVN higgsboundssvn
Event Timeline
Log In to Comment