C LAST UPDATE:
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C***************************
C**       LFN_PLFSUB      **
C***************************
C
C
      SUBROUTINE LFN_PLFSUB ( IPASS)
C
      INCLUDE 'implicit_none.finc'
C
C PURPOSE:     1)Improve the background by omitting pixels outside the
C                expected gaussian distribution about mean background.
C              2)Integrate intensities
C
C **** Service Routine Only ****
C
C Author: Q. Hao
C
C ARGUMENTS:
C
      INTEGER  IPASS
C      INTEGER  IPASS (R)   =1, reject pixels ( > 1.0(ip),2.0(film)*STD) 
C                               from background
C                           =2, integrate background subtracted intensities
C			    =3, profile fitting
C
C====== Parameter statements
C
      INCLUDE 'saint_limits.finc'
C
C====== GLOBALS: 
C
      INCLUDE 'lfn_saint.finc'
C
C====== LOCALS:
C
      REAL AVB
      REAL AVB2
      REAL PBCK
      REAL PEAK
      REAL SUMP
      REAL SUMP2
      REAL SUM2
      REAL STD
      REAL STDP2
      REAL RAST
      REAL G1OD
      REAL QFAC
      REAL SFAC
      CHARACTER*10 TYPSTR 
      INTEGER I
      INTEGER J
      INTEGER IMTYP
      INTEGER IFLAG
      INTEGER ND
      INTEGER N1OD
C
C     REAL AVB		!Average background
C     REAL AVB2		!Average squared background 
C     REAL PBCK		!Number of background pixels
C     REAL PEAK		!Number of peak pixels
C     REAL SUMP		!Sum of profile
C     REAL SUMP2	!Sum of squared profile 
C     REAL SUM2		!Sum of squared background subtracted profile 
C     REAL STDP2	!Squared standard deviation of peak distribution
C     REAL STD		!Standard deviation of background distribution
C     REAL RAST		!Raster size in microns
C     REAL G1OD		!Selwyn granularity at unit optical density
C     REAL QFAC		!Selwyn granularity factor
C     REAL SFAC         !Sd's rejection factor
C     CHARACTER*10 TYPSTR     !Parameter as a string 'film' or 'i2'
C     INTEGER I         !Temp/loop variable
C     INTEGER J         !Temp/loop variable
C     INTEGER IMTYP     !Code number 1=film, 2=ip
C     INTEGER ND        !Number of decimal places for print
C     INTEGER IFLAG     !Flag = -1 undefined value 
C                             =  0 default value
C                             =  1 value set explicitly
C                             =  2 value set globally  
C     INTEGER N1OD	!Pixel value for 1 OD
C
C------------------------------------------------------------------------------
      CALL LDM_GET_IMTYP (TYPSTR, IMTYP, IFLAG) 
C
C  Calculate background subtracted intensity SUBPR and average background AVB
C
      AVB = 0.0
      AVB2 = 0.0
      PBCK = 0.0
      PEAK = 0.0
      INTEN = 0.0
      SUM2 = 0.0
      SUMP = 0.0
      SUMP2 = 0.0
      STDP2 = 0.0
      DO 200 I = 1, NUM_IN_SPOT
         SUBPR(I) = PIXV(I) - IXP(I)*BACK(1) - IYP(I)*BACK(2) - BACK(3) 
         IF (SHAPE(I) .LT. 0) GO TO 200
         IF (SHAPE(I) .LE. 1) THEN
C
C====== This is a background pixel
C
             AVB = AVB + PIXV(I)
             AVB2 = AVB2 + PIXV(I)**2
             PBCK = PBCK + 1.0
         ENDIF
         IF (SHAPE(I) .GE. 10) THEN
C
C====== This is a peak pixel
C
             PEAK = PEAK + 1.0
             SUM2 = SUM2 + SUBPR(I)**2
             IF (IPASS.EQ.3) THEN
                INTEN = INTEN + SUBPR(I)*PROF_IND(I)
                SUMP = SUMP + PROF_IND(I)
                SUMP2 = SUMP2 + PROF_IND(I)**2
             ELSE
                INTEN = INTEN + SUBPR(I)
             ENDIF
         ENDIF
200   CONTINUE
      IF (IPASS.EQ.3) THEN
         IF (SUMP2.GT.0) THEN
            STDP2 = SUM2 - INTEN*INTEN/SUMP2
            INTEN = INTEN*SUMP/SUMP2
         END IF
      ELSE
         IF (PEAK.GT.0.0) STDP2 = SUM2 - INTEN*INTEN/PEAK
      ENDIF
      IF ( (IPASS .GE. 2) .AND. (NREJ .GT. 0) ) THEN 
C
C====== Correct sums
C
         DO 300 J = 1, NREJ
            I = IREJ(J)
            PBCK = PBCK - 1.0
            AVB = AVB - PIXV(I)
            AVB2 = AVB2 - PIXV(I)**2
300      CONTINUE
      ENDIF
      IF (PBCK.GT.0.0) THEN
         AVB = AVB/PBCK
         AVB2 = AVB2/PBCK
      END IF
      STD = SQRT(AVB2 - AVB**2)
      IF (IMTYP.EQ.1) THEN
C
C====== Get film parameters for Selwyn granularity calculation
C
         CALL LDM_GET_RAST (RAST, IFLAG, ND)
         CALL LDM_GET_G1OD (G1OD, IFLAG, ND)
         CALL LDM_GET_N1OD (N1OD, IFLAG)
         QFAC = N1OD*G1OD*G1OD/(2.0*RAST*RAST)
C
C====== Calculate standard deviation of intensity 
C
         SIGI = SQRT(STD*STD*PEAK + 4.0*QFAC*ABS(INTEN) +STDP2)
      ELSE
         SIGI = SQRT(STD*STD*PEAK + STDP2)
      ENDIF
      IF (IPASS .GE. 2) RETURN
C                       ======
      NREJ = 0
C
C====== Reject pixels ( > 1.0(ip),2.0(film)*STD) from background
C
      SFAC = 1.0
      IF (IMTYP.EQ.1) SFAC = 2.0
         DO 400 I = 1, NUM_IN_SPOT
            IF (SHAPE(I).LT.0) GO TO 400
            IF (SHAPE(I).LE.1) THEN
               IF (ABS(SUBPR(I)).GT.SFAC*STD) THEN
                  NREJ = NREJ + 1
                  IREJ(NREJ) = I
               ENDIF
            ENDIF
400      CONTINUE
      RETURN
      END
