;AnalyseGeant4Image
;
;Purpose: This script will attempt to find the cosmic rays in a monte carlo
;	  image and characterize them by their shape and intensity
;
;	  The count of a given pixel will be zero unless hit by a cosmic ray,
;	  but the keyword THRESHOLD is included for adaptation purposes.
;	  
;	  The script goes through the image and finds pixels that 
;  	  have a greater count then all others within a box of 
;	  SEARCHREGIONxSEARCHREGION centered around it.  It qualifies this
;	  as a hit and takes the 1D coordinate in the array Highest.
;
;	  For each coordinate, it considers a Viewing Region of VRxVR pixels 
;  	  centered around it.  The center coordinate of the event is calculated
;	  and the second moments are calculated with respect to this center.
;	  By diagonalizaing the second moments matrix, the frame in which 
; 	  the second moments are a minimum and maximum is found.  
;
;KEYWORDS:
;
;   DATE:  
;          The date of the night of observing (for real data)
;   SINGLEIMAGE:
;	   1=Look at one image; 0=Look at the difference between 2 images
;   PITCH;
;	   The pitch of the pixels in millimeters
;   THICKNESS:
;          The thickness of the CCD in millimeters
;   NEVENTS:
;          The number of events generated in the GEANT4 Simulation
;   PTYPE:
;	   The type of the particle in the primary event
;   OPTICSFLAG:
;          A flag incidicating whether optics were on or off in simulation
;   BORDERSIZE:
;          The size of the border around the image where events will not be 
;	   flagged in the initial search
;   MINTHRESHOLD:
;          The count value above which a pixel is considered to be part of an 
;	   event. This value is for the first row of the image (y=0)
;   THRESHOLDSLOPE:
;	   The threshold at a given y coordinate is calculated as 
;	   THRESHOLD=MINTHRESHOLD+THRESHOLDSLOPE*Y/NAXIS2
;	
;	   Where the slope is Delta_Threshold/Naxis2
;
Pro analysedarkImage,date=date,SingleImage=SingleImage,$
    Pitch=Pitch,Thickness=Thickness,NEvents=NEvents,$
    PType=Ptype,OpticsFlag=OpticsFlag,BorderSize=BorderSize,$
    Threshold=Threshold,SearchRegion=SearchRegion,$
    Sigma=Sigma,FileName=FileName,$
    ThresholdSlope=ThresholdSlope,$
    MinThreshold=MinThreshold

If (not keyword_set(SingleImage)) then SingleImage = 1
If (not keyword_set(Pitch)) then Pitch = 15.
If (not keyword_set(Thickness)) then Thickness = 25.
If (not keyword_set(NEvents)) then NEvents = 400000
If (not keyword_set(PType)) then Ptype = 'Muons'
If (not keyword_set(OpticsFlag)) then OpticsFlag = 'OpticsOff'
If (not keyword_set(MinThreshold)) then MinThreshold = 1500
If (not keyword_set(ThresholdSlope)) then ThresholdSlope = 0
If (not keyword_set(Sigma)) then Sigma = 1.1
If (not keyword_set(FileName)) then Filename = 'd1365.fits'

;PRIMITIVE SUBTRACTION AND FLAT FIELD
FileParts=strsplit(FileName,'.',/extract)
DarkDir='/nfs/slac/g/ki/ki08/lsst/lsstgeant4/data/'
FitsFileName=DarkDir+FileName
fits_read,FitsFileName,Im1,header1
Naxis=SXPAR(header1,'NAXIS')
Naxis1=SXPAR(header1,'NAXIS1')
Naxis2=SXPAR(header1,'NAXIS2')

;Hot Rows and Columns
Im1(989,*)=0
Im1(1152,*)=0
Im1(2040,*)=0

;Background Varies 
if FileName eq 'd1365.fits' then begin 
	MinThreshold = 20080
	ThresholdSlope = 0
endif else if FileName eq 'd3239.fits' then begin
        MinThreshold = 4280
	ThresholdSlope = 300
endif

TvFlag=2
if TvFlag eq 1 then begin
   set_plot,'x'
