;function returnEvents ;purpose: ; This function takes in an Naxis1xNaxis2 image, Im1, and returns ; all of the sets of contiguous pixels. ; ; For each set of contiugous pixels, it returns the second moments, ; ellipticity, length (as calculated by sqrt((MaxX-MinX)^2-(MaxY-MinY)^2), ; and a number of other paramters ;INPUTS: ; im1 An Naxis1xNaxis2 array with data (i.e. fits array) ; Naxis1 Number of Columns ; Naxis2 Number of Rows ; Threshold Value above which pixel is considered part of event ; Pitch The width/legnth of pixel in microns ; Thickness The thickness of CCD in microns ; ;RETURN VALUES: An Array of (11,Number of Events) quantities ; will be returned. To extract the values, use these ; indices. (t=thickness, L =track length) ; ; Lambda1Arr=ReturnArr(0,*) Max Second Moment ; Lambda2Arr=ReturnArr(1,*) Min Second Moment ; EArr=ReturnArr(2,*) Ellipticity ; FluxArr=ReturnArr(3,*) Total Counts ; TrackLength2DArr=ReturnArr(4,*) 2D Length in pixels ; TrackLength3DArr=ReturnArr(5,*) 3D Length in microns ; PerpCountsArr=ReturnArr(6,*) PerpCounts=Flux*t/L ; HitSizeArr=ReturnArr(7,*) Number of Pixels affected in Hit ; HitWidthArr=ReturnArr(8,*) Width = HitSize/L ; XExtentArr=ReturnArr(9,*) Max X value - Min X Value ; YExtentArr=ReturnArr(10,*) Max Y value - Min Y Value function returnevents,Im1,Naxis1,Naxis2,MinThreshold,ThresholdSlope,$ Sigma,Pitch,Thickness,TvFlagOriginal,ClassFlag TvFlag=TvFlagOriginal ;TVFlag = 1 XWindow =2 Postscript ImV=Im1 ;Array for viewing. Im will be erased as events are flagged ;The indices of the Array NoPixels=N_elements(Im1) PixInd=lindgen(NoPixels) PixX=(PixInd mod Naxis1) PixY=(PixInd/Naxis1) ;///////////////////////////////////////////////////////////////////// ;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 EventColor=[0.] ;Color of Event to Plot NumQuantities=13 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 Pix=0.,NoPixels-1 do begin ;Dark Images Have significant variation in dark current from top to bottom Threshold=MinThreshold+ThresholdSlope*float(CurY)/Naxis2 if Im1(Pix) gt Sigma*Threshold then begin ;Arrays to Hold Pixels ContX=[PixX(Pix)] ContY=[PixY(Pix)] ContInd=[PixInd(Pix)] ContCou=[Im1(Pix)] PossibleDelta=-1 PossibleWorm=-1 Im1(Pix)=0 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 eq 0 then begin LeftEdge=1 RightEdge=0 endif else if CurX eq Naxis1-1 then begin RightEdge=1 LeftEdge=0 endif else begin RightEdge=0 LeftEdge=0 endelse if CurY eq 0 then begin LowerEdge=1 UpperEdge=0 endif else if CurY eq Naxis2-1 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=Im1(CurX+1,CurY) if SRight gt Threshold Then begin ContCou=[ContCou,SRight] ContInd=[ContInd,CurInd+1] ContX=[ContX,CurX+1] ContY=[ContY,CurY] NumPixelsTotal=NumPixelsTotal+1 Im1(CurX+1,CurY)=0 endif Endif If (RightEdge ne 1) and (UpperEdge ne 1) then begin ;UP,RIGHT URight=Im1(CurX+1,CurY+1) if URight gt Threshold then begin ContCou=[ContCou,URight] ContInd=[ContInd,CurInd+Naxis1+1] ContX=[ContX,CurX+1] ContY=[ContY,CurY+1] NumPixelsTotal=NumPixelsTotal+1 Im1(CurX+1,CurY+1)=0 endif Endif If (UpperEdge ne 1) then begin ;UP SUp=Im1(CurX,CurY+1) if SUp gt Threshold then begin ContCou=[ContCou,SUp] ContInd=[ContInd,CurInd+Naxis1] ContX=[ContX,CurX] ContY=[ContY,CurY+1] NumPixelsTotal=NumPixelsTotal+1 Im1(CurX,CurY+1)=0 endif Endif If (LeftEdge ne 1) and (UpperEdge ne 1) then begin ;UP,LEFT ULeft=Im1(CurX-1,CurY+1) if ULeft gt Threshold then begin ContCou=[ContCou,ULeft] ContInd=[ContInd,CurInd+Naxis1-1] ContX=[ContX,CurX-1] ContY=[ContY,CurY+1] NumPixelsTotal=NumPixelsTotal+1 Im1(CurX-1,CurY+1)=0 endif Endif If (LeftEdge ne 1) then begin ;LEFT SLeft=Im1(CurX-1,CurY) if SLeft gt Threshold then begin ContCou=[ContCou,SLeft] ContInd=[ContInd,CurInd-1] ContX=[ContX,CurX-1] ContY=[ContY,CurY] NumPixelsTotal=NumPixelsTotal+1 Im1(CurX-1,CurY)=0 endif Endif If (LeftEdge ne 1) and (LowerEdge ne 1) then begin ;DOWN,LEFT DLeft=Im1(CurX-1,CurY-1) if DLeft gt Threshold then begin ContCou=[ContCou,DLeft] ContInd=[ContInd,CurInd-Naxis1-1] ContX=[ContX,CurX-1] ContY=[ContY,CurY-1] NumPixelsTotal=NumPixelsTotal+1 Im1(CurX-1,CurY-1)=0 endif Endif If (LowerEdge ne 1) then begin ;DOWN SDown=Im1(CurX,CurY-1) if SDown gt Threshold then begin ContCou=[ContCou,SDown] ContInd=[ContInd,CurInd-Naxis1] ContX=[ContX,CurX] ContY=[ContY,CurY-1] NumPixelsTotal=NumPixelsTotal+1 Im1(CurX,CurY-1)=0 endif Endif If (RightEdge ne 1) and (LowerEdge ne 1) then begin ;DOWN,RIGHT DRight=Im1(CurX+1,CurY-1) if DRight gt Threshold then begin ContCou=[ContCou,DRight] ContInd=[ContInd,CurInd-Naxis1+1] ContX=[ContX,CurX+1] ContY=[ContY,CurY-1] NumPixelsTotal=NumPixelsTotal+1 Im1(CurX+1,CurY-1)=0 endif Endif NumPixelsSearched=NumPixelsSearched+1 CurSubInd=CurSubInd+1 endwhile If N_elements(ContX) gt 1 then begin NoHits=NoHits+1. ;Event Has been Found Now Classify it ////////////////////////////////// Flux=total(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)] ;Find Breaks in the Event to see if it is a delta ray or worm if HitSizeArr(NoHits) gt 5 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 ClassFlag eq 2 then begin If (XExtentArr(NoHits) le 3) and (YExtentArr(NoHits) le 3) then begin NoSpots=NoSpots+1 EventColor=[EventColor,0] ;Tiny Spot print,"Classification: Spot" endif else if (EArr(noHits) lt 0.10) and $ (Lambda1Arr(noHits) gt 0.9) and (Lambda2Arr(noHits) gt 0.9) and $ (Lambda2Arr(noHits)/Lambda1Arr(noHits) gt 0.8) and $ (Lambda2Arr(noHits)/Lambda1Arr(noHits) lt 1.2) $ then begin NoFakeStars=NoFakeStars+1 EventColor=[EventColor,50.] print,"Classification: Fake Star" Tvflag=TvFlagOriginal ; 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) then begin ; NoSpots=NoSpots+1 ; print,"Classification: Fake Star" ; 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" ; endif else if ((Lambda1Arr(noHits) gt 2.2) and (Lambda2Arr(noHits) gt 2.2)) $ ; then begin ; NoFakeStars=NoFakeStars+1 ; print,"Classification: Fake Star" ; TvFlag=TvFlagOriginal ; endif ; endif else if abs(XExtentArr(NoHits)-YExtentArr(NoHits)) lt 2 and $ ; HitSizeArr(NoHits)/(XExtentArr(NoHits)*YExtentArr(NoHits)) gt .5 and $ ; (XExtentArr(NoHits) lt 5) and (YExtentArr(NoHits) lt 5) then begin ; NoSpots=NoSpots+1 ; print,"Classification: Spot" endif else if PossibleWorm(0) ne -1 then begin NoWorms=NoWorms+1 EventColor=[EventColor,100.] print,"Classification: Worm" endif else if PossibleDelta(0) ne -1 then begin NoDeltas=NoDeltas+1 EventColor=[EventColor,150.] print,"Classification: Delta" endif else if Lambda1Arr(Nohits) gt 2 and $ Lambda2Arr(NoHits)/Lambda1Arr(NoHits) lt 0.45 then begin NoMuons=NoMuons+1 EventColor=[EventColor,200.] print,"Classification: Muon" endif else begin NoUnidentified=NoUnidentified+1 EventColor=[EventColor,250.] print,"Classification: Unidentified" endelse EndIf if (TVFlag ge 0) and (TVFlag le 1) then begin erase !p.noerase=1 print,'Event Number = ',NoHits print,'Coordinates = ',CurX,' ',CurY print,'Ellipticity = ',Earr(noHits) print,'E Angle = ',PA 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) print,'TrackLength2D = ',TrackLength2DArr(NoHits) print,'Hit Width = ',HitWidthArr(NoHits) print,'PerpCounts = ',PerpCountsArr(NoHits) if TvFlag eq 1 then begin Device, Get_Screen_Size=ScreenSize window,xsize=.25*ScreenSize(0),ysize=.5*ScreenSize(1) tvscl,congrid($ (ImV(Min(AbsX):Max(AbsX),Min(AbsY):Max(AbsY))),$ 0.125*ScreenSize(0),0.5*ScreenSize(1),/center) endif else if TvFlag eq 2 then begin tvscl,congrid($ (ImV(Min(AbsX):Max(AbsX),Min(AbsY):Max(AbsY))),$ 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 HLineInd=where(contY eq 0) VLineInd=where(ContX eq 0) HLineContX=ContX(HLineInd) HLineContCou=ContCou(HLineInd) VLineContY=ContY(VLineInd) VLineContCou=ContCou(VLineInd) VSort=sort(VLineContY) VLineContY=VLineContY(VSort) VLineContCou=VLineContCou(VSort) plot,VLineContY,VLineContCou-Threshold,$ psym=2,title='Vertical Profile at x=0',$ position=[0.55,0.05,0.95,.5],/norm plot,HLineContX,HLineContCou-Threshold,$ psym=2,title='Horizontal Profile at y=0',$ position=[0.55,0.55,0.95,.95],/norm if TvFlag eq 1 then stop endif endif endif endfor ReturnArr=fltarr(NumQuantities,NoHits-1) ReturnArr(0,*)=Lambda1Arr(1:NoHits-1) ReturnArr(1,*)=Lambda2Arr(1:NoHits-1) ReturnArr(2,*)=EArr(1:NoHits-1) ReturnArr(3,*)=FluxArr(1:NoHits-1) ReturnArr(4,*)=TrackLength2DArr(1:NoHits-1) ReturnArr(5,*)=TrackLength3DArr(1:NoHits-1) ReturnArr(6,*)=PerpCountsArr(1:NoHits-1) ReturnArr(7,*)=HitSizeArr(1:NoHits-1) ReturnArr(8,*)=HitWidthArr(1:NoHits-1) ReturnArr(9,*)=XExtentArr(1:NoHits-1) ReturnArr(10,*)=YExtentArr(1:NoHits-1) ReturnArr(11,*)=EventColor(1:NoHits-1) ReturnArr(12,*)=Deltas(1:NoHits) ReturnArr(13,*)=Worms(1:NoHits) ReturnArr(14,0)=NoSpots ReturnArr(14,1)=NoDeltas ReturnArr(14,2)=NoWorms ReturnArr(14,3)=NoFakeStars ReturnArr(14,4)=NoStreaks ReturnArr(14,5)=NoUnidentified !p.noerase=0 return,ReturnArr end