*-- Author : Stephen J. Maxfield
SUBROUTINE FPOKEM(NRUN0)
**----------------------------------------------------------------------
*
PARAMETER (NBIN=200)
PARAMETER (NBINLR=40)
PARAMETER (RADEG=57.295779)
DIMENSION CVEC(6)
DIMENSION IOUT(4)
DIMENSION FOUT(5)
DIMENSION GOUT(8)
DIMENSION PDAT(5,5)
DIMENSION PAR(2), PMIN(2), PMAX(2), EPAR(2), COV(3)
DIMENSION XST(4), YST(4)
DIMENSION AVEC(8)
DATA NEMAX/20/
*--------------------------------------------------------------------
* Get number of events...
CALL SAREA('FTREC', 0)
CALL GHSTAT('HS', 1, 0, NENT, SUMW, RNEFF, XST, YST)
ITOTAN = NENT
CALL SAREA('FPOKE', 0)
* Book histograms...
CALL BVEC(100, 0, 6)
CALL STEXT(100, 4,'D-time (scaled) vs. Dist(predicted)')
* Lorentz Angle stuff...
CALL BVEC(7500, 0, 6)
CALL BHD(7501, 0, 100, -25., 25., 100, -25., 25.)
CALL STEXT(7501, 4,' Delta R-perp vs. Pred drift')
* Analyse the data.
* ------- --- ----
* Analysis of histogram results.
* Write(6,*) ' Fpkanl >> Begin peaks analysis...'
* Do peakparm analysis of the pred drift histograms...
NPEAK = 0
DO KBIN = 1, NBIN
DLO = -5. + (KBIN-1) * 0.05
DHI = DLO + 0.05
JHIS1 = 2000 + KBIN
JHIS2 = 3000 + KBIN
* Get average predicted drift distance in the slice...
CALL GHSTAT('HS', JHIS2, 0, NENT, SUMW, RNEFF, XST, YST)
DMAV = XST(3)
IF(NENT .GT. NEMAX) THEN
CALL HPEAK('HS',JHIS1, 0, NPK, PDAT)
IF (NPK .GE. 1) THEN
NPEAK = NPEAK + 1
* peak position and error on...
PPOS = PDAT(1,1)
PERR = ABS(PDAT(2,1))
PINT = ABS(PDAT(3,1))
* comment out next line for 'full width errors'
PERR = 2.0*(PERR / SQRT(PINT))
* Hence Drift distance vs. drift time:-
CVEC(1) = DMAV
CVEC(2) = PPOS
CVEC(3) = DMAV - DLO
CVEC(4) = DHI - DMAV
CVEC(5) = PERR
CVEC(6) = PERR
CALL SVEC(100, 0, CVEC)
ENDIF
ENDIF
* Now purge figures - no longer needed.
CALL PURGEF(JHIS1)
CALL PURGEF(JHIS2)
ENDDO
* Lorentz Angle stuff...
NPEAKL = 0
DO KBIN = 1, NBINLR
DLO = -5.0 + (KBIN-1)*0.25
DHI = DLO + 0.25
JHIS1 = 10000 + KBIN
JHIS2 = 11000 + KBIN
* Get average predicted drift distance in the slice...
CALL GHSTAT('HS', JHIS2, 0, NENT, SUMW, RNEFF, XST, YST)
DMAV = XST(3)
* Write(6,*) KBIN, NENT, DMAV
IF(NENT .GT. 20) THEN
CALL HPEAK('HS',JHIS1, 0, NPK, PDAT)
IF (NPK .GE. 1) THEN
NPEAKL = NPEAKL + 1
* peak position and error on...
PPOS = PDAT(1,1)
PERR = ABS(PDAT(2,1))
PINT = ABS(PDAT(3,1))
* comment out next line for 'full width errors'
PERR = 2.0*(PERR / SQRT(PINT))
* Hence Delta R vs. predicted drift...
CVEC(1) = DMAV
CVEC(2) = PPOS
CVEC(3) = DMAV - DLO
CVEC(4) = DHI - DMAV
CVEC(5) = PERR
CVEC(6) = PERR
CALL SVEC(7500, 0, CVEC)
ENDIF
ENDIF
ENDDO
* Extraction of calibration data. Not for online at moment...
* Now fit the drift stuff to two straight lines...
* Set x-ranges.
CALL SAREA('FPOKE', 0)
IF(NPEAK.GE.4) THEN
KKHIS = 100
XMI = 0.6
XMA = 3.5
PAR(1) = 0.0
PMIN(1)= -0.5
PMAX(1)= 0.5
PAR(2) = 1.0
PMIN(2)= 0.5
PMAX(2)= 1.5
CALL LKFPAR (2,PAR,PMIN,PMAX,IERR)
CALL LKFITY(KKHIS, 0, 1, XMI, XMA, 'P1', IERRP)
CALL LKFGET(IFLAG,CHIS,NPT,NPAR,PAR,EPAR,COV)
IF(IERRP .EQ. 0) THEN
SLPLUS = PAR(2)
DSPLUS = EPAR(2)
ALPLUS = PAR(1)
DAPLUS = EPAR(1)
VPLUS = 7.9/(0.192308*PAR(2))
DVPLUS = VPLUS * EPAR(2)/ PAR(2)
NPPLUS = NPT
CHPLUS = CHIS
ELSE
SLPLUS = 0.
DSPLUS = 0.
ALPLUS = 0.
DAPLUS = 0.
VPLUS = 0.
DVPLUS = 0.
NPPLUS = 0.
CHPLUS = 0.
ENDIF
PAR(1) = 0.0
PMIN(1)= -0.5
PMAX(1)= 0.5
PAR(2) = 1.0
PMIN(2)= 0.5
PMAX(2)= 1.5
CALL LKFPAR (2,PAR,PMIN,PMAX,IERR)
CALL LKFITY(KKHIS, 0, 2, -XMA, -XMI, 'P1', IERRM)
CALL LKFGET(IFLAG,CHIS,NPT,NPAR,PAR,EPAR,COV)
IF(IERRM .EQ. 0) THEN
SLMINU = PAR(2)
DSMINU = EPAR(2)
ALMINU = PAR(1)
DAMINU = EPAR(1)
VMINU = 7.9/(0.192308*PAR(2))
DVMINU = VPLUS * EPAR(2)/ PAR(2)
NPMINU = NPT
CHMINU = CHIS
ELSE
SLMINU = 0.
DSMINU = 0.
ALMINU = 0.
DAMINU = 0.
VMINU = 0.
DVMINU = 0.
NPMINU = 0.
CHMINU = 0.
ENDIF
* Compute mean of +/- sides and F0R8 value...
IF(DVPLUS .GT. 0.0) THEN
WTP = 1.0/(DVPLUS**2)
ELSE
WTP = 0.0
ENDIF
IF(DVMINU .GT. 0.0) THEN
WTM = 1.0/(DVMINU**2)
ELSE
WTM = 0.0
ENDIF
IF((WTP+WTM) .GT. 0.0) THEN
VMEAN = ( VPLUS*WTP + VMINU*WTM ) / (WTP + WTM)
DVMEAN = SQRT(1.0/(WTP+WTM))
ELSE
VMEAN = 0.0
DVMEAN = 0.0
ENDIF
IF(DAPLUS .GT. 0.0) THEN
WTP = 1.0/(DAPLUS**2)
ELSE
WTP = 0.0
ENDIF
IF(DAMINU .GT. 0.0) THEN
WTM = 1.0/(DAMINU**2)
ELSE
WTM = 0.0
ENDIF
IF((WTP+WTM) .GT. 0.0) THEN
ALMEAN = ( ABS(ALPLUS)*WTP + ABS(ALMINU)*WTM ) / (WTP + WTM)
DAMEAN = SQRT(1.0/(WTP+WTM))
ELSE
ALMEAN = 0.0
DAMEAN = 0.0
ENDIF
VF0R8 = VMEAN * 0.211475
* ( entry in f0r8 bank is 10**-4 times this)
DVF0R8 = DVMEAN * 0.211475
ELSE
SLPLUS = 0.
DSPLUS = 0.
ALPLUS = 0.
DAPLUS = 0.
VPLUS = 0.
DVPLUS = 0.
NPPLUS = 0.
CHPLUS = 0.
SLMINU = 0.
DSMINU = 0.
ALMINU = 0.
DAMINU = 0.
VMINU = 0.
DVMINU = 0.
NPMINU = 0.
CHMINU = 0.
VMEAN = 0.
DVMEAN = 0.
ALMEAN = 0.
DAMEAN = 0.
DVF0R8 = 0.
ENDIF
IF(NPEAKL.GE.4) THEN
* Fit the Lorentz Angle data...
*
* Set x-ranges...
XMI = -3.5
XMA = 3.5
PAR(1) = 0.0
PMIN(1)= -2.0
PMAX(1)= 2.0
PAR(2) = -1.0
PMIN(2)= -1.5
PMAX(2)= -0.2
CALL LKFPAR (2,PAR,PMIN,PMAX,IERR)
CALL LKFITY(7500, 0, 1, XMI, XMA, 'P1', IERRL)
CALL LKFGET(IFLAG,CHIS,NPT,NPAR,PAR,EPAR,COV)
IF(IERRL .EQ. 0) THEN
TLORR = PAR(2)
ERSLP = EPAR(2)
BLOR = PAR(1)
DBLOR = EPAR(1)
NPLOR = NPT
CHLOR = CHIS
ANGLOR = ATAN(TLORR)
DANG = ANGLOR*RADEG
DANGL = ERSLP / (1. + TLORR**2)
DDANG = DANGL*RADEG
ELSE
TLORR = 0.0
ERSLP = 0.0
BLOR = 0.0
DBLOR = 0.0
NPLOR = 0.0
CHLOR = 0.0
ANGLOR = 0.0
DANG = 0.0
DANGL = 0.0
DDANG = 0.0
ENDIF
* Set some attributes...
* CALL SATTR('VEC',7500,0,'full curv')
ELSE
TLORR = 0.0
ERSLP = 0.0
BLOR = 0.0
DBLOR = 0.0
NPLOR = 0.0
CHLOR = 0.0
ANGLOR = 0.0
DANG = 0.0
DANGL = 0.0
DDANG = 0.0
ENDIF
* Output to history n-tuple.
CALL SAREA('FPOKER',0)
* Fill Ntuples...
AVEC(1) = FLOAT(NRUN0)
AVEC(2) = FLOAT(ITOTAN)
AVEC(3) = VMEAN
AVEC(4) = DVMEAN
AVEC(5) = VF0R8
AVEC(6) = DVF0R8
AVEC(7) = ALMEAN
AVEC(8) = DAMEAN
CALL SVEC(1, 0, AVEC)
AVEC(1) = VPLUS
AVEC(2) = DVPLUS
AVEC(3) = SLPLUS
AVEC(4) = DSPLUS
AVEC(5) = ALPLUS
AVEC(6) = DAPLUS
AVEC(7) = FLOAT(NPPLUS)
AVEC(8) = CHPLUS
CALL SVEC(2, 0, AVEC)
AVEC(1) = VMINU
AVEC(2) = DVMINU
AVEC(3) = SLMINU
AVEC(4) = DSMINU
AVEC(5) = ALMINU
AVEC(6) = DAMINU
AVEC(7) = FLOAT(NPMINU)
AVEC(8) = CHMINU
CALL SVEC(3, 0, AVEC)
AVEC(1) = TLORR
AVEC(2) = ERSLP
AVEC(3) = BLOR
AVEC(4) = DBLOR
AVEC(5) = DANG
AVEC(6) = DANGL
AVEC(7) = FLOAT(NPLOR)
AVEC(8) = CHLOR
CALL SVEC(4, 0, AVEC)
*-------------------------------------------------------------
999 RETURN
END
*
Function Fwiebl(x)
Common/PawPar/ P(4)
* Normalised Wiebull Function
* P(1) = area
* P(2) = x-position of the maximum [1050.0]
* P(3) = decay parameter [2.5]
* P(4) = x-threshold [1000.0]
* The values in [square] brackets are for typical Planar back-edge DOS.
* JVM 24/4/95
Fwiebl = 0.0
z = x-P(4)
z0 = P(2)-P(4)
if (z.gt.0.) then
if (z0.ge.0.) then
if (P(3).ge.1.) then
v = (P(3)-1.)/(P(3)* z0**P(3))
y = v*P(1)*P(3)
y = y* z**(P(3)-1.)
Fwiebl = y*exp( -v* z**P(3) )
endif
endif
endif
RETURN
END