;^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009 ;classify_one_pixel ; ;PURPOSE: ; To classify the signal vs. time for a pixel that has been hit by a cosmic ray ;INPUTS: ; Row: int ; The row of the pixel ; Col: int ; The column of the pixel ; RawRamp: long int array ; The raw signal values for the pixel (Col, Row) ; CosRead: int ; The read before the cosmic ray hit ; BorderPix: int ; 0 - The pixel is not on the periphery of the event ; 1 - The pixel is on the border ; Temp: int ; The temperature that the exposures were taken at (In Kelvin) ; XCen: int ; The center x pixel of the event ; YCen: int ; The center y pixel of the event ; Temp: int ; The temperature ; VSUB: float ; The voltage of the backside ;KEYWORDS: ; DarkSub: int ; 0 - Don't bother with subtracting a dark ; 1 - Subtract the median dark ; ZDepth: float ; The depth of the event (based upon charge lost) ; pro Classify_One_Pixel, Col, Row, RawRamp, CosRead, BorderPix, $ Temp, VSUB, $ XCen, YCen, CosEDep, ChargeLost, ChargeRise, TvFlag, $ TvPPT, UpdateDB, Verbose, MySQLLun=MySQLLun, DarkSub = DarkSub, $ ZDepth=ZDepth ;Get the common block Common KeyParams, KeyStr Common func_cosmic_decay, ReadTimesAC, Decay If N_Elements(DarkSub) eq 0 then DarkSub = 0 If N_Elements(ZDepth) eq 0 then ZDepth = 0 if TvFlag eq 1 and TvPPT eq 1 then begin window, 1, xpos=0,ypos=500 !p.noerase=0 endif if UpdateDB eq 1 then begin SingleTableName = KeyStr.DetStr+'_'+KeyStr.ElecStr+'_'+'Sin_CRHits' strreplace, SingleTableName, '-','_' strreplace, SingleTableName, '-','_' endif ;Get read times ReadTimes = findgen(KeyStr.NAxis3)*KeyStr.FrameTime ;Get the dark ramp for this pixel if DarkSub eq 1 then begin Fits_Read_DataCube, KeyStr.DarkName , DarkRamp , DarkHeader, $ XStart= Col , XStop = Col, $ YStart= Row , YStop = Row EndIf Else if DarkSub eq 2 then Begin LineFit_xRG, ReadTimes, RawRamp, $ m_d, b_d, sigm_d, sigb_d, Error, $ CosEDep2, CosRead2, CosType2, $ Col= Col, Row= Row, ThresholdFactor=3 if CosRead2 eq 0 then begin LineFit_xRG, ReadTimes, RawRamp, $ m_d, b_d, sigm_d, sigb_d, Error, $ CosEDep2, CosRead2, CosType2, $ Col= Col, Row= Row, ThresholdFactor=6 endif DarkRamp = m_d*ReadTimes+b_d EndIf Else Begin DarkRamp = fltarr(N_Elements(RawRamp)) EndElse ;Generate Ramp-Dark MaxDark = max((DarkRamp(1:*)-float(DarkRamp(0)))) RampDS = long(RawRamp)-DarkRamp ;Get some values CosEDep = (RampDS(CosRead)-RampDS(CosRead-1))*KeyStr.CGain MaxRamp = max(RampDS) MaxRampRead = where(RampDS eq MaxRamp) MinRamp = min(RampDS) MinRampRead = where(RampDS eq MinRamp) ;AC=After Cosmic Hit RampAC = RampDS(CosRead:*) MaxAC = Max(RampAC) MaxReadACS = where(RampAC eq MaxAC) MaxReadAC = MaxReadACS+CosRead MinAC = Min(RampAC) MinReadACS = where(RampAC eq MinAC) MinReadAC = MinReadACS+CosRead ChargeLost = (RampAC(0) - MinAC)*KeyStr.CGain ;Check to see if signal rises or decays immediately after if N_Elements(RampAC) gt 2 then begin AveAftHit = mean(RampAC(1: (5 < (N_Elements(RampAC)-1)))) endif else begin AveAftHit = RampAC(0) endelse ;Check to see if it's a decay or not if RampAC(0) gt AveAftHit then begin ;DECAY -- CHECK TIME CONSTANT=========================================== ;Get the ReadTimes of the Decay Decay = RampDS(CosRead:MinReadAC(0))-long(MinAC) DecayLen = N_Elements(Decay) ReadTimesAC = findgen(DecayLen)*KeyStr.FrameTime ChargeRise = -1 rs=findgen(2) tries=20000 ;Attempt to fit rs=amoeba(1.0e-10, scale=[2.*ChargeLost, 300.], $ p0=[1.5*ChargeLost, 300.], $ function_value=fval,$ function_name='fit_cosmic_decay', ncalls=tries,$ nmax=100000) ;Fit Failure------------------------------------------------------------ if rs(0) eq -1 then begin FitErrTot = -1 FitErrStd = -1 rs=[-1,-1] if TvFlag eq 1 and TvPPT eq 1 then begin erase plot, ReadTimesAC, Decay, color=fsc_color('black'), $ background = fsc_color('white') endif ;Fit Success------------------------------------------------------------ endif else begin FitErrTot = total(Decay-rs(0)*exp(-ReadTimesAC/rs(1))) FitErrStd = stddev(Decay-rs(0)*exp(-ReadTimesAC/rs(1))) if TvFlag eq 1 and TvPPT eq 1 then begin erase !p.noerase=0 plot, [0],[0], /nodata, background=fsc_color('white') !p.noerase=1 plot, ReadTimesAC, Decay, color=fsc_color('black'), $ background = fsc_color('white'), position = [0.05, 0.05, .475, 0.95], $ /normal oplot,ReadTimesAC, rs(0)*exp(-ReadTimesAC/rs(1)), $ color=fsc_color('black') plot, ReadTimes, RampDS, $ color=fsc_color('black'), background=fsc_color('white'),$ position = [0.55, 0.53, .98, 0.95], $ /normal plot, ReadTimes, RawRamp, yrange=[min(RawRamp), max(RawRamp)], $ color=fsc_color('black'), background=fsc_color('white'),$ position = [0.55, 0.05, .98, 0.48], $ /normal !p.noerase=0 endif endelse endif else begin ;RISE -- CHECK GAIN IN SIGNAL========================================= ChargeRise = (RampDS(MaxReadAC(0))-long(RampDS(CosRead)))*KeyStr.CGain ChargeLost = -1 rs = [-1,-1] FitErrTot = -1 FitErrStd = -1 if TvFlag eq 1 and TvPPT eq 1 then begin plot, ReadTimes, RampDS(1:*)-long(RampDS(1)), color=fsc_color('black'), $ background = fsc_color('white') endif endelse If Verbose eq 2 or Verbose eq 3 then begin print,'X = '+Strtrim(Col,2) print,'Y = '+Strtrim(Row,2) print,'CosRead = '+Strtrim(CosRead,2) print,'CosEDep = '+Strtrim(CosEDep,2) print,'Q Lost = '+Strtrim(ChargeLost, 2) print,'Q Gained = '+Strtrim(ChargeRise, 2) print,'MinReadAC = '+Strtrim(MinReadAC,2) print,'MaxReadAC = '+Strtrim(MaxReadAC,2) print,'Fit A = '+Strtrim(rs(0),2) print,'Fit Tau = '+Strtrim(rs(1),2) print,'Fit Err. Tot = '+Strtrim(FitErrTot,2) print,'Fit Err. Std = '+Strtrim(FitErrStd,2) print,'Border Pixel = '+Strtrim(long(BorderPix), 2) endif if UpdateDB eq 1 then begin CosID = KeyStr.RawObjectKey+'_'+strtrim(XCen,2)+'_'+strtrim(YCen,2) MySQLString = 'INSERT INTO ' + SingleTableName+' VALUES('+$ "'"+strtrim(CosID,2)+"',"+$ strtrim(Col,2)+','+$ strtrim(Row,2)+','+$ strtrim(long(CosRead),2)+','+$ strtrim(long(CosEDep),2)+','+$ strtrim(ChargeLost,2)+','+$ strtrim(ChargeRise,2)+','+$ strtrim(MinReadAC,2)+','+$ strtrim(MaxReadAC,2)+','+$ strtrim(rs(0),2)+','+$ strtrim(rs(1),2)+','+$ strtrim(FitErrTot,2)+','+$ strtrim(FitErrStd,2)+','+$ strtrim(long(Temp),2)+","+$ strtrim(long(VSUB),2)+","+$ strtrim(MaxDark,2)+","+$ strtrim(ZDepth, 2)+","+$ strtrim(long(BorderPix),2)+","+$ "'"+strtrim(KeyStr.Date,2)+"'"+$ ");" mysqlcmd, MySQLLun, MySQLString, Answer, NLines, debug=0 endif if TvFlag eq 1 and TvPPT eq 1 then stop end ;fit_cosmic_decay ; ;PURPOSE: ; To fit the decay of the ramp after the cosmic ray has hit ; and the traps begin to do their thing ; ;INPUTS: ; x - The read times ; y - The signal value ; ;COEFF: ; Equation: Ae^-(tau*t) ; p(0): A : The overall trap number ; p(1): Tau : the time constant ; function fit_cosmic_decay, p common func_cosmic_decay, x, y CosDiff = total(abs(y-p(0)*exp(-x/p(1)))) return, CosDiff end ;classify_cosmics ;purpose: ; To find and classify the cosmic ray events from a datacube generated ; by a HyViSi device. If UpdateDB is set, the results will be written to ; a mysql database ; ; To avoid memory overflow in the computer, only half of the datacube is ; processed at a given time ; ; Take in Half the frame at a time and shift it upward, as shown ; ; 1 00000 2 00000 3 ----- ; 00000 ----- ----- ; ----- ----- 00000 ; ----- 00000 00000 ; ;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: ; DARKSUB: int ; 0 - Don't bother with subtracting a dark ; 1 - Subtract the dark ; 2 - Fit a line to the dark current and subtract this ; 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 ; TESTROW: int ; If this keyword is set, the cosmics will not be evaluated until the loop has ; reached this row. ; MUTHRESH: int ; The threshold to use for determining the depth of cosmic muons ;EXAMPLES: ; One file: ; classify_cosmics_xrg, '/nfs/slac/g/ki/ki04/lances/H1RG-022/LEACH/05May08/Dark/dark_130K_65_min_1_04.fits', updatedb=0, tvflag=1, tvppt=1, avgtemp=120 ; Multiple Files from 05May08: ; classify_cosmics_xrg, '/nfs/slac/g/ki/ki04/lances/H1RG-022/LEACH/05May08/Dark/DarkFiles.lst',/updatedb, /darksub ; Multiple File from H2RG-001 ; classify_cosmics_xrg, '/nfs/slac/g/ki/ki04/lances/H2RG-001/ASIC/08Dec16/Dark/DarkFiles_110K_170K_20V.txt', /updatedb, darksub=2 ; pro Classify_Cosmics_xRG, FitsFileNameOrList, TvFlag=TvFlag, $ UpdateDB = UpdateDB, AvgTemp = AvgTemp, $ TvPPT = TvPPT, Verbose=Verbose, DarkSub = DarkSub, $ TestRow = TestRow, MuThresh=MuThresh ;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(DarkSub) eq 0 then DarkSub = 0 If N_Elements(TestRow) eq 0 then TestRow = 0 If N_Elements(MuThresh) eq 0 then MuThresh = 400 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 if KeyStr.DetStr eq 'H2RG-001' then begin Fits_Read, KeyStr.BadMaskDir+'BadMask.fits', BadMask BadMask(where(BadMask ne 0)) = -1 endif else begin BadMask = bytarr(KeyStr.NAxis1, KeyStr.NAxis2) endelse ;Database Handling if UpdateDB eq 1 then begin openmysql, MySQLLun, 'CRHit_Decays', MySQLErr GlobalTableName = KeyStr.DetStr+'_'+KeyStr.ElecStr+'_'+'Glo_CRHits' SingleTableName = KeyStr.DetStr+'_'+KeyStr.ElecStr+'_'+'Sin_CRHits' strreplace, GlobalTableName, '-','_' strreplace, GlobalTableName, '-','_' strreplace, SingleTableName, '-','_' strreplace, SingleTableName, '-','_' endif For FileNum = 0, NumFitsFiles-1 do begin FitsFileName = FitsFiles(FileNum) if NumFitsFiles gt 1 then begin AvgTemp = long(strmid(FitsFiles(FileNum),stregex(FitsFiles(FileNum),'K_')-3,3)) ObjectText = Strtrim(AvgTemp,2)+'K' KeyStr.ObjectText = ObjectText endif else begin FitsFileName = FitsFileName(0) endelse ;Read the header first - once to get the detector temp and then ;run the file params again to get the proper dark FileParam = Return_File_Params_xRG(FitsFileName) SlopeStr = Return_Slope_Params_xRG(RawHeader) TimePoints = SlopeStr.ReadTimes Pitch = KeyStr.PixPitch Thickness = KeyStr.Thickness if not (KeyStr.DetStr eq 'H1RG-022' and KeyStr.Date eq '05May08') then begin AvgTemp = KeyStr.DetTemp ;Round the temperature to the tens place AvgTemp = 10*round(float(AvgTemp)/10) ObjectText = Strtrim(AvgTemp,2)+'K' KeyStr.ObjectText = ObjectText endif ;The indices of the Array NoPixels=float(KeyStr.NAxis1)*KeyStr.NAxis2 NoPixelsPerThird = NoPixels/3 PixInd=lindgen(NoPixels) PixX=(PixInd mod KeyStr.NAxis1) PixY=(PixInd/KeyStr.NAxis1) ;Make a mask of the pixels that have been contaminated or are reference Mask = lonarr(KeyStr.NAxis1, KeyStr.NAxis2) Mask(*) = 0 ;Apply a bad mask if (where(BadMask eq 255))[0] ne -1 then Mask(where(BadMask eq 255)) = -1 ;The threhsold is the far end of the read noise tail Threshold = 80 if KeyStr.Date eq '05May08' then Threshold=300 ViewingBorder = 40 ;///////////////////////////////////////////////////////////////////// ;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. ;First Get the Reference Pixels RefPix = Return_RefPix_xRG(FitsFileName, 0) FirstFrThird = long(TestRow)/(KeyStr.NAxis2/3) For FrThird=FirstFrThird, 2 do begin ;Read in Half of the raw/median image Fits_Read_DataCube, KeyStr.RawObjectPath , RawHalf , RawHeader, $ XStart= 0, XStop = KeyStr.Naxis1-1, $ YStart= FrThird*KeyStr.Naxis2/4, $ YStop = (FrThird+2)*KeyStr.Naxis2/4-1, $ Offset = Offset For PixPerThirdInd=0.,NoPixelsPerThird-1 do begin Pix = PixPerThirdInd+NoPixelsPerThird*FrThird ;Get coordinates and output X = PixX(Pix) SubY = PixY(Pix)-FrThird*KeyStr.Naxis2/4 FullY = PixY(Pix) If X eq 0 then print, 'Row: '+strtrim(FullY,2) If Mask(Pix) ne -1 then begin ;Try to Linefit RawRamp = long(RawHalf(X, SubY, *)) LineFit_xRG, TimePoints, RawRamp, $ m, b, sigm, sigb, Error, $ CosEDep, CosRead, CosType, $ Col= X, Row= Y, ThresholdFactor=6 EndIf else begin CosRead = 0 EndElse ;Only classify the pixel if CosRead > 0 ; CosRead = 0 or -1 means no cosmic found If CosRead gt 0 and X gt KeyStr.Ref/2 and FullY gt KeyStr.Ref/2 and $ FullY gt TestRow then begin ;Arrays to Hold Pixels ContX=[X] ContY=[FullY] ContInd=[PixInd(Pix)] ContCou=[Long(RawRamp(CosRead))-RawRamp(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) CurSubY = ContY(CurSubInd)-FrThird*KeyStr.Naxis2/4 CurFullY = ContY(CurSubInd) ;A very slow, but exact way to find CR ;Being Careful About Edges and Corners if CurX le KeyStr.Ref/2-1 then begin LeftEdge=1 RightEdge=0 endif else if CurX ge KeyStr.Naxis1-KeyStr.Ref/2-1 then begin RightEdge=1 LeftEdge=0 endif else begin RightEdge=0 LeftEdge=0 endelse if CurFullY le KeyStr.Ref/2-1 then begin LowerEdge=1 UpperEdge=0 endif else if CurSubY ge KeyStr.Naxis2/2-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=Long(RawHalf(CurX+1,CurSubY,CosRead))-RawHalf(CurX+1,CurSubY, CosRead-1) if (SRight) gt Threshold and Mask(CurX+1, CurFullY) ne -1 Then begin ContCou=[ContCou,SRight] ContInd=[ContInd,CurInd+1] ContX=[ContX,CurX+1] ContY=[ContY,CurFullY] NumPixelsTotal=NumPixelsTotal+1 Mask(CurX+1,CurFullY)=-1 endif Endif If (RightEdge ne 1) and (UpperEdge ne 1) then begin ;UP,RIGHT URight=Long(RawHalf(CurX+1,CurSubY+1,CosRead))-RawHalf(CurX+1,CurSubY+1, CosRead-1) if (URight) gt Threshold and Mask(CurX+1, CurFullY+1) ne -1 then begin ContCou=[ContCou,URight] ContInd=[ContInd,CurInd+Naxis1+1] ContX=[ContX,CurX+1] ContY=[ContY,CurFullY+1] NumPixelsTotal=NumPixelsTotal+1 Mask(CurX+1,CurFullY+1)=-1 endif Endif If (UpperEdge ne 1) then begin ;UP SUp=Long(RawHalf(CurX,CurSubY+1,CosRead))-RawHalf(CurX,CurSubY+1, CosRead-1) if (SUp) gt Threshold and Mask(CurX, CurFullY+1) ne -1 then begin ContCou=[ContCou,SUp] ContInd=[ContInd,CurInd+Naxis1] ContX=[ContX,CurX] ContY=[ContY,CurFullY+1] NumPixelsTotal=NumPixelsTotal+1 Mask(CurX,CurFullY+1)=-1 endif Endif If (LeftEdge ne 1) and (UpperEdge ne 1) then begin ;UP,LEFT ULeft=Long(RawHalf(CurX-1,CurSubY+1,CosRead))-RawHalf(CurX-1,CurSubY+1, CosRead-1) if abs(ULeft) gt Threshold and Mask(CurX-1, CurFullY+1) ne -1 then begin ContCou=[ContCou,ULeft] ContInd=[ContInd,CurInd+Naxis1-1] ContX=[ContX,CurX-1] ContY=[ContY,CurFullY+1] NumPixelsTotal=NumPixelsTotal+1 Mask(CurX-1,CurFullY+1)=-1 endif Endif If (LeftEdge ne 1) then begin ;LEFT SLeft=Long(RawHalf(CurX-1,CurSubY,CosRead))-RawHalf(CurX-1,CurSubY, CosRead-1) if (SLeft) gt Threshold and Mask(CurX-1, CurFullY) ne - 1 then begin ContCou=[ContCou,SLeft] ContInd=[ContInd,CurInd-1] ContX=[ContX,CurX-1] ContY=[ContY,CurFullY] NumPixelsTotal=NumPixelsTotal+1 Mask(CurX-1,CurFullY)=-1 endif Endif If (LeftEdge ne 1) and (LowerEdge ne 1) then begin ;DOWN,LEFT DLeft=Long(RawHalf(CurX-1,CurSubY-1,CosRead))-RawHalf(CurX-1,CurSubY-1, CosRead-1) if (DLeft) gt Threshold and Mask(CurX-1, CurFullY-1) ne -1 then begin ContCou=[ContCou,DLeft] ContInd=[ContInd,CurInd-Naxis1-1] ContX=[ContX,CurX-1] ContY=[ContY,CurFullY-1] NumPixelsTotal=NumPixelsTotal+1 Mask(CurX-1,CurFullY-1)=-1 endif Endif If (LowerEdge ne 1) then begin ;DOWN SDown=Long(RawHalf(CurX,CurSubY-1,CosRead))-RawHalf(CurX,CurSubY-1, CosRead-1) if (SDown) gt Threshold and Mask(CurX, CurFullY-1) ne -1 then begin ContCou=[ContCou,SDown] ContInd=[ContInd,CurInd-Naxis1] ContX=[ContX,CurX] ContY=[ContY,CurFullY-1] NumPixelsTotal=NumPixelsTotal+1 Mask(CurX,CurFullY-1)=-1 endif Endif If (RightEdge ne 1) and (LowerEdge ne 1) then begin ;DOWN,RIGHT DRight=Long(RawHalf(CurX+1,CurSubY-1,CosRead))-RawHalf(CurX+1,CurSubY-1, CosRead-1) if (DRight) gt Threshold and Mask(CurX+1, CurFullY-1) ne -1 then begin ContCou=[ContCou,DRight] ContInd=[ContInd,CurInd-Naxis1+1] ContX=[ContX,CurX+1] ContY=[ContY,CurFullY-1] NumPixelsTotal=NumPixelsTotal+1 Mask(CurX+1,CurFullY-1)=-1 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)+1] YExtentArr=[YExtentArr,Max(ContY)-Min(ContY)+1] 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 ;CLASSIFICATION=================================================================== ;Set TVFlag equal to 1 if you want to see the event If (HitSizeArr(NoHits) lt 5) then begin ;Tiny Spot print,"Classification: Spot" Type = 'Spot' endif else if (EArr(noHits) lt 0.08) and $ (XExtentArr(NoHits) gt 8) and (YExtentArr(NoHits) gt 8) and $ (HitSizeArr(NoHits)/(XExtentArr(NoHits)*YExtentArr(NoHits))) gt .64 and $ PossibleDelta(0) eq -1 and PossibleWorm(0) eq -1 $ then begin ;Big Fake Star print,"Classification: Big 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.20 then begin ;Medium Fake Star print,"Classification: Big Spot" Type='Spot' 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 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 and $ PossibleDelta(0) eq -1 and PossibleWorm(0) eq -1 $ then begin 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 print,"Classification: Spot" Type='Spot' endif else if PossibleWorm(0) ne -1 then begin print,"Classification: Worm" Type='Worm' endif else if PossibleDelta(0) ne -1 then begin print,"Classification: Delta" Type='Delta' endif else if Lambda1Arr(Nohits) gt 2 and $ Lambda2Arr(NoHits)/Lambda1Arr(NoHits) lt 0.45 then begin print,"Classification: Muon" Type='Muon' endif else begin print,"Classification: Unidentified" Type='UnIdent' endelse ;Allow for one more check on the muons that might have been delta or worm flagged if ((XExtentArr(NoHits) gt 5) or (YExtentArr(NoHits) gt 5)) and $ (Lambda2Arr(NoHits)/Lambda1Arr(NoHits) lt 0.1) then begin NoMuons=NoMuons+1 print,"Classification: Muon" Type='Muon' endif ;ADD TO TOTALS if Type eq 'FakeStar' then begin NoFakeStars=NoFakeStars+1 endif else if Type eq 'Spot' then begin NoSpots=NoSpots+1 endif else if Type eq 'Worm' then begin NoWorms=NoWorms+1 endif else if Type eq 'Delta' then begin NoDeltas=NoDeltas+1 endif else if Type eq 'Muon' then begin NoMuons=NoMuons+1 endif else if Type eq 'UnIdent' then begin NoUnidentified=NoUnidentified+1 endif ;============================================================================== ;************************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,' ',FullY 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) endif ThisViewingBorderX = (ViewingBorder) > XExtentArr(NoHits) XLeft = (ContXAvg-ThisViewingBorderX/2) > $ KeyStr.Ref/2 XRight = (ContXAvg+ThisViewingBorderX/2) < $ (KeyStr.NAxis1-KeyStr.Ref/2-1) ThisViewingBorderY = (ViewingBorder) > YExtentArr(NoHits) YBottom = (ContYAvg-FrThird*KeyStr.Naxis2/4-ThisViewingBorderY/2) > $ KeyStr.Ref/2 YTop = (ContYAvg-FrThird*KeyStr.Naxis2/4+ThisViewingBorderY/2) < $ (KeyStr.NAxis2-KeyStr.Ref/2-1) TvIm = float(RawHalf(XLeft:XRight, YBottom:YTop, CosRead))-$ RawHalf(XLeft:XRight, YBottom:YTop, CosRead-1) if tvflag eq 1 then begin window,0,xsize=1000,ysize=768 plot, [0],[0], /nodata, background=fsc_color('white') MinTVIm = mean(TvIm)-2*Stddev(TvIm) MaxTVIm = mean(TvIm)+2*Stddev(TvIm) tvimage,congrid(bytscl(TvIm, min=MinTVIm, max=MaxTVIm), $ 400,400,/center), position=[0.01,0.01,0.5,0.5], $ background=fsc_color('white') if Stregex(Type, 'Muon', /boolean) then begin window, 1 linefit_xrg, xordv, yordv, Mu_m, Mu_b, sigm, sigb, Error, MuED, MR XRat = (400./(XRight-XLeft)) YRat = (400./(YTop-YBottom)) AspR = YRat/XRat Xvals = indgen((max(absxordh)-min(absxordh))*XRat) XMinIm = XRat*(min(absxordh)-xleft) YMinIm = YRat*(min(absyordh)-ybottom) tvimage,congrid(bytscl(TvIm), 400,400,/center), $ background=fsc_color('white') !p.noerase=1 plot, Xvals+XMinIm, Mu_m*AspR*(Xvals)+YMinIm, $ color=fsc_color('red'), $ /normal, xrange=[0,400], yrange=[0,400], xstyle=5, ystyle=5 wset, 0 endif !p.noerase=1 endif else if tvflag eq 2 then begin MinTVIm = mean(TvIm)-2*Stddev(TvIm) MaxTVIm = mean(TvIm)+2*Stddev(TvIm) tvimage,congrid(bytscl(TvIm, min=MinTVIm, max=MaxTVIm),$ 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 ;Figure out the reads that are valuable to plot to see ;a dip after the cosmic hit FirstRead = 0 > (CosRead - 10) LastRead = (CosRead+80) < KeyStr.NAxis3-1 ;Now to try to sort the other pixels in order of intensity Incr = reverse(sort(ContCou)) XInc = ContX(Incr)+ContXAvg YInc = ContY(Incr)+ContYAvg-FrThird*KeyStr.Naxis2/4 ContCouInc = ContCou(Incr) MinRawRamp = Min(RawHalf(XInc(0),YInc(0),FirstRead:LastRead)) MaxRawRamp = Max(RawHalf(XInc(0),YInc(0),FirstRead:LastRead)) RawBias = RawRamp(FirstRead) plot, RawHalf(XInc(0),YInc(0),FirstRead:LastRead) - $ long(RawHalf(XInc(0),YInc(0),FirstRead)), $ yrange=[MinRawRamp-long(RawHalf(XInc(0),YInc(0),FirstRead)),$ long(MaxRawRamp)-long(RawHalf(XInc(0),YInc(0),FirstRead))], $ position = [0.05,0.60,0.5,0.95], $ color=fsc_color('black'), background=fsc_color('white') MaxI = 9 < N_Elements(Incr) for i=0, MaxI-1 do begin oplot, long(RawHalf(XInc(i), YInc(i), FirstRead:LastRead))-$ RawHalf(XInc(i),YInc(i),FirstRead), $ color=ColorNum.ColorNums(i) endfor HLineInd=where(ContY eq 0) if HLineInd(0) eq -1 then HLineInd = 0 VLineInd=where(ContX eq 0) if VLineInd(0) eq -1 then VLineInd = 0 HLineContX=ContX(HLineInd) HLineContCou=ContCou(HLineInd) VLineContY=ContY(VLineInd) VLineContCou=ContCou(VLineInd) VSort=sort(VLineContY) VLineContY=VLineContY(VSort) VLineContCou=VLineContCou(VSort) if Not (Stregex(Type, 'Muon', /boolean)) then begin plot,VLineContY,VLineContCou-Threshold,$ psym=2,title='Vertical Profile at x=0',$ position=[0.55,0.05,0.95,.5],/norm, $ color=fsc_color('black'), background=fsc_color('white') plot,HLineContX,HLineContCou-Threshold,$ psym=2,title='Horizontal Profile at y=0',$ position=[0.55,0.55,0.95,.95],/norm, $ color=fsc_color('black'), background=fsc_color('white') endif if TvFlag eq 1 then stop endif ;****MUON DEPTH DETERMINATION*************************************** ;Go through XVals or YVals depending on which is bigger and make an estimate of depth If stregex(Type, 'Muon', /boolean) then begin linefit_xrg, xordv, yordv, Mu_m, Mu_b, sigm, sigb, Error, MuED, MR ;Muon track is horizontal or 45 degrees if XExtentArr(NoHits) ge YExtentArr(NoHits) then begin MuXOrd = AbsXordH MuYOrd = AbsYordH ThisOrd = MuXOrd endif else begin ;Muon track is vertical or > 45 degrees MuXOrd = AbsXordV MuYOrd = AbsYordV ThisOrd = MuYOrd endelse ;Create Array for Charge Lost Array ChargeLostArr = fltarr(HitSizeArr(NoHits)) ChargeGainArr = fltarr(HitSizeArr(NoHits)) ChargeDepoArr = fltarr(HitSizeArr(NoHits)) For i = 0, HitSizeArr(NoHits)-1 do begin BorderPix = 0 ThisX = MuXOrd(i) ThisSubY = MuYOrd(i)-FrThird*KeyStr.Naxis2/4 ThisFullY = MuYOrd(i) RawRamp = long(RawHalf(ThisX, ThisSubY, *)) Classify_One_Pixel, ThisX, ThisFullY, RawRamp, CosRead, BorderPix, $ AvgTemp, KeyStr.VSUB, $ XCen, YCen, CosEDep, ChargeLost, ChargeRise, TvFlag, TvPPT, 0, Verbose, $ MySQLLun=MySQLLun, DarkSub=DarkSub ChargeLostArr(i) = ChargeLost ChargeGainArr(i) = ChargeRise ChargeDepoArr(i) = CosEDep EndFor If TvFlag gt 0 then begin plot, ChargeDepoArr(where(ChargeDepoArr gt MuThresh)), $ ChargeLostArr(where(ChargeDepoArr gt MuThresh)), $ psym=2,title='Charge Lost Vs. Deposit',$ position=[0.55,0.05,0.95,.5],/norm, $ color=fsc_color('black'), background=fsc_color('white') plot, ThisOrd(where(chargedepoarr gt MuThresh)), $ chargelostarr(where(chargedepoarr gt MuThresh))/chargedepoarr(where(chargedepoarr gt MuThresh)),$ psym=2,title='Charge Ratio Vs Position',$ position=[0.55,0.55,0.95,.95],/norm, $ color=fsc_color('black'), background=fsc_color('white') stop EndIf ;Now figure out where muon entered and left detector based on charge lost ;It is very conclusive that muon loses more charge near top of detector than ;closer to pixel based upon 5.0 Volt VSUB data ChargePoints = where(chargedepoarr gt MuThresh) if ChargePoints(0) ne -1 then begin MuFit=linfit(ThisOrd(ChargePoints), $ ChargeLostArr(ChargePoints)/ChargeDepoArr(ChargePoints)) endif else begin MuFit = [0,1] endelse MuSlope = MuFit(1) ;MuSlope is the slope of the counts vs. position; Mu_m is the slope of the track if (MuSlope ge 0 and (Mu_M ge 0 and Mu_M lt 1)) or $ (MuSlope lt 0 and (Mu_M ge 1)) then begin XFrontSide = min(MuXOrd) YFrontSide = min(MuYOrd) XBackSide = max(MuXOrd) YBackSide = max(MuYOrd) endif else if (MuSlope ge 0 and (Mu_M lt 0 and Mu_M gt -1)) or $ (MuSlope lt 0 and (Mu_M le -1)) then begin XFrontSide = min(MuXOrd) YFrontSide = max(MuYOrd) XBackSide = max(MuXOrd) YBackSide = min(MuYOrd) endif else if (MuSlope lt 0 and (Mu_M ge 0 and Mu_M lt 1)) or $ (MuSlope gt 0 and (Mu_M ge 1)) then begin XFrontSide = max(MuXOrd) YFrontSide = max(MuYOrd) XBackSide = min(MuXOrd) YBackSide = min(MuYOrd) endif else if (MuSlope lt 0 and (Mu_M lt 0 and Mu_M gt -1)) or $ (MuSlope gt 0 and (Mu_M le -1)) then begin XFrontSide = max(MuXOrd) YFrontSide = min(MuYOrd) XBackSide = min(MuXOrd) YBackSide = max(MuYOrd) endif TrackLength = sqrt((YBackSide-YFrontSide)^2+(XBackSide-XFrontSide)^2) 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 ;Hold on to the largest rise and the largest decay MaxChargeRise = -1 MaxChargeLost = -1 MaxCosEDep = -1 TotChargeRise = -1 TotChargeLost = -1 TotCosEDep = -1 ;Classify Individual pixels and fit the decay For i = 0, HitSizeArr(NoHits)-1 do begin BorderPix = BorderMask(i) ThisX = AbsXordV(i) ThisSubY = AbsYordV(i)-FrThird*KeyStr.Naxis2/4 ThisFullY = AbsYordV(i) ;If it's a muon, calculate distance along track If Stregex(Type, 'Muon', /boolean) then begin Z = (sqrt((ThisX-XBackSide)^2+(ThisFullY-YBackSide)^2)/TrackLength)*100 endif else begin Z = 0 endelse RawRamp = long(RawHalf(ThisX, ThisSubY, *)) Classify_One_Pixel, ThisX, ThisFullY, RawRamp, CosRead, BorderPix, $ AvgTemp, KeyStr.VSUB, $ XCen, YCen, CosEDep, ChargeLost, ChargeRise, TvFlag, TvPPT, UpdateDB, Verbose, $ MySQLLun=MySQLLun, DarkSub=DarkSub, ZDepth=Z MaxChargeRise = MaxChargeRise > ChargeRise MaxChargeLost = MaxChargeLost > ChargeLost MaxCosEDep = MaxCosEDep > CosEDep TotChargeRise = TotChargeRise + (ChargeRise > 0) TotChargeLost = TotChargeLost + (ChargeLost > 0) TotCosEDep = TotCosEDep + (CosEDep > 0) EndFor ;ooooooooooooooooooooooooooooooooooooooooooooDATABASE If UpdateDB eq 1 then begin CosID = KeyStr.RawObjectKey+'_'+strtrim(XCen,2)+'_'+strtrim(YCen,2) GloCmd = 'INSERT INTO ' + GlobalTableName+' VALUES('+$ "'"+strtrim(CosID,2)+"',"+$ strtrim(XCen,2)+','+$ strtrim(YCen,2)+','+$ strtrim(CosRead,2)+','+$ strtrim(NumPix,2)+','+$ strtrim(long(MaxCosEDep),2)+','+$ strtrim(long(MaxChargeLost),2)+','+$ strtrim(long(MaxChargeRise),2)+','+$ strtrim(long(TotCosEDep),2)+','+$ strtrim(long(TotChargeLost),2)+','+$ strtrim(long(TotChargeRise),2)+','+$ strtrim(XExtentArr(NoHits),2)+','+$ strtrim(YExtentArr(NoHits),2)+','+$ strtrim(TrackLength2DArr(NoHits),2)+','+$ strtrim(PerpCountsArr(NoHits),2)+','+$ strtrim(string(EArr(NoHits),format='(F10.3)'),2)+','+$ strtrim(string(Lambda1Arr(NoHits),format='(F10.2)'),2)+','+$ strtrim(string(Lambda2Arr(NoHits),format='(F10.2)'),2)+','+$ "'"+strtrim(Type, 2)+"',"+$ "'"+strtrim(KeyStr.ThisFilter,2)+"',"+$ strtrim(AvgTemp,2)+","+$ "'"+strtrim(KeyStr.Date,2)+"'"+$ ");" mysqlcmd, MySQLLun, GloCmd, Answer, NLines EndIf ;oooooooooooooooooooooooooooooooooooooooooooooooooooooo endif endif endfor delvarx, RawHalf endfor endfor !p.noerase=0 end