;^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009 ;classify_fe55_xrg ; ;purpose: ; To find and classify the Fe55 events in ; ;KEYWORDS: ; UPDATEDB: int ; 0 - Don't update the mysql database ; 1 - Update the mysql database ; ;INPUTS: ; FitsFileNameOrList: str ; This can be the full path to a fits file or a .lst file that ; contains the names of many fits files. ; '****.fits' or '****.lst' ;KEYWORDS: ; UPDATEDB: int ; 0 - Don't update the mysql database ; 1 - Update the mysql database ; OBJECTTEMP: int. ; For the date of 05May08, this will pick the proper dark and also ; serve as the temperature entered into the database. The values ; are [100,110,120,.....210] ; TVPPT: int ; 0 - only plot the ramps for the full cosmic ray after all pixels have been assessed ; 1 - plot each pixel as it is evaluated ; VERBOSE: ; 0 - Don't output any text ; 1 - Output text for a global event ; 2 - Output text for every pixel ; 3 - Output all text ; SUBTRACTREFMEAN: ; 0 - Don't subtract the mean difference between successive frames ; 1 - Subtract the mean difference of the reference pixels between ; successive frames ;EXAMPLES: ; One file: ; classify_fe55_xrg,'/nfs/slac/g/ki/ki03/lances/H1RG-022/ASIC/08Jun29/FE55_Gain4_8ADC_H1RG_SIPIN_OneWindow1000_Reads_Jun30_2008_11_10_44.fits', /tvflag,/verbose ; classify_fe55_xrg,'/nfs/slac/g/ki/ki03/lances/H2RG-001/ASIC/08Dec04/Fe55/FE55_Gain4_8ADC_H2RG_SIPIN_OneWindow10_Reads_Dec04_2008_10_58_35.fits',/tvflag,/verbose ; ; Multiple Files: ; classify_fe55_xrg, '/nfs/slac/g/ki/ki03/lances/H2RG-001/ASIC/08Dec04/Fe55/Fe55_Gain4_Files.lst',/updatedb ; classify_fe55_xrg, '/nfs/slac/g/ki/ki03/lances/H2RG-32-147/ASIC/07Dec12/Fe55/Fe55Gain4.lst', /updatedb ; classify_fe55_xrg, '/nfs/slac/g/ki/ki03/lances/H2RG-32-147/ASIC/07Dec12/Fe55/Fe55Gain1_FullFrame.lst', /updatedb ; pro Classify_Fe55_xRG, FitsFileNameOrList, TvFlag=TvFlag, $ UpdateDB = UpdateDB, AvgTemp = AvgTemp, $ TvPPT = TvPPT, Verbose=Verbose, SubtractRefMean = SubtractRefMean ;Obtain the paramters and put them in the structures Common ColorStr, ColNameStr, ColNumStr, ColorTri, ColorNumTri If N_Elements(AvgTemp) eq 0 then AvgTemp = '' If N_Elements(UpdateDB) eq 0 then UpdateDB = 0 If N_Elements(TvFlag) eq 0 then TvFlag = 0 If N_Elements(TvPPT) eq 0 then TvPPT = 0 If N_Elements(Verbose) eq 0 then Verbose = 0 If N_Elements(SubtractRefMean) eq 0 then SubtractRefMean = 1 TvFlagOriginal = TvFlag ObjectText = Strtrim(AvgTemp,2)+'K' if stregex(FitsFileNameOrList, '.fits', /boolean) then begin FitsFiles = [FitsFileNameOrList] NumFitsFiles = 1 FitsFileName = FitsFileNameOrList endif else begin readcol, FitsFileNameOrList, FitsFiles, format='A' NumFitsFiles = N_Elements(FitsFiles) FitsFileName = FitsFiles(0) endelse ;Place the AvgTemp into the object text to get @KeywordStruct_xRG.pro @PlotSettings_xRG.pro ;Database Handling if UpdateDB eq 1 then begin openmysql, MySQLLun, 'FE55_HitSize', MySQLErr endif For FileNum = 0, NumFitsFiles-1 do begin FitsFileName = FitsFiles(FileNum) if NumFitsFiles gt 1 then begin if KeyStr.Date eq '07Dec12' then begin Fits_Read, FitsFileName, d, h, /header_only AvgTemp = long(strtrim(SxPar(h, 'DETTEMP'),2)) endif else begin AvgTemp = long(strmid(FitsFiles(FileNum),stregex(FitsFiles(FileNum),'K_')-3,3)) ObjectText = Strtrim(AvgTemp,2)+'K' KeyStr.ObjectText = ObjectText endelse endif else begin FitsFileName = FitsFileName(0) endelse ;Read the header first FileParam = Return_File_Params_xRG(FitsFileName) SlopeStr = Return_Slope_Params_xRG(RawHeader) TimePoints = SlopeStr.ReadTimes Pitch = KeyStr.PixPitch Thickness = KeyStr.Thickness AvgTemp = round(KeyStr.DetTemp) ThisVSUB = round(KeyStr.VSUB) ;The table name might change depending on the file being used if KeyStr.WindowMode eq 1 then begin TableName = KeyStr.DetStr+'_'+KeyStr.ElecStr+'_Window_Gain'+$ strtrim(long(KeyStr.Gain),2)+'_Fe55' endif else begin TableName = KeyStr.DetStr+'_'+KeyStr.ElecStr+'_FullFrame_Gain'+$ strtrim(long(KeyStr.Gain),2)+'_Fe55' endelse strreplace, TableName, '-','_' strreplace, TableName, '-','_' ;See if there were channels averaged if stregex(KeyStr.RawObjectKey, 'ADC', /boolean) then begin ChAvg = long(strmid(keystr.rawobjectkey,stregex(keystr.rawobjectkey,'ADC')-1,1)) endif else begin ChAvg = 1 endelse ;The indices of the Array NoPixels = float(KeyStr.NAxis1)*KeyStr.NAxis2 PixInd = lindgen(NoPixels) PixX = (PixInd mod KeyStr.NAxis1) PixY = (PixInd/KeyStr.NAxis1) ;//////////////////////////////////////////////////////////////////////////////// ;-----------------------------------------------------OPERATING PARAMETERS ;Determine borders ;FULL FRAME ooooooooooooooooooooooooooooooooooooooooooooooooooo if KeyStr.WindowMode eq 0 then begin LeftBor = KeyStr.Ref/2-1 RightBor = KeyStr.NAxis1-KeyStr.Ref/2-1 BotBor = KeyStr.Ref/2-1 TopBor = KeyStr.NAxis2-KeyStr.Ref/2-1 UpTheRamp = 1 NumRamps = 1 LastRead = KeyStr.NAxis3-1 if KeyStr.Date eq '07Dec12' then begin ReadNoise = 15.0 NeiThreshold = ReadNoise HitThreshold = 5*ReadNoise endif endif else begin ;WINDOW MODE ooooooooooooooooooooooooooooooooooooooooooooooooooo LeftBor = (KeyStr.XWindowStart(0) > (KeyStr.Ref/2-1)) RightBor = (KeyStr.XWindowStop(0) < (KeyStr.Hxrg*1024-KeyStr.Ref/2-1)) BotBor = (KeyStr.YWindowStart(0) > (KeyStr.Ref/2-1)) TopBor = (KeyStr.YWindowStop(0) < (KeyStr.Hxrg*1024-KeyStr.Ref/2-1)) ;Determine whether we are up the ramp or not if stregex(KeyStr.MCDFile, 'UpTheRamp', /boolean) then begin NumRamps = 1 LastRead = KeyStr.NAxis3-1 UpTheRamp = 1 ReadNoise = 14.0 ; ADU of CDS read noise if KeyStr.Gain eq 1 then ReadNoise = 50 HitThreshold = KeyStr.Gain*ReadNoise ; Threshold to flag an Fe55 event (max pixel) NeiThreshold = 3*ReadNoise ; Threshold to include a neighboring pixel if KeyStr.Date eq '07Dec12' then begin NumRamps = 1 LeftBor = 0 RightBor = KeyStr.NAxis1-1 BotBor = 0 TopBor = KeyStr.NAxis2-1 ReadNoise = 14.0 ; ADU of CDS read noise if KeyStr.Gain eq 1 then ReadNoise = 50 HitThreshold = KeyStr.Gain*ReadNoise ; Threshold to flag an Fe55 event (max pixel) NeiThreshold = 3*ReadNoise ; Threshold to include a neighboring pixel endif endif else if stregex(KeyStr.DetStr, 'H2RG-001', /boolean) then begin if KeyStr.Date eq '08Dec04' or KeyStr.Date eq '08Dec05' or KeyStr.Date eq '08Dec06' or $ KeyStr.Date eq '08Dec07' or KeyStr.Date eq '08Dec08' or KeyStr.Date eq '08Dec10' or $ KeyStr.Date eq '08Dec10' or KeyStr.Date eq '08Dec11' then begin NumRamps = KeyStr.NAxis3/(KeyStr.NReads*KeyStr.NGroups) LastRead = KeyStr.NReads-1 UpTheRamp = 1 ReadNoise = 8.0 ; ADU of CDS read noise LeftBor = 4 RightBor = KeyStr.NAxis1-1 BotBor = 1 TopBor = KeyStr.NAxis2-1 HitThreshold = 100 ; Threshold to flag an Fe55 event (max pixel) NeiThreshold = 3*ReadNoise ; Threshold to include a neighboring pixel endif endif else begin LastRead = KeyStr.NAxis3/2 UpTheRamp = 0 NumRamps = 1 ReadNoise = 14.0 ; ADU of CDS read noise if KeyStr.Gain eq 1 then ReadNoise = 10 HitThreshold = 5*ReadNoise ; Threshold to flag an Fe55 event (max pixel) NeiThreshold = 2*ReadNoise ; Threshold to include a neighboring pixel endelse endelse ;//////////////////////////////////////////////////////////////////////////////////////////// ;-------------------------------------------------------------------------------------------- ;The threhsold is the far end of the read noise tail NegReadNoise = -ReadNoise CosRead = 1 ; Always 1st read - 0th read ViewingBorder = 6 ;///////////////////////////////////////////////////////////////////// ;DATA CONTAINERS NumHits=0. Lambda1Arr=[0.] ;The long second moment Lambda2Arr=[0.] ;The short second moment EArr=[0.] ;Scalar Ellipticity FluxArr=[0.] ;Total Counts in a region TrackLength2DArr=[0.] ;Length of the track in Pixels 2D TrackLength3DArr=[0.] ;Lenght of the track in physical units PerpCountsArr=[0.] ;Counts scaled by length HitSizeArr=[0.] ;Total number of pixels affected HitWidthArr=[0.] ;HitWidth as judged by Length XExtentArr=[0.] ;Extent of Hit in X pixels YExtentArr=[0.] ;Extent of Hit in Y pixels NumQuantities=12 NoSpots=0 ;Number of "Spots" NoMuons=0 ;Number of "Muons" NoWorms=0 ;Number of "Worms" NoDeltas=0 ;Number of "Delta Rays" NoFakeStars=0 ;Number of Objects that look like stars NoUnidentified=0 ;Number of objects not identified as above ;////////////////////////////////////////////////////////////////////// CurX=0 CurY=0 NoHits=0. For RampNum = 0, NumRamps-1 do begin For ReadNum = 0, LastRead-1 do begin ;Make a mask of the pixels that have been contaminated or are reference Mask = lonarr(KeyStr.NAxis1, KeyStr.NAxis2) Mask(*) = 0 if UpTheRamp eq 0 then begin ;Read in Half of the raw/median image Fits_Read_DataCube, KeyStr.RawObjectPath , RawHalf , RawHeader, $ ZStart = 2*ReadNum, ZStop=2*ReadNum+1,$ Offset = Offset endif else begin Fits_Read_DataCube, KeyStr.RawObjectPath , RawHalf , RawHeader, $ ZStart = RampNum*KeyStr.NReads+ReadNum, ZStop=RampNum*KeyStr.NReads+ReadNum+1,$ Offset = Offset endelse ;Allow the subtraction of the reference pixels (e.g. if we reset preamp once per frame) If SubtractRefMean eq 1 then begin RefMean = mean(RawHalf(0:KeyStr.Ref/2-1,*,1)-long(RawHalf(0:KeyStr.Ref/2-1,*,0))) EndIf else begin RefMean = 0 EndElse ;Loop through all of the pixels in the difference frame For Pix=0., NoPixels-1 do begin ;Get coordinates and output X = PixX(Pix) Y = PixY(Pix) If X eq 0 then print, 'Row: '+strtrim(Y,2) If Mask(Pix) ne -1 then begin ;Try to Linefit HitVal = RawHalf(X,Y,CosRead)-long(RawHalf(X,Y,CosRead-1))-RefMean EndIf else begin HitVal = 0 EndElse If (HitVal gt HitThreshold) and (X gt LeftBor) and (Y gt BotBor) then begin ;Arrays to Hold Pixels NegativePix = 0 ContX=[X] ContY=[Y] ContInd=[PixInd(Pix)] ContCou=[RawHalf(X,Y,CosRead)-long(RawHalf(X,Y,CosRead-1))] PossibleDelta=-1 PossibleWorm=-1 Mask(Pix)=-1 CurSubInd=0 ;SubIndex NumPixelsSearched=0 ;NumPixels already searched NumPixelsTotal=1 ;Total NumPixels While NumPixelsSearched lt NumPixelsTotal do begin CurInd = ContInd(CurSubInd) CurX = ContX(CurSubInd) CurY = ContY(CurSubInd) ;A very slow, but exact way to find CR ;Being Careful About Edges and Corners if CurX le LeftBor then begin LeftEdge=1 RightEdge=0 endif else if CurX ge RightBor then begin RightEdge=1 LeftEdge=0 endif else begin RightEdge=0 LeftEdge=0 endelse if CurY le BotBor then begin LowerEdge=1 UpperEdge=0 endif else if CurY ge TopBor then begin UpperEdge=1 LowerEdge=0 endif else begin LowerEdge=0 UpperEdge=0 endelse ;Look Through Adjacent Pixels. Set them to zero after they've been tacked If (RightEdge ne 1) then begin ;RIGHT SRight=Long(RawHalf(CurX+1,CurY,CosRead))-RawHalf(CurX+1,CurY, CosRead-1) if (SRight gt NeiThreshold) and (Mask(CurX+1, CurY) ne -1) Then begin ContCou=[ContCou,SRight] ContInd=[ContInd,CurInd+1] ContX=[ContX,CurX+1] ContY=[ContY,CurY] NumPixelsTotal=NumPixelsTotal+1 Mask(CurX+1,CurY)=-1 endif else if (SRight lt NegReadNoise) and (Mask(CurX+1, CurY) ne -1) then begin NegativePix = 1 endif EndIf If (RightEdge ne 1) and (UpperEdge ne 1) then begin ;UP,RIGHT URight=Long(RawHalf(CurX+1,CurY+1,CosRead))-RawHalf(CurX+1,CurY+1, CosRead-1) if (URight gt NeiThreshold) and (Mask(CurX+1, CurY+1) ne -1) then begin ContCou=[ContCou,URight] ContInd=[ContInd,CurInd+Naxis1+1] ContX=[ContX,CurX+1] ContY=[ContY,CurY+1] NumPixelsTotal=NumPixelsTotal+1 Mask(CurX+1,CurY+1)=-1 endif else if (URight lt NegReadNoise) and (Mask(CurX+1, CurY+1) ne -1) then begin NegativePix = 1 endif Endif If (UpperEdge ne 1) then begin ;UP SUp=Long(RawHalf(CurX,CurY+1,CosRead))-RawHalf(CurX,CurY+1, CosRead-1) if (SUp gt NeiThreshold) and (Mask(CurX, CurY+1) ne -1) then begin ContCou=[ContCou,SUp] ContInd=[ContInd,CurInd+Naxis1] ContX=[ContX,CurX] ContY=[ContY,CurY+1] NumPixelsTotal=NumPixelsTotal+1 Mask(CurX,CurY+1)=-1 endif else if (SUp lt NegReadNoise) and (Mask(CurX, CurY+1) ne -1) then begin NegativePix = 1 endif Endif If (LeftEdge ne 1) and (UpperEdge ne 1) then begin ;UP,LEFT ULeft=Long(RawHalf(CurX-1,CurY+1,CosRead))-RawHalf(CurX-1,CurY+1, CosRead-1) if (ULeft gt NeiThreshold) and (Mask(CurX-1, CurY+1) ne -1) then begin ContCou=[ContCou,ULeft] ContInd=[ContInd,CurInd+Naxis1-1] ContX=[ContX,CurX-1] ContY=[ContY,CurY+1] NumPixelsTotal=NumPixelsTotal+1 Mask(CurX-1,CurY+1)=-1 endif else if (ULeft lt NegReadNoise) and (Mask(CurX-1, CurY+1) ne -1) then begin NegativePix = 1 endif Endif If (LeftEdge ne 1) then begin ;LEFT SLeft=Long(RawHalf(CurX-1,CurY,CosRead))-RawHalf(CurX-1,CurY, CosRead-1) if (SLeft gt NeiThreshold) and (Mask(CurX-1, CurY) ne -1) then begin ContCou=[ContCou,SLeft] ContInd=[ContInd,CurInd-1] ContX=[ContX,CurX-1] ContY=[ContY,CurY] NumPixelsTotal=NumPixelsTotal+1 Mask(CurX-1,CurY)=-1 endif else if (SLeft lt NegReadNoise) and (Mask(CurX-1, CurY) ne -1) then begin NegativePix = 1 endif Endif If (LeftEdge ne 1) and (LowerEdge ne 1) then begin ;DOWN,LEFT DLeft=Long(RawHalf(CurX-1,CurY-1,CosRead))-RawHalf(CurX-1,CurY-1, CosRead-1) if (DLeft gt NeiThreshold) and (Mask(CurX-1, CurY-1) ne -1) then begin ContCou=[ContCou,DLeft] ContInd=[ContInd,CurInd-Naxis1-1] ContX=[ContX,CurX-1] ContY=[ContY,CurY-1] NumPixelsTotal=NumPixelsTotal+1 Mask(CurX-1,CurY-1)=-1 endif else if (Dleft lt NegReadNoise) and (Mask(CurX-1, CurY-1) ne -1) then begin NegativePix = 1 endif Endif If (LowerEdge ne 1) then begin ;DOWN SDown=Long(RawHalf(CurX,CurY-1,CosRead))-RawHalf(CurX,CurY-1, CosRead-1) if (SDown gt NeiThreshold) and (Mask(CurX, CurY-1) ne -1) then begin ContCou=[ContCou,SDown] ContInd=[ContInd,CurInd-Naxis1] ContX=[ContX,CurX] ContY=[ContY,CurY-1] NumPixelsTotal=NumPixelsTotal+1 Mask(CurX,CurY-1)=-1 endif else if (SDown lt NegReadNoise) and (Mask(CurX, CurY-1) ne -1) then begin NegativePix = 1 endif Endif If (RightEdge ne 1) and (LowerEdge ne 1) then begin ;DOWN,RIGHT DRight=Long(RawHalf(CurX+1,CurY-1,CosRead))-RawHalf(CurX+1,CurY-1, CosRead-1) if (DRight gt NeiThreshold) and (Mask(CurX+1, CurY-1) ne -1) then begin ContCou=[ContCou,DRight] ContInd=[ContInd,CurInd-Naxis1+1] ContX=[ContX,CurX+1] ContY=[ContY,CurY-1] NumPixelsTotal=NumPixelsTotal+1 Mask(CurX+1,CurY-1)=-1 endif else if (DRight lt NegReadNoise) and (Mask(CurX+1, CurY-1) ne -1) then begin NegativePix = 1 endif Endif NumPixelsSearched=NumPixelsSearched+1 CurSubInd=CurSubInd+1 Endwhile If N_elements(ContX) gt 0 then begin NoHits=NoHits+1. ;Event Has been Found Now Classify it ////////////////////////////////// Flux=total(ContCou) MaxHit = Max(ContCou) ContXAvg=fix(total(ContCou*ContX)/Flux) ;First Moment in X ContYAvg=fix(total(ContCou*ContY)/Flux) ;First Moment in Y AbsX=ContX ;Maintain Absloute Coordinates AbsY=ContY ContX=ContX-ContXAvg ContY=ContY-ContYAvg FluxArr=[FluxArr,Flux] ContCountNorm=float(ContCou)/Flux ;Calculate Second Moments and find them in rotated frame I11=total(ContX^2*ContCou)/Flux I22=total(ContY^2*ContCou)/Flux I12=total(ContX*ContY*ContCou)/Flux E1=(I11-I22)/(I11+I22) E2=(2*I12)/(I11+I22) PA=.5*atan(E2,E1) XMat=[[I11,I12],[I12,I22]] Evals=Hqr(Elmhes(XMat),/Double) Lambda1Arr=[Lambda1Arr,Real_part(max(Evals))] Lambda2Arr=[Lambda2Arr,Real_part(min(Evals))] EArr=[EArr,E1^2+E2^2] ;Get Track Length and Number of Perpendicular Counts XExtentArr=[XExtentArr,Max(ContX)-Min(ContX)] YExtentArr=[YExtentArr,Max(ContY)-Min(ContY)] TrackLength3DArr=[TrackLength3DArr,$ float(sqrt((float(Pitch)*XExtentArr(NoHits))^2+$ (float(Pitch)*YExtentArr(NoHits))^2+$ float(Thickness)^2))] TrackLength2DArr=[TrackLength2DArr,sqrt((XExtentArr(NoHits)+1)^2+$ (YExtentArr(NoHits)+1)^2)] PerpCountsArr=[PerpCountsArr,FluxArr(NoHits)*Thickness/TrackLength3DArr(NoHits)] HitSizeArr=[HitSizeArr,N_elements(ContCou)] HitWidthArr=[HitWidthArr,HitSizeArr(NoHits)/TrackLength2DArr(NoHits)] NumPix = HitSizeArr(NoHits) ;=================================================================== ;Determine pixels on boundary ;Go up vertically and find left and right borders IndV=sort(ContInd) AbsXordV=AbsX(IndV) AbsYordV=AbsY(IndV) UniqueY=uniq(AbsYordV) NumUniqueY = N_Elements(UniqueY) UniqueY=[-1, UniqueY] BorderMask = BytArr(HitSizeArr(NoHits)) For i = 0, NumUniqueY-1 do begin ThisY = AbsYordV(UniqueY(i+1)) TheseX = AbsXordV(UniqueY(i)+1:UniqueY(i+1)) if i eq 0 then begin ;Bottom Boundary BorderMask(UniqueY(i)+1:UniqueY(i+1))=1 endif else if i eq NumUniqueY-1 then begin ;Top Boundary BorderMask(UniqueY(i)+1:UniqueY(i+1))=1 endif else begin ;Middle Pixels BorderMask(UniqueY(i)+1) = 1 BorderMask(UniqueY(i+1)) = 1 endelse if Verbose eq 3 then print, BorderMask(UniqueY(i)+1:UniqueY(i+1)) EndFor ;Go left to right and find top and bottom borders ContIndHor = ContInd(sort(ContInd)) IndH = sort(AbsX*Naxis1+AbsY) AbsXordH = AbsX(IndH) AbsYordH = AbsY(IndH) ContIndordH = ContInd(IndH) UniqueX = uniq(AbsXordH) NumUniqueX = N_Elements(UniqueX) UniqueX = [-1, UniqueX] BorderMask = BorderMask(sort(IndV)) For i = 0, NumUniqueX-1 do begin ThisX = AbsXordH(UniqueX(i+1)) TheseY = AbsYordH(UniqueX(i)+1:UniqueX(i+1)) if i eq 0 then begin ;Bottom Boundary BorderMask(UniqueX(i)+1:UniqueX(i+1))=1 endif else if i eq NumUniqueX-1 then begin ;Top Boundary BorderMask(UniqueX(i)+1:UniqueX(i+1))=1 endif else begin ;Middle Pixels BorderMask(UniqueX(i)+1) = 1 BorderMask(UniqueX(i+1)) = 1 endelse if Verbose eq 3 then print, BorderMask(UniqueX(i)+1:UniqueX(i+1)) EndFor ;The border mask should contain 1s for pixels on the border and 0s for interior BorderMask = BorderMask(IndV) ;==================================================================== ;Find Breaks in the Event to see if it is a delta ray or worm if HitSizeArr(NoHits) gt 4 then begin IndH=sort(ContInd) IndV=AbsX*Naxis1+AbsY XordH=ContX(IndH) YordH=ContY(IndH) XordV=ContX(sort(IndV)) YordV=ContY(sort(IndV)) XBreakH=XordH(1:HitSizeArr(NoHits)-1)-XordH(0:HitSizeArr(NoHits)-2) YBreakH=YordH(1:HitSizeArr(NoHits)-1)-YordH(0:HitSizeArr(NoHits)-2) XBreakV=XordV(1:HitSizeArr(NoHits)-1)-XordV(0:HitSizeArr(NoHits)-2) YBreakV=YordV(1:HitSizeArr(NoHits)-1)-YordV(0:HitSizeArr(NoHits)-2) PossibleDelta=where((((abs(XBreakH) gt 1) and (abs(YBreakH) eq 0)) or $ ((abs(XBreakH) eq 0) and (abs(YBreakH) gt 1))) or $ (((abs(XBreakV) gt 1) and (abs(YBreakV) eq 0)) or $ ((abs(XBreakV) eq 0) and (abs(YBreakV) gt 1)))) PossibleWorm=where((((abs(XBreakH) gt 3) and (abs(YBreakH) eq 0)) or $ ((abs(XBreakH) eq 0) and (abs(YBreakH) gt 3))) or $ (((abs(XBreakV) gt 3) and (abs(YBreakV) eq 0)) or $ ((abs(XBreakV) eq 0) and (abs(YBreakV) gt 3)))) endif ;Set TVFlag equal to 1 if you want to see the event If (HitSizeArr(NoHits) lt 5) then begin NoSpots=NoSpots+1 ;Tiny Spot print,"Classification: Spot" Type = 'Spot' endif else if (EArr(noHits) lt 0.08) and $ ;Big Fake Star (XExtentArr(NoHits) gt 8) and (YExtentArr(NoHits) gt 8) and $ (HitSizeArr(NoHits)/(XExtentArr(NoHits)*YExtentArr(NoHits))) gt .64 $ then begin NoFakeStars=NoFakeStars+1 print,"Classification: Fake Star" Tvflag=TvFlagOriginal Type = 'FakeStar' endif else if abs(XExtentArr(NoHits)-YExtentArr(NoHits)) lt 2 and $ HitSizeArr(NoHits)/(XExtentArr(NoHits)*YExtentArr(NoHits)) gt .64 and $ (XExtentArr(NoHits) ge 4) and (YExtentArr(NoHits) ge 4) and $ EArr(NoHits) lt 0.2 then begin ;Medium Fake Star NoSpots=NoSpots+1 print,"Classification: Fake Star" Type='FakeStar' endif else if((EArr(noHits) lt 0.03) and $ (Lambda2Arr(NoHits)/Lambda1Arr(NoHits) gt 0.8) $ and PossibleDelta(0) eq -1) then begin if ((Lambda1Arr(noHits) lt 2.2) and (Lambda2Arr(noHits) lt 2.2)) then begin NoSpots=NoSpots+1 print,"Classification: Spot" Type='Spot' endif else if ((Lambda1Arr(noHits) gt 2.2) and (Lambda2Arr(noHits) gt 2.2)) $ and EArr(NoHits) lt 0.2 then begin NoFakeStars=NoFakeStars+1 print,"Classification: Fake Star" TvFlag=TvFlagOriginal Type='FakeStar' endif endif else if abs(XExtentArr(NoHits)-YExtentArr(NoHits)) lt 2 and $ HitSizeArr(NoHits)/(XExtentArr(NoHits)*YExtentArr(NoHits)) gt .5 and $ (XExtentArr(NoHits) lt 4) and (YExtentArr(NoHits) lt 4) then begin NoSpots=NoSpots+1 print,"Classification: Spot" Type='Spot' endif else if PossibleWorm(0) ne -1 then begin NoWorms=NoWorms+1 print,"Classification: Worm" Type='Worm' endif else if PossibleDelta(0) ne -1 then begin NoDeltas=NoDeltas+1 print,"Classification: Delta" Type='Delta' endif else if Lambda1Arr(Nohits) gt 2 and $ Lambda2Arr(NoHits)/Lambda1Arr(NoHits) lt 0.45 then begin NoMuons=NoMuons+1 print,"Classification: Muon" Type='Muon' endif else begin NoUnidentified=NoUnidentified+1 print,"Classification: Unidentified" Type='UnIdent' endelse ;************************DISPLAY TO XWIN******************************* if TVFlag ne 0 then begin erase !p.noerase=0 if Verbose eq 1 or Verbose eq 3 then begin print,'Event Number = ',NoHits print,'Coordinates = ',CurX,' ',CurY print,'Max = ',MaxHit print,'Total = ',Flux print,'Ellipticity = ',Earr(noHits) print,'XExtent = ',XExtentArr(NoHits) print,'YExtent = ',YExtentArr(NoHits) print,'Lambda1 = ',Lambda1Arr(NoHits) print,'Lambda2 = ',Lambda2Arr(NoHits) print,'Lambda Ratio = ',Lambda2Arr(NoHits)/Lambda1Arr(NoHits) print,'Hit Size = ',HitSizeArr(NoHits) endif ThisViewingBorderX = (ViewingBorder) > XExtentArr(NoHits) XLeft = (ContXAvg-ThisViewingBorderX/2) > $ LeftBor XRight = (ContXAvg+ThisViewingBorderX/2) < $ RightBor ThisViewingBorderY = (ViewingBorder) > YExtentArr(NoHits) YBottom = (ContYAvg-BotBor-ThisViewingBorderY/2) > $ BotBor YTop = (ContYAvg+ThisViewingBorderY/2) < $ TopBor TvIm = float(RawHalf(XLeft:XRight, YBottom:YTop, CosRead))-$ RawHalf(XLeft:XRight, YBottom:YTop, CosRead-1) if TvFlag eq 1 then begin window,xsize=512,ysize=512 plot, [0],[0], /nodata, background=fsc_color('white') tvimage,congrid(bytscl(TvIm), 400,400,/center), $ background=fsc_color('white') !p.noerase=1 endif else if tvflag eq 2 then begin tvimage,congrid(bytscl(TvIm),512,512,/center),$ 0.2,0.4,$ xsize=2.5,ysize=2.5,/inches !p.charsize=0.6 xyouts,0.05,0.94,'Ellipticity = '+$ strtrim(Earr(noHits),2),/norm xyouts,0.05,0.90,'E Angle = '+$ strtrim(PA,2),/norm xyouts,0.05,0.86,'XExtent = '+$ strtrim(XExtentArr(NoHits),2),/norm xyouts,0.05,0.82,'YExtent = '+$ strtrim(YExtentArr(NoHits),2),/norm xyouts,0.05,0.78,'Lambda1 = '+$ strtrim(Lambda1Arr(NoHits),2),/norm xyouts,0.05,0.74,'Lambda2 = '+$ strtrim(Lambda2Arr(NoHits),2),/norm xyouts,0.05,0.70,'Lambda Ratio = '+$ strtrim(Lambda2Arr(NoHits)/Lambda1Arr(NoHits),2),/norm xyouts,0.05,0.66,'Hit Size = '+$ strtrim(HitSizeArr(NoHits),2),/norm xyouts,0.05,0.62,'TrackLength2D = '+$ strtrim(TrackLength2DArr(NoHits),2),/norm xyouts,0.05,0.58,'Hit Width = '+$ strtrim(HitWidthArr(NoHits),2),/norm xyouts,0.05,0.54,'PerpCounts = '+$ strtrim(PerpCountsArr(NoHits),2),/norm endif if TvFlag eq 1 then stop endif ;**************************************************************** ;==================================================================== ;Get the center of the event to use in the database entry name XCen = (Max(AbsX)+Min(AbsX))/2 YCen = (Max(AbsY)+Min(AbsY))/2 ;Now to try to sort the other pixels in order of intensity Incr = reverse(sort(ContCou)) ContCouSort = ContCou(Incr) AbsXSort = AbsX(Incr) AbsYSort = AbsY(Incr) SortedCounts = intarr(9) SortedX = intarr(9) SortedY = intarr(9) FinalCount = N_Elements(Incr) < 9 for Count=0, FinalCount-1 do begin SortedCounts(Count) = ContCouSort(Count) SortedX(Count) = AbsXSort(Count) SortedY(Count) = AbsYSort(Count) endfor ;If there is only one or two pixels, gather the outer pixel values ;to estimate IPC X = SortedX(0) Y = SortedY(0) If (Y gt (BotBor+1)) and (Y lt (TopBor-1)) and $ (X gt (LeftBor+1)) and (X lt (RightBor-1)) then begin ;Max pixel value is not on edge SortedCounts(0)=Long(RawHalf(X+0,Y+0, CosRead))-RawHalf(X+0,Y+0, CosRead-1) SortedCounts(1)=Long(RawHalf(X+0,Y+1, CosRead))-RawHalf(X+0,Y+1, CosRead-1) SortedCounts(2)=Long(RawHalf(X+0,Y-1, CosRead))-RawHalf(X+0,Y-1, CosRead-1) SortedCounts(3)=Long(RawHalf(X+1,Y+0, CosRead))-RawHalf(X+1,Y+0, CosRead-1) SortedCounts(4)=Long(RawHalf(X-1,Y+0, CosRead))-RawHalf(X-1,Y+0, CosRead-1) SortedCounts(5)=Long(RawHalf(X-1,Y-1, CosRead))-RawHalf(X-1,Y-1, CosRead-1) SortedCounts(6)=Long(RawHalf(X-1,Y+1, CosRead))-RawHalf(X-1,Y+1, CosRead-1) SortedCounts(7)=Long(RawHalf(X+1,Y-1, CosRead))-RawHalf(X+1,Y-1, CosRead-1) SortedCounts(8)=Long(RawHalf(X+1,Y+1, CosRead))-RawHalf(X+1,Y+1, CosRead-1) Edge = 0 endif else begin Edge = 1 endelse if (where(SortedCounts < 0))[0] ne -1 then NegativeBorder = 1 else NegativeBorder=0 ;ooooooooooooooooooooooooooooooooooooooooooooDATABASE If UpdateDB eq 1 and Edge eq 0 then begin GloCmd = 'INSERT INTO ' + TableName+' VALUES('+$ strtrim(RampNum*KeyStr.NReads+ReadNum,2)+','+$ strtrim(XCen,2)+','+$ strtrim(YCen,2)+','+$ strtrim(ChAvg,2)+','+$ "'"+strtrim(KeyStr.RawObjectKey,2)+"',"+$ strtrim(long(AvgTemp),2)+','+$ strtrim(long(ThisVSUB),2)+','+$ strtrim(long(NumPix),2)+','+$ strtrim(long(Flux),2)+','+$ strtrim(long(MaxHit),2)+','+$ strtrim(SortedCounts(0),2)+','+$ strtrim(SortedCounts(1),2)+','+$ strtrim(SortedCounts(2),2)+','+$ strtrim(SortedCounts(3),2)+','+$ strtrim(SortedCounts(4),2)+','+$ strtrim(SortedCounts(5),2)+','+$ strtrim(SortedCounts(6),2)+','+$ strtrim(SortedCounts(7),2)+','+$ strtrim(SortedCounts(8),2)+','+$ strtrim(NegativePix,2)+','+$ strtrim(NegativeBorder,2)+$ ");" mysqlcmd, MySQLLun, GloCmd, Answer, NLines EndIf ;oooooooooooooooooooooooooooooooooooooooooooooooooooooo endif endif endfor delvarx, RawHalf endfor endfor endfor !p.noerase=0 end