*-- Author : John V. Morris
SUBROUTINE FVFITP(NC,IDF)
************************************************************************
* *
* Fit the Class I plot IDF to determine Drift Velocity and resolution *
* JVM 24/9/91 *
* *
* Scale errors on fitted parameters by SQRT(CHISQ/NDF) for bad Chisq!? *
* JVM 29/7/92 not very sanitary, but takes care of some systematics *
************************************************************************
*KEEP,FCOMF.
PARAMETER( NP = 8 ) ! number of fit parameters
COMMON/FCOMF/PAR(NP),STAG
*KEEP,FVOUTP.
PARAMETER( MXC=40 ) ! maximum stored results
COMMON/FVOUTP/VEL(MXC),EVEL(MXC),RMIC(MXC),EMIC(MXC),IDVF(MXC),
& SEGS(MXC),ESEG(MXC),ASYM(MXC),EASYM(MXC),
& SEGN(MXC),ESEGN(MXC),SEGF(MXC),ESEGF(MXC)
*KEEP,FHLUN.
COMMON/FHLUN/ LUNH, LUNS, LMES, IL4L5
*KEND.
CHARACTER*30 TITL
EXTERNAL FM2SOV
DATA IC/ 2/
DIMENSION P(NP),ST0(NP),PMI0(NP),PMA0(NP),SIG(NP),COV(NP*(NP+1)/2)
DIMENSION ST(NP) ,PMI(NP) ,PMA(NP)
DATA P/0.0,0.0,0.0,300.,300.,40.0,4.0,0.0/
DATA ST0/10.,1.0,1.0, 30., 30., 4., 0.4,0.5/
DATA PMI0/0.0,-99.,-99.,0., 0.,10.0,1.0,-10./
DATA PMA0/999.,99.,99.,9999.,9999.,80.0,20.,10./
* is there anything to fit ??
IF( NC.GT.MXC )THEN
* WRITE(LUNH,1002)IDF,MXC
1002 FORMAT(1X,'*** FVFITP ID ',I10,' Max calls exceeded ',I10)
RETURN
ENDIF
VEL(NC) = 0.0
EVEL(NC)= 0.0
RMIC(NC)= 0.0
EMIC(NC)= 0.0
IDVF(NC)= IDF
SEGS(NC)= 0.0
ESEG(NC)= 0.0
SEGN(NC)= 0.0
ESEGN(NC)= 0.0
SEGF(NC)= 0.0
ESEGF(NC)= 0.0
IF( HSUM(IDF).LT.300. )THEN
* WRITE(LUNH,1000)IDF,HSUM(IDF)
1000 FORMAT(1X,'*** FVFITP ID ',I10,' Low contents',F10.1)
RETURN
ENDIF
P(6) = 25.0 ! start at standard Full Field Drift V
P(4) = HMAX(IDF)
P(5) = P(4)
P(1) = HMIN(IDF)
PMA0(1) = HMAX(IDF)
PMA0(4) = 2.0*HMAX(IDF)
PMA0(5) = PMA0(4)
CALL UCOPY( P,PAR,NP )
CALL UCOPY( ST0,ST,NP )
CALL UCOPY( PMI0,PMI,NP )
CALL UCOPY( PMA0,PMA,NP )
CALL HFIT(IDF,FM2SOV,NP,PAR,CHISQ,IC,SIG,COV,ST,PMI,PMA)
* compute scaling factor for errors .. ?? ..
CALL HGIVE(IDF,TITL,NBIN,XMIN,XMAX,NY,YMI,YMA,NWT,LOC) ! get bin s
SCER = SQRT( CHISQ/FLOAT(NBIN-NP) )
IF( SCER.LT.1. ) SCER = 1.
* compute resolution in microns and print results
VEL(NC) = PAR(6)
EVEL(NC) = SIG(6)*SCER
RMIC(NC) = PAR(6)*PAR(7)
EMIC(NC) = ( SIG(6)*PAR(7) )**2 + ( SIG(7)*PAR(6) )**2 +
& 2.0*COV(27)*PAR(6)*PAR(7)
IF(EMIC(NC).GT.0.0) THEN
EMIC(NC) = SQRT(EMIC(NC))*SCER
ELSE
EMIC(NC) = -1.0
ENDIF
ASYM(NC) = PAR(8)
EASYM(NC) = SIG(8)
* WRITE(LUNH,9600)PAR(6),SIG(6),PAR(7),SIG(7),RMIC(NC),EMIC(NC)
* & ,ASYM(NC),EASYM(NC)
*9600 FORMAT(/,1X,'Drift Velocity =',F10.2,' +/-',F6.2,' micr/nsec',
* & /,1X,'Resolution =',F10.2,' +/-',F6.2,' nsec ',
* & /,1X,' =',F10.2,' +/-',F6.2,' microns ',
* & /,1X,'Asymmetry =',F10.2,' +/-',F6.2,' nsec ')
*
* recompute RMIC in nsecs for compatibility with FSUMP
RMIC(NC) = PAR(7)
EMIC(NC) = SIG(7)*SCER
* compute number of fitted 4-hit track segments and error
BIN = (XMAX-XMIN)/FLOAT(NBIN)
* BIN = 3.0
CON = 2.8024956/BIN ! Sqrt(10pi)/2*Bin size
* note than this constant includes the conversion factor of
* sqrt(5)/2 which converts PAR(7) to the actual RMS of the peaks.
SEGS(NC) = (PAR(4)+PAR(5))*PAR(7)*CON ! number of 4-hit segs
DXDP4 = PAR(7)*CON
DXDP5 = DXDP4
DXDP7 = (PAR(4)+PAR(5))*CON
ESEG(NC)=(DXDP4*SIG(4))**2+(DXDP5*SIG(5))**2+(DXDP7*SIG(7))**2
& +2.0*DXDP4*DXDP5*COV(23) +2.0*DXDP4*DXDP7*COV(25)
& +2.0*DXDP5*DXDP7*COV(29)
IF(ESEG(NC).GT.0.0) THEN
ESEG(NC)= SQRT( ESEG(NC) )*SCER
ELSE
ESEG(NC)= -1.0
ENDIF
* WRITE(LUNH,9601)SEGS(NC),ESEG(NC)
9601 FORMAT(/,1X,'4-hit Segments =',F10.1,' +/-',F6.1,' ')
SEGF(NC) = PAR(4)*PAR(7)*CON ! number of 4-hit segs
DXDP4 = PAR(7)*CON
DXDP7 = PAR(4)*CON
ESEGF(NC)=(DXDP4*SIG(4))**2+(DXDP7*SIG(7))**2
& +2.0*DXDP4*DXDP7*COV(25)
IF(ESEGF(NC).GT.0.0) THEN
ESEGF(NC)= SQRT( ESEGF(NC) )*SCER
ELSE
ESEGF(NC)= -1.0
ENDIF
* WRITE(LUNH,9602)SEGF(NC),ESEGF(NC)
9602 FORMAT(/,1X,'4-hit Far Seg =',F10.1,' +/-',F6.1,' ')
SEGN(NC) = PAR(5)*PAR(7)*CON ! number of 4-hit segs
DXDP5 = PAR(7)*CON
DXDP7 = PAR(5)*CON
ESEGN(NC)=(DXDP5*SIG(5))**2+(DXDP7*SIG(7))**2
& +2.0*DXDP5*DXDP7*COV(29)
IF(ESEGN(NC).GT.0.0) THEN
ESEGN(NC)= SQRT( ESEGN(NC) )*SCER
ELSE
ESEGN(NC)= -1.0
ENDIF
* WRITE(LUNH,9603)SEGN(NC),ESEGN(NC)
9603 FORMAT(/,1X,'4-hit Near Seg =',F10.1,' +/-',F6.1,' ')
RETURN
END