Page MenuHomeHEPForge
Files F7879424

No OneTemporary
Actions

View Options

Index: trunk/HiggsBounds_KW/HiggsBounds_subroutines.F90
===================================================================
--- trunk/HiggsBounds_KW/HiggsBounds_subroutines.F90 (revision 500)
+++ trunk/HiggsBounds_KW/HiggsBounds_subroutines.F90 (revision 501)
@@ -1,1719 +1,1804 @@
! This file is part of HiggsBounds
! -KW
!************************************************************
subroutine initialize_HiggsBounds(nHiggsneut,nHiggsplus,whichanalyses_in)
! This the first Higgsbounds subroutine that should be called
! by the user.
! It calls subroutines to read in the tables of Standard Model data,
! read in the tables of LEP, Tevatron and LHC data,
! set up lists of processes which should be checked against
! the experimental results, allocate arrays etc
! Arguments (input):
! * nHiggs= number of neutral Higgs in the model
! (see subroutine check_nH_nHplus in input.f90 for more details)
! * nHiggsplus= number of singly,positively charged Higgs in the model
! (see subroutine check_nH_nHplus in input.f90 for more details)
! * whichanalyses_in= which combination of experimental results to use
! (see subroutine check_whichanalyses in input.f90 for more details)
!************************************************************
use usefulbits, only : np,Hneut,Hplus,Chineut,Chiplus,debug,inputmethod, &
& inputsub,theo,whichanalyses,HiggsBounds_info,just_after_run,&
& file_id_debug1,file_id_debug2,allocate_if_stats_required,run_HB_classic
use input, only : setup_input,check_number_of_particles,check_whichanalyses
use S95tables, only : setup_S95tables,S95_t2
use likelihoods, only : setup_likelihoods
use theory_BRfunctions, only : setup_BRSM
use channels, only : setup_channels
use output, only : setup_output
#ifdef enableCHISQ
use S95tables_type3, only : clsb_t3,fillt3needs_M2_gt_2M1
#endif
#if defined(NAGf90Fortran)
use F90_UNIX_IO, only : flush
#endif
!#define FORFITTINO
implicit none
!--------------------------------------input
integer,intent(in) :: nHiggsneut
! integer,intent(in),optional :: nHiggsplus
! character(LEN=5),intent(in),optional :: whichanalyses_in
integer,intent(in) :: nHiggsplus
character(LEN=5),intent(in) :: whichanalyses_in
!-----------------------------------internal
integer :: i
logical :: messages
!-------------------------------------------
! if((.not.present(nHiggsplus)).or.(.not.present(whichanalyses_in)))then
!Actually, this doesn't work as I wanted it to
!because if initialize_HiggsBounds is called in the old way, the program
!usually just crashes..... but leaving it in for now, in case
!some compilers accept it
! call attempting_to_use_an_old_HB_version('init')
! endif
#ifdef FORFITTINO
write(*,*)'The arguments passed to initialize_HiggsBounds are:'
write(*,*)'nHiggsneut=',nHiggsneut
write(*,*)'nHiggsplus=',nHiggsplus
write(*,*)'whichanalyses_in=','~'//trim(adjustl(whichanalyses_in))//'~'
#endif
#ifdef DEBUGGING
debug=.True.
#else
debug=.False.
#endif
messages=debug.or.(inputmethod=='datfile')
! inputmethod='subrout' !('datfile' or 'website' are also possible, but not here)
np(Hneut)=nHiggsneut
np(Hplus)=nHiggsplus
np(Chineut)=0! do not change this without contacting us first!
np(Chiplus)=0! do not change this without contacting us first!
whichanalyses=whichanalyses_in
if(inputmethod=='subrout') then
if(allocated(theo))then
stop 'subroutine HiggsBounds_initialize has already been called once'
endif
if(messages)write(*,*)'doing other preliminary tasks...' ; call flush(6)
call setup_input
allocate(inputsub( 4 )) !(1)np(Hneut)>0 (2)np(Hplus)>0 (3)np(Chineut)>0 (4)np(Chineut)>0 and np(Chiplus)>0
! | np
! |Hneu Hcha Chineut Chiplus
! | ==0 ==0 ==0 ==0
inputsub(1)%desc='HiggsBounds_neutral_input_*'
inputsub(1)%req=req( 0, 1, 1, 1)
inputsub(2)%desc='HiggsBounds_charged_input'
inputsub(2)%req=req( 1, 0, 1, 1)
inputsub(3)%desc='SUSYBounds_neutralinoonly_input'
inputsub(3)%req=req( 1, 1, 0, 1)
inputsub(4)%desc='SUSYBounds_neutralinochargino_input'
inputsub(4)%req=req( 1, 1, 0, 0)
do i=1,ubound(inputsub,dim=1)
inputsub(i)%stat=0
enddo
endif
#ifndef WEBVERSION
if(inputmethod.ne.'datfile') call HiggsBounds_info
if (run_HB_classic.EQV..True.) then
PRINT *, "run_HB_classic=True - HiggsBounds is running in classic mode"
endif
#endif
if(messages)write(*,*)'reading in Standard Model tables...' ; call flush(6)
call setup_BRSM
if(messages)write(*,*)'reading in S95tables...' ; call flush(6)
call setup_S95tables
if(messages)write(*,*)'reading in likelihoods...' ; call flush(6)
call setup_likelihoods
!if(debug)write(*,*)'doing other preliminary tasks...' ; call flush(6)
!call setup_input
if(messages)then
open(file_id_debug2,file='debug_predratio.txt')
open(file_id_debug1,file='debug_channels.txt')
endif
if(messages)write(*,*)'sorting out processes to be checked...'; call flush(6)
call setup_channels
if(messages)write(*,*)'preparing output arrays...' ; call flush(6)
call setup_output
#ifdef enableCHISQ
if(allocated(allocate_if_stats_required))then
call fillt3needs_M2_gt_2M1(clsb_t3,S95_t2)
endif
#endif
just_after_run=.False.
contains
! | np
! |Hneu Hcha Chineut Chiplus
! | ==0 ==0 ==0 ==0
function req(Hneu,Hcha, Chneu, Chcha)
integer, intent(in) ::Hneu,Hcha, Chneu, Chcha
integer :: req
req=1
if(np(Hneut)==0) req= Hneu * req
if(np(Hplus)==0) req= Hcha * req
if(np(Chineut)==0)req= Chneu * req
if(np(Chiplus)==0)req= Chcha * req
end function req
end subroutine initialize_HiggsBounds
!************************************************************
!************************************************************
! Version of initialize_HiggsBounds which takes an integer as
! the third argument. More useful for library linking to
! non-Fortran codes.
subroutine initialize_HiggsBounds_int(nHn,nHp,flag)
implicit none
integer nHn,nHp,flag
interface
subroutine initialize_HiggsBounds(nHiggsneut, nHiggsplus, whichanalyses_in)
integer,intent(in) :: nHiggsneut
integer,intent(in),optional :: nHiggsplus
character(LEN=5),intent(in),optional :: whichanalyses_in
end subroutine initialize_HiggsBounds
end interface
IF (flag.EQ.1) then
call initialize_HiggsBounds(nHn,nHp, "onlyL")
elseif (flag.EQ.2) then
call initialize_HiggsBounds(nHn,nHp, "onlyH")
elseif (flag.EQ.3) then
call initialize_HiggsBounds(nHn,nHp, "LandH")
elseif (flag.EQ.4) then
call initialize_HiggsBounds(nHn,nHp, "onlyP")
else
stop "Illegal value for flag in call to initialize_HB"
endif
end subroutine
!************************************************************
!************************************************************
subroutine attempting_to_use_an_old_HB_version(subroutineid)
use usefulbits, only : vers
character(len=4),intent(in) :: subroutineid
select case(subroutineid)
case('init')
write(*,*)'The subroutine initialize_HiggsBounds has been called with the'
write(*,*)'wrong number of arguments. It should be called as:'
write(*,*)'initialize_HiggsBounds(nHiggsneut,nHiggsplus,whichanalyses)'
write(*,*)
write(*,*)'Note that in early versions of HiggsBounds (HB 1.*.*)'
write(*,*)'this subroutine was called as:'
write(*,*)'initialize_HiggsBounds(nHiggsneut,whichanalyses)'
write(*,*)
case('effC','part','hadr')
write(*,*)'The subroutine run_HiggsBounds_'//subroutineid//' has been discontinued in this'
write(*,*)'version of HiggsBounds.'
case default
stop'wrong input to subroutine attempting_to_use_an_old_HB_version'
end select
write(*,*)'If you have code written for use with HB 1.*.*, you have two choices:'
write(*,*)
write(*,*)' (1) You can edit your code, such that it works with this'
write(*,*)' version of HiggsBounds (HB'//trim(adjustl(vers))//').'
write(*,*)' This has the advantage that you can test your model against many, many'
write(*,*)' more Higgs search limits , including charged Higgs search limits.'
write(*,*)' See the updated manual for more information.'
write(*,*)
write(*,*)' (2) You can download the most recent HB 1.*.* from the HiggsBounds'
write(*,*)' website. This contains the LEP Higgs search limits which are'
write(*,*)' generally the most useful when constraining new physics models.'
write(*,*)' We will continue to support this code.'
stop'Incorrect call to a HiggsBounds subroutine.'
end subroutine attempting_to_use_an_old_HB_version
!************************************************************
subroutine HiggsBounds_input_SLHA(infile)
! This subroutine can be called by the user after subroutine initialize_HiggsBounds
! has been called.
! Arguments (input): SLHA filename
!************************************************************
use usefulbits, only : whichinput,inputsub,infile1,theo,g2,just_after_run, &
& np,Hneut,Hplus
use extra_bits_for_SLHA
#if defined(NAGf90Fortran)
use F90_UNIX_IO, only : flush
#endif
implicit none
!----------------------------------------input
character(len=100),intent(in) :: infile
!--------------------------------------internal
integer :: n
!----------------------------------------------
whichinput='SLHA'
if(np(Hneut).gt.0)inputsub(Hneut)%stat=inputsub(Hneut)%stat+1
if(np(Hplus).gt.0)inputsub(Hplus)%stat=inputsub(Hplus)%stat+1
! note: can't be used for charginos or neutralinos yet
n=1
if(.not.allocated(theo))then
stop 'subroutine HiggsBounds_initialize must be called first'
endif
infile1=infile
call getSLHAdata(theo(n),g2(n),infile1)
just_after_run=.False.
end subroutine HiggsBounds_input_SLHA
!************************************************************
subroutine HiggsBounds_neutral_input_effC(Mh,GammaTotal_hj, &
& g2hjss_s,g2hjss_p,g2hjcc_s,g2hjcc_p, &
& g2hjbb_s,g2hjbb_p,g2hjtoptop_s,g2hjtoptop_p, &
& g2hjmumu_s,g2hjmumu_p, &
& g2hjtautau_s,g2hjtautau_p, &
& g2hjWW,g2hjZZ,g2hjZga, &
& g2hjgaga,g2hjgg,g2hjggZ,g2hjhiZ_nHbynH, &
& BR_hjinvisible,BR_hjhihi_nHbynH )
! This subroutine can be called by the user after subroutine initialize_HiggsBounds
! has been called.
! Arguments (input): theoretical predictions (see manual for definitions)
!************************************************************
use usefulbits, only : theo,np,Hneut,g2,whichinput,inputsub,just_after_run
#if defined(NAGf90Fortran)
use F90_UNIX_IO, only : flush
#endif
implicit none
!----------------------------------------input
double precision,intent(in) :: Mh( np(Hneut) ),GammaTotal_hj( np(Hneut) ), &
& g2hjss_s( np(Hneut) ),g2hjss_p( np(Hneut) ),g2hjcc_s( np(Hneut) ),g2hjcc_p( np(Hneut) ), &
& g2hjbb_s( np(Hneut) ),g2hjbb_p( np(Hneut) ),g2hjtoptop_s( np(Hneut) ),g2hjtoptop_p( np(Hneut) ),&
& g2hjmumu_s( np(Hneut) ),g2hjmumu_p( np(Hneut) ), &
& g2hjtautau_s( np(Hneut) ),g2hjtautau_p( np(Hneut) ), &
& g2hjWW( np(Hneut) ),g2hjZZ( np(Hneut) ),g2hjZga( np(Hneut) ), &
& g2hjgaga( np(Hneut) ),g2hjgg( np(Hneut) ),g2hjggZ( np(Hneut) ),g2hjhiZ_nHbynH(np(Hneut),np(Hneut)),&
& BR_hjinvisible( np(Hneut) ),BR_hjhihi_nHbynH(np(Hneut),np(Hneut))
!--------------------------------------internal
integer :: n
integer :: subtype
!----------------------------------------------
whichinput='effC'
subtype=1
n=1
inputsub(subtype)%stat=inputsub(subtype)%stat+1
if(.not.allocated(theo))then
stop 'subroutine HiggsBounds_initialize must be called first'
endif
if(np(Hneut).eq.0)then
write(*,*)'subroutine HiggsBounds_neutral_input_effC should'
write(*,*)'only be called if np(Hneut)>0'
stop 'error in subroutine HiggsBounds_neutral_input_effC'
endif
theo(n)%particle(Hneut)%M = Mh
theo(n)%particle(Hneut)%Mc = Mh
theo(n)%particle(Hneut)%GammaTot= GammaTotal_hj
g2(n)%hjss_s = g2hjss_s
g2(n)%hjss_p = g2hjss_p
g2(n)%hjcc_s = g2hjcc_s
g2(n)%hjcc_p = g2hjcc_p
g2(n)%hjbb_s = g2hjbb_s
g2(n)%hjbb_p = g2hjbb_p
g2(n)%hjtoptop_s = g2hjtoptop_s
g2(n)%hjtoptop_p = g2hjtoptop_p
g2(n)%hjmumu_s = g2hjmumu_s
g2(n)%hjmumu_p = g2hjmumu_p
g2(n)%hjtautau_s = g2hjtautau_s
g2(n)%hjtautau_p = g2hjtautau_p
g2(n)%hjWW = g2hjWW
g2(n)%hjZZ = g2hjZZ
g2(n)%hjZga = g2hjZga
g2(n)%hjgaga = g2hjgaga
g2(n)%hjgg = g2hjgg
g2(n)%hjggZ = g2hjggZ
g2(n)%hjhiZ = g2hjhiZ_nHbynH
theo(n)%BR_hjinvisible = BR_hjinvisible
theo(n)%BR_hjhihi = BR_hjhihi_nHbynH
just_after_run=.False.
end subroutine HiggsBounds_neutral_input_effC
!************************************************************
subroutine HiggsBounds_neutral_input_part(Mh,GammaTotal_hj,CP_value, &
& CS_lep_hjZ_ratio, &
& CS_lep_bbhj_ratio,CS_lep_tautauhj_ratio, &
& CS_lep_hjhi_ratio_nHbynH, &
& CS_gg_hj_ratio,CS_bb_hj_ratio, &
& CS_bg_hjb_ratio, &
& CS_ud_hjWp_ratio,CS_cs_hjWp_ratio, &
& CS_ud_hjWm_ratio,CS_cs_hjWm_ratio, &
& CS_gg_hjZ_ratio, &
& CS_dd_hjZ_ratio,CS_uu_hjZ_ratio, &
& CS_ss_hjZ_ratio,CS_cc_hjZ_ratio, &
& CS_bb_hjZ_ratio, &
& CS_tev_vbf_ratio,CS_tev_tthj_ratio, &
& CS_lhc7_vbf_ratio,CS_lhc7_tthj_ratio, &
& CS_lhc8_vbf_ratio,CS_lhc8_tthj_ratio, &
& BR_hjss,BR_hjcc, &
& BR_hjbb,BR_hjmumu,BR_hjtautau, &
& BR_hjWW,BR_hjZZ,BR_hjZga, BR_hjgaga,BR_hjgg, &
& BR_hjinvisible,BR_hjhihi_nHbynH )
! This subroutine can be called by the user after subroutine initialize_HiggsBounds
! has been called.
! (see manual for full description)
!************************************************************
use usefulbits, only : theo,np,Hneut,partR,whichinput,inputsub,just_after_run
#if defined(NAGf90Fortran)
use F90_UNIX_IO, only : flush
#endif
implicit none
!----------------------------------------input
double precision,intent(in) :: Mh( np(Hneut) ),GammaTotal_hj( np(Hneut) )
integer,intent(in) ::CP_value( np(Hneut) )
double precision,intent(in) :: CS_lep_hjZ_ratio( np(Hneut) ), &
& CS_lep_bbhj_ratio( np(Hneut) ),CS_lep_tautauhj_ratio( np(Hneut) ), &
& CS_lep_hjhi_ratio_nHbynH(np(Hneut),np(Hneut)), &
& CS_gg_hj_ratio( np(Hneut) ),CS_bb_hj_ratio( np(Hneut) ), &
& CS_bg_hjb_ratio( np(Hneut) ), &
& CS_ud_hjWp_ratio( np(Hneut) ),CS_cs_hjWp_ratio( np(Hneut) ), &
& CS_ud_hjWm_ratio( np(Hneut) ),CS_cs_hjWm_ratio( np(Hneut) ), &
& CS_gg_hjZ_ratio( np(Hneut) ), &
& CS_dd_hjZ_ratio( np(Hneut) ),CS_uu_hjZ_ratio( np(Hneut) ), &
& CS_ss_hjZ_ratio( np(Hneut) ),CS_cc_hjZ_ratio( np(Hneut) ), &
& CS_bb_hjZ_ratio( np(Hneut) ), &
& CS_tev_vbf_ratio( np(Hneut) ),CS_tev_tthj_ratio( np(Hneut) ), &
& CS_lhc7_vbf_ratio( np(Hneut) ),CS_lhc7_tthj_ratio( np(Hneut) ), &
& CS_lhc8_vbf_ratio( np(Hneut) ),CS_lhc8_tthj_ratio( np(Hneut) ), &
& BR_hjss( np(Hneut) ),BR_hjcc( np(Hneut) ), &
& BR_hjbb( np(Hneut) ),BR_hjmumu( np(Hneut) ),BR_hjtautau( np(Hneut) ), &
& BR_hjWW( np(Hneut) ),BR_hjZZ( np(Hneut) ),BR_hjZga( np(Hneut) ), &
& BR_hjgaga( np(Hneut) ),BR_hjgg( np(Hneut) ), &
& BR_hjinvisible( np(Hneut) ),BR_hjhihi_nHbynH(np(Hneut),np(Hneut))
!---------------------------------------internal
integer :: n
integer :: subtype
!-----------------------------------------------
whichinput='part'
subtype=1
n=1
inputsub(subtype)%stat=inputsub(subtype)%stat+1
if(.not.allocated(theo))then
stop 'subroutine HiggsBounds_initialize must be called first'
endif
if(np(Hneut).eq.0)then
write(*,*)'subroutine HiggsBounds_neutral_input_part should'
write(*,*)'only be called if np(Hneut)>0'
stop 'error in subroutine HiggsBounds_neutral_input_part'
endif
theo(n)%particle(Hneut)%M = Mh
theo(n)%particle(Hneut)%Mc = Mh
theo(n)%particle(Hneut)%GammaTot= GammaTotal_hj
theo(n)%CP_value = CP_value
theo(n)%lep%XS_hjZ_ratio = CS_lep_hjZ_ratio
theo(n)%lep%XS_bbhj_ratio = CS_lep_bbhj_ratio
theo(n)%lep%XS_tautauhj_ratio = CS_lep_tautauhj_ratio
theo(n)%lep%XS_hjhi_ratio = CS_lep_hjhi_ratio_nHbynH
partR(n)%gg_hj = CS_gg_hj_ratio
partR(n)%qq_hj(5,:) = CS_bb_hj_ratio
partR(n)%bg_hjb = CS_bg_hjb_ratio
partR(n)%qq_hjWp(1,:) = CS_ud_hjWp_ratio
partR(n)%qq_hjWp(2,:) = CS_cs_hjWp_ratio
partR(n)%qq_hjWm(1,:) = CS_ud_hjWm_ratio
partR(n)%qq_hjWm(2,:) = CS_cs_hjWm_ratio
partR(n)%gg_hjZ(:) = CS_gg_hjZ_ratio
partR(n)%qq_hjZ(1,:) = CS_dd_hjZ_ratio
partR(n)%qq_hjZ(2,:) = CS_uu_hjZ_ratio
partR(n)%qq_hjZ(3,:) = CS_ss_hjZ_ratio
partR(n)%qq_hjZ(4,:) = CS_cc_hjZ_ratio
partR(n)%qq_hjZ(5,:) = CS_bb_hjZ_ratio
theo(n)%tev%XS_vbf_ratio = CS_tev_vbf_ratio
theo(n)%tev%XS_tthj_ratio = CS_tev_tthj_ratio
theo(n)%lhc7%XS_vbf_ratio = CS_lhc7_vbf_ratio
theo(n)%lhc7%XS_tthj_ratio= CS_lhc7_tthj_ratio
theo(n)%lhc8%XS_vbf_ratio = CS_lhc8_vbf_ratio
theo(n)%lhc8%XS_tthj_ratio= CS_lhc8_tthj_ratio
theo(n)%BR_hjss = BR_hjss
theo(n)%BR_hjcc = BR_hjcc
theo(n)%BR_hjbb = BR_hjbb
theo(n)%BR_hjmumu = BR_hjmumu
theo(n)%BR_hjtautau = BR_hjtautau
theo(n)%BR_hjWW = BR_hjWW
theo(n)%BR_hjZZ = BR_hjZZ
theo(n)%BR_hjZga = BR_hjZga
theo(n)%BR_hjgaga = BR_hjgaga
theo(n)%BR_hjgg = BR_hjgg
theo(n)%BR_hjinvisible = BR_hjinvisible
theo(n)%BR_hjhihi = BR_hjhihi_nHbynH
just_after_run=.False.
end subroutine HiggsBounds_neutral_input_part
!************************************************************
subroutine HiggsBounds_neutral_input_hadr(Mh,GammaTotal_hj,CP_value, &
& CS_lep_hjZ_ratio, &
& CS_lep_bbhj_ratio,CS_lep_tautauhj_ratio, &
& CS_lep_hjhi_ratio_nHbynH, &
& CS_tev_hj_ratio ,CS_tev_hjb_ratio, &
& CS_tev_hjW_ratio,CS_tev_hjZ_ratio, &
& CS_tev_vbf_ratio,CS_tev_tthj_ratio, &
& CS_lhc7_hj_ratio ,CS_lhc7_hjb_ratio, &
& CS_lhc7_hjW_ratio,CS_lhc7_hjZ_ratio, &
& CS_lhc7_vbf_ratio,CS_lhc7_tthj_ratio, &
& CS_lhc8_hj_ratio ,CS_lhc8_hjb_ratio, &
& CS_lhc8_hjW_ratio,CS_lhc8_hjZ_ratio, &
& CS_lhc8_vbf_ratio,CS_lhc8_tthj_ratio, &
& BR_hjss,BR_hjcc, &
& BR_hjbb, &
& BR_hjmumu, &
& BR_hjtautau, &
& BR_hjWW,BR_hjZZ,BR_hjZga,BR_hjgaga, &
& BR_hjgg, BR_hjinvisible, &
& BR_hjhihi_nHbynH )
! This subroutine can be called by the user after subroutine initialize_HiggsBounds
! has been called.
! (see manual for full description)
!************************************************************
use usefulbits, only : theo,np,Hneut,whichinput,inputsub,just_after_run
#if defined(NAGf90Fortran)
use F90_UNIX_IO, only : flush
#endif
implicit none
!----------------------------------------input
double precision,intent(in) :: Mh( np(Hneut) ),GammaTotal_hj( np(Hneut) )
integer,intent(in) :: CP_value( np(Hneut) )
double precision,intent(in) :: CS_lep_hjZ_ratio( np(Hneut) ), &
& CS_lep_bbhj_ratio( np(Hneut) ),CS_lep_tautauhj_ratio( np(Hneut) ), &
& CS_lep_hjhi_ratio_nHbynH(np(Hneut),np(Hneut)), &
& CS_tev_hj_ratio( np(Hneut) ) ,CS_tev_hjb_ratio( np(Hneut) ), &
& CS_tev_hjW_ratio( np(Hneut) ) ,CS_tev_hjZ_ratio( np(Hneut) ), &
& CS_tev_vbf_ratio( np(Hneut) ) ,CS_tev_tthj_ratio( np(Hneut)), &
& CS_lhc7_hj_ratio( np(Hneut) ),CS_lhc7_hjb_ratio( np(Hneut) ), &
& CS_lhc7_hjW_ratio( np(Hneut) ),CS_lhc7_hjZ_ratio( np(Hneut) ), &
& CS_lhc7_vbf_ratio( np(Hneut) ),CS_lhc7_tthj_ratio( np(Hneut)), &
& CS_lhc8_hj_ratio( np(Hneut) ),CS_lhc8_hjb_ratio( np(Hneut) ), &
& CS_lhc8_hjW_ratio( np(Hneut) ),CS_lhc8_hjZ_ratio( np(Hneut) ), &
& CS_lhc8_vbf_ratio( np(Hneut) ),CS_lhc8_tthj_ratio( np(Hneut)), &
& BR_hjss( np(Hneut) ),BR_hjcc( np(Hneut) ), &
& BR_hjbb( np(Hneut) ), &
& BR_hjmumu( np(Hneut) ),BR_hjtautau( np(Hneut) ), &
& BR_hjWW( np(Hneut) ),BR_hjZZ( np(Hneut) ), &
& BR_hjZga( np(Hneut) ),BR_hjgaga( np(Hneut) ), &
& BR_hjgg( np(Hneut) ), BR_hjinvisible( np(Hneut) ), &
& BR_hjhihi_nHbynH(np(Hneut),np(Hneut))
!-------------------------------------internal
integer :: n
integer :: subtype
!---------------------------------------------
whichinput='hadr'
subtype=1
n=1
inputsub(subtype)%stat=inputsub(subtype)%stat+1
if(.not.allocated(theo))then
stop 'subroutine HiggsBounds_initialize must be called first'
endif
if(np(Hneut).eq.0)then
write(*,*)'subroutine HiggsBounds_neutral_input_hadr should'
write(*,*)'only be called if np(Hneut)>0'
stop 'error in subroutine HiggsBounds_neutral_input_hadr'
endif
! write(*,*) "DEBUG HB: before hadronic input. Mass is ",theo(n)%particle(Hneut)%M
theo(n)%particle(Hneut)%M = Mh
theo(n)%particle(Hneut)%Mc = Mh
theo(n)%particle(Hneut)%GammaTot= GammaTotal_hj
theo(n)%CP_value = CP_value
theo(n)%lep%XS_hjZ_ratio = CS_lep_hjZ_ratio
theo(n)%lep%XS_bbhj_ratio = CS_lep_bbhj_ratio
theo(n)%lep%XS_tautauhj_ratio = CS_lep_tautauhj_ratio
theo(n)%lep%XS_hjhi_ratio = CS_lep_hjhi_ratio_nHbynH
theo(n)%tev%XS_hj_ratio = CS_tev_hj_ratio
theo(n)%tev%XS_hjb_ratio = CS_tev_hjb_ratio
theo(n)%tev%XS_hjW_ratio = CS_tev_hjW_ratio
theo(n)%tev%XS_hjZ_ratio = CS_tev_hjZ_ratio
theo(n)%tev%XS_vbf_ratio = CS_tev_vbf_ratio
theo(n)%tev%XS_tthj_ratio = CS_tev_tthj_ratio
theo(n)%lhc7%XS_hj_ratio = CS_lhc7_hj_ratio
theo(n)%lhc7%XS_hjb_ratio = CS_lhc7_hjb_ratio
theo(n)%lhc7%XS_hjW_ratio = CS_lhc7_hjW_ratio
theo(n)%lhc7%XS_hjZ_ratio = CS_lhc7_hjZ_ratio
theo(n)%lhc7%XS_vbf_ratio = CS_lhc7_vbf_ratio
theo(n)%lhc7%XS_tthj_ratio = CS_lhc7_tthj_ratio
theo(n)%lhc8%XS_hj_ratio = CS_lhc8_hj_ratio
theo(n)%lhc8%XS_hjb_ratio = CS_lhc8_hjb_ratio
theo(n)%lhc8%XS_hjW_ratio = CS_lhc8_hjW_ratio
theo(n)%lhc8%XS_hjZ_ratio = CS_lhc8_hjZ_ratio
theo(n)%lhc8%XS_vbf_ratio = CS_lhc8_vbf_ratio
theo(n)%lhc8%XS_tthj_ratio = CS_lhc8_tthj_ratio
theo(n)%BR_hjss = BR_hjss
theo(n)%BR_hjcc = BR_hjcc
theo(n)%BR_hjbb = BR_hjbb
theo(n)%BR_hjmumu = BR_hjmumu
theo(n)%BR_hjtautau = BR_hjtautau
theo(n)%BR_hjWW = BR_hjWW
theo(n)%BR_hjZZ = BR_hjZZ
theo(n)%BR_hjZga = BR_hjZga
theo(n)%BR_hjgaga = BR_hjgaga
theo(n)%BR_hjgg = BR_hjgg
theo(n)%BR_hjinvisible = BR_hjinvisible
theo(n)%BR_hjhihi = BR_hjhihi_nHbynH
just_after_run=.False.
! write(*,*) "DEBUG HB: filled hadronic input. Mass is ",theo(n)%particle(Hneut)%M
end subroutine HiggsBounds_neutral_input_hadr
!************************************************************
subroutine HiggsBounds_charged_input(Mhplus,GammaTotal_Hpj, &
& CS_lep_HpjHmj_ratio, &
& BR_tWpb,BR_tHpjb, &
& BR_Hpjcs,BR_Hpjcb,BR_Hpjtaunu)
! This subroutine can be called by the user after subroutine initialize_HiggsBounds
! has been called.
! Arguments (input): theoretical predictions (see manual for definitions)
!************************************************************
use usefulbits, only : theo,np,Hplus,inputsub,just_after_run
#if defined(NAGf90Fortran)
use F90_UNIX_IO, only : flush
#endif
implicit none
!----------------------------------------input
double precision,intent(in) :: Mhplus( np(Hplus) ),GammaTotal_Hpj( np(Hplus) ), &
& CS_lep_HpjHmj_ratio( np(Hplus) ), &
& BR_tWpb,BR_tHpjb( np(Hplus) ), &
& BR_Hpjcs( np(Hplus) ),BR_Hpjcb( np(Hplus) ),BR_Hpjtaunu( np(Hplus) )
!--------------------------------------internal
integer :: n
integer :: subtype
!----------------------------------------------
n=1
subtype=2
inputsub(subtype)%stat=inputsub(subtype)%stat+1
if(.not.allocated(theo))then
stop 'subroutine HiggsBounds_initialize must be called first'
endif
if(np(Hplus).eq.0)then
write(*,*)'subroutine HiggsBounds_charged_input should'
write(*,*)'only be called if np(Hplus)>0'
stop 'error in subroutine HiggsBounds_charged_input'
endif
theo(n)%particle(Hplus)%M = Mhplus
theo(n)%particle(Hplus)%Mc = Mhplus
theo(n)%particle(Hplus)%GammaTot= GammaTotal_Hpj
theo(n)%lep%XS_HpjHmj_ratio = CS_lep_HpjHmj_ratio
theo(n)%BR_tWpb = BR_tWpb
theo(n)%BR_tHpjb = BR_tHpjb
theo(n)%BR_Hpjcs = BR_Hpjcs
theo(n)%BR_Hpjcb = BR_Hpjcb
theo(n)%BR_Hpjtaunu = BR_Hpjtaunu
just_after_run=.False.
end subroutine HiggsBounds_charged_input
!************************************************************
subroutine HiggsBounds_set_mass_uncertainties(dMhneut, dMhch)
! Assigns the mass uncertainties in the subroutine version.
!
use usefulbits, only : theo,np,Hneut,Hplus
implicit none
double precision, intent(in) :: dMhneut(np(Hneut))
double precision, intent(in) :: dMhch(np(Hplus))
theo(1)%particle(Hneut)%dMh = dMhneut
theo(1)%particle(Hplus)%dMh = dMhch
end subroutine HiggsBounds_set_mass_uncertainties
!************************************************************
subroutine get_mass_variation_param(n)
use usefulbits, only : theo,np,Hneut,Hplus,diffMhneut,diffMhch,ndmh,dmhsteps,small_mh
implicit none
integer, intent(in) :: n
double precision :: dMhneut(np(Hneut))
double precision :: dMhch(np(Hplus))
integer :: km(np(Hneut)+np(Hplus))
integer :: dm(dmhsteps**(np(Hneut)+np(Hplus)),np(Hneut)+np(Hplus))
integer i,j,k,kp
if(np(Hneut).gt.0) dMhneut = theo(n)%particle(Hneut)%dMh
if(np(Hplus).gt.0) dMhch = theo(n)%particle(Hplus)%dMh
if (modulo(dmhsteps,2).NE.1) then
stop 'Wrong number of steps in set_mass_uncertainty: must be odd (>=3)'
endif
ndmh = 0
do i=1,np(Hneut)
IF (dMhneut(i).GT.small_mh) THEN
ndmh = ndmh + 1
ENDIF
km(i)=-(dmhsteps-1)/2
enddo
do i=1,np(Hplus)
IF (dMhch(i).GT.small_mh) ndmh = ndmh + 1
km(i+np(Hneut))=-(dmhsteps-1)/2
enddo
IF (ndmh.EQ.0) THEN
RETURN
ENDIF
! print *, "Number of mass uncertainties: ", ndmh
if(allocated(diffMhneut)) deallocate(diffMhneut)
if(allocated(diffMhch)) deallocate(diffMhch)
allocate(diffMhneut(dmhsteps**(np(Hneut)+np(Hplus)),np(Hneut)))
allocate(diffMhch(dmhsteps**(np(Hneut)+np(Hplus)),np(Hplus)))
k = 1
do i=1,dmhsteps**ndmh
do j=1,ndmh
dm(i,j) = km(j)
enddo
km(k) = km(k)+1
do j=2,ndmh
IF (modulo(i,dmhsteps**(j-1)).EQ.0) THEN
km(j) = km(j)+1
km(j-1) = -1
ENDIF
ENDDO
enddo
do i=1,dmhsteps**ndmh
k=1
do j=1,np(Hneut)
IF (dMhneut(j).GT.small_mh) THEN
diffMhneut(i,j)=theo(n)%particle(Hneut)%M(j)+dm(i,k)*dMhneut(k)/((dmhsteps-1)/2)
k = k +1
ELSE
diffMhneut(i,j)=theo(n)%particle(Hneut)%M(j)
ENDIF
enddo
kp = k
do j=1,np(Hplus)
IF (dMhch(j).GT.small_mh) THEN
diffMhch(i,j)=theo(n)%particle(Hplus)%M(j)+dm(i,k)*dMhch(k-(kp-1))/((dmhsteps-1)/2)
k = k +1
ELSE
diffMhch(i,j)=theo(n)%particle(Hplus)%M(j)
ENDIF
enddo
! print *, i, (diffMhneut(i,j),j=1,np(Hneut)),(diffMhch(i,j),j=1,np(Hplus))
enddo
end subroutine get_mass_variation_param
subroutine SUSYBounds_neutralinoonly_input(MN,GammaTotal_N, &
& CS_NjNi, &
& BR_NjqqNi,BR_NjZNi &
& )
! This subroutine can be called by the user after subroutine initialize_HiggsBounds
! has been called.
! Arguments (input): theoretical predictions (see manual for definitions)
!************************************************************
use usefulbits, only : theo,np,Chineut,inputsub,just_after_run
#if defined(NAGf90Fortran)
use F90_UNIX_IO, only : flush
#endif
implicit none
!----------------------------------------input
double precision,intent(in) :: MN( np(Chineut) ),GammaTotal_N( np(Chineut) ) , &
& CS_NjNi( np(Chineut),np(Chineut) ), &
& BR_NjqqNi( np(Chineut),np(Chineut) ),BR_NjZNi( np(Chineut),np(Chineut) )
!--------------------------------------internal
integer :: n
integer :: subtype
!----------------------------------------------
n=1
subtype=3
inputsub(subtype)%stat=inputsub(subtype)%stat+1
if(.not.allocated(theo))then
stop 'subroutine HiggsBounds_initialize must be called first'
endif
if(np(Chineut).eq.0)then
write(*,*)'subroutine SUSYBounds_neutralinoonly_input should'
write(*,*)'only be called if np(Chineut)>0'
stop'error in SUSYBounds_neutralinoonly_input'
endif
theo(n)%particle(Chineut)%M = MN
theo(n)%particle(Chineut)%GammaTot= GammaTotal_N
theo(n)%lep%XS_NjNi = CS_NjNi
theo(n)%BR_NjqqNi = BR_NjqqNi
theo(n)%BR_NjZNi = BR_NjZNi
just_after_run=.False.
end subroutine SUSYBounds_neutralinoonly_input
!************************************************************
subroutine SUSYBounds_neutralinochargino_input(MC,GammaTotal_C, &
& CS_CpjCmj, &
& BR_CjqqNi, &
& BR_CjlnuNi, &
& BR_CjWNi &
& )
! This subroutine can be called by the user after subroutine initialize_HiggsBounds
! has been called.
! Arguments (input): theoretical predictions (see manual for definitions)
!************************************************************
use usefulbits, only : theo,np,Chineut,Chiplus,inputsub,just_after_run
#if defined(NAGf90Fortran)
use F90_UNIX_IO, only : flush
#endif
implicit none
!----------------------------------------input
double precision,intent(in) :: MC( np(Chiplus) ),GammaTotal_C( np(Chiplus) ), &
& CS_CpjCmj( np(Chiplus) ), &
& BR_CjqqNi( np(Chiplus),np(Chineut) ), &
& BR_CjlnuNi( np(Chiplus),np(Chineut) ), &
& BR_CjWNi( np(Chiplus),np(Chineut) )
!--------------------------------------internal
integer :: n
integer :: subtype
!----------------------------------------------
n=1
subtype=4
inputsub(subtype)%stat=inputsub(subtype)%stat+1
if(.not.allocated(theo))then
stop 'subroutine HiggsBounds_initialize must be called first'
endif
if((np(Chineut).eq.0).or.(np(Chiplus).eq.0))then
write(*,*)'subroutine SUSYBounds_neutralinochargino_input should'
write(*,*)'only be called if np(Chineut)>0 and np(Chiplus)>0'
stop 'error in subroutine SUSYBounds_neutralinochargino_input'
endif
theo(n)%particle(Chineut)%M = MC
theo(n)%particle(Chineut)%GammaTot= GammaTotal_C
theo(n)%lep%XS_CpjCmj = CS_CpjCmj
theo(n)%BR_CjqqNi = BR_CjqqNi
theo(n)%BR_CjlnuNi = BR_CjlnuNi
theo(n)%BR_CjWNi = BR_CjWNi
just_after_run=.False.
end subroutine SUSYBounds_neutralinochargino_input
!************************************************************
subroutine run_HiggsBounds(HBresult, chan, obsratio, ncombined)
! This subroutine can be called by the user after HiggsBounds_initialize has been called.
! The input routines, where required, should be called before calling run_HiggsBounds.
! It takes theoretical predictions for a particular parameter point
! in the model and calls subroutines which compare these predictions
! to the experimental limits
! Arguments (output):
! * HBresult = 1 if point is unexcluded, 0 if excluded, -1 if parameter point is invalid
! * chan = number of channel predicted to have the highest statistical sensitivity, as defined in Key.dat
! * obsratio = ratio of the theoretical rate to the observed limit for this channel
! * ncombined = number of Higgs combined in order to calculate this obsratio
! (see manual for more precise definitions))
! (TS 30/01/2012): Note, that if many data points are tested at the same time (as for
! inputmethod==datfiles), this subroutine only returns the results of
! the last datapoint. The full results are saved in fullHBres.
use usefulbits, only : np, Hneut, Hplus, run_HB_classic
implicit none
integer HBresult, chan, ncombined
double precision obsratio
integer hbres(0:np(Hneut)+np(Hplus)), hbchan(0:np(Hneut)+np(Hplus)), hbcomb(0:np(Hneut)+np(Hplus))
double precision hbobs(0:np(Hneut)+np(Hplus))
! Check if we are using the old 'classic' method
if (run_HB_classic.EQV..True.) then
call run_HiggsBounds_classic(HBresult,chan,obsratio,ncombined)
return
endif
! Call the new ('full') method
call run_HiggsBounds_full(hbres, hbchan, hbobs, hbcomb)
! Combined results are contained in the zero elements of result arrays
HBresult = hbres(0)
chan = hbchan(0)
obsratio = hbobs(0)
ncombined = hbcomb(0)
end subroutine run_HiggsBounds
!************************************************************
subroutine run_HiggsBounds_single(h, HBresult, chan, obsratio, ncombined)
! This subroutine can be used to get the exclusion results
! for a single Higgs boson (specified by the index h).
!
! To obtain individual results from more than one Higgs boson, it
! is more efficient to use run_HiggsBounds_full rather than this method.
use usefulbits, only : np, Hneut, Hplus
implicit none
integer, intent(in) :: h
integer, intent(out) :: HBresult, chan, ncombined
double precision, intent(out) :: obsratio
integer hbres(0:np(Hneut)+np(Hplus)), hbchan(0:np(Hneut)+np(Hplus)), hbcomb(0:np(Hneut)+np(Hplus))
double precision hbobs(0:np(Hneut)+np(Hplus))
IF (h.LT.0) stop "Illegal number of Higgs boson: h < 0"
if (h.GT.np(Hneut)+np(Hplus)) stop "Illegal number of Higgs boson"
call run_HiggsBounds_full(hbres, hbchan, hbobs, hbcomb)
HBresult = hbres(h)
chan = hbchan(h)
obsratio = hbobs(h)
ncombined = hbcomb(h)
end subroutine run_HiggsBounds_single
!************************************************************
subroutine run_HiggsBounds_full( HBresult,chan, &
& obsratio, ncombined )
! This subroutine can be called by the user after HiggsBounds_initialize has been called.
! The input routines, where required, should be called before calling run_HiggsBounds.
! It takes theoretical predictions for a particular parameter point
! in the model and calls subroutines which compare these predictions
! to the experimental limits.
!
! The results are given as (n+1)-component arrays (starting from 0),
! where n is the total number of Higgs bosons in the model (neutral+charged).
! The zeroth component gives the combined results (equivalent to run_HiggsBounds).
!
! Arguments (output):
! * HBresult = 1 if point is unexcluded, 0 if excluded, -1 if parameter point is invalid
! * chan = number of channel predicted to have the highest statistical sensitivity, as defined in Key.dat
! * obsratio = ratio of the theoretical rate to the observed limit for this channel
! * ncombined = number of Higgs combined in order to calculate this obsratio
! (see manual for more precise definitions))
use usefulbits, only : theo,res,inputsub,just_after_run,ndmh,debug, &
& np,Hneut,Hplus,dmhsteps,ndat,fullHBres,small_mh
use channels, only : check_channels
!use input, only : test_input
use theo_manip, only : complete_theo, recalculate_theo_for_datapoint
#if defined(NAGf90Fortran)
use F90_UNIX_IO, only : flush
#endif
implicit none
!----------------------------------------output
integer,intent(out):: HBresult(0:np(Hneut)+np(Hplus))
integer,intent(out):: chan(0:np(Hneut)+np(Hplus))
integer,intent(out):: ncombined(0:np(Hneut)+np(Hplus))
double precision,intent(out) :: obsratio(0:np(Hneut)+np(Hplus))
double precision :: Mhneut(np(Hneut))
double precision :: Mhch(np(Hplus))
!-------------------------------------internal
integer :: n,i,j,ind,part
!---------------------------------------------
! print *, "Running HiggsBounds in Normal Mode (most sensitive limit considered for each Higgs boson)"
if (lbound(HBresult,dim=1).NE.0) stop "run_HiggsBounds_full: Array HBresult must begin with element 0"
if (ubound(HBresult,dim=1).NE.(np(Hneut)+np(Hplus))) then
stop "run_HiggsBounds_full: Upper limit of array HBresult must be equal to number of Higgses"
endif
if (lbound(chan,dim=1).NE.0) stop "run_HiggsBounds_full: Array chan must begin with element 0"
if (ubound(chan,dim=1).NE.(np(Hneut)+np(Hplus))) then
stop "run_HiggsBounds_full: Upper limit of array chan must be equal to number of Higgses"
endif
if (lbound(obsratio,dim=1).NE.0) stop "run_HiggsBounds_full: Array obsratio must begin with element 0"
if (ubound(obsratio,dim=1).NE.(np(Hneut)+np(Hplus))) then
stop "run_HiggsBounds_full: Upper limit of array obsratio must be equal to number of Higgses"
endif
if (lbound(ncombined,dim=1).NE.0) stop "run_HiggsBounds_full: Array ncombined must begin with element 0"
if (ubound(ncombined,dim=1).NE.(np(Hneut)+np(Hplus))) then
stop "run_HiggsBounds_full: Upper limit of array ncombined must be equal to number of Higgses"
endif
if(.not.allocated(theo))then
stop 'subroutine HiggsBounds_initialize must be called first'
endif
do i=1,ubound(inputsub,dim=1)
if( inputsub(i)%req .ne. inputsub(i)%stat )then
write(*,*)'subroutine '//trim(adjustl(inputsub(i)%desc))
write(*,*)'should be called once and only once before each call to'
write(*,*)'subroutine run_HiggsBounds.'
stop'error in subroutine run_HiggsBounds'
endif
inputsub(i)%stat=0!now we have used this input, set back to zero
enddo
call complete_theo
do n=1,ndat
! if(debug) then
! write(*,*) "DEBUG BRs: ", theo(n)%BR_hjWW, theo(n)%BR_hjZZ, theo(n)%BR_hjgaga
! endif
theo(n)%particle(Hneut)%Mc = theo(n)%particle(Hneut)%M
theo(n)%particle(Hplus)%Mc = theo(n)%particle(Hplus)%M
call get_mass_variation_param(n)
do i=0,ubound(Hbresult,dim=1)
obsratio(i) = -999d0
HBresult(i) = 1
chan(i) = -999
ncombined(i) = -999
enddo
! Do we have mass uncertainties to take care off
IF (ndmh.GT.0) THEN
! print *, "Running HiggsBounds with Higgs mass uncertainties"
! write(*,*) theo(n)%particle(Hplus)%dM
if(np(Hneut).ne.0) Mhneut = theo(n)%particle(Hneut)%M
if(np(Hplus).ne.0) Mhch = theo(n)%particle(Hplus)%M
! Loop over all Higgses
do i=1,np(Hneut)+np(Hplus)
obsratio(i) = 1.D23
IF (i.LE.np(Hneut)) THEN
ind = i
part = Hneut
ELSE
ind = i-np(Hneut)
part = Hplus
ENDIF
! Check for mass steps for this particular Higgs boson
IF(theo(n)%particle(part)%dMh(ind).GT.small_mh) THEN
! theo(n)%particle(part)%M(ind)=theo(n)%particle(part)%M(ind) &
! & -(dmhsteps-1)/2*theo(n)%particle(part)%dMh(ind)
theo(n)%particle(part)%M(ind)=theo(n)%particle(part)%M(ind) &
& -theo(n)%particle(part)%dMh(ind)
do j=1,dmhsteps
! print *, theo(n)%particle(Hneut)%M, theo(n)%particle(Hplus)%M
call recalculate_theo_for_datapoint(n)
call check_channels(theo(n),res(n),i)
IF (res(n)%obsratio(1).LT.obsratio(i)) THEN
HBresult(i) = res(n)%allowed95(1)
chan(i) = res(n)%chan(1)
obsratio(i) = res(n)%obsratio(1)
ncombined(i) = res(n)%ncombined(1)
ENDIF
! print *, i,theo(n)%particle(part)%M(ind),HBresult(i),chan(i),obsratio(i),ncombined(i)
theo(n)%particle(part)%M(ind)= theo(n)%particle(part)%M(ind) &
& +theo(n)%particle(part)%dMh(ind)/(dmhsteps-1)*2
enddo
else
call recalculate_theo_for_datapoint(n)
call check_channels(theo(n),res(n),i)
HBresult(i) = res(n)%allowed95(1)
chan(i) = res(n)%chan(1)
obsratio(i) = res(n)%obsratio(1)
ncombined(i) = res(n)%ncombined(1)
endif
! Logical OR between exclusions (one Higgs excluded = combined exclusion)
HBresult(0) = HBresult(0) * HBresult(i)
! Save the data for the Higgs that has the highest ratio of theory/obs
IF (obsratio(i).GT.obsratio(0)) THEN
chan(0) = chan(i)
obsratio(0) = obsratio(i)
ncombined(0) = ncombined(i)
ENDIF
theo(n)%particle(Hneut)%M = Mhneut
theo(n)%particle(Hplus)%M = Mhch
enddo
! return
ELSE
! print *, "Running HiggsBounds without Higgs mass uncertainties"
call recalculate_theo_for_datapoint(n)
! write(*,*) "Higgses = " , np(Hneut)+np(Hplus)
do i=1,np(Hneut)+np(Hplus)
call check_channels(theo(n),res(n),i)
HBresult(i) = res(n)%allowed95(1)
chan(i) = res(n)%chan(1)
obsratio(i) = res(n)%obsratio(1)
ncombined(i) = res(n)%ncombined(1)
HBresult(0) = HBresult(0) * res(n)%allowed95(1)
IF (obsratio(i).GT.obsratio(0)) THEN
! write(*,*) "hello: ", n, i
chan(0) = res(n)%chan(1)
obsratio(0) = res(n)%obsratio(1)
ncombined(0) = res(n)%ncombined(1)
ENDIF
! IF (i.LE.np(Hneut)) THEN
! print *, i,theo(n)%particle(Hneut)%M(i),HBresult(i),chan(i),obsratio(i),ncombined(i),HBresult(0), obsratio(0)
! ELSE
! print *, i,theo(n)%particle(Hplus)%M(i-np(Hneut)),HBresult(i),chan(i),obsratio(i),ncombined(i),HBresult(0), obsratio(0)
! endif
enddo
ENDIF
fullHBres(n)%allowed95=HBresult(0)
fullHBres(n)%chan=chan(0)
fullHBres(n)%obsratio=obsratio(0)
fullHBres(n)%ncombined=ncombined(0)
enddo
just_after_run=.True.
end subroutine run_HiggsBounds_full
!************************************************************
subroutine run_HiggsBounds_classic( HBresult,chan,obsratio,ncombined)
! This subroutine can be called by the user after HiggsBounds_initialize has been called.
! The input routines, where required, should be called before calling run_HiggsBounds.
! It takes theoretical predictions for a particular parameter point
! in the model and calls subroutines which compare these predictions
! to the experimental limits
! Arguments (output):
! * HBresult = 1 if point is unexcluded, 0 if excluded, -1 if parameter point is invalid
! * chan = number of channel predicted to have the highest statistical sensitivity, as defined in Key.dat
! * obsratio = ratio of the theoretical rate to the observed limit for this channel
! * ncombined = number of Higgs combined in order to calculate this obsratio
! (see manual for more precise definitions))
use usefulbits, only : theo,res,debug,inputsub,just_after_run,ndmh,diffmhneut,diffmhch, &
np,Hneut,Hplus,full_dmth_variation,dmhsteps, ndat,fullHBres
use channels, only : check_channels
!use input, only : test_input
use theo_manip, only : complete_theo, recalculate_theo_for_datapoint
#if defined(NAGf90Fortran)
use F90_UNIX_IO, only : flush
#endif
implicit none
!----------------------------------------output
integer,intent(out):: HBresult,chan,ncombined
double precision,intent(out) :: obsratio
double precision :: Mhneut(np(Hneut))
double precision :: Mhch(np(Hplus))
!-------------------------------------internal
integer :: n,i
integer :: HBresult_tmp,chan_tmp,ncombined_tmp
double precision :: obsratio_tmp
!---------------------------------------------
! n=1
! print *, "Running HiggsBounds in Classic Mode (globally most sensitive limit only)"
if(.not.allocated(theo))then
stop 'subroutine HiggsBounds_initialize must be called first'
endif
do i=1,ubound(inputsub,dim=1)
if( inputsub(i)%req .ne. inputsub(i)%stat )then
write(*,*)'subroutine '//trim(adjustl(inputsub(i)%desc))
write(*,*)'should be called once and only once before each call to'
write(*,*)'subroutine run_HiggsBounds.'
stop'error in subroutine run_HiggsBounds'
endif
inputsub(i)%stat=0!now we have used this input, set back to zero
enddo
call complete_theo
do n=1,ndat
theo(n)%particle(Hneut)%Mc = theo(n)%particle(Hneut)%M
call get_mass_variation_param(n)
IF (ndmh.GT.0) THEN
if(np(Hneut).ne.0) Mhneut = theo(n)%particle(Hneut)%M
if(np(Hplus).ne.0) Mhch = theo(n)%particle(Hplus)%M
obsratio_tmp = 10.0E6 ! Set to very large initial value
do i=1,dmhsteps**ndmh
theo(n)%particle(Hneut)%M = diffMhneut(i,:)
theo(n)%particle(Hplus)%M = diffMhch(i,:)
if(debug)write(*,*)'manipulating input...' ; call flush(6)
call recalculate_theo_for_datapoint(n)
if(debug)write(*,*)'compare each data point to the experimental bounds...' ; call flush(6)
call check_channels(theo(n),res(n),0)
HBresult = res(n)%allowed95(1)
chan = res(n)%chan(1)
obsratio = res(n)%obsratio(1)
ncombined = res(n)%ncombined(1)
! print *, HBresult, chan, obsratio, ncombined
IF (.NOT.full_dmth_variation) THEN
IF (HBresult.EQ.1) THEN
! theo(n)%particle(Hneut)%M = Mhneut
! theo(n)%particle(Hplus)%M = Mhch
just_after_run=.True.
exit
ENDIF
ELSE
IF (obsratio.lt.obsratio_tmp) THEN
HBresult_tmp = HBresult
chan_tmp = chan
obsratio_tmp = obsratio
ncombined_tmp = ncombined
ENDIF
ENDIF
enddo
IF (full_dmth_variation) THEN
HBresult = HBresult_tmp
chan = chan_tmp
obsratio = obsratio_tmp
ncombined = ncombined
! theo(n)%particle(Hneut)%M = Mhneut
! theo(n)%particle(Hplus)%M = Mhch
just_after_run=.True.
! return
ENDIF
theo(n)%particle(Hneut)%M = Mhneut
theo(n)%particle(Hplus)%M = Mhch
call recalculate_theo_for_datapoint(n)
call check_channels(theo(n),res(n),0)
ELSE
if(debug)write(*,*)'manipulating input...' ; call flush(6)
call recalculate_theo_for_datapoint(n)
if(debug)write(*,*)'compare each data point to the experimental bounds...' ; call flush(6)
call check_channels(theo(n),res(n),0)
HBresult = res(n)%allowed95(1)
chan = res(n)%chan(1)
obsratio = res(n)%obsratio(1)
ncombined = res(n)%ncombined(1)
just_after_run=.True.
ENDIF
fullHBres(n)%allowed95=HBresult
fullHBres(n)%chan=chan
fullHBres(n)%obsratio=obsratio
fullHBres(n)%ncombined=ncombined
enddo
just_after_run=.True.
end subroutine run_HiggsBounds_classic
!************************************************************
-subroutine HiggsBounds_get_chisq(analysisID, M_av, Hindex, nc, llh)
+subroutine HiggsBounds_get_chisq(analysisID, Hindex, M_av, nc, llh)
!************************************************************
use usefulbits, only : theo,np,Hneut,Hplus
use theo_manip, only : complete_theo
use likelihoods, only : get_likelihood, calcpredratio_llh
implicit none
integer, intent(in) :: analysisID
integer, intent(out) :: Hindex, nc
double precision, intent(out) :: llh, M_av
- integer :: c,i
+ integer :: c,i,cbin
double precision, allocatable :: expllh(:)
! double precision :: fact
double precision, allocatable :: mass(:) ! predratio(:)
integer, allocatable :: nclist(:)
! call complete_theo
! allocate(predratio(np(Hneut)))
! predratio = 0.0D0
! write(*,*) "Calling HiggsBounds_get_chisq..."
allocate(expllh(np(Hneut)),mass(np(Hneut)),nclist(np(Hneut)))
expllh = 0.0D0
select case(analysisID)
case(3316)
c=1
case default
stop 'Unknown analysisID in subroutine HiggsBounds_get_chisq!'
- end select
-
-
+ end select
! Determine most sensitive combination
do i=1,np(Hneut)
-! Change this here
-! call calcpredratio_llh(c,i,theo(1),mass(i),fact,nclist(i),predratio(i))
-! to determining the most sensitive combination judged by maximal expected -2lnL:
- call get_likelihood(analysisID, i, theo(1), expllh(i), mass(i), nclist(i), 'pred')
+ call get_likelihood(analysisID, i, theo(1), expllh(i), mass(i), nclist(i), cbin, 'pred')
enddo
-! write(*,*) 'HB_get_chisq, expllh = ',expllh
-! write(*,*) 'HB_get_chisq, mass = ',mass
-! write(*,*) 'HB_get_chisq, nclist = ',nclist
-
Hindex = maxloc(expllh,dim=1)
- call get_likelihood(analysisID, Hindex, theo(1), llh, M_av, nc, 'obs')
-! nc = nclist(Hindex)
-! M_av = mass(Hindex)
+ call get_likelihood(analysisID, Hindex, theo(1), llh, M_av, nc, cbin, 'obs')
deallocate(mass,nclist,expllh) !predratio
end subroutine HiggsBounds_get_chisq
!************************************************************
-subroutine HiggsBounds_get_maximal_chisq(analysisID, M_av, Hindex, nc, llh)
+subroutine HiggsBounds_get_combined_chisq(analysisID, llh)
+!************************************************************
+ use usefulbits, only : theo,np,Hneut,Hplus, vsmall
+
+ integer, intent(in) :: analysisID
+ double precision, intent(out) :: llh
+
+ double precision :: M_av, llh_tmp
+ integer :: i, j, nc, cbin, Hindex, cbin_end, cbin_in
+
+ cbin_end = 0
+ do i= 1,np(Hneut)
+ cbin_end = cbin_end + 2**(i-1)
+ enddo
+
+ llh = -1.0D0
+ cbin_in = 0
+ llh_tmp = 0.0D0
+
+ do while(cbin_in.lt.cbin_end)
+ call HiggsBounds_get_maximal_chisq_for_comb(analysisID, 'obs', cbin_in, Hindex, cbin, nc, M_av, llh)
+ if(llh.ge.0.0D0) then
+ llh_tmp = llh_tmp + llh
+ else
+ exit
+ endif
+ cbin_in = cbin_in + cbin
+! write(*,*) "Found combination ",cbin, " with llh = ", llh, ". cbin_in = ", cbin_in
+ enddo
+
+ if(llh_tmp.gt.0.0D0) then
+ llh = llh_tmp
+ endif
+! allocate(obsllh(np(Hneut)),mass(np(Hneut)),nclist(np(Hneut)),added(np(Hneut)),cbin(np(Hneut)))
+!
+!obsllh = 0.0D0
+! mass = 0.0D0
+! llh = 0.0D0
+! added = .False.
+!
+! do i=1,np(Hneut)
+! call HiggsBounds_get_chisq_for_Higgs(analysisID, i, mass(i), nclist(i), cbin(i), obsllh(i))
+! enddo
+!
+! do i=1,np(Hneut)
+! if(.not.added(i)) then
+! do j=1,np(Hneut)
+! if(.not.added(j)) then
+! if(abs(obsllh(i)-obsllh(j)).le.vsmall) then
+! added(j) = .True.
+! endif
+! endif
+! enddo
+! if(obsllh(i).ge.0.0D0) then
+! llh = llh + obsllh(i)
+! endif
+! endif
+! enddo
+!
+! deallocate(obsllh,mass,nclist,added,cbin)
+
+end subroutine HiggsBounds_get_combined_chisq
+!************************************************************
+subroutine HiggsBounds_get_chisq_for_Higgs(analysisID, cbin_in, Hindex, nc, cbin, M_av, llh)
!************************************************************
use usefulbits, only : theo,np,Hneut,Hplus
use theo_manip, only : complete_theo
use likelihoods, only : get_likelihood, calcpredratio_llh
implicit none
+ integer, intent(in) :: analysisID,Hindex
+ integer, intent(out) :: nc, cbin
+ double precision, intent(out) :: llh, M_av
+ integer, intent(in) :: cbin_in
+ integer :: c,i
+
+! write(*,*) "Calling HiggsBounds_get_chisq_for_Higgs..."
+ select case(analysisID)
+ case(3316)
+ c=1
+ case default
+ stop 'Unknown analysisID in subroutine HiggsBounds_get_chisq!'
+ end select
+
+ call get_likelihood(analysisID, Hindex, theo(1), llh, M_av, nc, cbin,'obs',cbin_in)
+
+! write(*,*) "h(",Hindex,") llh: ", llh, M_av, nc, cbin
+
+end subroutine HiggsBounds_get_chisq_for_Higgs
+!************************************************************
+subroutine HiggsBounds_get_maximal_chisq(analysisID, Hindex, nc, M_av, llh)
+! Wrapper subroutine for HiggsBounds_get_maximal_chisq_for_comb considering
+! all neutral Higgs bosons
+!************************************************************
+
+ use usefulbits, only : theo,np,Hneut,Hplus
+
integer, intent(in) :: analysisID
integer, intent(out) :: Hindex, nc
double precision, intent(out) :: llh, M_av
+
+ integer :: cbin
+
+ call HiggsBounds_get_maximal_chisq_for_comb(analysisID, 'obs', 0, Hindex, cbin, nc, M_av, llh)
+
+end subroutine HiggsBounds_get_maximal_chisq
+!************************************************************
+subroutine HiggsBounds_get_maximal_chisq_for_comb(analysisID, obspred, cbin_in, Hindex, cbin, nc, M_av, llh)
+!************************************************************
+ use usefulbits, only : theo,np,Hneut,Hplus
+ use theo_manip, only : complete_theo
+ use likelihoods, only : get_likelihood, calcpredratio_llh
+ implicit none
+
+ integer, intent(in) :: analysisID, cbin_in
+ integer, intent(out) :: Hindex, nc, cbin
+ double precision, intent(out) :: llh, M_av
+ character(LEN=*), intent(in) :: obspred
integer :: c,i
double precision, allocatable :: obsllh(:)
! double precision :: fact
double precision, allocatable :: mass(:) ! predratio(:)
- integer, allocatable :: nclist(:)
+ integer, allocatable :: nclist(:), cbinlist(:)
! call complete_theo
! allocate(predratio(np(Hneut)))
! predratio = 0.0D0
-! write(*,*) "Calling HiggsBounds_get_chisq..."
+! write(*,*) "Calling HiggsBounds_get_maximal_chisq_for_comb with dataset: ", obspred
- allocate(obsllh(np(Hneut)),mass(np(Hneut)),nclist(np(Hneut)))
+ allocate(obsllh(np(Hneut)),mass(np(Hneut)),nclist(np(Hneut)),cbinlist(np(Hneut)))
obsllh = 0.0D0
select case(analysisID)
case(3316)
c=1
case default
- stop 'Unknown analysisID in subroutine HiggsBounds_get_chisq!'
+ stop 'Unknown analysisID in subroutine HiggsBounds_get_maximal_chisq!'
end select
-
-
! Determine most sensitive combination
do i=1,np(Hneut)
-! Change this here
-! call calcpredratio_llh(c,i,theo(1),mass(i),fact,nclist(i),predratio(i))
-! to determining the most sensitive combination judged by maximal expected -2lnL:
- call get_likelihood(analysisID, i, theo(1), obsllh(i), mass(i), nclist(i), 'obs')
+ call get_likelihood(analysisID, i, theo(1), obsllh(i), mass(i), nclist(i),cbinlist(i),obspred, cbin_in)
enddo
-! write(*,*) 'HB_get_chisq, expllh = ',expllh
-! write(*,*) 'HB_get_chisq, mass = ',mass
-! write(*,*) 'HB_get_chisq, nclist = ',nclist
-
Hindex = maxloc(obsllh,dim=1)
llh = obsllh(Hindex)
M_av = mass(Hindex)
nc = nclist(Hindex)
+ cbin = cbinlist(Hindex)
-! call get_likelihood(analysisID, Hindex, theo(1), llh, M_av, nc, 'obs')
-! nc = nclist(Hindex)
-! M_av = mass(Hindex)
-
- deallocate(mass,nclist,obsllh) !predratio
+ deallocate(mass,nclist,obsllh,cbinlist) !predratio
-end subroutine HiggsBounds_get_maximal_chisq
+end subroutine HiggsBounds_get_maximal_chisq_for_comb
!************************************************************
subroutine HiggsBounds_SLHA_output
!**** ********************************************************
use usefulbits, only : whichinput,just_after_run
use output, only : do_output
if(.not.just_after_run)then
stop'subroutine run_HiggsBounds should be called before subroutine HiggsBounds_SLHA_output'
endif
select case(whichinput)
case('SLHA')
call do_output
case default
stop'The subroutine HiggsBounds_SLHA_output should only be used when whichinput=SLHA'
end select
end subroutine HiggsBounds_SLHA_output
#ifdef enableCHISQ
!************************************************************
subroutine initialize_HiggsBounds_chisqtables
! use S95tables, only : S95_t2
use S95tables_type3
use usefulbits, only : allocate_if_stats_required,theo
implicit none
if(allocated(theo))then
stop 'subroutine initialize_HiggsBounds_chisqtables should be called before subroutine HiggsBounds_initialize'
elseif(allocated(clsb_t3))then
stop 'subroutine initialize_HiggsBounds_chisqtables has already been called once'
endif
allocate(clsb_t3(ntable3))
call initializetables_type3_blank(clsb_t3)
call initializetables3(clsb_t3)
call readclsbfiles_binary
if(allocated(allocate_if_stats_required))then
stop'error in subroutine initialize_HiggsBounds_chisqtables'
else
allocate(allocate_if_stats_required(1))
endif
end subroutine initialize_HiggsBounds_chisqtables
!************************************************************
subroutine finish_HiggsBounds_chisqtables
!************************************************************
use S95tables_type3
use usefulbits, only : allocate_if_stats_required
implicit none
integer :: x
if(.not.allocated(clsb_t3))then
stop 'initialize_HiggsBounds_chisqtables should be called first'
endif
do x=lbound(clsb_t3,dim=1),ubound(clsb_t3,dim=1)
deallocate(clsb_t3(x)%dat)
enddo
deallocate(filename)
deallocate(clsb_t3)
deallocate(allocate_if_stats_required)
end subroutine finish_HiggsBounds_chisqtables
!************************************************************
subroutine HB_calc_stats(theory_uncertainty_1s,chisq_withouttheory,chisq_withtheory,chan2)
!************************************************************
! this is in the middle of development! DO NOT USE!
use usefulbits, only : res,theo,pr,just_after_run,vsmall
use interpolate
use S95tables_type1
use S95tables_type3
use S95tables
use extra_bits_for_chisquared
implicit none
integer,intent(out)::chan2
integer :: x,c,z,y
integer :: id
double precision, intent(in) :: theory_uncertainty_1s
double precision :: chisq_withouttheory,chisq_withtheory
double precision :: low_chisq,sigma
x=1
low_chisq=1.0D-2
if(.not.allocated(theo))then
stop 'subroutine HiggsBounds_initialize must be called first'
elseif(.not.allocated(clsb_t3))then
stop 'subroutine initialize_HiggsBounds_chisqtables must be called first'
elseif(.not.just_after_run)then
stop 'subroutine run_HiggsBounds must be called first'
endif
sigma=theory_uncertainty_1s
if(sigma.lt.vsmall)then
write(*,*)'Warning: will not calculate chi^2 with theory uncertainty'
endif
chisq_withtheory = -2.0D0
chisq_withouttheory = -2.0D0
z=2;
c= res(x)%chan(z)
chan2=c
if(res(x)%allowed95(z).eq.-1)then! labels an unphysical parameter point
chisq_withtheory =-1.0D0
chisq_withouttheory =-1.0D0
elseif( c.gt.0 )then ! labels a physical parameter point and a real channel
id=S95_t1_or_S95_t2_idfromelementnumber(pr(c)%ttype,pr(c)%tlist)
y=clsb_t3elementnumber_from_S95table(pr(c)%ttype,id)
if(y.gt.0)then
!------------------------------
call get_chisq(sigma,res(x)%axis_i(z),res(x)%axis_j(z),res(x)%sfactor(z), &
& y,chisq_withouttheory,chisq_withtheory)
!-------------------------------
else
write(*,*)'hello y=',y
stop'problem here with y'
endif
else
chisq_withtheory =0.0D0
chisq_withouttheory =0.0D0
endif
end subroutine HB_calc_stats
#endif
!************************************************************
subroutine finish_HiggsBounds
! This subroutine needs to be called right at the end, to close files
! and deallocate arrays
!************************************************************
use usefulbits, only : deallocate_usefulbits,debug,theo,debug,inputsub, &
& file_id_debug1,file_id_debug2
use S95tables, only : deallocate_S95tables
use theory_BRfunctions, only : deallocate_BRSM
#if defined(NAGf90Fortran)
use F90_UNIX_IO, only : flush
#endif
if(debug)then
close(file_id_debug2)
close(file_id_debug1)
endif
if(.not.allocated(theo))then
stop 'HiggsBounds_initialize should be called first'
endif
if(debug)write(*,*)'finishing off...' ; call flush(6)
call deallocate_BRSM
call deallocate_S95tables
call deallocate_usefulbits
if(debug)write(*,*)'finished' ; call flush(6)
if(allocated(inputsub)) deallocate(inputsub)
end subroutine finish_HiggsBounds
!
!!************************************************************
!
!subroutine run_HiggsBounds_effC(nHdummy,Mh,GammaTotal, &
! & g2hjbb,g2hjtautau,g2hjWW,g2hjZZ, &
! & g2hjgaga,g2hjgg,g2hjhiZ_nHbynH, &
! & BR_hjhihi_nHbynH, &
! & HBresult,chan, &
! & obsratio, ncombined )
!! Obsolete subroutine
!!************************************************************
!
! implicit none
!
! !----------------------------------------input
! integer,intent(in) :: nHdummy
! double precision,intent(in) :: Mh(nHdummy),GammaTotal(nHdummy),&
! & g2hjbb(nHdummy),g2hjtautau(nHdummy), &
! & g2hjWW(nHdummy),g2hjZZ(nHdummy), &
! & g2hjgaga(nHdummy),g2hjgg(nHdummy), &
! & g2hjhiZ_nHbynH(nHdummy,nHdummy), &
! & BR_hjhihi_nHbynH(nHdummy,nHdummy)
! !----------------------------------------output
! integer :: HBresult,chan,ncombined
! double precision :: obsratio
! !----------------------------------------------
!
! call attempting_to_use_an_old_HB_version('effC')
!
!end subroutine run_HiggsBounds_effC
!!************************************************************
!subroutine run_HiggsBounds_part(nHdummy,Mh, &
! & CS_lep_hjZ_ratio, &
! & CS_lep_hjhi_ratio_nHbynH, &
! & CS_tev_gg_hj_ratio,CS_tev_bb_hj_ratio, &
! & CS_tev_bg_hjb_ratio, &
! & CS_tev_ud_hjWp_ratio,CS_tev_cs_hjWp_ratio, &
! & CS_tev_ud_hjWm_ratio,CS_tev_cs_hjWm_ratio, &
! & CS_tev_dd_hjZ_ratio,CS_tev_uu_hjZ_ratio, &
! & CS_tev_ss_hjZ_ratio,CS_tev_cc_hjZ_ratio, &
! & CS_tev_bb_hjZ_ratio, &
! & CS_tev_pp_vbf_ratio, &
! & BR_hjbb,BR_hjtautau, &
! & BR_hjWW,BR_hjgaga, &
! & BR_hjhihi_nHbynH, &
! & HBresult,chan, &
! & obsratio, ncombined )
!! Obsolete subroutine
!!************************************************************
!
! implicit none
!
! !----------------------------------------input
! integer , intent(in) :: nHdummy
! double precision,intent(in) ::Mh(nHdummy), &
! & CS_lep_hjZ_ratio(nHdummy), &
! & CS_lep_hjhi_ratio_nHbynH(nHdummy,nHdummy), &
! & CS_tev_gg_hj_ratio(nHdummy),CS_tev_bb_hj_ratio(nHdummy), &
! & CS_tev_bg_hjb_ratio(nHdummy), &
! & CS_tev_ud_hjWp_ratio(nHdummy),CS_tev_cs_hjWp_ratio(nHdummy),&
! & CS_tev_ud_hjWm_ratio(nHdummy),CS_tev_cs_hjWm_ratio(nHdummy),&
! & CS_tev_dd_hjZ_ratio(nHdummy),CS_tev_uu_hjZ_ratio(nHdummy), &
! & CS_tev_ss_hjZ_ratio(nHdummy),CS_tev_cc_hjZ_ratio(nHdummy), &
! & CS_tev_bb_hjZ_ratio(nHdummy), &
! & CS_tev_pp_vbf_ratio(nHdummy), &
! & BR_hjbb(nHdummy),BR_hjtautau(nHdummy), &
! & BR_hjWW(nHdummy),BR_hjgaga(nHdummy), &
! & BR_hjhihi_nHbynH(nHdummy,nHdummy)
! !---------------------------------------output
! integer :: HBresult,chan,ncombined
! double precision :: obsratio
! !-----------------------------------------------
!
! call attempting_to_use_an_old_HB_version('part')
!
!end subroutine run_HiggsBounds_part
!!************************************************************
!subroutine run_HiggsBounds_hadr(nHdummy,Mh, &
! & CS_lep_hjZ_ratio,CS_lep_hjhi_ratio_nHbynH, &
! & CS_tev_pp_hj_ratio, CS_tev_pp_hjb_ratio, &
! & CS_tev_pp_hjW_ratio,CS_tev_pp_hjZ_ratio, &
! & CS_tev_pp_vbf_ratio, &
! & BR_hjbb,BR_hjtautau, &
! & BR_hjWW,BR_hjgaga, &
! & BR_hjhihi_nHbynH, &
! & HBresult,chan, &
! & obsratio, ncombined )
!! Obsolete subroutine
!!************************************************************
!
! implicit none
! !----------------------------------------input
! integer,intent(in) :: nHdummy
! double precision,intent(in) :: Mh(nHdummy), &
! & CS_lep_hjZ_ratio(nHdummy), &
! & CS_lep_hjhi_ratio_nHbynH(nHdummy,nHdummy),&
! & CS_tev_pp_hj_ratio(nHdummy), &
! & CS_tev_pp_hjb_ratio(nHdummy), &
! & CS_tev_pp_hjW_ratio(nHdummy), &
! & CS_tev_pp_hjZ_ratio(nHdummy), &
! & CS_tev_pp_vbf_ratio(nHdummy), &
! & BR_hjbb(nHdummy),BR_hjtautau(nHdummy), &
! & BR_hjWW(nHdummy),BR_hjgaga(nHdummy), &
! & BR_hjhihi_nHbynH(nHdummy,nHdummy)
! !---------------------------------------output
! integer :: HBresult,chan,ncombined
! double precision :: obsratio
! !---------------------------------------------
!
! call attempting_to_use_an_old_HB_version('hadr')
!
!
!end subroutine run_HiggsBounds_hadr
!!************************************************************
Index: trunk/HiggsBounds_KW/likelihoods.F90
===================================================================
--- trunk/HiggsBounds_KW/likelihoods.F90 (revision 500)
+++ trunk/HiggsBounds_KW/likelihoods.F90 (revision 501)
@@ -1,995 +1,1056 @@
! This file is part of HiggsBounds
! -TS (23/09/2014)
!******************************************************************
module likelihoods
!******************************************************************
! Some options:
logical :: rescale = .False.
!--
type likelihood2D
integer :: analysisID
integer :: particle_x
character(LEN=45) :: label
character(LEN=3) :: expt
double precision :: lumi, energy
character(LEN=100) :: description
double precision :: mass_sep
double precision :: mass
double precision :: xstep
double precision :: ystep
double precision :: xmin
double precision :: xmax
double precision :: ymin
double precision :: ymax
integer :: size
integer :: xsize
integer :: ysize
double precision, dimension(200,200) :: llhgrid_pred, llhgrid_obs
double precision, dimension(1:3) :: bestfit_pred, bestfit_obs
end type
type threetuple
double precision, dimension(1:3) :: val
end type
integer, parameter :: nllhs = 1
! If more likelihoods become available (nllhs > 1), have to generalize this
! to a list of likelihoods
! CMS-3316 related stuff
integer, parameter :: n2Dslices = 38
type(likelihood2D), dimension(1:n2Dslices) :: CMS_3316_data
contains
!--------------------------------------------------------------
subroutine setup_likelihoods
!--------------------------------------------------------------
use usefulbits, only : file_id_common4, Hneut
use store_pathname, only : pathname
implicit none
integer :: ios,i,j,status
! for CMS-3316 likelihoods
integer :: analysisID = 3316
character(LEN=100) :: label = '[hep-ex] arXiv:1408.3316 (CMS)'
character(LEN=100) :: description = '(pp) -> h -> tautau, using -2ln(L) reconstruction'
character(LEN=100) :: stem = 'Expt_tables/CMStables/CMS_tautau_llh_1408.3316/'
character(LEN=3) :: expt = 'CMS'
integer :: particle_type = Hneut
double precision :: lumi = 19.7D0
double precision :: energy = 8.0D0
double precision :: mass_sep = 0.135D0
double precision, dimension(1:n2Dslices) :: masses, xmin, xmax, ymin, ymax, xstep, ystep
integer, dimension(1:n2Dslices) :: xsize, ysize
masses = (/ 90.0D0, 100.0D0, 110.0D0, 120.0D0, 125.0D0, 130.0D0, 140.0D0, &
& 150.0D0, 160.0D0, 170.0D0, 180.0D0, 190.0D0, 200.0D0, 210.0D0,&
& 220.0D0, 230.0D0, 240.0D0, 250.0D0, 275.0D0, 300.0D0, 325.0D0,&
& 350.0D0, 375.0D0, 400.0D0, 425.0D0, 450.0D0, 475.0D0, 500.0D0,&
& 550.0D0, 600.0D0, 650.0D0, 700.0D0, 750.0D0, 800.0D0, 850.0D0,&
& 900.0D0, 950.0D0, 1000.0D0 /)
xmin = (/ 0.1750000, 0.1000000, 0.0300000, 0.0250000, 0.0250000, 0.0150000, 0.0075000,&
& 0.0050000, 0.0037500, 0.0025000, 0.0025000, 0.0025000, 0.0020000, 0.0020000,&
& 0.0020000, 0.0020000, 0.0015000, 0.0012500, 0.0010000, 0.0005000, 0.0004500,&
& 0.0003750, 0.0003750, 0.0003750, 0.0003250, 0.0002500, 0.0002000, 0.0001250,&
& 0.0001000, 0.0000750, 0.0000700, 0.0000625, 0.0000575, 0.0000500, 0.0000500,&
& 0.0000375, 0.0000375, 0.0000375/)
xmax = (/69.8250000,39.9000000,11.9700000, 9.9750000, 9.9750000, 5.9850000, 2.9925000,&
& 1.9950000, 1.4962500, 0.9975000, 0.9975000, 0.9975000, 0.7980000, 0.7980000,&
& 0.7980000, 0.7980000, 0.5985000, 0.4987500, 0.3990000, 0.1995000, 0.1795500,&
& 0.1496250, 0.1496250, 0.1496250, 0.1296750, 0.0997500, 0.0798000, 0.0498750,&
& 0.0399000, 0.0299250, 0.0279300, 0.0249375, 0.0229425, 0.0199500, 0.0199500,&
& 0.0149625, 0.0149625, 0.0149625/)
ymin = (/ 0.0325000, 0.0200000, 0.0175000, 0.0087500, 0.0087500, 0.0062500, 0.0062500,&
& 0.0050000, 0.0037500, 0.0030000, 0.0025000, 0.0025000, 0.0022500, 0.0020000,&
& 0.0017500, 0.0015000, 0.0015000, 0.0010000, 0.0008750, 0.0005000, 0.0004500,&
& 0.0003000, 0.0003000, 0.0003000, 0.0002500, 0.0002000, 0.0001750, 0.0001500,&
& 0.0001250, 0.0000750, 0.0000700, 0.0000625, 0.0000625, 0.0000625, 0.0000625,&
& 0.0000625, 0.0000625, 0.0000625/)
ymax = (/12.9675000, 7.9800000, 6.9825000, 3.4912500, 3.4912500, 2.4937500, 2.4937500,&
& 1.9950000, 1.4962500, 1.1970000, 0.9975000, 0.9975000, 0.8977500, 0.7980000,&
& 0.6982500, 0.5985000, 0.5985000, 0.3990000, 0.3491250, 0.1995000, 0.1795500,&
& 0.1197000, 0.1197000, 0.1197000, 0.0997500, 0.0798000, 0.0698250, 0.0598500,&
& 0.0498750, 0.0299250, 0.0279300, 0.0249375, 0.0249375, 0.0249375, 0.0249375,&
& 0.0249375, 0.0249375, 0.0249375/)
xstep = (/0.350, 0.200, 0.060, 0.050, 0.050, 0.030, 0.015,&
& 0.010, 0.0075, 0.0050, 0.0050, 0.0050, 0.0040, 0.0040,&
& 0.0040, 0.0040, 0.0030, 0.0025, 0.0020, 0.0010, 0.00090,&
& 0.00075, 0.00075, 0.00075, 0.00065, 0.00050, 0.00040, 0.00025,&
& 0.00020, 0.00015, 0.00014, 0.000125, 0.000115, 0.00010, 0.00010,&
& 0.000075, 0.000075, 0.000075/)
ystep = (/0.065, 0.040, 0.035, 0.0175, 0.0175, 0.0125, 0.0125,&
& 0.010, 0.0075, 0.0060, 0.0050, 0.0050, 0.0045, 0.0040,&
& 0.0035, 0.0030, 0.0030, 0.0020, 0.00175, 0.00100, 0.00090,&
& 0.00060, 0.00060, 0.00060, 0.00050, 0.00040, 0.00035, 0.00030,&
& 0.00025, 0.00015, 0.00014, 0.000125, 0.000125, 0.000125, 0.000125,&
& 0.000125, 0.000125, 0.000125/)
xsize = (/ 200, 200, 200, 200, 200, 200, 200,&
& 200, 200, 200, 200, 200, 200, 200,&
& 200, 200, 200, 200, 200, 200, 200,&
& 200, 200, 200, 200, 200, 200, 200,&
& 200, 200, 200, 200, 200, 200, 200,&
& 200, 200, 200/)
ysize = (/ 200, 200, 200, 200, 200, 200, 200,&
& 200, 200, 200, 200, 200, 200, 200,&
& 200, 200, 200, 200, 200, 200, 200,&
& 200, 200, 200, 200, 200, 200, 200,&
& 200, 200, 200, 200, 200, 200, 200,&
& 200, 200, 200/)
open(file_id_common4,file = trim(adjustl(pathname))// &
& '/Expt_tables/CMS_tautau_llh_1408.3316.binary',form='unformatted')
read(file_id_common4,iostat=ios) CMS_3316_data
if(ios.ne.0)then
do j=1,n2Dslices
call write_metadata(analysisID,label,description,expt,lumi,energy,particle_type,mass_sep,masses(j), &
& xmin(j),xmax(j),xstep(j),xsize(j),ymin(j),ymax(j),ystep(j),ysize(j),CMS_3316_data(j))
call read2Ddata(masses(j),CMS_3316_data(j),stem,'obs',status)
if(status.ne.0) stop 'Error in setup_likelihoods: Data files not found!'
call read2Ddata(masses(j),CMS_3316_data(j),stem,'pred',status)
if(status.ne.0) stop 'Error in setup_likelihoods: Data files not found!'
enddo
! write(*,*) 'Creating binaries...'
rewind(file_id_common4)
write(file_id_common4) CMS_3316_data
endif
close(file_id_common4)
!--------------------------------------------------------------
end subroutine setup_likelihoods
!--------------------------------------------------------------
subroutine write_metadata(analysisID,label,description,expt,lumi,energy,particle_type,&
& mass_sep,mass,xmin,xmax,xstep,xsize,ymin,ymax,ystep,ysize,llh2D)
!--------------------------------------------------------------
use usefulbits, only : small
implicit none
integer, intent(in) :: analysisID,particle_type
character(LEN=*),intent(in) :: label, description,expt
double precision, intent(in) :: mass_sep,mass,xmin,xmax,ymin,ymax,xstep,ystep,lumi,energy
integer, intent(in) :: xsize, ysize
type(likelihood2D), intent(inout) :: llh2D
llh2D%analysisID = analysisID
llh2D%label = trim(adjustl(label))
llh2D%description = trim(adjustl(description))
llh2D%expt = expt
llh2D%energy = energy
llh2D%lumi = lumi
llh2D%particle_x = particle_type
llh2D%mass_sep = mass_sep
llh2D%mass = mass
llh2D%xmin = xmin
llh2D%xmax = xmax
llh2D%xstep = xstep
llh2D%xsize = xsize
llh2D%ymin = ymin
llh2D%ymax = ymax
llh2D%ystep = ystep
llh2D%ysize = ysize
! if(abs((xmax-xmin)/xstep+1.-xsize).gt.small) then
! write(*,*) 'Warning: xsize does not match: ',(xmax-xmin)/xstep+1,' vs. ', xsize
! endif
! if(abs((ymax-ymin)/ystep+1.-ysize).gt.small) then
! write(*,*) 'Warning: ysize does not match: ',(ymax-ymin)/ystep+1,' vs. ', ysize
! endif
!--------------------------------------------------------------
end subroutine write_metadata
!--------------------------------------------------------------
subroutine read2Ddata(mass,llh2D,stem,obspred,status)
!--------------------------------------------------------------
use usefulbits, only : small,file_id_common2
use store_pathname, only : pathname
implicit none
double precision, intent(in) :: mass
type(likelihood2D), intent(inout) :: llh2D
character(LEN=*), intent(in) :: stem,obspred
integer, intent(out) :: status
integer :: size,i,j,k,remove,posx,posy
double precision, dimension(1:3) :: values
type(threetuple), allocatable :: dummylikelihood(:), likelihood(:)
double precision, dimension(200,200) :: llhgrid
double precision, dimension(1:3) :: bestfit
character(LEN=4) :: intstring
character(LEN=100) :: filename
double precision :: modx, mody
write(intstring,"(I4)") int(mass)
if(trim(adjustl(obspred)).eq.'obs') then
filename = 'L_data_SMHb_'//trim(adjustl(intstring))//'.out'
+! filename = 'L_data_b_'//trim(adjustl(intstring))//'.out'
else if(trim(adjustl(obspred)).eq.'pred') then
filename = 'L_asimovSMH_SMHb_'//trim(adjustl(intstring))//'.out'
+! filename = 'L_asimov_b_'//trim(adjustl(intstring))//'.out'
else
stop 'Error in read2Ddata: obspred unknown.'
endif
llh2D%mass = mass
call get_length_of_file(trim(adjustl(pathname))//trim(adjustl(stem))// &
& trim(adjustl(filename)),size,status)
if(status.ne.0) return
open(file_id_common2,file=trim(adjustl(pathname))//trim(adjustl(stem))// &
& trim(adjustl(filename)))
allocate(dummylikelihood(size))
k=0
remove=0
do i=1, size
read(file_id_common2,*) values
if(values(3).eq.0.0D0) then
bestfit = values
remove = remove + 1
else
k=k+1
dummylikelihood(k)%val = values
endif
enddo
close(file_id_common2)
llh2D%size = k
allocate(likelihood(k))
do i=1,k
likelihood(i)%val = dummylikelihood(i)%val
enddo
deallocate(dummylikelihood)
! write(*,*) 'size = ',llh2D%size
llhgrid = 0.0D0
do i=lbound(likelihood,dim=1),ubound(likelihood,dim=1)
posx = nint((likelihood(i)%val(1)-llh2D%xmin)/llh2D%xstep)+1
posy = nint((likelihood(i)%val(2)-llh2D%ymin)/llh2D%ystep)+1
llhgrid(posx,posy) = likelihood(i)%val(3)
enddo
! k=0
! do i=1,llh2D%xsize
! do j=1, llh2D%ysize
! if(llhgrid(i,j).lt.small) then
! k=k+1
!! write(*,*) i,j,llh2D%llhgrid(i,j)
! endif
! enddo
! enddo
! write(*,*) k, 'grid points are unfilled!'
if(trim(adjustl(obspred)).eq.'obs') then
llh2D%llhgrid_obs = llhgrid
llh2D%bestfit_obs = bestfit
else if(trim(adjustl(obspred)).eq.'pred') then
llh2D%llhgrid_pred = llhgrid
llh2D%bestfit_pred = bestfit
endif
deallocate(likelihood)
!--------------------------------------------------------------
end subroutine read2Ddata
!--------------------------------------------------------------
- subroutine get_likelihood(analysisID, jj, t, llh, M_av, nc, obspred)
+ subroutine get_likelihood(analysisID, jj, t, llh, M_av, nc, cbin, obspred, cbin_in)
!--------------------------------------------------------------
use usefulbits, only : dataset
implicit none
type(dataset), intent(in) :: t
integer, intent(in) :: analysisID, jj
character(LEN=*), intent(in) :: obspred
double precision, intent(out) :: llh,M_av
- integer, intent(out) :: nc
+ integer, intent(out) :: nc, cbin
+ integer, optional, intent(in) :: cbin_in
double precision :: cfact1,cfact2
llh = -1.0D0
-
+
select case(analysisID)
case(3316)
- call calcfact_llh(1,jj,t,cfact1,cfact2,M_av,nc)
- llh = get_llh_CMS_tautau_3316(M_av, cfact1, cfact2, obspred)
+ if(present(cbin_in)) then
+! write(*,*) "Calling get_likelihood with cbin_in = ", cbin_in, ", obspred = ", obspred
+ call calcfact_llh(1,jj,t,cfact1,cfact2,M_av,nc,cbin,cbin_in)
+ else
+ call calcfact_llh(1,jj,t,cfact1,cfact2,M_av,nc,cbin)
+ endif
+ if(cbin.ne.0) then
+ llh = get_llh_CMS_tautau_3316(M_av, cfact1, cfact2, obspred)
+ endif
+! write(*,*) jj, cfact1, cfact2, M_av, nc, cbin, "llh = ", llh
+! if(nc.eq.3) then
+! write(*,*) "nc=3: ", jj, cfact1, cfact2, M_av, llh
+! endif
case default
stop 'Error: Unknown analysis ID to retrieve chi^2 value.'
end select
end subroutine get_likelihood
- !--------------------------------------------------------------
+ !--------------------------------------------------------------
subroutine calcpredratio_llh(c,jj,t,M_av,fact,nc,predratio)
!--------------------------------------------------------------
use usefulbits, only : dataset,small,vsmall
implicit none
!prsep(h,n)%tlist,prsep(h,n)%findj,t,axis_i(n),ncomb(n),predratio(n)
type(dataset), intent(in) :: t
integer, intent(in) :: c,jj
double precision, intent(out) :: M_av, predratio, fact
integer, intent(out) :: nc
double precision :: cfact1,cfact2
double precision :: sf, expllh, interval,expllh0
double precision :: diff(2)
- integer :: kk
+ integer :: kk, cbin
double precision sf_low,sf_high,llh_low,llh_high, sf_max
- call calcfact_llh(c,jj,t,cfact1,cfact2,M_av,nc)
+ call calcfact_llh(c,jj,t,cfact1,cfact2,M_av,nc,cbin)
fact = cfact1 + cfact2
expllh = get_llh_CMS_tautau_3316(M_av, cfact1, cfact2, 'pred')
expllh0 = expllh
! write(*,*) 'c = ',c,' jj = ',jj,'nc = ',nc,' M_av = ', M_av, 'ggH = ',cfact1, 'bbH = ',cfact2
if(cfact1.lt.small.and.cfact2.lt.small) then
! write(*,*) 'Very small cross sections. CMS tautau is not competitive, expected llh =',expllh
predratio = -1.0D0
return
endif
! write(*,*) "Find predicted ratio..."
if(expllh.lt.0.0D0) then
predratio = -1.0D0
return
! NEW
endif
! sf = 1.0D0
! kk = 0
! interval = 0.0005D0
! interval = 0.50D0
! diff = 0.0D0
! TS implementation
! do while(abs(expllh-5.99D0).gt.0.01D0)
! kk = kk + 1
! if(abs(diff(2)-diff(1)).le.vsmall) then
! interval = 0.5D0*sf
! else
! interval = min(abs(rand(0)*interval/(diff(2)-diff(1))),0.5D0*sf)
! endif
!
! sf = sf + sign(1.0D0,5.99D0-expllh)*interval
! if (sf.le.0.0D0) exit
!
! diff(2)=abs(expllh-5.99D0)
!
! expllh = get_llh_CMS_tautau_3316(M_av, sf*cfact1, sf*cfact2, 'pred')
!
! diff(1)=abs(expllh-5.99D0)
!
!! write(*,*) kk, sf, interval, diff(1), diff(2), expllh
! enddo
! OS implementation of bisection
! Maximum scale factor corresponds to min_predratio=1/sf_max
sf_max = 1.d0
llh_high = get_llh_CMS_tautau_3316(M_av, sf_max*cfact1, sf_max*cfact2, 'pred')
do while(llh_high-5.99d0.lt.0d0)
sf_max = sf_max*10.d0
llh_high = get_llh_CMS_tautau_3316(M_av, sf_max*cfact1, sf_max*cfact2, 'pred')
enddo
sf_low = 0.d0
sf_high = sf_max
kk=0
sf = (sf_high+sf_low)/2.
llh_low = get_llh_CMS_tautau_3316(M_av, sf_low*cfact1, sf_low*cfact2, 'pred')
llh_high = get_llh_CMS_tautau_3316(M_av, sf_high*cfact1, sf_high*cfact2, 'pred')
! Check that root exists in starting interval
if ((llh_low-5.99)*(llh_high-5.99).ge.0) then
sf = sf_max
else
expllh = get_llh_CMS_tautau_3316(M_av, sf*cfact1, sf*cfact2, 'pred')
! write(*,*) kk, sf_low, sf_high, sf, expllh
do while(abs(expllh-5.99D0).gt.0.01D0.and.kk.lt.100)
! Five digit precision on scale factor in case of oscillating solutions
if ((sf_high-sf_low).LT.1D-5) exit
kk = kk + 1
if((expllh-5.99d0).lt.0) sf_low = sf
if((expllh-5.99d0).gt.0) sf_high= sf
sf = (sf_high+sf_low)/2.
expllh = get_llh_CMS_tautau_3316(M_av, sf*cfact1, sf*cfact2, 'pred')
! write(*,*) kk, sf_low, sf_high, sf, expllh
enddo
endif
if( kk.ge.100) write(*,*) 'Warning: Maximum number of iterations reached with no convergence: predratio might be unreliable.'
! write(*,*) kk, sf_low, sf_high, sf, expllh
! write(*,*) "Ended predratio finder for h",jj," (",nc," combined, average mass = "&
! & ,M_av,") with ",kk,"steps: scalefactor = ", sf,", original -2ln(L) = ",expllh0
! --
! elseif(expllh.gt.5.99D0) then
! do while(expllh.gt.5.99D0)
!
!! Adjust scanning intervals to distance
! kk = kk + 1
! if(abs(expllh-5.99D0) > 2.5D0) then
! interval = 0.05D0
! elseif(abs(expllh-5.99D0) > 0.5D0) then
! interval = 0.01D0
! elseif(abs(expllh-5.99D0) > 0.25D0) then
! interval = 0.001D0
! elseif(abs(expllh-5.99D0) > 0.05D0) then
! interval = 0.0005D0
! endif
!
! sf = sf - interval
! if (sf.le.0.0D0) exit
!
! expllh = get_llh_CMS_tautau_3316(M_av, sf*cfact1, sf*cfact2, 'pred')
!
! write(*,*) kk, sf, interval, expllh
! enddo
! elseif(expllh.lt.5.99D0) then
! do while(expllh.lt.5.99D0)
!
!! Adjust scanning intervals to distance
! kk = kk + 1
! if(abs(expllh-5.99D0) > 2.5D0) then
! interval = 0.05D0
! elseif(abs(expllh-5.99D0) > 0.5D0) then
! interval = 0.01D0
! elseif(abs(expllh-5.99D0) > 0.25D0) then
! interval = 0.001D0
! elseif(abs(expllh-5.99D0) > 0.05D0) then
! interval = 0.0005D0
! endif
!
! sf = sf + interval
! expllh = get_llh_CMS_tautau_3316(M_av, sf*cfact1, sf*cfact2, 'pred')
!
! write(*,*) kk, sf, interval, expllh
! enddo
! endif
if(sf.le.0.0D0) sf = 1.0D-10
predratio = 1./sf
end subroutine calcpredratio_llh
!--------------------------------------------------------------
subroutine check_against_bound_llh(c,jj,t,M_av,nc,obsratio)
!--------------------------------------------------------------
use usefulbits, only : dataset,vsmall
implicit none
type(dataset), intent(in) :: t
integer, intent(in) :: c,jj , nc
double precision, intent(in) :: M_av
double precision, intent(out) :: obsratio
double precision :: cfact1,cfact2
double precision :: sf, obsllh,M_av2, interval,obsllh0
double precision :: diff(2)
- integer :: kk,nc2
+ integer :: kk,nc2,cbin
double precision sf_low,sf_high,llh_low,llh_high,sf_max
- call calcfact_llh(c,jj,t,cfact1,cfact2,M_av2,nc2)
+ call calcfact_llh(c,jj,t,cfact1,cfact2,M_av2,nc2,cbin)
if(abs(M_av2-M_av).gt.vsmall.or.nc.ne.nc2) then
write(*,*) M_av2, M_av, nc, nc2
stop 'Error in subroutine check_against_bound_llh !'
endif
obsllh = get_llh_CMS_tautau_3316(M_av, cfact1, cfact2, 'obs')
obsllh0 = obsllh
! write(*,*) "Find observed ratio..."
if(obsllh.lt.0.0D0) then
obsratio = -1.0D0
return
endif
sf = 1.0D0
kk = 0
interval = 0.50D0
diff = 0.0D0
!! TS implementation
! do while(abs(obsllh-5.99D0).gt.0.05D0)
! kk = kk + 1
! if(abs(diff(2)-diff(1)).le.vsmall) then
! interval = 0.5D0*sf
! else
! interval = min(abs(rand(0)*interval/(diff(2)-diff(1))),0.5D0*sf)
! endif
!
! sf = sf + sign(1.0D0,5.99D0-obsllh)*interval
! if (sf.le.0.0D0) exit
!
! diff(2)=abs(obsllh-5.99D0)
!
! obsllh = get_llh_CMS_tautau_3316(M_av, sf*cfact1, sf*cfact2, 'obs')
!
! diff(1)=abs(obsllh-5.99D0)
!
! write(*,*) kk, sf, interval, diff(1), diff(2), obsllh
!
! enddo
!! OS implementation of bisection
sf_max = 1.d0
llh_high = get_llh_CMS_tautau_3316(M_av, sf_max*cfact1, sf_max*cfact2, 'obs')
do while(llh_high-5.99d0.lt.0d0)
sf_max = sf_max*10.d0
llh_high = get_llh_CMS_tautau_3316(M_av, sf_max*cfact1, sf_max*cfact2, 'obs')
enddo
sf_low = 0.d0
sf_high = sf_max
kk=0
sf = (sf_high+sf_low)/2.
llh_low = get_llh_CMS_tautau_3316(M_av, sf_low*cfact1, sf_low*cfact2, 'obs')
llh_high = get_llh_CMS_tautau_3316(M_av, sf_high*cfact1, sf_high*cfact2, 'obs')
! Check that root exists in starting interval
if ((llh_low-5.99d0)*(llh_high-5.99d0).ge.0) then
sf = sf_max
else
obsllh = get_llh_CMS_tautau_3316(M_av, sf*cfact1, sf*cfact2, 'obs')
! write(*,*) kk, sf_low, sf_high, sf, obsllh
do while(abs(obsllh-5.99D0).gt.0.01D0.and.kk.lt.100)
! Five digit precision on scale factor in case of oscillating solutions
if ((sf_high-sf_low).LT.1D-5) exit
kk = kk + 1
if((obsllh-5.99d0).lt.0) sf_low = sf
if((obsllh-5.99d0).gt.0) sf_high= sf
sf = (sf_high+sf_low)/2.
obsllh = get_llh_CMS_tautau_3316(M_av, sf*cfact1, sf*cfact2, 'obs')
! write(*,*) kk, sf_low, sf_high, sf, obsllh
enddo
endif
if( kk.ge.100) write(*,*) 'Warning: Maximum number of iterations reached with no convergence: obsratio might be unreliable.'
! write(*,*) "Ended obsratio finder with ",kk,"steps: ", sf,obsllh, ", original -2ln(L) = ",obsllh0
! elseif(obsllh.gt.5.99D0) then
! do while(obsllh.gt.5.99D0)
!
!! Adjust scanning intervals to distance
! kk = kk + 1
! if(abs(obsllh-5.99D0) > 3.0D0) then
! interval = 0.05D0
! elseif(abs(obsllh-5.99D0) > 0.5D0) then
! interval = 0.01D0
! elseif(abs(obsllh-5.99D0) > 0.25D0) then
! interval = 0.001D0
! elseif(abs(obsllh-5.99D0) > 0.05D0) then
! interval = 0.0005D0
! endif
!
! sf = sf - interval
! if(sf.le.0.0D0) exit
! obsllh = get_llh_CMS_tautau_3316(M_av, sf*cfact1, sf*cfact2, 'obs')
!
! write(*,*) kk, sf, interval, obsllh
! enddo
! elseif(obsllh.lt.5.99D0) then
! do while(obsllh.lt.5.99D0)
!
!! Adjust scanning intervals to distance
! kk = kk + 1
! if(abs(obsllh-5.99D0) > 3.0D0) then
! interval = 0.05D0
! elseif(abs(obsllh-5.99D0) > 0.5D0) then
! interval = 0.01D0
! elseif(abs(obsllh-5.99D0) > 0.25D0) then
! interval = 0.001D0
! elseif(abs(obsllh-5.99D0) > 0.05D0) then
! interval = 0.0005D0
! endif
!
! sf = sf + interval
! obsllh = get_llh_CMS_tautau_3316(M_av, sf*cfact1, sf*cfact2, 'obs')
!
! write(*,*) kk, sf, interval, obsllh
! enddo
! endif
if(sf.le.0.0D0) sf = 1.0D-10
obsratio = 1./sf
end subroutine check_against_bound_llh
-!--------------------------------------------------------------
- subroutine calcfact_llh(c,jj,t,cfact1,cfact2,M_av,nc)
+!--------------------------------------------------------------
+ subroutine calcfact_llh(c,jj,t,cfact1,cfact2,M_av,nc,cbin,cbin_in)
!--------------------------------------------------------------
! This routine calculates the predicted rates for each process
! Need to generalize this if more likelihoods become available.
use usefulbits, only : dataset, np, div, vvsmall
use theory_BRfunctions
use theory_XS_SM_functions
implicit none
!--------------------------------------input
type(dataset), intent(in) :: t
integer, intent(in) :: c,jj
+ integer, optional, intent(in) :: cbin_in
!-----------------------------------output
double precision, intent(out) :: cfact1, cfact2 ,M_av
- integer, intent(out) :: nc
+ integer, intent(out) :: nc,cbin
!-----------------------------------internal
double precision,allocatable :: mass(:),fact1(:),fact2(:)
integer :: npart,f,j !number of particles
double precision :: numerator, denominator
+ logical, allocatable :: Havail(:)
cfact1=0.0D0
cfact2=0.0D0
if(c.eq.1) then ! this is the CMS tautau likelihood
npart=np( CMS_3316_data(1)%particle_x )
allocate(mass(npart),fact1(npart),fact2(npart))
mass(:)=t%particle( CMS_3316_data(1)%particle_x )%M(:)
fact1= 0.0D0
fact2= 0.0D0
+
+ allocate(Havail(npart))
+
+ if(present(cbin_in)) then
+ call available_Higgses(Havail,npart,cbin_in)
+ else
+ Havail = .True.
+ endif
+
! write(*,*) 'calcfact_llh: mass = ', mass
+ cbin = 0
- do j=1,npart
- if((abs(mass(jj)-mass(j)).le.CMS_3316_data(1)%mass_sep*mass(jj)))then ! &
-! & .and.(mass(jj).le.mass(j)))then
- if((t%gg_hj_ratio(j).ge.0.0D0).and.(t%bb_hj_ratio(j).ge.0.0D0)) then
+ !write(*,*) "Calling calcfact_llh..."
+
+ do j=1,npart
+ if(Havail(j)) then
+ if(abs(mass(jj)-mass(j)).le.CMS_3316_data(1)%mass_sep*mass(jj))then
+! if((abs(mass(jj)-mass(j)).le.CMS_3316_data(1)%mass_sep*100.0D0)) then
+! if(abs(mass(jj)-mass(j)).le.CMS_3316_data(1)%mass_sep*max(mass(jj),mass(j)))then
+ if((t%gg_hj_ratio(j).ge.0.0D0).and.(t%bb_hj_ratio(j).ge.0.0D0)) then
! write(*,*) "Using partonic input for CMS tautau analysis..."
- fact1(j)=t%gg_hj_ratio(j)*XS_lhc8_gg_H_SM(mass(j))*t%BR_hjtautau(j)
- fact2(j)=t%bb_hj_ratio(j)*XS_lhc8_bb_H_SM(mass(j))*t%BR_hjtautau(j)
- else
+ fact1(j)=t%gg_hj_ratio(j)*XS_lhc8_gg_H_SM(mass(j))*t%BR_hjtautau(j)
+ fact2(j)=t%bb_hj_ratio(j)*XS_lhc8_bb_H_SM(mass(j))*t%BR_hjtautau(j)
+ else
! write(*,*) "Using hadronic input for CMS tautau analysis..."
- fact1(j)=t%lhc8%XS_hj_ratio(j)*XS_lhc8_gg_H_SM(mass(j))*t%BR_hjtautau(j)
- fact2(j)=t%lhc8%XS_hjb_ratio(j)*XS_lhc8_bb_H_SM(mass(j))*t%BR_hjtautau(j)
+ fact1(j)=t%lhc8%XS_hj_ratio(j)*XS_lhc8_gg_H_SM(mass(j))*t%BR_hjtautau(j)
+ fact2(j)=t%lhc8%XS_hjb_ratio(j)*XS_lhc8_bb_H_SM(mass(j))*t%BR_hjtautau(j)
+ endif
+ cbin = cbin + 2**(j-1)
+! write(*,*) cbin
endif
- endif
+ endif
enddo
! write(*,*) 'calcfact_llh: fact1 = ', fact1
! write(*,*) 'calcfact_llh: fact2 = ', fact2
if(fact1(jj).le.vvsmall.and.fact2(jj).le.vvsmall)then
!Higgs jj doesn't contribute - wait until another call of this subroutine before
!looking at nearby masses
M_av = mass(jj)
nc=0
cfact1=0.0D0
cfact2=0.0D0
else
!find M_av weighted by the rates (only using higgs which have non-zero fact):
f=0
numerator=0.0D0
denominator=0.0D0
do j=1,npart
if( fact1(j).gt.vvsmall.or.fact2(j).gt.vvsmall )then
f=f+1
numerator=numerator+mass(j)*(fact1(j)+fact2(j))
denominator=denominator+(fact1(j)+fact2(j))
endif
enddo
if(denominator.ne.0.0D0) then
M_av = numerator/denominator
else
M_av = 0.0D0
write(*,*) 'Warning: could not find average mass for CMS tautau analysis!'
endif
nc=f !f will always be > 0 because we've already made sure that fact(jj)>0.0D0
cfact1=sum(fact1)
cfact2=sum(fact2)
endif
- deallocate(mass,fact1,fact2)
+ deallocate(mass,fact1,fact2)
+ deallocate(Havail)
endif
end subroutine calcfact_llh
!--------------------------------------------------------------
function get_llh_from_interpolation(llh2D, xval0, yval0,obspred)
! This routine returns the - 2 Delta ln(llh) from grid interpolation
!--------------------------------------------------------------
use usefulbits, only : small
type(likelihood2D), intent(in) :: llh2D
double precision, intent(in) :: xval0, yval0
character(LEN=*), intent(in) :: obspred
double precision :: get_llh_from_interpolation
integer :: posx1, posx2, posy1, posy2, posxlow, posxup
double precision :: xval, yval, xval1, xval2, yval1, yval2, fvaly1, fvaly2, fval
integer :: i, max_expansion
logical :: lowerxfound, upperxfound
double precision, dimension(200,200) :: llhgrid
if(trim(adjustl(obspred)).eq.'obs') then
llhgrid = llh2D%llhgrid_obs
else if(trim(adjustl(obspred)).eq.'pred') then
llhgrid = llh2D%llhgrid_pred
else
stop 'Error in get_llh_from_interpolation: obspred unknown.'
endif
! This number is the maximal number of tries, that the algorithm steps further up or
! downwards the grid in case of non-existent grid points.
max_expansion = 3
xval = xval0
yval = yval0
if(xval.lt.llh2D%xmin) then
! write(*,*) 'Warning: ggH rate outside grid -- value too small.'
xval = llh2D%xmin
else if (xval.gt.llh2D%xmax) then
! write(*,*) 'Warning: ggH rate outside grid -- value too large.'
xval = llh2D%xmax-small
endif
if(yval.lt.llh2D%ymin) then
! write(*,*) 'Warning: bbH rate outside grid -- value too small.'
yval = llh2D%ymin
else if (yval.gt.llh2D%ymax) then
! write(*,*) 'Warning: ggH rate outside grid -- value too large.'
yval = llh2D%ymax -small
endif
! Get coordinates of neighboring points
posx1 = (xval-llh2D%xmin) / llh2D%xstep + 1
if(posx1.lt.llh2D%xsize) then
posx2 = (xval-llh2D%xmin) / llh2D%xstep + 2
else
write(*,*) "Warning: on the edge of the grid in x-position!"
posx1 = (xval-llh2D%xmin) / llh2D%xstep
posx2 = (xval-llh2D%xmin) / llh2D%xstep + 1
endif
posy1 = (yval-llh2D%ymin) / llh2D%ystep + 1
if(posy1.lt.llh2D%ysize) then
posy2 = (yval-llh2D%ymin) / llh2D%ystep + 2
else
write(*,*) "Warning: on the edge of the grid in y-position!"
posy1 = (yval-llh2D%ymin) / llh2D%ystep
posy2 = (yval-llh2D%ymin) / llh2D%ystep + 1
endif
! Get x,y values at these points
xval1 = (posx1 - 1)*llh2D%xstep + llh2D%xmin
xval2 = (posx2 - 1)*llh2D%xstep + llh2D%xmin
yval1 = (posy1 - 1)*llh2D%ystep + llh2D%ymin
yval2 = (posy2 - 1)*llh2D%ystep + llh2D%ymin
! Do bilinear interpolation and retrieve likelihood value
if(llhgrid(posx1,posy1).ne.0.0D0.and.llhgrid(posx2,posy1).ne.0.0D0) then
fvaly1 = (xval2 - xval) / (xval2 - xval1) * llhgrid(posx1,posy1) + &
& (xval - xval1) / (xval2 - xval1) * llhgrid(posx2,posy1)
else if (llhgrid(posx1,posy1).ne.0.0D0) then
fvaly1 = llhgrid(posx1,posy1)
else if (llhgrid(posx2,posy1).ne.0.0D0) then
fvaly1 = llhgrid(posx2,posy1)
else
! write(*,*) 'Warning in interpolating grid: Both gridpoints in x-direction missing!'
! write(*,*) '(', xval1, yval1, '), (', xval2, yval1,')'
lowerxfound = .False.
upperxfound = .False.
do i=1,max_expansion
if(posx1-i.gt.0) then
if (llhgrid(posx1-i,posy1).ne.0.0D0) then
xval1 = (posx1 - i - 1)*llh2D%xstep + llh2D%xmin
lowerxfound = .True.
posxlow = posx1-i
endif
endif
if(posx2+i.le.ubound(llhgrid,dim=1)) then
if (llhgrid(posx2+i,posy1).ne.0.0D0) then
xval2 = (posx2 + i - 1)*llh2D%xstep + llh2D%xmin
upperxfound = .True.
posxup = posx2+i
endif
endif
if(lowerxfound.or.upperxfound) exit
enddo
if(lowerxfound.and.upperxfound) then
fvaly1 = (xval2 - xval) / (xval2 - xval1) * llhgrid(posxlow,posy1) + &
& (xval - xval1) / (xval2 - xval1) * llhgrid(posxup,posy1)
else if(lowerxfound) then
fvaly1 = llhgrid(posxlow,posy1)
else if(upperxfound) then
fvaly1 = llhgrid(posxup,posy1)
else
write(*,*) 'Found no surrounding points!'
fvaly1 = 99999.0D0
endif
endif
if(llhgrid(posx1,posy2).ne.0.0D0.and.llhgrid(posx2,posy2).ne.0.0D0) then
fvaly2 = (xval2 - xval) / (xval2 - xval1) * llhgrid(posx1,posy2) + &
& (xval - xval1) / (xval2 - xval1) * llhgrid(posx2,posy2)
else if (llhgrid(posx1,posy2).ne.0.0D0) then
fvaly2 = llhgrid(posx1,posy2)
else if (llhgrid(posx2,posy2).ne.0.0D0) then
fvaly2 = llhgrid(posx2,posy2)
else
! write(*,*) 'Warning in interpolating grid: Both gridpoints in x-direction missing!'
! write(*,*) '(', xval1, yval2, '), (', xval2, yval2,')'
lowerxfound = .False.
upperxfound = .False.
do i=1,max_expansion
if(posx1-i.gt.0) then
if (llhgrid(posx1-i,posy2).ne.0.0D0) then
xval1 = (posx1 - i - 1)*llh2D%xstep + llh2D%xmin
lowerxfound = .True.
posxlow = posx1-i
endif
endif
if(posx2+i.le.ubound(llhgrid,dim=1)) then
if (llhgrid(posx2+i,posy2).ne.0.0D0) then
xval2 = (posx2 + i - 1)*llh2D%xstep + llh2D%xmin
upperxfound = .True.
posxup = posx2+i
endif
endif
! write(*,*) i, lowerxfound, upperxfound
if(lowerxfound.or.upperxfound) exit
enddo
if(lowerxfound.and.upperxfound) then
fvaly2 = (xval2 - xval) / (xval2 - xval1) * llhgrid(posxlow,posy2) + &
& (xval - xval1) / (xval2 - xval1) * llhgrid(posxup,posy2)
else if(lowerxfound) then
fvaly2 = llhgrid(posxlow,posy2)
else if(upperxfound) then
fvaly2 = llhgrid(posxup,posy2)
else
! write(*,*) 'Found no surrounding points!'
fvaly2 = 99999.0D0
endif
endif
fval = (yval2 - yval) / (yval2 - yval1) * fvaly1 + &
& (yval - yval1) / (yval2 - yval1) * fvaly2
! write(*,*) 'The interpolated value is ', fval
! we want to obtain - 2 * Delta log(likelihood) \simeq chi^2
get_llh_from_interpolation = 2*fval
return
end function get_llh_from_interpolation
!--------------------------------------------------------------
function get_llh_CMS_tautau_3316(mpred, rate1, rate2, obspred)
!--------------------------------------------------------------
use theory_XS_SM_functions, only : XS_lhc8_gg_H_SM, XS_lhc8_bb_H_SM
use theory_BRfunctions, only : BRSM_Htautau
implicit none
double precision, intent(in) :: mpred, rate1, rate2
character(LEN=*), intent(in) :: obspred
double precision :: get_llh_CMS_tautau_3316
double precision :: ggH1, bbH1, llh1, ggH2, bbH2, llh2
integer :: i,j
if (mpred.lt.CMS_3316_data(1)%mass) then
! write(*,*) 'Warning: predicted mass below lowest mass value of CMS analysis.'
get_llh_CMS_tautau_3316 = -1.0D0
else if (mpred.gt.CMS_3316_data(n2Dslices)%mass) then
! write(*,*) 'Warning: predicted mass above highest mass value of CMS analysis.'
get_llh_CMS_tautau_3316 = -1.0D0
else
do i=lbound(CMS_3316_data,dim=1),ubound(CMS_3316_data,dim=1)-1
if(((mpred - CMS_3316_data(i)%mass).ge.0.0D0).and.(mpred-CMS_3316_data(i+1)%mass.lt.0.0D0)) then
! Rescale to grid mass values by SM rate ratios
if(rescale) then
ggH1 = rate1 * XS_lhc8_gg_H_SM(CMS_3316_data(i)%mass) * BRSM_Htautau(CMS_3316_data(i)%mass) / &
& ( XS_lhc8_gg_H_SM(mpred) * BRSM_Htautau(mpred) )
bbH1 = rate2 * XS_lhc8_bb_H_SM(CMS_3316_data(i)%mass) * BRSM_Htautau(CMS_3316_data(i)%mass) / &
& ( XS_lhc8_bb_H_SM(mpred) * BRSM_Htautau(mpred) )
else
ggH1 = rate1
bbH1 = rate2
endif
llh1 = get_llh_from_interpolation(CMS_3316_data(i),ggH1,bbH1,obspred)
if(rescale) then
ggH2 = rate1*XS_lhc8_gg_H_SM(CMS_3316_data(i+1)%mass)*BRSM_Htautau(CMS_3316_data(i+1)%mass) / &
& ( XS_lhc8_gg_H_SM(mpred) * BRSM_Htautau(mpred) )
bbH2 = rate2*XS_lhc8_bb_H_SM(CMS_3316_data(i+1)%mass)*BRSM_Htautau(CMS_3316_data(i+1)%mass) / &
& ( XS_lhc8_bb_H_SM(mpred) * BRSM_Htautau(mpred) )
else
ggH2 = rate1
bbH2 = rate2
endif
llh2 = get_llh_from_interpolation(CMS_3316_data(i+1),ggH2,bbH2,obspred)
get_llh_CMS_tautau_3316 = llh1 + (mpred - CMS_3316_data(i)%mass)/(CMS_3316_data(i+1)%mass-CMS_3316_data(i)%mass) * &
& (llh2 - llh1)
exit
else if(mpred-CMS_3316_data(i+1)%mass.eq.0.0D0) then
get_llh_CMS_tautau_3316 = get_llh_from_interpolation(CMS_3316_data(i+1),rate1,rate2,obspred)
endif
enddo
endif
return
!--------------------------------------------------------------
end function get_llh_CMS_tautau_3316
!--------------------------------------------------------------
subroutine outputproc_llh(proc,file_id_common,descrip)
!--------------------------------------------------------------
use usefulbits, only : listprocesses
type(listprocesses),intent(in) :: proc
integer, intent(in) :: file_id_common
character(LEN=200), intent(out) :: descrip
character(LEN=45) :: label
integer :: i,j
character(LEN=1) :: jj
i=proc%findi
j=proc%findj
write(jj,'(I1)')j
label='('//trim(CMS_3316_data(1)%label)//')'
descrip = '(pp)->h'//jj//'->tautau, using -2ln(L) reconstruction '//label
end subroutine outputproc_llh
!--------------------------------------------------------------
subroutine get_length_of_file(filename, n, status)
!--------------------------------------------------------------
use usefulbits, only : file_id_common3
implicit none
character(LEN=*),intent(in) :: filename
integer, intent(out) :: n, status
integer :: error
open(file_id_common3,file=trim(adjustl(filename)),action='read',status='old',iostat=status)
if(status.ne.0) then
write(*,*) 'Bad status', status, 'with the following file:'
write(*,*) trim(adjustl(filename))
endif
n = 0
do
read(file_id_common3,*, iostat = error)
if (error == -1) exit
n = n + 1
enddo
close(file_id_common3)
end subroutine get_length_of_file
!--------------------------------------------------------------
-
+ subroutine available_Higgses(Havail,nH,cbin_in)
+ implicit none
+
+ integer, intent(in) :: nH, cbin_in
+ logical, intent(inout) :: Havail(:)
+
+ integer :: c,i,r
+
+ if(size(Havail,dim=1).ne.nH) then
+ write(*,*) "Warning in subroutine available_Higgses: sizes do not match!"
+ Havail = .True.
+ else
+ c = cbin_in
+ do i=1,nH
+ r = c / (2**(nH-i))
+ if(r.eq.1) then
+ Havail(nH-i+1) = .False.
+ elseif (r.eq.0) then
+ Havail(nH-i+1) = .True.
+ else
+ stop "Error in subroutine available_Higgses: not valid binary!"
+ endif
+ c = mod(c,2**(nH-i))
+ enddo
+! write(*,*) "cbin = ",cbin_in, " Havail = ",Havail
+ endif
+
+ end subroutine available_Higgses
!******************************************************************
end module likelihoods
!******************************************************************
\ No newline at end of file

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 8:21 PM (1 d, 5 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806009
Default Alt Text
(106 KB)

Event Timeline