;AnalyseGeant4Image
;
;Purpose: This script will attempt to find the cosmic rays in a monte carlo
;	  image and characterize them by their shape and intensity


Pro tvwormsgeant4image,date=date,SingleImage=SingleImage,$
    Pitch=Pitch,Thickness=Thickness,NEvents=NEvents,$
    PType=Ptype,OpticsFlag=OpticsFlag,BorderSize=BorderSize

If (not keyword_set(SingleImage)) then SingleImage = 1
If (not keyword_set(Pitch)) then Pitch = 10
If (not keyword_set(Thickness)) then Thickness = 100
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(BorderSize)) then BorderSize = 8
;PRIMITIVE SUBTRACTION AND FLAT FIELD
GEANT4Dir='/nfs/slac/g/ki/ki08/lsst/lsstgeant4/'
FitsFileName=GEANT4Dir+strtrim(Pitch,2)+'x'+strtrim(Pitch,2)+'x'+$
	     strtrim(Thickness,2)+'Microns'+strtrim(NEvents,2)+$
	     PType+OpticsFlag+'.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)

SearchRegion=5
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)
Lambda2Arr=fltarr(NumHits)
FluxArr=fltarr(NumHits)
TrackLengthArr=fltarr(NumHits)
PerpCountsArr=fltarr(NumHits)

VR = 7;	Viewing Region
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)
  XMat=[[I11,I12],[I12,I22]]
  Evals=Hqr(Elmhes(XMat),/Double)
  Lambda1Arr(Cosmics)=Real_part(max(Evals))
  Lambda2Arr(Cosmics)=Real_part(min(Evals))
  ;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)
  TrackLengthArr(Cosmics)=sqrt((Pitch*(max(XsCosmic)-min(XsCosmic)))^2+$	
		   	       (Pitch*(max(YsCosmic)-min(YsCosmic)))^2+$
				Thickness^2)
  PerpCountsArr(Cosmics)=FluxArr(Cosmics)*Thickness/TrackLengthArr(Cosmics)

  if (lambda1arr(cosmics) gt 2) and (lambda2arr(cosmics) gt 1.5) then begin
       window,xsize=512,ysize=512
       tvscl,congrid(IRegion,512,512,/center)
       stop
  endif
endfor
stop

end