;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 CosmicSecondMoments,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 = 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(BorderSize)) then BorderSize = 5
;PRIMITIVE SUBTRACTION AND FLAT FIELD

VR = 5;	Viewing Region
  ;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)
  Cosmic=findgen(2*VR+1)
  IRegion=fltarr(2*VR+1,2*VR+1)
  ;IRegion(Cosmic,Cosmic)=10;(Cosmic-5)^2
  ;IRegion(5,Cosmic)=Cosmic
  ;IRegion(Cosmic,5)=10
  IRegion(Cosmic,Cosmic/2+2)=10
  IRegionNorm=float(IRegion)/total(IRegion)
  print,total(XRegion*IRegionNorm),total(YRegion*IRegionNorm)

  XAvgReg=fix(total(XRegion*IRegionNorm))
  YAvgReg=fix(total(YRegion*IRegionNorm))
  IRegionNorm=float(IRegion)/total(IRegion)
  print,total(XRegion*IRegionNorm),total(YRegion*IRegionNorm)
  X11=total(XRegion^2*IRegionNorm)
  X22=total(YRegion^2*IregionNorm)
  X12=total(XRegion*YRegion*IRegionNorm)
  XMat=[[X11,X12],[X12,X22]]
  Evals=Hqr(Elmhes(XMat),/Double)
  window,xsize=868,ysize=868
  tvscl,congrid(IRegion,768,768,/center)
  stop
stop

end