;BoxAnalyseGeant4Image
;
;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
;   VR:
;          The size of the region in which the second moments are calculated
;   THRESHOLD:
;          The count value that constitutes a hit rather than noise.
;   SEARCHREGION:
;	   The width and length of the box in pixels where a maximum
;	   count is found.

Pro analysegeant4Image,date=date,SingleImage=SingleImage,$
    Pitch=Pitch,Thickness=Thickness,NEvents=NEvents,$
    PType=Ptype,OpticsFlag=OpticsFlag,BorderSize=BorderSize,$
    Threshold=Threshold,SearchRegion=SearchRegion
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(Threshold)) then Threshold = 0
If (not keyword_set(SearchRegion)) then SearchRegion = 5
If (not keyword_set(VR)) then VR = 7
If (not keyword_set(BorderSize)) then BorderSize = VR+1

;PRIMITIVE SUBTRACTION AND FLAT FIELD
GEANT4Dir='/nfs/slac/g/ki/ki08/lsst/lsstgeant4/'
ParamsString=strtrim(Pitch,2)+'x'+strtrim(Pitch,2)+'x'+$
             strtrim(Thickness,2)+'Microns'+strtrim(NEvents,2)+$
             PType+OpticsFlag
FitsFileName=GEANT4Dir+ParamsString+'.fits'
fits_read,FitsFileName,Im1,header1
Naxis=SXPAR(header1,'NAXIS')
Naxis1=SXPAR(header1,'NAXIS1')
Naxis2=SXPAR(header1,'NAXIS2')
Im1Hist=histogram(Im1,Nbins=100,max=10000,locations=locIm1)

Im1MB=Im1(BorderSize-1:Naxis1-1-Bordersize,BorderSize-1:Naxis2-1-BorderSize)

Pass = 0
for i=-SearchRegion,SearchRegion do begin
  for j=-SearchRegion,SearchRegion do begin
    if not ((i eq 0) and (j eq 0)) then begin 
      FirstCut=where(Im1MB gt Im1(Bordersize-1-i:Naxis2-1-BorderSize-i,$
		            BorderSize-1-j:Naxis2-1-BorderSize-j))
      if Pass eq 0 then begin
        Highest=FirstCut
      endif else begin
        Highest=[Highest,Firstcut]
        Highest=Highest(sort(Highest))
        Dups=where(Highest(0:N_elements(Highest)-2) eq Highest(1:N_elements(Highest)-1))
        Highest=Highest(Dups)
      endelse
      Pass=Pass+1
    endif
  endfor
endfor

;/////////////////////////////////////////////////////////////////////
;FILL ARRAYS WITH DATA FROM COSMIC RAYS
NumHits=N_elements(Highest)
Lambda1Arr=fltarr(NumHits)        ;The long second moment
Lambda2Arr=fltarr(NumHits)	  ;The short second moment
EArr=fltarr(NumHits)	 	  ;Scalar Ellipticity
FluxArr=fltarr(NumHits)		  ;Total Counts in a region
TrackLength2DArr=fltarr(NumHits)  ;Length of the track in Pixels 2D
TrackLength3DArr=fltarr(NumHits)  ;Lenght of the track in physical units
PerpCountsArr=fltarr(NumHits)	  ;Counts scaled by length
HitSizeArr=fltarr(NumHits)	  ;Total number of pixels affected
HitWidthArr=fltarr(NumHits)	  ;HitWidth as judged by Length
XExtentArr=Uintarr(NumHits)	  ;Extent of Hit in X pixels
YExtentArr=Uintarr(NumHits)	  ;Extent of Hit in Y pixels
FillFactorArr=fltarr(NumHits)     ;Amount of space filled in hit

NoSpots=0			  ;Number of "Spots" 
NoMuons=0			  ;Number of "Muons"
NoWorms=0			  ;Number of "Worms"
NoPossStars=0			  ;Number of "Possible Star" Candidates
NoDeltas=0			  ;Number of "Delta Rays"
;//////////////////////////////////////////////////////////////////////

