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