;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 finds regions of contiguous pixels and classifies them using 
;	  the function ReturnEvents.pro
;
;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
;   THRESHOLD:
;          The count value that constitutes a hit rather than noise.
;
;///////////////////////////////////////////////////////////////////////
;\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

Pro analysegeant4Image,date=date,SingleImage=SingleImage,$
    Pitch=Pitch,Thickness=Thickness,NEvents=NEvents,$
    PType=Ptype,OpticsFlag=OpticsFlag,$
    Sigma=Sigma,FileName=FileName,$
    ThresholdSlope=ThresholdSlope,$
    MinThreshold=MinThreshold,TvFlag=TvFlag

If (not keyword_set(SingleImage)) then SingleImage = 1
If (not keyword_set(Pitch)) then Pitch = 10
If (not keyword_set(Thickness)) then Thickness = 200
If (not keyword_set(NEvents)) then NEvents = 1600
If (not keyword_set(PType)) then Ptype = 'Muons'
If (not keyword_set(OpticsFlag)) then OpticsFlag = 'OpticsOn'
If (not keyword_set(MinThreshold)) then MinThreshold = 0
If (not keyword_set(ThresholdSlope)) then ThresholdSlope = 0
If (not keyword_set(Sigma)) then Sigma = 1.1
If (not keyword_set(TvFlag)) then TvFlag = 1

;PRIMITIVE SUBTRACTION AND FLAT FIELD
GEANT4Dir='~/Desktop/Cosmic_Rays/SingleEData/'
GEANT4DataDir=GEANT4Dir
ParamsString=strtrim(Pitch,2)+'x'+strtrim(Pitch,2)+'x'+$
             strtrim(Thickness,2)+'Microns'+strtrim(NEvents,2)+$
             PType+OpticsFlag
FitsFileName=GEANT4DataDir+ParamsString+'.fits'
fits_read,FitsFileName,Im1,header1
Naxis=SXPAR(header1,'NAXIS')
Naxis1=SXPAR(header1,'NAXIS1')
Naxis2=SXPAR(header1,'NAXIS2')

if TvFlag eq 1 then begin
   set_plot,'x'
endif else if TVFlag eq 2 then begin
   set_plot,'ps'
   device,filename=Geant4DataDir+ParamsString+'.ps'
   loadct,0	                   ;load color table 39
   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=1.0
   !P.CHARTHICK=2
   !P.THICK=4
endif

!p.multi=[0,1,1]
tv,congrid(Im1,1200,1200,/center)
ReturnArr=ReturnEvents(Im1,Naxis1,Naxis2,MinThreshold,ThresholdSlope,$
                       Sigma,Pitch,Thickness,TvFlag)

;////////////////////////////////////////////////////////////////////////
;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
loadct,39
!p.multi=[0,1,1]
bar_plot,returnArr(11,0:5),$
         barnames=['Spots','Deltas','Worms','"Stars"','Muons','UnIdentified'],$
	 title='Distribution of Events for Simulated Data'
xyouts,0.5,0.90,'Thickness  ='+Strtrim(Thickness,2)+' Microns',/norm
xyouts,0.5,0.85,'Pitch	    ='+Strtrim(Pitch,2)+' Microns',/norm

!p.multi=[0,2,2]
HitSizeHist=histogram(HitSizeArr,Nbins=25,max=25,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^(-4.),psym=2,$
		      xtitle='Lambda 2',ytitle='Perp Counts'
if tvflag eq 2 then device,/close
stop

end