Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F11222317
add_M0_and_Mexcl_to_clsb_tables.F90
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
11 KB
Subscribers
None
add_M0_and_Mexcl_to_clsb_tables.F90
View Options
!-------------------------------------------------------------------
program add_M0_and_Mexcl_to_clsb_tables
!-------------------------------------------------------------------
! This program find M0 and Mexcl for each point in the chi^2 tables.
! It follows the procedure outlined in section 1.1.3 of description_chisquaredextenstion.ps
! It uses files containing chi^2 values, as generated by transform_clsb_tables.f90
! with the setting MexclM0=.False.
! It creates files containing
! massj massi log10sf Mexcl M0
! You should then run transform_clsb_tables.f90 with the setting MexclM0=.True.
! to create files which include columns of chi^2, Mexcl and M0, which are suitable
! for distribution with the HiggsBounds package.
! Fixme: currently reads in whatever readclsbfiles is set to read-can change by #define oldformat). But: should make
! readclsbfiles more general, so can give filename as input and it will count number f columns itself.
! Then can use it in transform_clsb_usefulbits.f90 as well.
! It is compiled and run using:
! gfortran -L../HiggsBounds_KW/ -I../HiggsBounds_KW/ add_M0_and_Mexcl_to_clsb_tables.F90 -o add_M0_and_Mexcl_to_clsb_tables.exe -lHB
! ./add_M0_and_Mexcl_to_clsb_tables.exe
!-------------------------------------------------------------------
!#define testnewmethod
use S95tables_type2
use S95tables_type3
use usefulbits, only : vsmall
use interpolate
use usefulbits, only : pr,vsmall
use S95tables
use theory_BRfunctions, only : setup_BRSM
implicit none
!----------------------------------
integer :: i,j,k,x,ji,z,n,zi,b,nsf,y,jj
integer :: itot,findmass,tablenum,n_xy_combinations
integer :: clsbtable_selection_id
double precision :: M0,Mexcl,log10sf,Mi_av,Mj_av,M_dummy,Mselect,Mi_target,Mi_middle
double precision :: m_interpol(n_points_max),m_temp
integer :: n_points
double precision :: fact,massi_in,massj_in,xmin1,test_var
integer :: zi_selection(2),ftype_selection(2)
double precision,allocatable :: store_results(:,:)
integer,allocatable :: relevent_points(:,:)
character(LEN=10) :: createtime
character(LEN=8) :: createdate
type(table1) :: f_t1_for_findM0
type(table2) :: slices_t2(2)
double precision :: sepmultfactor
double precision :: valueoutsidetable
integer :: vmassm1orm2
double precision :: m1_at_ref_point_1,m2_at_ref_point_1,m1_at_ref_point_2,m2_at_ref_point_2
!----------------------------------
allocate(clsb_t3(ntable3))
call setup_BRSM !needed for setup_S95tables
call setup_S95tables
call initializetables_type2_blank(clsbslices_t2)
call initializetables_type3_blank(clsb_t3)
call initializetables3(clsb_t3)
call fillt3needs_M2_gt_2M1(clsb_t3,S95_t2)
tablenum=t3elementnumberfromid(clsb_t3,1900)
do z=lbound(clsb_t3,dim=1),ubound(clsb_t3,dim=1)
!do z=tablenum,tablenum
!write(*,*)'take this out!!'
!This is where the results will be outputted to
open(18,file='/home/Karina/Work/HiggsBounds/clsb_tables/csboutput_trans_test/'// &
& trim(adjustl(filename(z)))//'_Mexcl_M0.txt')
clsbtable_selection_id = clsb_t3(z)%id
write(*,*)'hello, about to read file'
call readclsbfiles(clsbtable_selection_id)
write(*,*)'hello, read file'
select case(clsb_t3(z)%type_of_assoc_table)
case(1)
f_t1_for_findM0%xmin = clsb_t3(z)%xmin1
f_t1_for_findM0%xmax = clsb_t3(z)%xmax1
f_t1_for_findM0%sep = clsb_t3(z)%sep1
f_t1_for_findM0%nx = clsb_t3(z)%nx1
allocate(f_t1_for_findM0%dat(f_t1_for_findM0%nx,1))
allocate(store_results(f_t1_for_findM0%nx,5))
Mj_av=clsb_t3(z)%xmin2
do zi=1,clsb_t3(z)%nz
!do zi=91,91
!write(*,*)'take this out (2)!!'
!picks out a 1D table from the 3D table, parallel to z-axis of table
call fill_t1dat_from_t3(f_t1_for_findM0,clsb_t3(z),zi,clsb_t3_dat_3rdcol_chisq)
!write(*,*)'hello',f_t1_for_findM0%dat
log10sf=clsb_t3(z)%zmin+clsb_t3(z)%zsep*dble(zi-1)
Mexcl=-1.0D0
M0=-1.0D0
do findmass=1,2
massi_in=f_t1_for_findM0%xmin
massj_in=Mj_av
select case(findmass)
case(1) !Mexcl
fact=2.7D0
case(2) !M0
fact=1.0D-2
case default
stop'uh oh problem 445'
end select
!finds all the points where the chisq equals 'fact':
call interpolate_1D_inv(fact,f_t1_for_findM0,1,massi_in,m_interpol,M_dummy,n_points)
!old method (have checked it gives same result):
!call find_M0(3,z,clsb_t3_dat_3rdcol_chisq,fact,massi_in,massj_in,log10sf,m_interpol,M_dummy,n_points)
if(n_points.gt.n_points_max)then
stop'problem: n_points.gt.n_points_max'
endif
Mi_middle=0.5D0*(f_t1_for_findM0%xmin+f_t1_for_findM0%xmax)
do i=1,f_t1_for_findM0%nx
!do i=1158,1158
!write(*,*)'take this out (3)!!'
Mi_av=f_t1_for_findM0%xmin+f_t1_for_findM0%sep*dble(i-1)
select case(findmass)
case(1) !Mexcl
Mi_target=Mi_av
case(2) !M0
if(store_results(i,4).gt.vsmall)then !recall store_results(i,4) contains Mexcl
Mi_target=store_results(i,4)
else
Mi_target=Mi_middle
endif
case default
stop'uh oh problem 446'
end select
Mselect=-1.0D8
do b=1,ubound(m_interpol,dim=1)
if(m_interpol(b).gt.0.0D0)then
if(abs(m_interpol(b)-Mi_target).lt.abs(Mselect-Mi_target))then
Mselect=m_interpol(b) !want Mselect to be as near to Mi_target as possible
endif
endif
enddo
store_results(i,1)=Mj_av
store_results(i,2)=Mi_av
store_results(i,3)=log10sf
select case(findmass)
case(1) !Mexcl
store_results(i,4)=Mselect
case(2) !M0
store_results(i,5)=Mselect
case default
stop'uh oh problem 446'
end select
enddo
enddo
do i=1,f_t1_for_findM0%nx
write(18,*)(store_results(i,n),n=1,5)
if((store_results(i,4).gt.0.0D0).and.(store_results(i,5).lt.0.0D0))then
write(*,*)Mj_av,Mi_av,log10sf,store_results(i,4),store_results(i,5)
stop'there has been a problem(1)'
endif
enddo
enddo
deallocate(f_t1_for_findM0%dat)
close(18)
case(2)
call date_and_time(createdate,createtime)
write(*,*)'about to find Mexcl', &
& trim(createtime(1:2)//':'//createtime(3:4))//':'//createtime(5:6)//':'//createtime(7:8)//':'//createtime(9:10);call flush(6)
n_xy_combinations=clsb_t3(z)%nx2*(clsb_t3(z)%nx1+int(clsb_t3(z)%nx1/clsb_t3(z)%nx2))/2
allocate(store_results(n_xy_combinations,5))
#ifdef testnewmethod
Mi_middle=0.5D0*(clsb_t3(z)%xmin2+clsb_t3(z)%xmax2)
#else
Mi_middle=0.5D0*(clsb_t3(z)%xmin1+clsb_t3(z)%xmax1)
#endif
!write(*,*)'take this out! zi'
!do zi=38,38
do zi=1,clsb_t3(z)%nz
zi_selection=zi
ftype_selection=clsb_t3_dat_3rdcol_chisq
!call fill_clsbslices_t2_from_slices_of_t3(clsb_t3(z),zi_selection,ftype_selection)
call fill_slices_t2_from_slices_of_t3(clsb_t3(z),zi_selection,ftype_selection,slices_t2)
y=1
ji=0
log10sf=clsb_t3(z)%zmin+clsb_t3(z)%zsep*dble(zi-1)
!write(*,*)'take this out! j'
!do j=19,19
do j=1,clsb_t3(z)%nx2
Mj_av=clsb_t3(z)%xmin2+clsb_t3(z)%sep2*dble(j-1)
itot=int((clsb_t3(z)%xmin2+clsb_t3(z)%sep2*(j-1))*(clsb_t3(z)%xmax1/clsb_t3(z)%xmax2)/clsb_t3(z)%sep2)
!write(*,*)'take this out! i'
!do i=2,2
do i=1,itot
ji=ji+1
Mi_av=clsb_t3(z)%xmin1+clsb_t3(z)%sep1*dble(i-1)
store_results(ji,1)=Mj_av
store_results(ji,2)=Mi_av
store_results(ji,3)=log10sf
sepmultfactor=0.1D0
valueoutsidetable=-1.0D0
#ifdef testnewmethod
select case(clsb_t3(z)%type_of_assoc_table)
case(1)
vmassm1orm2=1
xlower=clsbslices_t2(y)%xmin1
xhigher=clsbslices_t2(y)%xmax1
m1_at_ref_point_1=xlower
m2_at_ref_point_1=clsbslices_t2(y)%xmin2
m1_at_ref_point_2=xhigher
m2_at_ref_point_2=clsbslices_t2(y)%xmin2
case(2)
vmassm1orm2=2
xlower=clsbslices_t2(y)%xmin2
xhigher=clsbslices_t2(y)%xmax2
m1_at_ref_point_1=Mi_av
m2_at_ref_point_1=xlower
m1_at_ref_point_2=Mi_av
m2_at_ref_point_2=xhigher
!call f_at_const_m1(clsbslices_t2(y),clsbslices_t2(y)%xmin2,clsbslices_t2(y)%xmax2,1,Mi_av,Mj_av,f_t1_for_findM0, &
! & -1.0D0,clsb_t3(z)%sep2)
case default
stop'problem here (31)'
end select
#else
vmassm1orm2=1
xlower=clsbslices_t2(y)%xmin1
xhigher=clsbslices_t2(y)%xmax1
m1_at_ref_point_1=0.0D0
m2_at_ref_point_1=0.0D0
m1_at_ref_point_2=Mi_av
m2_at_ref_point_2=Mj_av
!call f_at_const_massratio(clsbslices_t2(y),clsbslices_t2(y)%xmin1,clsbslices_t2(y)%xmax1,1,Mi_av,Mj_av,f_t1_for_findM0, &
! & -1.0D0,0.1D0)
#endif
call f_from_t2(slices_t2(y),m1_at_ref_point_1,m2_at_ref_point_1,m1_at_ref_point_2,m2_at_ref_point_2, &
& vmassm1orm2,xlower,xhigher,sepmultfactor,1, &
& f_t1_for_findM0,valueoutsidetable)
do findmass=1,2
select case(findmass)
case(1) !Mexcl
fact=2.7D0
Mi_target=Mi_av
case(2) !M0
fact=1.0D-2
if(store_results(ji,4).gt.vsmall)then!store_results(i,4) contains Mexcl
Mi_target=store_results(ji,4)
else
Mi_target=Mi_middle
endif
case default
stop'uh oh problem 447'
end select
call interpolate_1D_inv(fact,f_t1_for_findM0,1,Mi_target,m_interpol,Mselect,n_points)
select case(findmass)
case(1) !Mexcl
store_results(ji,4)=Mselect
case(2) !M0
store_results(ji,5)=Mselect
if((store_results(ji,4).gt.0.0D0).and.(store_results(ji,5).lt.0.0D0))then
!do x=1,ubound(f_t1_for_findM0%dat,dim=1)
! write(*,*)f_t1_for_findM0%xmin+f_t1_for_findM0%sep*dble(x-1),f_t1_for_findM0%dat(x,1)
!enddo
write(*,*)(store_results(ji,n),n=1,5)
write(*,*)'j,i,zi',j,i,zi
#ifdef testnewmethod
store_results(ji,5)=clsbslices_t2(y)%xmax2;
write(*,*)'store_results(ji,5) is now ',store_results(ji,5)
write(*,*)'where chi^2=',f_t1_for_findM0%dat(f_t1_for_findM0%nx-1,1)
write(*,*)'there has been a problem(2)'
!stop'there has been a problem(2)'
#else
stop'there has been a problem(2)'
#endif
endif
case default
stop'uh oh problem 448'
end select
enddo
deallocate(f_t1_for_findM0%dat)
enddo
enddo
if(ji.ne.n_xy_combinations)stop'problem here'
do x=1,n_xy_combinations
write(18,*)(store_results(x,n),n=1,5)
enddo
!do x=1,ubound(clsbslices_t2,dim=1)
! deallocate(clsbslices_t2(x)%dat)
!enddo
do x=1,ubound(slices_t2,dim=1)
deallocate(slices_t2(x)%dat)
enddo
enddo
call date_and_time(createdate,createtime)
write(*,*)'found Mexcl', &
& trim(createtime(1:2)//':'//createtime(3:4))//':'//createtime(5:6)//':'//createtime(7:8)//':'//createtime(9:10);call flush(6)
case default
stop'problem in program add_M0_and_Mexcl_to_clsb_tables'
end select
deallocate(store_results)
enddo
stop'stop here for the moment (just before deallocate statements in program add_M0_and_Mexcl_to_clsb_tables)'
do x=lbound(clsb_t3,dim=1),ubound(clsb_t3,dim=1)
deallocate(clsb_t3(x)%dat)
enddo
deallocate(filename)
deallocate(clsb_t3)
end program add_M0_and_Mexcl_to_clsb_tables
File Metadata
Details
Attached
Mime Type
text/plain
Expires
Wed, May 14, 11:45 AM (5 h, 17 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5111506
Default Alt Text
add_M0_and_Mexcl_to_clsb_tables.F90 (11 KB)
Attached To
rHIGGSBOUNDSSVN higgsboundssvn
Event Timeline
Log In to Comment