*file member=acf_CorChiSq library=snoman language=fortran77 date=24:Jul:2002
<<< Subroutine acf_CorChiSq(nHitPMT,FirstRTCPMT,LastRTCPMT,rEvt,udotr,
+ chisq,Probchisq,aLnDetCM,deltaCDF,DISCARD)
* ACF: Calculate correlation chi-square
* Contact: W. Heintzelman, Penn
* Verified:
* Refereed:
* Parameters:-
* ==========
* nHitPMT in Number of hit PMTs used to calculate correl. fn. w/o RTC
* FirstRTCPMT in Index of first PMT used in RTC case
* LastRTCPMT in Index of last PMT used in RTC case
* rEvt in Fitted event radius
* udotr in Value of u.r
* chisq out Correlation chi-square
* Probchisq out Probability of exceeding chi-square value
* aLnDetCM out Natural log of the determinant of the covariance matrix
* deltaCDF out Signed difference between cumulative event correlation
* function and that of the corresponding standard function
* that has the largest absolute value.
IMPLICIT NONE
* Common Block Access:-
* ===================
! To access titles file data for angular correlation function
COMMON/tacfBankPtrs/ LAtacf(8)
INTEGER LAtacf,Ltacf,LDtacf
EQUIVALENCE (LAtacf(6),Ltacf), (LAtacf(7),LDtacf)
* Specification:-
* =============
* Unpack the angular correlation function for this event
* from HBOOK histogram 99, and calculate the chi-square for its fit
* to the corresponding standard correlation function. Also calculate
* the probability of exceeding the observed chi-square value, and the
* maximum difference between the cumulative event correlation
* function and that of the corresponding standard function.
* Revision History:-
* ================
* 4.02 W. Heintzelman First version
* W. Heintzelman Second version
* N. West Convert " -> '.
include 'id_errors.inc'
INCLUDE 'mt.inc'
INCLUDE 'zunit.inc'
* include 'su_mnemonics.inc' (This statement is in acf_evalconst.inc)
include 'acf_evalconst.inc'
* Argument Declarations:-
* =====================
integer nHitPMT,FirstRTCPMT,LastRTCPMT
real rEvt,udotr
real chisq(0:1,MaxEvtType),Probchisq(0:1,MaxEvtType)
+ ,aLnDetCM(0:1,MaxEvtType),deltaCDF(0:1,MaxEvtType)
logical DISCARD
*ENDHEADER
* Program Notes:-
* =============
* Local Variable Declarations:-
* ===========================
integer ireq,nhdr,nwords,jstdfn,iwarnct
integer i,j,k,l,jhist,nstdcases,maxstdcases,nHits,iCase
parameter (maxstdcases=400)
integer lRTC(maxstdcases),number(maxstdcases),leastdiff
integer iEType(maxstdcases)
integer lowbound(maxstdcases),highbound(maxstdcases)
real rmin(maxstdcases),rmax(maxstdcases)
real udotrmin(maxstdcases),udotrmax(maxstdcases)
integer nc(0:1,MaxEvtType)
integer lowest(0:1,MaxEvtType), highest(0:1,MaxEvtType)
real temp(nBina-1),rEvtmod
real fn(nBina),difvec(nBina),binprob(nBina)
logical hexist,startflag,LateInitFlag
character*7 cwwo(0:1)
* real cminv(nBina-1,nBina-1)
* integer jj,kk,imax
integer inttype(2)
character*8 txttype
equivalence (txttype,inttype(1))
data startflag,LateInitFlag/.true.,.false./
data cwwo/'without','with'/
data iwarnct/0/
* External Functions:-
* ===========================
real prob,acf_cdfdif,vdot
SAVE
! Initialize now if call to acf_ini was not made.
if(startflag) then ! Initialization not yet done
LateInitFlag = .true.
write(iqlog,*) 'Initialization in ChiSqCalc at first call.'
goto 1000
endif
100 continue
! write(*,*) 'Entering chisqcalc.'
! write(*,*) ' rEvt =',rEvt
! write(*,*) ' udotr =',udotr
rEvtmod = min(rEvt,600.)
if (rEvt.gt.600.) then
iwarnct = iwarnct + 1
if(iwarnct.le.5)
+ write(iqlog,'(''acf_corchisq: Warning. Event radius of '',
+ f6.1,'' exceeds 600 cm, the maximum for which standard '',
+ ''functions were calculated.'')') rEvt
if(iwarnct.eq.5) write(iqlog,*)
+ ' No further event radius warnings will be printed.'
call ztell(IDZT_ACF_STD_FN_WARN, 0)
endif
do l=0,1
jhist = 99 + l
if ( hexist(jhist) ) then
call HUnpak(jhist,fn,' ',0) ! Copy event correl function into fn
! Convert bin values from correlation function normalization to
! bin probability normalization
do i=1,nBina
binprob(i) = (1./float(NBina))*fn(i)
enddo
if(l.eq.0) then
nHits = nHitPMT
else
nHits= LastRTCPMT - FirstRTCPMT + 1
endif
do k=1,NEvtType
! Find the standard set of bin probabilities that best
! corresponds to this event:
iCase = 0
leastdiff = 99999
i=0
dowhile(i.lt.nstdCases .and. leastdiff.gt.0)
i = i+1
if(l.eq.lRTC(i) .and. k.eq.iEType(i) .and.
+ rEvtmod.ge.rmin(i) .and. rEvtmod.le.rmax(i) .and.
+ udotr.ge.udotrmin(i).and. udotr.le.udotrmax(i)) then
if(nHits.ge.lowbound(i) .and.
+ nHits.le.highbound(i)) then
iCase=i
leastdiff = 0
endif
if(nHits.lt.lowbound(i) .and.
+ (lowbound(i)-nHits).lt.leastdiff ) then
iCase = i
leastdiff = lowbound(i)-nHits
endif
if(nHits.gt.highbound(i) .and.
+ (nHits-highbound(i).lt.leastdiff) ) then
iCase = i
leastdiff = nHits-highbound(i)
endif
endif
enddo
* write(iqlog,*)
* + 'acf_corchisq: l,k,rEvt,udotr,nHits,iCase ... ='
* write(iqlog,*)
* + l,k,rEvt,udotr,nHits,iCase,lRTC(iCase),
* + iEType(iCase),rmin(iCase),rmax(iCase),
* + udotrmin(iCase),udotrmax(icase),
* + lowbound(iCase),highbound(iCase)
if(iCase.eq.0) then
write(iqlog,*)
+ 'acf_chiSqCalc: No corresponding standard case found.'
write(iqlog,*)'l, nHits, EvtType, rEvtmod, udotr =',
+ l,nHits,k,rEvtmod,udotr
DISCARD = .true.
return
endif
* write(iqlog,*) 'Case:',iCase
* write(iqlog,'(''Event correl fn:''/(10f8.5))' ) fn
! Calculate vector of bin deviations from the standard function:
jstdfn = LDtacf + nhdr + nwords*(iCase-1) + 10
! rcons(jstdfn) is the first word, less one, of the standard
! correlation function for this case.
do i=1,nBina
difvec(i) = binprob(i) - rcons(i+jstdfn)
enddo
* write(iqlog,'(''Standard correl fn:''/(10f8.5))' )
* + (rcons(i+jstdfn),i=1,nBina)
* Extract Ln of determinant of covariance matrix from titile bank:
aLnDetCM(l,k) = rcons( jstdfn + NBina + 1 )
j = jstdfn + NBina + 2
* write(iqlog,'
* + (''First entries of cov matrix inv. in triangular form:''
* + /(10e10.3))' ) (rcons(i+j),i=0,nBina)
* write(iqlog,*)
! rcons(j) is the first word of the lower left triangle of the
! inverse of the covariance matrix for this case.
! Multiple difvec by the matrix.
call TRSA(rcons(j),difvec,temp,nBina-1,1)
! Expand inverse of cov matrix from triangular to square form:
* call SymMatrixExpand(cminv,cminvll(1,iCase),NBina-1)
* ! Calculate the chi-square statistic:
* call matrixmult(cminv,difvec,temp,nBina-1,nBina-1,1)
chisq(l,k) = vdot(difvec,temp,nBina-1)
* write(iqlog,*) 'Event l,k,chisq =',l,k,chisq(l,k)
! Determine value of cum chi-square distn fn at calculated
! chi-square:
! Probchisq(l,k) = gammp( 0.5*float(NBina-1) , 0.5*chisq(l))
! Determine value of 1 minus cum chi-square distn fn at calculated
! chi-square:
Probchisq(l,k) = prob(chisq(l,k),NBina-1)
! Determine and accumulate in a histogram the maximum difference
! (over the bins) between the cdf for the event and the cdf of the
! standard case (acf_cdfdif is an external function):
* write(iqlog,'(2i5/(10f10.7))') l,nhits,(rcons(i+jstdfn),i=1,nBina)
deltaCDF(l,k) = acf_cdfdif(binprob,rcons(1+jstdfn),nBina)
enddo
endif
enddo
* write(iqlog,*) 'Leaving chisqcalc.'
return
*----------------------------------------
<<< Entry CorChiSqInit
! Initialization performed either by call to this entry point or by
! tranfer here after first call to main entry point, acf_CorChiSq.
1000 continue
startflag = .false.
write(iqlog,'('
+ // '''Using Angular Correlation Code 30 Nov 2001 Version.'')')
! Set up access to TACF titles bank:
! Declare the size of the area to hold user data to be zero
! to prevent user data from being copied to the COMMON.
LAtacf(8) = 0
IREQ = 2 !Request level = 2, compulsory.
CALL MT_REQUEST_TITLES('TACF',2,IREQ,LAtacf,KSU_CLN)
! write(iqlog,*) 'Pointer to acf titles bank:',LDtacf
! Extract some data items from TACF titles bank:
IF (Ltacf .eq. 0) THEN
write(iqprnt,*)
+ 'acf_corchisq: Fatal error. TACF titles bank not found.'
write(iqlog,*)
+ 'acf_corchisq: Fatal error. TACF titles bank not found.'
stop
endif
nstdCases = icons(LDtacf+1)
nhdr = 1 ! number of words of my header data
write(iqlog,*)
+ 'acf_corchisq: Number of standard cases read =',nstdCases
if(nstdcases.gt.maxstdcases) then
write(iqprnt,*) 'acf_corchisq: Fatal error. '
write(iqlog,*) 'acf_corchisq: Fatal error. '
+ ,'Not able to store all standard cases.'
stop
endif
nwords = 10 + NBina + 1 + NBina*(NBina-1)/2 ! Data storage per case
* write(iqlog,*) 'i,lRTC,lowbound,highbound,number:'
do i=1,NEvtType
do j=0,1
nc(j,i)=0
lowest(j,i)=9999
highest(j,i)=0
enddo
enddo
do j=1,nstdCases
i = LDtacf +nhdr + (j-1)*nwords
inttype(1) = icons(i+1)
inttype(2) = icons(i+2)
iEType(j) = 0
do k=1,NEvtType ! txttype is equivalenced to inttype
if(txttype.eq.EvtType(k)) iEType(j)=k
enddo
if(iEType(j).eq.0) then
write(iqprnt,*) 'acf_corchisq: Fatal error. '
write(iqlog,*) 'acf_corchisq: Fatal error. '
+ ,'Unrecognized event type -- ',txttype
stop
endif
lRTC(j) = icons(i+3)
lowbound(j) = icons(i+4)
highbound(j) = icons(i+5)
number(j) = icons(i+6)
rmin(j) = rcons(i+7)
rmax(j) = rcons(i+8)
udotrmin(j) = rcons(i+9)
udotrmax(j) = rcons(i+10)
* write(iqlog,*)
* + ,lRTC(j),lowbound(j),highbound(j),number(j)
if(udotrmin(j).lt.-1. .or. udotrmax(j).gt.1.) then
write(iqprnt,*) 'acf_corchisq: Fatal error. '
write(iqlog,*) 'acf_corchisq: Fatal error. '
+ ,'Illegal udotr input.', udotrmin(j),udotrmax(j)
stop
endif
if(rmin(j).lt.0. .or. rmax(j).lt.rmin(j)) then
write(iqprnt,*) 'acf_corchisq: Fatal error. '
write(iqlog,*) 'acf_corchisq: Fatal error. '
+ ,'Illegal rmin/rmax input.', rmin(j),rmax(j)
stop
endif
if((lRTC(j).ne.0.and.lRTC(j).ne.1).or.lowbound(j).lt.1)then
write(iqprnt,*) 'acf_corchisq: Fatal error. '
write(iqlog,*) 'acf_corchisq: Fatal error. '
+ ,'Illegal input. lRTC, lowbound = ',lRTC(j),lowbound(j)
stop
endif
nc(lRTC(j),iEType(j)) = nc(lRTC(j),iEType(j)) + 1
lowest(lRTC(j),iEType(j)) =
+ min(lowest(lRTC(j),iEType(j)),lowbound(j))
highest(lRTC(j),iEType(j)) =
+ max(highest(lRTC(j),iEType(j)),highbound(j))
enddo
do i=1,NEvtType
write(iqlog,*) ' ',EvtType(i),' event type:'
do j=0,1
write(iqlog,
+ '(5x,i4,'' standard cases, NHit range '''
+ // ',i4,'' to '',i4,'' '',a7,'' residual time cut'')' )
+ nc(j,i),lowest(j,i),highest(j,i),cwwo(j)
enddo
enddo
* write(iqlog,*)'List of standard Cases'
* write(iqlog,*)'EType RTC lowPMT hiPMT rmin rmax udrmin udrmax'
* do j=1,nstdCases
* write(iqlog,'(2i6,2i4,2f8.2,2f10.6)')
* + iEType(j), lRTC(j), lowbound(j), highbound(j),
* + rmin(j), rmax(j), udotrmin(j), udotrmax(j)
* enddo
if(LateInitFlag) go to 100
return
end
*endfile member=acf_CorChiSq