Page MenuHomeHEPForge
Files F11222495

datatables_old.f90
No OneTemporary
Actions

datatables_old.f90
View Options

!******************************************************************
module datatables
!******************************************************************
use usefulbits_HS, only : mutable, mutables, mupeak, peaks, Exptdir
use usefulbits, only : np
implicit none
integer :: ntable1,ntable2
contains
!************************************************************
subroutine initializemutables_blank(tables)
!***********************************************************
! still leaves dat unallocated
integer:: i
type(mutable) :: tables(:)
do i=lbound(tables,dim=1),ubound(tables,dim=1)
tables(i)%id = -1
tables(i)%nx = -1
tables(i)%particle_x = -1
tables(i)%label = ''
tables(i)%desc = ''
tables(i)%expt = ''
tables(i)%lumi = -1.0D0
tables(i)%dlumi = -1.0D0
tables(i)%energy = -1.0D0
tables(i)%Nc = -1
tables(i)%xmax = -1.0D0
tables(i)%xmin = -1.0D0
tables(i)%sep = -1.0D0
tables(i)%deltam = -1.0D0
tables(i)%deltax = -1.0D0
tables(i)%mhchisq = 0
enddo
end subroutine initializemutables_blank
subroutine initializemutables(tables)
use store_pathname_HS
use usefulbits, only: np,Hneut,Hplus,file_id_common3
implicit none
character(LEN=150) :: filename
character(LEN=80) :: datafile(500)
integer n_datafiles
!--------------------------------------input
type(mutable),allocatable :: tables(:)
!-----------------------------------internal
logical :: newtables
integer :: x, n, xbeg, xend, i
character(LEN=pathname_length+150) :: fullfilename
character(LEN=200) :: comment
character(LEN=1) :: firstchar
integer :: col
integer :: ios
integer, allocatable :: skip(:)
!-------------------------------------------
! PRINT *, "Path: ", pathname_HS
OPEN(file_id_common3, file=trim(adjustl(pathname_HS))//"Expt_tables/analyses.txt", form='formatted')
READ(file_id_common3,'(A)',iostat=ios) datafile(1)
if(ios.ne.0)then
CLOSE(file_id_common3)
CALL system('ls -1 -p '//trim(adjustl(pathname_HS))// &
& 'Expt_tables/'//trim(adjustl(Exptdir))//' > '//trim(adjustl(pathname_HS))//'Expt_tables/analyses.txt')
OPEN(file_id_common3, file=trim(adjustl(pathname_HS))//"Expt_tables/analyses.txt", form='formatted')
else
rewind(file_id_common3)
! call get_environment_variable("PWD",fullfilename)
! write(*,*) "Reading "//trim(adjustl(fullfilename))//"/HiggsSignals_implemented_analyses.txt"
endif
! print *, "Using automagic reading of datafiles"
print *, "Reading in the following datafiles from analysis-set "//trim(adjustl(Exptdir))//":"
n = 0
n_datafiles = 0
DO
n = n+1
READ(file_id_common3,'(A)', END=99) datafile(n)
PRINT *, n, datafile(n)
ENDDO
99 n_datafiles = n - 1
CLOSE(file_id_common3)
allocate(tables(n_datafiles),skip(n_datafiles))
call initializemutables_blank(tables)
! print *, "Done initializing data files"
DO n=1,n_datafiles
skip(n)=11
OPEN(file_id_common3, file=trim(adjustl(pathname_HS)) //'Expt_tables/'// &
& trim(adjustl(Exptdir))//'/' // datafile(n))
DO
READ(file_id_common3,'(A)') comment
comment = trim(adjustl(comment))
write(firstchar,'(A1)') comment
if(firstchar.ne.'#') then
exit
else
!! write(*,*) "Skipping line with comment in datafile ",datafile(n)
skip(n)=skip(n)+1
endif
ENDDO
backspace(file_id_common3)
READ(file_id_common3,*) tables(n)%id
print *, tables(n)%id
READ(file_id_common3,'(A)') tables(n)%label
READ(file_id_common3,*) tables(n)%collider, tables(n)%collaboration, tables(n)%expt
READ(file_id_common3,'(A)') tables(n)%desc
READ(file_id_common3,*) tables(n)%energy, tables(n)%lumi, tables(n)%dlumi
READ(file_id_common3,*) tables(n)%particle_x, tables(n)%mhchisq
READ(file_id_common3,*) tables(n)%deltam
READ(file_id_common3,*) tables(n)%xmin, tables(n)%xmax, tables(n)%sep
READ(file_id_common3,*) tables(n)%Nc, tables(n)%eff_ref_mass
allocate(tables(n)%channel_id(tables(n)%Nc))
READ(file_id_common3,*) (tables(n)%channel_id(i),i=1,tables(n)%Nc)
allocate(tables(n)%channel_eff(tables(n)%Nc))
if(tables(n)%eff_ref_mass.ge.0D0) then
READ(file_id_common3,*) (tables(n)%channel_eff(i),i=1,tables(n)%Nc)
else
do i=1,tables(n)%Nc
tables(n)%channel_eff(i)=1.0D0
enddo
READ(file_id_common3,*)
endif
allocate(tables(n)%channel_description(tables(n)%Nc,2))
CLOSE(file_id_common3)
ENDDO
col = 4
do x=1,n_datafiles
IF (tables(x)%deltam.LE.tables(x)%sep) THEN
print *, "Warning: Mass resolution for " // trim(tables(x)%label) // &
" very small - may lead to unreliable results"
endif
tables(x)%nx=nint((tables(x)%xmax-tables(x)%xmin)/tables(x)%sep)+1
allocate(tables(x)%mu(tables(x)%nx,col-1))
allocate(tables(x)%sig(tables(x)%nx,4))
allocate(tables(x)%p0(tables(x)%nx))
allocate(tables(x)%channel_mu(tables(x)%Nc,np(tables(x)%particle_x)))
allocate(tables(x)%channel_systSM(tables(x)%Nc,np(tables(x)%particle_x)))
allocate(tables(x)%channel_syst(tables(x)%Nc,np(tables(x)%particle_x)))
allocate(tables(x)%channel_w(tables(x)%Nc,np(tables(x)%particle_x)))
enddo
!! print *, "Allocated data storage"
open(file_id_common3,file = trim(adjustl(pathname_HS))//'Expt_tables/' // &
& 'mutables.binary',form='unformatted')
read(file_id_common3,iostat=ios) tables(1)%mu, tables(1)%p0
if(ios.eq.0)then
do x=2,n_datafiles
read(file_id_common3) tables(x)%mu, tables(x)%p0
enddo
else
! print *, "No data, going to write"
rewind(file_id_common3)
do x=1, n_datafiles
fullfilename=trim(adjustl(pathname_HS))//'Expt_tables/'//trim(adjustl(Exptdir))//'/' &
//trim(datafile(x))
call read_mutable(tables(x),skip(x),col,fullfilename)
write(file_id_common3) tables(x)%mu, tables(x)%p0
enddo
endif
close(file_id_common3)
deallocate(skip)
end subroutine initializemutables
!************************************************************
subroutine read_mutable(t1,skip,col,fullfilename)
!************************************************************
!--------------------------------------input
type(mutable) :: t1
integer :: skip,col
character(LEN=*) :: fullfilename
!-----------------------------------internal
integer :: i,n , file_id_1
double precision :: xdummy,xdummy_store
!-------------------------------------------
file_id_1 = 666
t1%mu=-1.0D0
t1%p0=-1.0D0
open(file_id_1, file=(trim(fullfilename)))
do i=1,skip
read(file_id_1,*) !skip lines
enddo
xdummy_store = t1%xmin-t1%sep
do i=1,t1%nx
read(file_id_1,*)xdummy,(t1%mu(i,n),n=1,3),t1%p0(i)
! checks that x are evenly spaced as expected
if((abs(xdummy-xdummy_store-t1%sep).gt.1.0D-7) &
& .or.(abs(xdummy-(t1%xmin+dble(i-1)*t1%sep)).gt.1.0D-7))then
!! write(*,*)i,t1%id,xdummy,t1%xmin+dble(i-1)*t1%sep
stop 'error in read_mutable (a1)'
endif
xdummy_store=xdummy
enddo
if(abs(xdummy-t1%xmax).gt.1.0D-7)stop 'error in read_mutable (a2)'
close(file_id_1)
end subroutine read_mutable
!**********************************************************
subroutine setup_mutables
!**********************************************************
!-n.b. Have to take "whichanalyses" into account. BETTER: Here, we read in all available
!- analyses. In a following subroutine we can then make a selection of the relevant
!- analyses we want to look at.
use usefulbits, only : debug,np,not_a_particle
implicit none
!-----------------------------------internal
integer :: i,c, nmutables
call initializemutables(mutables)
! checks that none of the id's are repeated in S95_t1
do i=lbound(mutables,dim=1),ubound(mutables,dim=1)
if( count(mutables%id.eq.mutables(i)%id).ne.1)then
write(*,*)'the id',mutables(i)%id,'is repeated'
stop 'error in setup_mutables (c1)'
endif
enddo
!check to make sure that mutables(i)%particle_x are all particles
do i=lbound(mutables,dim=1),ubound(mutables,dim=1)
if(mutables(i)%particle_x.eq.not_a_particle)then
write(*,*)mutables(i)%id,'particle_x=not_a_particle.'
stop 'error in setup_mutables (d1)'
endif
enddo
do i=lbound(mutables,dim=1),ubound(mutables,dim=1)
! call mutosigma(mutables(i))
call p0tosigma(mutables(i))
enddo
end subroutine setup_mutables
!**********************************************************
subroutine mutosigma(t)
implicit none
integer i
double precision s2_low, s2_up, s2_mean
double precision s1_low, s1_up, s1_mean
type(mutable) :: t
do i=lbound(t%mu,dim=1),ubound(t%mu,dim=1)
if ((t%mu(i,2).ne.0d0).AND.(t%mu(i,1).ne.0d0)) then
s2_low = t%mu(i,2)/(t%mu(i,2)-t%mu(i,1))
s2_up = t%mu(i,2)/(t%mu(i,3)-t%mu(i,2))
s2_mean = t%mu(i,2)/(0.5d0*(t%mu(i,3)-t%mu(i,1)))
s1_low = sigma1s(s2_low)
s1_up = sigma1s(s2_up)
s1_mean = sigma1s(s2_mean)
t%sig(i,1) = s1_low
t%sig(i,2) = s1_up
t%sig(i,3) = s1_mean
else
t%sig(i,1) = 0d0
t%sig(i,2) = 0d0
t%sig(i,3) = 0d0
endif
enddo
end subroutine
!**********************************************************
subroutine p0tosigma(t)
!**********************************************************
implicit none
integer i
double precision s2_low, s2_up, s2_mean
double precision s1_low, s1_up, s1_mean
double precision sigma
type(mutable) :: t
do i=lbound(t%mu,dim=1),ubound(t%mu,dim=1)
sigma = DSQRT(2.D0)*erfinv(1D0-2.D0*t%p0(i))
! write(*,*) i, t%p0(i), sigma
t%sig(i,1) = sigma
t%sig(i,2) = sigma
t%sig(i,3) = sigma
enddo
end subroutine
!**********************************************************
! Finds the peaks with significance above the 'speak' value,
! using the mean significance of upper and lower uncertainty bands
subroutine findpeaks(t, speak, peaks)
!**********************************************************
type(mutable) :: t
type(mupeak),allocatable :: peaks(:)
double precision mass, speak, sigmax, mumax
integer :: i,imin,imax,npeaks,newp, p
imin = lbound(t%sig,dim=1)
imax = ubound(t%sig,dim=1)
! print *, ""
! print *, "Searching for peaks in: ", t%label
! PRINT *, "Peak limit is sigma > ", speak
! First survey the number of peaks present
npeaks = 0
newp = 1
do i=imin,imax
IF (t%sig(i,3).GE.speak) THEN
IF (newp.eq.1) THEN
newp = 0
npeaks = npeaks + 1
ENDIF
ELSE
IF (newp.eq.0) THEN
newp = 1
ENDIF
ENDIF
enddo
! PRINT *, "Total number of peaks found: ", npeaks
allocate(peaks(1:npeaks))
do i=lbound(peaks,dim=1),ubound(peaks,dim=1)
call initpeak_blank(peaks(i))
enddo
newp = 1
p = 1
sigmax = -1d0
do i=imin,imax
mass = t%xmin+(i-1)*t%sep
IF (t%sig(i,3).GE.speak) THEN
! print *, i,mass,t%sig(i,1), t%sig(i,2), t%sig(i,3)
IF (newp.eq.1) THEN
newp = 0
peaks(p)%ilow = i
peaks(p)%mlow = mass
ENDIF
IF (t%sig(i,3).GT.sigmax) THEN
sigmax = t%sig(i,3)
peaks(p)%sig = sigmax
peaks(p)%dmulow = abs(t%mu(i,2)-t%mu(i,1))
peaks(p)%mu = t%mu(i,2)
peaks(p)%dmuup = abs(t%mu(i,3)-t%mu(i,2))
peaks(p)%ipeak = i
peaks(p)%mpeak = mass
peaks(p)%dmlow = peaks(p)%mpeak-peaks(p)%mlow
ENDIF
ELSE
IF (newp.eq.0) THEN
newp = 1
peaks(p)%iup = i-1
peaks(p)%mup = mass-t%sep
peaks(p)%dm = peaks(p)%mup-peaks(p)%mlow
peaks(p)%dmup = peaks(p)%mup-peaks(p)%mpeak
p = p + 1
sigmax = -1d0
mumax = -1d0
ENDIF
ENDIF
enddo
IF (newp.eq.0) THEN
peaks(p)%iup = imax
peaks(p)%mup = mass
peaks(p)%dm = peaks(p)%mup-peaks(p)%mlow
peaks(p)%dmup = peaks(p)%mup-peaks(p)%mpeak
ENDIF
do p=lbound(peaks,dim=1), ubound(peaks,dim=1)
IF (peaks(p)%dm.LT.t%deltam) THEN
peaks(p)%dmlow = peaks(p)%dmlow+(t%deltam-peaks(p)%dm)/2D0
peaks(p)%dmup = peaks(p)%dmup+(t%deltam-peaks(p)%dm)/2D0
peaks(p)%dm = t%deltam
peaks(p)%mlow = peaks(p)%mpeak-peaks(p)%dmlow
peaks(p)%mup = peaks(p)%mpeak+peaks(p)%dmup
ENDIF
!-(TS 06/07/2012) Have to set mass uncertainty of peak always equal to analysis'
peaks(p)%dm = t%deltam
peaks(p)%dlumi = t%dlumi
allocate(peaks(p)%Higgs_comb(np(t%particle_x)))
allocate(peaks(p)%Higgses(np(t%particle_x)))
peaks(p)%Nc = t%Nc
! have to allocate channel info and set them initially to zero!
allocate(peaks(p)%channel_id(peaks(p)%Nc))
peaks(p)%channel_id(:) = t%channel_id(:)
allocate(peaks(p)%channel_mu(peaks(p)%Nc))
allocate(peaks(p)%channel_w(peaks(p)%Nc))
allocate(peaks(p)%channel_systSM(peaks(p)%Nc))
allocate(peaks(p)%channel_syst(peaks(p)%Nc))
do i=1,peaks(p)%Nc
peaks(p)%channel_mu(i)=0.0D0
peaks(p)%channel_w(i)=0.0D0
peaks(p)%channel_systSM(i)=0.0D0
peaks(p)%channel_syst(i)=0.0D0
enddo
enddo
do i=lbound(peaks,dim=1), ubound(peaks,dim=1)
! WRITE(*, '(4I8,6F10.2,4F12.3)') i, peaks(i)%ilow, peaks(i)%ipeak, peaks(i)%iup, peaks(i)%mlow, &
! peaks(i)%mpeak, peaks(i)%mup, peaks(i)%dmlow, peaks(i)%dm, peaks(i)%dmup, &
! peaks(i)%sig, peaks(i)%mu, peaks(i)%dmulow, peaks(i)%dmuup
!! WRITE(*, '(2I8,5F10.3,1I10,5X,A25,1X,A)') i, npeaks, peaks(i)%mpeak,&
!! & peaks(i)%dm, peaks(i)%mu, peaks(i)%dmulow, peaks(i)%dmuup, t%id, t%label, &
!! & trim(t%desc)
enddo
end subroutine
subroutine initpeak_blank(p)
implicit none
type(mupeak) :: p
p%ilow = -1
p%iup = -1
p%ipeak = -1
p%mpeak = -1.0D0
p%mlow = -1.0D0
p%mup = -1.0D0
p%dmlow = -1.0D0
p%dmup = -1.0D0
p%dm = -1.0D0
p%sig = -1.0D0
p%mu = -1.0D0
p%dmuup = -1.0D0
p%dmulow = -1.0D0
end subroutine
function erfinv(x)
implicit none
double precision erfinv
double precision x,w
double precision a, Pi
parameter (a = 8887D0/63473D0)
parameter (Pi = 3.141592D0)
w = DLOG(1D0-x**2)
erfinv = DSQRT(-2D0/(Pi*a)-w/2D0+DSQRT((2D0/(Pi*a)+w/2D0)**2-w/a))
end function
function sigma1s(s2s)
implicit none
double precision sigma1s
double precision s2s,y
IF (s2s.LT.0D0) THEN
sigma1s = s2s
RETURN
ENDIF
IF (s2s.GT.5D0) THEN
sigma1s = s2s
RETURN
ENDIF
y = 1D0-0.5d0*(1d0-DERF(s2s/DSQRT(2d0)))
sigma1s = DSQRT(2D0)*erfinv(y)
end function
subroutine deallocate_mutables
implicit none
integer :: x
do x=lbound(mutables,dim=1),ubound(mutables,dim=1)
deallocate(mutables(x)%mu)
deallocate(mutables(x)%sig)
deallocate(mutables(x)%p0)
deallocate(mutables(x)%channel_mu)
deallocate(mutables(x)%channel_systSM)
deallocate(mutables(x)%channel_syst)
deallocate(mutables(x)%channel_w)
deallocate(mutables(x)%channel_eff)
deallocate(mutables(x)%channel_description)
if(allocated(mutables(x)%Toys_mhobs)) deallocate(mutables(x)%Toys_mhobs)
if(allocated(mutables(x)%Toys_muobs)) deallocate(mutables(x)%Toys_muobs)
enddo
deallocate(mutables)
end subroutine deallocate_mutables
subroutine deallocate_peaks
implicit none
integer :: x
do x=lbound(peaks,dim=1),ubound(peaks,dim=1)
deallocate(peaks(x)%Higgs_comb)
deallocate(peaks(x)%Higgses)
deallocate(peaks(x)%channel_id)
deallocate(peaks(x)%channel_mu)
deallocate(peaks(x)%channel_w)
deallocate(peaks(x)%channel_systSM)
deallocate(peaks(x)%channel_syst)
enddo
deallocate(peaks)
end subroutine deallocate_peaks
subroutine remove_binaries
use store_pathname_HS
implicit none
call system('rm -f '//trim(adjustl(pathname_HS))//'Expt_tables/mutables.binary')
end subroutine remove_binaries
end module datatables
!************************************************************
!************************************************************
! subroutine initializemutables(tables)
!***********************************************************
! use usefulbits, only: Hneut,Hplus,file_id_common2
! implicit none
!
! !--------------------------------------input
! type(mutable) :: tables(:)
! !-----------------------------------internal
! integer :: x,xbeg,xend
! character(len=100),allocatable :: filename(:)
! character(LEN=150) :: fullfilename
! integer :: col
! integer :: ios
! !-------------------------------------------
!
! xbeg=lbound(tables,dim=1)
! xend=ubound(tables,dim=1)
!
! allocate(filename(xbeg:xend))
! x=xbeg-1
!
! !instead, could read in the values of xmin,xmax,sep from the
! !files, but it's kinda nice having them all here to refer to
!
! x=x+1
! tables(x)%id=12021414
! tables(x)%particle_x=Hneut
! tables(x)%expt='ATL'
! tables(x)%label='[hep-ex] arXiv:1202.1414 (ATLAS)'
! tables(x)%energy=7.0D0
! tables(x)%lumi=4.9D0
! tables(x)%SMlike=1
! tables(x)%xmin=110.0D0
! tables(x)%xmax=150.0D0
! tables(x)%sep=0.1D0
! tables(x)%deltax=0.0D0
! filename(x)='12021414'
!
! ! checks we've filled the whole array
! if(x.ne.xend)then
! stop'error in initializetables1 (a)'
! endif
!
! ! do loop to read in mu plot tables
! col=4
! do x=xbeg,xend
! tables(x)%nx=nint((tables(x)%xmax-tables(x)%xmin)/tables(x)%sep)+1
! allocate(tables(x)%mu(tables(x)%nx,col-1))
! allocate(tables(x)%sig(tables(x)%nx,4))
! enddo
!
! open(file_id_common2,file = 'Expt_tables/' // &
! & 'mutables.binary',form='unformatted')
!
! read(file_id_common2,iostat=ios)tables(xbeg)%mu
!
! if(ios.eq.0)then
!
! do x=xbeg+1,xend
! read(file_id_common2)tables(x)%mu
! enddo
!
! else
! rewind(file_id_common2)
! do x=xbeg,xend
! fullfilename='Expt_tables/' &
! & //trim(adjustl(tables(x)%expt))//'tables/' &
! & //trim(filename(x))//'.txt'
!
!
! call read_mutable(tables(x),2,col,fullfilename)
! write(file_id_common2)tables(x)%mu
! enddo
! endif
!
! close(file_id_common2)
!
! deallocate(filename)
!
! end subroutine initializemutables

File Metadata

Mime Type
text/plain
Expires
Wed, May 14, 11:59 AM (1 h, 23 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
5078139
Default Alt Text
datatables_old.f90 (18 KB)

Event Timeline