endif else if TVFlag eq 2 then begin
   set_plot,'ps'
   device,filename=DarkDir+FileParts(0)+'.ps'
   loadct,0			   ;BlackAndWhite for Original Intensity
   device,/color                   ;allow color on the postscript
   device,ysize=7.5,/inches        ;Height of plot in y
   device,xsize=7.0,/inches        ;Width of plot in x
   device,yoffset=0.50,/inches      ;Y position of lower left corner
   device,xoffset=0.00,/inches
   white='FFFFFF'x
   black='000000'x
   !P.CHARSIZE=0.8
   !P.CHARTHICK=2
   !P.THICK=4
endif

!p.multi=[0,1,1]
tv,congrid(Im1,1024,1024,/center)
ReturnArr=ReturnEvents(Im1,Naxis1,Naxis2,MinThreshold,ThresholdSlope,$
	               Sigma,Pitch,Thickness)
;////////////////////////////////////////////////////////////////////////
;HISTOGRAMS AND SUCH
Lambda1Arr=ReturnArr(0,*)
Lambda2Arr=ReturnArr(1,*)
EArr=ReturnArr(2,*)
FluxArr=ReturnArr(3,*)
TrackLength2DArr=ReturnArr(4,*)
TrackLength3DArr=ReturnArr(5,*)
PerpCountsArr=ReturnArr(6,*)
HitSizeArr=ReturnArr(7,*)
HitWidthArr=ReturnArr(8,*)
XExtentArr=ReturnArr(9,*)
YExtentArr=ReturnArr(10,*)

NoSpots=ReturnArr(11,0)
NoDeltas=ReturnArr(11,1)
NoWorms=ReturnArr(11,2)
NoFakeStars=ReturnArr(11,3)
NoMuons=ReturnArr(11,4)
NoUnidentified=ReturnArr(11,5)
NoHits=N_elements(XExtentArr)

print,"Number of Spots     = ",NoSpots
print,"Number of Deltas    = ",NoDeltas
print,"Number of Muons     = ",NoMuons
print,"Number of Worms     = ",NoWorms
print,"Number of FakeStars = ",NoFakeStars
print,"Number of UCR       = ",NoUnidentified
print,"Total               = ",NoHits

;////////////////////////////////////////////////////////////////////////
;HISTOGRAMS AND SUCH
!p.multi=[0,1,1]
loadct,39
bar_plot,returnArr(11,0:5),$
         barnames=['Spots','Deltas','Worms','"Stars"','Muons','UnIdentified'],$
         title='Distribution of Events for Dark Images From Lick CCD'
xyouts,0.5,0.90,'Thickness     =  200 Microns',/norm
xyouts,0.5,0.85,'Pitch         =   15 Microns',/norm
xyouts,0.5,0.80,'Exposure Time = 1000 Seconds',/norm
!p.multi=[0,2,2]
HitSizeHist=histogram(HitSizeArr,Nbins=40,locations=locHitSize)
plot,locHitSize,HitSizeHist,psym=10,title='Real Data Num pixels hit'
TrackLength3DHist=histogram(TrackLength3DArr,Nbins=10,$
	                  locations=locTrack3DLength)
plot,locTrack3DLength,TrackLength3DHist,psym=10,title='TrackLength'
TrackLength2DHist=histogram(TrackLength2DArr,Nbins=15,$
                          locations=locTrack2DLength)
plot,locTrack2DLength,TrackLength2DHist,psym=10,title='TrackLength 2D'
HitWidthHist=histogram(HitWidthArr,Nbins=15,$
                          locations=locHitWidth)
plot,locHitWidth,HitWidthHist,psym=10,title='Track Width'
lambda1Hist=histogram(Lambda1Arr,Nbins=10,$
	    	      locations=locLambda1)
plot,locLambda1,Lambda1Hist,psym=10,title='Real Data Lambda 1'
lambda2Hist=histogram(Lambda2Arr,Nbins=10,$
                      locations=locLambda2)
plot,locLambda2,Lambda2Hist,psym=10,title='Lambda 2'
lambda2Hist=histogram(Lambda2Arr,Nbins=10,$
		      locations=locLambda2)
plot,Lambda1Arr,Lambda2Arr,psym=2,xtitle='Lambda 1',ytitle='Lambda 2'
plot,lambda2arr,perpcountsarr*10^(-5.),psym=2,$
		      xtitle='Lambda 2',ytitle='Perp Counts'

if tvflag eq 2 then device,/close
stop

end