      SUBROUTINE LECSCORRECT
C
C----------------------------------------------------------------------
C
C.IDENTIFICATION: Subroutine LECSCORRECT
C.LIBRARY:        LECSLIB
C.AUTHOR:         L.Chiappetti - IFCTR Milano
C.VERSIONS:       0.0 - 24 Dec 98 - based on MECSCORRECT 0.5 style
C.PURPOSE:        perform corrections for LECS events
C.METHOD:         all LECS corrections + selection
C
C.SYNTAX:         called only by CORRECT
C.PARAMETERS:     none
C.RESTRICTIONS:
C.NOTES:          
C.FILES:
C.REFERENCES
C
C----------------------------------------------------------------------
C
C     positional correction (x',y') = f(x,y,e)
C     x'y' replace xy in place
C     dependency on energy e currently dummy
C
      INTEGER  IX,IY,IDUM,IE,K,IV,IB,INEW
      REAL     RR,RAN1,CX,DX,CY,DY,CE,CORFACT,EXTRPD
      DOUBLE PRECISION E,PROBENERGY,MEANPICHAN,MEANPI,REFPICHAN
      LOGICAL  BIT_GET
      INCLUDE  'hcommon.inc'
      INCLUDE  'bincommon.inc'
      INCLUDE  'accumcommon.inc'
      INCLUDE  'mecommon.inc'
      INCLUDE  'lecommon.inc'
C
C      extract from event
       IX=ACCUMCOMMON_J(ACCUMCOMMON_CORRINDX(1))
       IY=ACCUMCOMMON_J(ACCUMCOMMON_CORRINDX(2))
       IE=ACCUMCOMMON_J(ACCUMCOMMON_CORRINDX(3))
       IV=ACCUMCOMMON_J(ACCUMCOMMON_CORRINDX(4))
       IB=ACCUMCOMMON_J(ACCUMCOMMON_CORRINDX(5))
C
C      mandatory corrections (subtract 3 from X,Y,PHA,VETO)
       IX=MAX(IX-3,0)
       IY=MAX(IY-3,0)
       IE=MAX(IE-3,0)
       IV=MAX(IV-3,0)
       IF(.NOT.ACCUMCOMMON_CORRECT)GOTO 999
C
C      prepare gain vs time correction factor
C
       CORFACT=GAIN(1+IOF4)
       IF(GAIN_N.EQ.1)THEN
C         case of no correction or of a single gain bin
          CONTINUE
       ELSE
C         verify if we are still in the current gain window
          IF(ACCUMCOMMON_T.GT.T_NEXT.AND.GAIN_PTR.NE.GAIN_N)THEN
 10         CONTINUE
C           advance gain window pointer
            GAIN_PTR=GAIN_PTR+1
*           write(*,*)'advance ',GAIN_PTR,ACCUMCOMMON_T,T_NEXT
            IF(GAIN_PTR.LE.GAIN_N)THEN
              T_HERE=T_NEXT
              G_HERE=G_NEXT
              T_NEXT=GAIN_T(GAIN_PTR+IOF5)
              G_NEXT=GAIN(GAIN_PTR+IOF4)
C             iterate in case some window is skipped
              IF(ACCUMCOMMON_T.GT.T_NEXT)GOTO 10
            ENDIF
          ENDIF
C         in any case interpolate
          CORFACT=EXTRPD(G_HERE,G_NEXT,T_HERE,T_NEXT,ACCUMCOMMON_T)
*         write(*,*)' corfact ',ACCUMCOMMON_T,corfact
       ENDIF
C
C      perform correction here
C
C      gain correction : edet = f(eraw,xraw,yraw)
C      put gain back to centre of detector
C
         RR=RAN1(IDUM)
         CE=IE+RR-0.5
*        IE=NINT(CE/MAPE(IX+1,IY+1))
         K=(IY)*256+IX+1
*        IE=NINT(CE/MAPE(K+IOF1))
C####    this will be final
         IE=NINT(CE/MAPE(K+IOF1)*CORFACT)
C####    this is provisional
C        to have LECSGAINACCUM in SECRET mode bring back the gain to 1.0 !
*        IE=NINT(CE*CORFACT)
C-->     events rescaled out of range are assigned MAX coordinate
C-->     NB this is different from IFCAI which uses EITHER zero OR max 
C-->     and counted too ? so far not 
         IF(IE.LT.0.OR.IE.GT.1023)THEN
*           NREJE=NREJE+1
            IE=1023
         ENDIF
C
C      positional correction (xdet,ydet) = f (xraw,yraw,)
C####  no energy dependency for LECS
*      ekev=EBOUND(IE+1)
C
       DX=IX-AFOV
       DY=IY-AFOV
       IF(DX*DX+DY*DY.LE.10000) THEN
