Program CalcOptLim3
C Get data from a file named "ULinput". ULinput in the present version of the
C Optimum Interval (OI) code contains one DM mass at a time (as opposed to all
C DM masses at once). The input information is listed in a single column:
C * 1st entry in the column: probed parameter (e.g. the coupling of DM to a SM
C particle, or in general a point on the ordinate of the limit curve). An initial
C assumption is to be given here, which will be used as starting point for the
C calculation.
C * 2nd entry: model assumption (e.g. the DM mass, or in general a point on the
C abscissa of the limit curve).
C * All following entries: all events in order of increasing energy. For each event,
C the probability is given for the event having lower energy. So, for example,
C suppose the distribution of events is expected to be uniform from 5 to 100 keV,
C and suppose there are 2 events, one at 10 keV and one at 50 keV. Then given an
C event, it has probability (10-5)/(100-5) = 0.05263 of being below 10 keV and
C probability (50-5)/(100-5) = 0.47368 of being below 50 keV.
C Append the results per model assumption (e.g. the DM mass) to the end of a file
C named "ULoutput". The results are listed in rows:
C * 1st entry in the row: the value of the resulting parameter limit at 90% CL.
C * 2nd entry: the model assumption (e.g. the DM mass).
C * 3rd entry: status flag from the OI code. See comments in UpperLimNew3.f for details.
C * 4th and 5th entry: the values of I in FC(I) for the endpoints of the OI.
C * 6th and 7th entry: the lower and upper bound in case of an excluded range around
C the 90% CL limit returned. If the status flag in the 3rd entry is >= 256, more than
C one solution exists and the sixth and seventh entry define the range of solutions.
C Compile with
C f77 -o CalcOptLim3 CalcOptLim3.f UpperLimNew3.f y_vs_CLf2.f CMaxinfNew2.f ConfLev2.f ConfLevNew2.f Cinf2.f CERN_Stuff2.f
C or "f77" -> "gfortran -frecord-marker=4".
Implicit None
Integer*8 If,N,Iflag,EP
Real*8 sigma,mu0,sigma0,muB,CL,Mw,fc0
C Must increase length FC if more than 500000 events!!!
Real*8 UpperLim,FC(0:500000),FB(0:500000),El
Common/UpperLimcom/EP(2),El(2)
Open(unit=25,form='FORMATTED',file='ULoutput', access='APPEND')
Call FillFC(FC,N,mu0,sigma0,Mw)
C Now fill the background vectors
C Call FillFB(FB,muB)
CL=0.9D0 ! 90% Confidence level
If=1 ! fmin = 0.
muB=0.D0 ! Don't subtract background
C fc0=UpperLim(CL,If,N,FC,muB,FC,Iflag,Ist,Ien)
fc0=UpperLim(CL,If,N,FC,muB,FB,Iflag)
sigma=(sigma0/mu0)*fc0
El(1)=(sigma0/mu0)*El(1)
El(2)=(sigma0/mu0)*El(2)
C sigma = fc0
C sigma is now the upper limit cross section. Iflag is a status; 0 is good.
Write(25,10) sigma,Mw,Iflag,EP(1),EP(2),El(1),El(2)
C Write(25,20) FC
10 Format(E13.6,E13.6,I5,I16,I16,E13.6,E13.6)
C 20 Format(E13.6)
Stop
End
C
Subroutine FillFC(FC,N,mu0,sigma0,Mw)
C FillFC reads from a table "ULinput" (see header of this file for its format).
Implicit None
Real*8 mu0,Mw,sigma0
Real*8 FC(0:500000) ! FC(I) contains the probability for the I'th event.
Integer*8 N ! The number of events read in
Open(unit=50,form='FORMATTED',status='OLD',file='ULinput')
Read(50,30) sigma0
Read(50,20) Mw
Read(50,30) mu0
N=0
10 continue
N=N+1
Read(50,30,END=100) FC(N)
20 Format(F10.6)
Go to 10
30 Format(E12.6)
Go to 10
100 N = N-1
Return
End
Subroutine FillFB(FB,muB)
C Same as FillFC, but reads from a table "BGinput" for the background vectors. If no
C backgrounds are subtracted, BGinput simply contains one entry: 0.0e+00.
Implicit None
Real*8 muB
Real*8 FB(0:500000) ! FB(I) contains the probability for the I'th event.
Integer*8 N ! The number of events read in
Open(unit=50,form='FORMATTED',status='OLD',file='BGinput')
Read(50,30) muB
N=0
10 continue
N=N+1
Read(50,30,END=100) FB(N)
20 Format(F10.6)
Go to 10
30 Format(E12.6)
Go to 10
100 N = N-1
Return
End