for Cosmics=0, NumHits-1 do begin 
  ;Return Coordinate in overall Image
  XCosmic=(Highest(Cosmics) mod (Naxis1+1-2*BorderSize))+BorderSize-1
  YCosmic=(Highest(Cosmics)/(Naxis2+1-2*BorderSize))+BorderSize-1
  ;Create Analysis Region
  XRegion=rebin(findgen(2*VR+1)-VR,2*VR+1,2*VR+1)
  YRegion=rebin(transpose(findgen(2*VR+1)-VR),2*VR+1,2*VR+1)
  IRegion=Im1(XCosmic-VR:XCosmic+VR,YCosmic-VR:YCosmic+VR)
  IRegionNorm=float(IRegion)/total(IRegion)
  ;Get Average X and Y Positions.  Redefine Analysis Region around them.
  XAvgReg=fix(total(XRegion*IRegionNorm))
  YAvgReg=fix(total(YRegion*IRegionNorm))
  IRegion=Im1(XCosmic-VR+XAvgReg:XCosmic+VR+XAvgReg,$
	      YCosmic-VR+YAvgReg:YCosmic+VR+YAvgReg)
  FluxArr(Cosmics)=total(IRegion)
  IRegionNorm=float(IRegion)/FluxArr(Cosmics)
  ;Calculate Second Moments and find them in rotated frame
  I11=total(XRegion^2*IRegionNorm)
  I22=total(YRegion^2*IregionNorm)
  I12=total(XRegion*YRegion*IRegionNorm)
  E1=(I11-I22)/(I11+I22)
  E2=(2*I12)/(I11+I22)
  XMat=[[I11,I12],[I12,I22]]
  Evals=Hqr(Elmhes(XMat),/Double)
  Lambda1Arr(Cosmics)=Real_part(max(Evals))
  Lambda2Arr(Cosmics)=Real_part(min(Evals))
  EArr(Cosmics)=E1^2+E2^2
  ;Get Track Length and Number of Perpendicular Counts
  XsCosmic=where(IRegionNorm ne 0) mod (2*VR+1)
  YsCosmic=where(IRegionNorm ne 0)/ (2*VR+1)
  XExtentArr(Cosmics)=Max(XsCosmic)-Min(XsCosmic)
  YExtentArr(Cosmics)=Max(YsCosmic)-Min(YsCosmic)
  TrackLength3DArr(Cosmics)=sqrt(Pitch*XExtentArr(Cosmics)^2+$	
		   	         Pitch*YExtentArr(Cosmics)^2+$
				Thickness^2)
  TrackLength2DArr(Cosmics)=sqrt(XExtentArr(Cosmics)^2+$
                                 YExtentArr(Cosmics)^2)+1
  PerpCountsArr(Cosmics)=FluxArr(Cosmics)*Thickness/TrackLength3DArr(Cosmics)
  HitPixels=where(IRegionNorm gt Threshold)
  HitSizeArr(Cosmics)=N_elements(HitPixels)
  HitWidthArr(Cosmics)=HitSizeArr(Cosmics)/TrackLength2DArr(Cosmics)
  FillFactorArr(Cosmics)=HitSizeArr(Cosmics)/$
			 float((1+XExtentArr(Cosmics))*(1+YExtentArr(Cosmics)))
  ;Try to find gaps in the pixel arrays HS=Horizontally Sorted VS=Vertically
  PossibleBreaks=-1
  PossibleDeltas=-1
  DefiniteBreaks=-1

  if N_elements(HitPixels) gt 2 then begin
   HitPixels=HitPixels(sort(HitPixels))
   CosmicsV=(XsCosmic-1)*(2*VR+1)+YsCosmic
   XsCosmicHS=XsCosmic(sort(HitPixels))
   YsCosmicHS=YsCosmic(sort(HitPixels))
   XsCosmicVS=XsCosmic(sort(CosmicsV))
   YsCosmicVS=YsCosmic(sort(CosmicsV))
   XsLines=XsCosmic(sort(XsCosmic))
   YsLines=YsCosmic(sort(YsCosmic))
   XsLineBreaks=XsLines(1:HitSizeArr(Cosmics)-1)-$
                XsLines(0:HitSizeArr(Cosmics)-2)
   YsLineBreaks=YsLines(1:HitSizeArr(Cosmics)-2)-$
                YsLines(0:HitSizeArr(Cosmics)-1)
   XsBreakH=XsCosmicHS(1:HitSizeArr(Cosmics)-1)-$
           XsCosmicHS(0:HitSizeArr(Cosmics)-2)
   YsBreakH=YsCosmicHS(1:HitSizeArr(Cosmics)-1)-$
           YsCosmicHS(0:HitSizeArr(Cosmics)-2)
   XsBreakV=XsCosmicVS(1:HitSizeArr(Cosmics)-1)-$
           XsCosmicVS(0:HitSizeArr(Cosmics)-2)
   YsBreakV=YsCosmicVS(1:HitSizeArr(Cosmics)-1)-$
           YsCosmicVS(0:HitSizeArr(Cosmics)-2)
   DefiniteBreaks=where((XsLineBreaks gt 1) or (YsLineBreaks gt 1))
   PossibleBreaksH=where((abs(XsBreakH) gt 1) and (abs(YsBreakH) gt 1))
   PossibleBreaksV=where((abs(XsBreakV) gt 1) and (abs(YsBreakV) gt 1))
   PossibleDeltas=where((((abs(XsBreakH) gt 1) and (abs(YsBreakH) eq 0)) or $
			 ((abs(XsBreakH) eq 0) and (abs(YsBreakH) gt 1))) or $
			(((abs(XsBreakV) gt 1) and (abs(YsBreakV) eq 0)) or $
                         ((abs(XsBreakV) eq 0) and (abs(YsBreakV) gt 1))))

  endif 
  ;Diagnose Type of Hit: Spot, Worm, or Delta
  if DefiniteBreaks(0) ne -1 then begin
     print,"Double Hit in Viewing Region"
     window,xsize=512,ysize=512
     tvscl,congrid(alog(IRegion),512,512,/center)
     print,XExtentArr(Cosmics),YExtentArr(Cosmics),FillFactorArr(Cosmics)
     print,Earr(Cosmics)
     PossibleDeltas=-1
     stop
 
  endif else if PossibleDeltas(0) ne -1 then begin
     NoDeltas=NoDeltas+1
     window,xsize=512,ysize=512
     tvscl,congrid(alog(IRegion),512,512,/center)
     print,XExtentArr(Cosmics),YExtentArr(Cosmics),FillFactorArr(Cosmics)
     print,Earr(Cosmics)
     PossibleDeltas=-1
     stop
  endif else if TrackLength2DArr(Cosmics) eq 1 then begin
     NoSpots=NoSpots+1
  endif else if ((FillFactorArr(Cosmics) gt .5) and $
	         (EArr(Cosmics) lt .15)) then begin
     NoWorms=NoWorms+1
  endif 

endfor

;////////////////////////////////////////////////////////////////////////
;HISTOGRAMS AND SUCH

!p.multi=[0,2,2]
HitSizeHist=histogram(HitSizeArr,Nbins=25,max=25,locations=locHitSize)
plot,locHitSize,HitSizeHist,psym=10,title='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'
stop
lambda1Hist=histogram(Lambda1Arr,Nbins=10,$
	    	      locations=locLambda1)
plot,locLambda1,Lambda1Hist,psym=10,title='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,yrange=[0,2],xrange=[0,2],$
		      xtitle='Lambda 2',ytitle='Perp Counts'
stop

end