C####    LECS linearization defined only inside FOV (100 pix radius)
         RR=RAN1(IDUM)
         CX=IX+RR-0.5-AA1
         RR=RAN1(IDUM)
         CY=IY+RR-0.5-AA2
         DX=AS(0)+ AS(1)*CX+ AS(2)*CY+ AS(3)*CX*CY+ AS(4)*CX*CX+ 
     .      AS(5)*CY*CY+ AS(6)*CX*CX*CY+ AS(7)*CX*CY*CY+
     .      AS(8)*CX*CX*CX+ AS(9)*CY*CY*CY
         DY=BS(0)+ BS(1)*CX+ BS(2)*CY+ BS(3)*CX*CY+ BS(4)*CX*CX+ 
     .      BS(5)*CY*CY+ BS(6)*CX*CX*CY+ BS(7)*CX*CY*CY+
     .      BS(8)*CX*CX*CX+ BS(9)*CY*CY*CY
       ELSE
C####    plain geometric scale outside FOV
         DX=(IX+RR-0.5-AFOV)*OSIZE
         DY=(IY+RR-0.5-BFOV)*OSIZE 
       ENDIF
C
C      mm to nexpixels
C####  contrary to SAXDAS do NOT rotate by misalignments
C      in order to preserve compatibility with XAS MECS !
       IX=NINT((DX+HALF)/SIZE)
       IY=NINT((DY+HALF)/SIZE)
C-->   events rescaled out of range are assigned zero coordinate 
C-->   and counted too ? so far not 
C      attention !  range is 0 to NPIX-1 not 1 to NPIX !
       IF(IX.LT.0.OR.IX.Ge.NPIX)THEN
*           NREJX=NREJX+1
            IX=0
       ENDIF
       IF(IY.LT.0.OR.IY.Ge.NPIX)THEN
*           NREJY=NREJY+1
            IY=0
       ENDIF
C
C      if event is outside of spatial region flag it for rejection
       IF(SPATIAL_REGION)THEN
       IF(.NOT.BIT_GET(IX+1,IY+1))IX=-1
       ENDIF
C
C####  LECS specific corrections
C
C       BL-Veto filtering
C       (actually VETO filtering occurs outside via normal range selection)
C       if event BL is outside energy-dependent limits flag it for rejection
C       (and do not process such event any further)
C       lookup in which BL box the event energy falls
        IF(IB.LT.BLLO(IE+1).OR.IB.GT.BLHI(IE+1))THEN
          IB=-1  
          GOTO 999
        ENDIF
C
C       PI-PIC correction (BL dependent correction of energy)
        IF(PICENABLED)THEN
C         random energy corresponding to gain
          E =ProbEnergy(MAPY(1+IOF3),IE)
C         correction only above 2.9 keV
          IF(E.GE.2.9)THEN
            MeanPI = MeanPIChan(E,IB,RefPIChan)
C           this is corrected energy
            INEW = RefPIChan + IE - MeanPI
*           write(*,*)ie,e,inew
C           do not correct if it goes out of range
            IF(INEW.GE.0.AND.INEW.LE.1023)IE=INEW
         ENDIF
        ENDIF
C
C      reinsert in event
 999   CONTINUE
       ACCUMCOMMON_J(ACCUMCOMMON_CORRINDX(1))=IX
       ACCUMCOMMON_J(ACCUMCOMMON_CORRINDX(2))=IY
       ACCUMCOMMON_J(ACCUMCOMMON_CORRINDX(3))=IE
       ACCUMCOMMON_J(ACCUMCOMMON_CORRINDX(4))=IV
       ACCUMCOMMON_J(ACCUMCOMMON_CORRINDX(5))=IB

       RETURN
       END
C
       FUNCTION PROBENERGY(MAP,IE)
C      probability map is actually dynamically allocated
       INTEGER MAP(1024,1024),IDUM,IE,N
       REAL RR,RAN1
       INCLUDE  'lecommon.inc'
       DOUBLE PRECISION PROBENERGY
       RR=RAN1(IDUM)*LESCALE
       DO N=1,1023
       IF(MAP(IE+1,N).GE.RR)THEN
         PROBENERGY=LEEMIN+(N+0.5)*EBINSIZE
         GOTO 10
       ENDIF
       ENDDO
       write(*,*)' error on probability for channel',IE
 10    CONTINUE
       RETURN
       END
C
C      aux function for PI-BL correction
C
       FUNCTION  MeanPIChan(E,IB,RefPIChan)
       INCLUDE  'lecommon.inc'
       INTEGER IB
       DOUBLE PRECISION MeanPIChan,E,GPPB,BRKPI,ASYMPI,LINPI,CONST
       DOUBLE PRECISION REFPICHAN
       IF(E.LE.xeledge)THEN
          GPPB=gradloa0+gradloa1*E
       ELSE
          GPPB=gradhia0+gradhia1*E
       ENDIF
       BRKPI=pibrka0+pibrka1*E
       ASYMPI=piasyma0+piasyma1*E
       REFPICHAN=pirefa0+pirefa1*E
       LINPI=BRKPI-brkbl*GPPB
       CONST=(brkbl-offsetbl)*(ASYMPI-BRKPI)
       IF(DBLE(IB).LT.brkbl)THEN
          MeanPIChan=LINPI+GPPB*IB
       ELSE
          MeanPIChan=ASYMPI-CONST/(IB-offsetbl)
       ENDIF
       RETURN
       END
