*file member=acf_AngCorFn library=snoman language=fortran77 date=21:Mar:2002
<<< subroutine acf_AngCorFn
+ (nHitPMT,PMTpos,FirstRTCPMT,LastRTCPMT,iFitIndx,histflag,etaAvg)
* ACF: Calculate individual event angular correlation functions;
* Accumulate distributions of pair angles and angular correlation
* functions for total sample;
* Normalize and store distributions in histograms at termination.
* Contact: W. Heintzelman, Penn
* Verified:
* Refereed:
* Parameters:-
* ==========
* nHitPMT in Number of unique PMTs that were hit
* PMTpos in Locations of hit PMTs relative to fitted event vertex
* FirstRTCPMT in Index of the first PMT in the residual time window
* LastRTCPMT in Index of the last PMT in the residual time window
* iFitIndx in Index of fitter used for event location
* (0 is an allowed value)
* histflag in .true. if total sample histograms are to be output
* etaAvg(2) out Average value of the pair angle, without & with RTC
* ioffset in Offset to numbering of output histograms.
* Input parameter at call to entry AngCorTrm.
* Common Block Access:-
* ===================
* COMMON /ZUNIT/ IQPRNT,IQLOG
* Specification:-
* =============
* Calculate the angular correlation function from locations of the hit
* PMTs relative to the fitted event vertex, and store the function
* in an HBOOK histogram, both without and with the time exclustion
* window, and with the time exclusion window with nBinc bins.
* Incidently, calculate the average pair angle, etaAvg, referred to
* by others as theta_ij.
* Program Notes:-
* =============
* In array PMTpos, the PMTs are listed in order by the time of the hit.
* Revision History:-
* ================
* 4.02 W. Heintzelman First version
* W. Heintzelman Add histogram 101, which is used in acf_fnip
* 21 Mar 2002
IMPLICIT NONE
INCLUDE 'zunit.inc'
include 'constants.inc'
include 'acf_evalconst.inc'
* Argument Declarations:-
* =====================
integer nHitPMT,FirstRTCPMT,LastRTCPMT,iFitIndx,ioffset
real PMTpos(3,maxhits),etaAvg(0:1)
logical histflag
*ENDHEADER
* Local Variable Declarations:-
* ===========================
integer i,j,l,nAngCorEvts(0:NFitters),last,npair(0:1),nHits
real radj(3,maxhits),rmag(maxhits),Err(NBinc),temp(NBinc)
real sinEta,cosEta,eta,delcosa,delcosb,delcosc
real wta(0:1),wtb(0:1),wtc(0:1),ff,cv(3),etasum(0:1)
logical startflag
character*40 histname
real vdot ! External function
data startflag/.true./
SAVE
if(startflag) then
startflag = .false.
! Histograms used to accumulate correlation function for each
! event:
call hbook1(99,
+ 'Event Correlation Fn vs. Cos(eta)',NBina,-1.,1.,0.)
! For results with Residual Time Cut:
call hbook1(100,
+ 'Event Correlation Fn vs. Cos(eta)(RTC)',NBina,-1.,1.,0.)
! Following histogram used by acf_fnip
call hbook1(101,
+ 'Event Correlation Fn vs. Cos(eta)(RTC)',NBinc,-1.,1.,0.)
do i=0,NFitters
! Initialize histogram storage to accumulate overall number of
! photon pairs in each bin:
call HistInit(4*i+1,0.,180.,NBinb)
call HistInit(4*i+2,0.,180.,NBinb)
! Initialize histogram storage to accumulate overall correlation
! functions:
call HistInit(4*i+3,-1.,1.,NBinc)
call HistInit(4*i+4,-1.,1.,NBinc)
! Initialize counter of number of events accumulated for the
! overall correlation function. (Needed so it can be scaled
! correctly at the end of the run, by call to AngCorTrm)
nAngCorEvts(i) = 0
enddo
delcosa = 2./float(NBina)
delcosb = 2./float(NBinb)
delcosc = 2./float(NBinc)
endif
* write(iqlog,*)'Entering acf_angcorfn --'
* write(iqlog,*)'FirstRTCPMT,LastRTCPMT = ',FirstRTCPMT,LastRTCPMT
call HReset(99,' ') ! Reset inidividual event histograms
call HReset(100,' ')
call HReset(101,' ')
if(histflag) nAngCorEvts(iFitIndx) = nAngCorEvts(iFitIndx) + 1
do j=1,nHitPMT !Calculate normalized PMT relative position vectors
rmag(j) = sqrt( vdot(PMTpos(1,j),PMTpos(1,j),3) )
if(rmag(j).ne.0.) then
call vscale(PMTpos(1,j), 1./rmag(j), radj(1,j) ,3)
else
! write(*,*)'Event',evtno,', PMT',j
write(iqlog,*)' Event vertex coincides with PMT.'
endif
enddo
* write(iqlog,'(''angcorfn: nHitPMT='',i6/(i6,3f10.3))')
* + nHitPMT,(j,(PMTpos(i,j),i=1,3),j=1,nHitPMT)
* write(iqlog,*)
do l=0,1
if(l.eq.0) then
nHits = nHitPMT
else
nHits = lastRTCPMT - FirstRTCPMT +1
endif
ff = 4./(float(nHits)*float(nHits-1))
wta(l) = ff/delcosa
wtb(l) = ff/delcosb
wtc(l) = ff/delcosc
enddo
! Initialize variables for computing average pair angle
do l=0,1
npair(l)=0
etasum(l) = 0.
enddo
do j=1,nHitPMT
! if(mod(j,100).eq.0) write(*,*) j,i
do i=j+1,nHitPMT
! write(*,*) 'Correlfn: at B. j,i =',j,i
cosEta = vdot(radj(1,j),radj(1,i),3)
if(abs(cosEta).lt.0.8) then
sinEta = sqrt(1.-cosEta**2)
eta = raddeg*acos(cosEta)
else
call cross(radj(1,j),radj(1,i),cv)
sinEta = sqrt(cv(1)**2 + cv(2)**2 +cv(3)**2)
sinEta = amin1( sinEta , 1.0)
eta = raddeg*asin(sinEta)
if(cosEta.lt.0.) eta = 180. - eta
! Maximum allowed value of 179.99 in the following
! statement so that values of 180.0 get put into the
! last bin rather than into overflow.
if(eta.gt.179.99) then
eta = 179.99
endif
! Calculate cosEta more accurately:
cosEta = sign( sqrt(1.-sinEta**2) , cosEta)
endif
if(rmag(i).gt.0. .and. rmag(j).gt.0.) then
last=0
! Test if both PMTs are in the residual time window:
if(j.ge.FirstRTCPMT .and. i.le.LastRTCPMT) then
call hfill(101,cosEta,0.,wtc(1))
last=1
endif
do l=0,last
call hfill(99+l,cosEta,0.,wta(l))
if(histflag) then
call Histaccum(4*iFitIndx+1+l,eta,1.)
call Histaccum(4*iFitIndx+3+l,cosEta,wtc(l))
endif
etasum(l) = etasum(l)+eta
npair(l)= npair(l)+1
* if(l.eq.1) write(iqlog,
* + '(''Cor Fn entry: '',2i3,7f10.7)')
* + j,i,coseta,(radj(k,j),k=1,3),(radj(k,i),k=1,3)
enddo
endif
enddo
enddo
do l=0,1
etaAvg(l)=etasum(l)/float(npair(l))
enddo
return
*------------------------------------------------------------
! Run termination
<<< entry AngCorTrm(ioffset)
if(.not.startflag) then
do l=0,1
! Delete histograms 99 and 100 so they don't get written out:
call HDelet(99+l)
enddo
if(ioffset.lt.0) then
write(iqlog,*)'Error: File number offset less than zero.'
write(iqlog,*)' ACF Histograms will not be written.'
return
endif
* write(iqlog,'(5x,''Events accumulated for various fitters '',
* + 8i9)') nAngCorEvts
write(iqlog,*)
do j=0,NFitters
if(nAngCorEvts(j).gt.0) then
write(iqlog,'(5x,'
+ // '''AngCorTrm: Events accumulated for fitter '','
+ // 'i1,'' ='',i9)') j,nAngCorEvts(j)
ff=1./nAngCorEvts(j)
! Initialize histograms to output overall distributions:
write(histname,
+ '(''PMT Pairs vs. Eta(deg), Fitter'',i2)') j
call hbook1(ioffset+4*j+1,histname,NBinb,0.,180.,0.)
write(histname,
+ '(''PMT Pairs vs. Eta(deg)(RTC), Fitter'',i2)') j
call hbook1(ioffset+4*j+2,histname,NBinb,0.,180.,0.)
write(histname,
+ '(''Correl Fn vs. Cos(eta), Fitter'',i2)') j
call hbook1(ioffset+4*j+3,histname,nbinc,-1.,1.,0.)
write(histname,
+ '(''Correl Fn vs. Cos(eta) (RTC), Fitter'',i2)') j
call hbook1(ioffset+4*j+4,histname,nbinc,-1.,1.,0.)
do l=0,1
! Retrieve and store pair angle distribution in an HBOOK
! histogram:
call HistRet(4*j+1+l,temp,Err)
call HPak(ioffset+4*j+1+l,temp)
call HPakE(ioffset+4*j+1+l,Err)
! Retrieve combined-event correlation function, normalize
! and store in HBOOK histogram:
call HistRet(4*j+3+l,temp,Err)
do i=1,NBinc
temp(i) = ff*temp(i)
Err(i) = ff*Err(i)
enddo
call HPak(ioffset+4*j+3+l,temp)
call HPakE(ioffset+4*j+3+l,Err)
enddo
endif
enddo
startflag=.true. ! Avoid re-execution if this routine is
! mistakenly called more than once at termination
endif
return
end
*endfile member=acf_AngCorFn