;---------------------------------------------------------------------------
; Document name: hsi_modul_pattern__define
; Created by:    Andre Csillaghy, May 1999
;
; Last Modified: Tue Jul 31 17:18:02 2001 (csillag@sunlight)
;---------------------------------------------------------------------------
;
;+
; PROJECT:
;       HESSI
;
; NAME:
;       Hsi_Mpat__define
;
; PURPOSE:
;       Defines the multi energy image to count response
;
; CATEGORY:
;       HESSI Utilities
;
; CONSTRUCTION:
;      Based on hsi_square_mpat. Provides image response matrix on
;		getdata for selected energy band (eb_index control).  Reprocesses
;	   on changes in image or calib_eventlist parameters.  Builds up vectors needed
;	  to differentiate the response at different energies and stores these.
;
;
; SEE ALSO:
;      HESSI Utility Reference http://hessi.ssl.berkeley.edu/software/reference.html
;
; HISTORY:
;
;
;----------------------------------------------------------

FUNCTION HSI_Mpat::INIT, SOURCE = source, _EXTRA=_extra

IF NOT Obj_Valid( source ) THEN source = HSI_modul_pattern() ; has calib_eventlist and modul_pattern so it can do everything
ohmp = Source->Get(/obj, class='hsi_modul_pattern')
control = hsi_mpat_control()
CheckVar, info, {hsi_mpat_info}

ret=self->Framework::INIT( SOURCE = ohmp, CONTROL=control, INFO=info )
Self->Set, _extra=_extra
;strategy_available =  ['HSI_ANNSEC_PATTERN', 'HSI_VISMOD_PATTERN']
;ret=self->Strategy_Holder::INIT( strategy_available, $
;                                 CONTROL=hsi_mpat_control(), $
;                                 INFO={hsi_mpat_info}, $
;                                 SOURCE=source, $
;                                 _EXTRA=_extra )

;self->SetStrategy, 'HSI_ANNSEC_PATTERN'
;self->Set, IMG_STRATEGY_AVAILABLE = strategy_available, /THIS_CLASS_ONLY

RETURN, ret

END
;---------------------------------------------------------------------------

PRO HSI_MPAT::Set, $
    sp_energy_binning=sp_energy_binning, $
    im_energy_binning=im_energy_binning, $
    mpat_energy_binning = mpat_energy_binning, $


    _EXTRA=_extra, NOT_FOUND=NOT_found, DONE=done, fw_set_id=fw_set_id

;Make sure a2d_index_mask doesn't cause an update without a real change

; acs 2007-11-22 fw_set_ids is there to make it faster (doesnt go
; several times in the same class)


;Kluge to fix energy binning problem.
;Because we haven't separated making the energy bins from the detector calibration  that
;is only needed to accumulate, we have to put in a dummy time if one isn't found
;There won't be an absolute time if obs_time_interval or filename aren't set
;If exist(energy_band) or exist(eb_index) or exist(sp_energy_binning) then begin
;   need_dummy_time = (SIZE(Self->Get(/filename),/tname) ne 'STRING') and $
;   (total(abs(self->Get(/obs_time_interval))) eq 0)
;
;   ;If NEED_DUMMY_TIME then $
;   ; Self->Framework::Set, obs_time_interval = '20-feb-2002 11:'+['06','07']
;   If exist(eb_index ) then self->framework::set, eb_index=eb_index>0, /no_update
;   If exist(energy_band) and not exist( sp_energy_binning) then $
;     sp_energy_binning = energy_band
;
; 8-aug-03, ras - reconciling eb_index, sp_energy_binning, energy_band
; Since sp_energy_binning completely controls energy binning,
; just as sp_energy_binning does for spectrum.  This is all we need to set
; energy_band and eb_index are control parameters of the calib_eventlist
; and are only used to select the desired energy_band
;

if exist( im_energy_binning) or exist(sp_energy_binning) or exist(mpat_energy_binning) then begin
	;Set the energy binning in binned eventlist
	if exist(im_energy_binning) then mpat_energy_binning = im_energy_binning
	if exist(sp_energy_binning) then mpat_energy_binning = sp_energy_binning
	Self->framework::Set, mpat_energy_binning = mpat_energy_binning
	;Has mpat_energy_binning changed?
	this_need_update = Self->Get(/NEED_UPDATE)
	;iF it has make sure binned eventlist has it, too
	obe = Self->Get(/obj, class='hsi_binned_eventlist')
	obe_need_update = obe->Need_update()
	obe->Set, im_energy_binning =  mpat_energy_binning, /no_update
	obe->Set, need_update = obe_need_update or this_need_update

	endif
;if n_elements(sp_energy_binning) gt 2 then help, /tra


If Keyword_set(_extra) then Self->framework::set, _extra=_extra,$
 NOT_FOUND=NOT_found, DONE=done, fw_set_id = fw_set_id




END
;---------------------------------------------------------------------------
;
;Newly added 4-mar-2005 in case the Ut_binned_eventlist is requested.
;After process it's been deleted and will only be created on request.
;Function HSI_MPAT::Get, $
;
;         FOUND=found, NOT_FOUND=NOT_found, DONE=done, $
;         _EXTRA=_extra
;
;;We need to run control2binning for ut_binned_eventlist to
;;be certain ut_binned_eventlist computed with the correct parameters
;;control2binning functions like a little object that doesn't store it's data products
;
;
;
;return, self->hsi_spectrogram::get( FOUND=found, NOT_FOUND=NOT_found, DONE=done, $
;         _extra=_extra,ut_binned_eventlist=ut_binned_eventlist )
;end
;----------------------------------------------------------
;---------------------------------------------------------------------------
;
FUNCTION HSI_mpat::GetData, $

            ;this_harmonic_index removed, ras, 20-apr-2007
            _EXTRA=_extra


@hsi_insert_catch
;Here we manage the energy binning
;
If keyword_set( _EXTRA) then Self->Set, _Extra = _extra

data=self->Framework::GetData( _Extra = _extra )
;
;IF N_Elements( this_det_index ) NE 0 THEN BEGIN
;
;
;
;    IF Ptr_Valid( data[ this_det_index ] ) THEN BEGIN
;        RETURN, *data[ this_det_index ]
;    ENDIF ELSE BEGIN
;        RETURN, -1
;    ENDELSE
;
;ENDIF ELSE BEGIN
;
;    RETURN, data
;
;ENDELSE
;
return, data
END
;
;
;

;-------------------------------------------------------------------------
PRO hsi_mpat::Process, $
             _EXTRA=_extra


If Keyword_set(_Extra) then Self->Set, _Extra=_extra

energy_bins = get_edges( Self->Get(/im_energy_binning), /edges_2)
nchan = (size(/dim, energy_bins))[1] ;dim is 2 x nchan

for ichan = 0, nchan-1 do begin
	Self->set, eb_index = ichan

	hsi_square_mpat, self->Get(/source), mpat, $
	cc_i, $;cc_i includes all you see below. It carries all the time bin dependence
;	scale = scale, (includes det_eff.rel[k]*time_bin[k]*flux_var[i]*livetime[i]*duration,
;		k index is for grid, i index is for the bin for that grid.
;	calib_str = replicate( {gridtran: 0.0, gridmod: 0.0, scale: 0.0}, n_elements(old))
;	calib_str.gridtran = old.gridtran
;	calib_str.gridmod = old.gridtran*old.modamp[0] ;in case of harmonics
;	calib_str.scale    = scale
	nbins = nbins, $
	/calib_str_def, $
	effarea=effarea_i,$
	/modpat_only, $
	no_modpat = ichan gt 1 ;don't change mpat after first creation, for speed
	if ichan eq 0 then begin ; define data structures
		calib_str = reproduce( cc_i, nchan)
		effarea =  rebin(effarea_i+fltarr(1), nchan) ;fltarr(nchan)
		endif else begin
		calib_str[0,ichan] = cc_i
		help, effarea_i
		effarea[ichan] = effarea_i
		endelse
	endfor
source = Self->Get(/source)
source->Set, need_update=0
self->Set, calib_str = calib_str, nbins=nbins, effarea=effarea
self->SetData,  reform( [[mpat], [mpat^2]], [(size(/dim, mpat))[0:1],2])


END
function hsi_mpat::profile, img, idet=idet,  photon=photon, _extra=_extra
;+
;Name: HSI_MPAT::PROFILE
;
;Purpose: This function computes the expected counts profile from img#mpat where
;	mpat is the count bin x pixel modulation pattern matrix
;
;Inputs:
;	Img: Image. May be 1D or 2D.  Full vector or only partial for one of included dets(grids)
;Optional Keywords:
;	IDET: Idet references only the detectors selected by det_index_mask, not all 9
;		Use bins for this detector(grid), if not used then all bins are used
;	EB_index: Energy bin index corresponding to energy channels
;-
Self->Set, _extra=_extra
;if Self->need_update() then Self->process
eb_index = Self->Get(/eb_index)

info     = Self->Get(/info, /this_class_only)
nbins     = info.nbins
sel      = exist( idet) ? nbins[idet:idet+1] : [0, last_item(nbins)]
;mpt      = (*self.data)[sel[0]:sel[1]-1,*,0] ;mpt is mpt and mpt^2 concatenated
;Now we check the input and extract the image of interest and make sure the image is in 1D shape
im_dim   = (size(/dim, *self.data))[1]
this_img = reform(/over, img, im_dim, n_elements(img)/im_dim)
this_img = n_dimensions( this_img ) eq 2 ? this_img[*, eb_index] : this_img
cc       = info.calib_str[ sel[0]:sel[1]-1, eb_index]
effarea  = info.effarea[ eb_index]
;profile  = ( cc.gridtran * total(this_img) + ( mpt # this_img ) * cc.gridmod )
profile  = cc.scale * ( cc.gridtran * total(this_img) + ( (*self.data)[sel[0]:sel[1]-1,*,0] # this_img ) * cc.gridmod )
profile  = keyword_set(photon) ?  profile * effarea : profile
;       Qx = fltarr(npix*npix,n_det)
;        for k=0,n_det-1 do begin
;            for j=0,el_n-1 do begin
;                Qx[*,k] += Q[i,j,k] * xim[*,j]  ;;; Q is drm * brm,Q = fltarr(ph_n,el_n,n_det)
;                ;Qx gives the contribution in a pixel from each energy band, need to use upper matrix only
;            endfor
;            HQx[ index[k]:index[k+1]-1 ,i] = H_eps[ index[k]:index[k+1]-1 ,*,i] # Qx[*,k]


;cc.gridtran*total(im)*cc.scale+((mpt#im[*])*cc.gridmod)*cc.scale

return, profile
end
;;;;;;
function hsi_mpat::need_update
return, Self->Get(/need_update) or Self.source[0]->get(/need_update) or $
 (Self->Get(/obj, class='hsi_calib_eventlist'))->Get(/need_update)
end
;;;;;
function hsi_mpat::bproj, prf, $
	idet=idet,  $
	squared=squared, $
	shape=shape, $
	noscale = noscale, $
	_extra=_extra
;+
;Name: HSI_MPAT::BPROJ
;
;Purpose: This function computes the backprojection for a count profile, ie, prf#mpat where
;	mpat is the count bin x pixel modulation pattern matrix
;
;Inputs:
;	Prf: Counts or Expected counts. May be 1D or 2D.  Full vector or only partial for one of included dets(grids)
;Optional Keywords:
;	IDET: Idet references only the detectors selected by det_index_mask, not all 9
;		Use bins for this detector(grid), if not used then all bins are used and Prf must agree in length in the first dimension
;	EB_index: Energy bin index corresponding to energy channels
;	SQUARED: If set, use the equivalent of the MPAT^2, used for stopping criterion in EM, i.e. prf#mpat^2
;-


;
;  ;;; 3. H^T # z
;  ;;; y indica solo la struttura di z (che deve essere la stessa)
;  HT_z = z # sqm
;;	all_1 = fltarr(n_elements(l)) + 1.
;;	ht_one = all_1#sqm
;
;  ;;; 4. new estimate
;  factor = (f_div( HT_z , ht_one ))^alpha
;  x = x * factor
; p1s = counts * cc.scale
; bp1 = total(p1s*cc.gridtran)+ (p1s*cc.gridmod)#mpt ;sqm


;  ;;; z sono fotoni -observed counts array, from each energy
;    zex = f_div(y_eps,HQx)
;
;
;    ;;; HT z sono fotoni
;    U = fltarr(npix*npix,ph_n,n_det) ;because forward prob is slightly different for each detector
;    for i=0,ct_n-1 do begin
;        for k=0,n_det-1 do begin ;backprojection
;            U[*,i,k] = zex[index[k]:index[k+1]-1,i] # H_eps[ index[k]:index[k+1]-1 ,*,i]
;        endfor
;    endforSelf->Set, _extra=_extra
;
Self->Set, _extra=_extra
;if Self->need_update() then Self->process

eb_index = Self->Get(/eb_index)
isq      = keyword_set( squared ) ;need the squared modpat for stopping term
;mpt      = *self.data ;mpt is mpt and mpt^2 concatanated
info     = Self->Get(/info, /this_class_only)
nbins    = info.nbins
sz_prf   = size(/struct, prf)
isdet    = exist( idet )
sel      = isdet ? nbins[idet:idet+1] : [0, last_item(nbins)]
noscale  = keyword_set( noscale ) ; if set don't mult by cc.scale
;We reference these directly because otherwise they take a suprising amount of time to create over and over
;mpt      = (*self.data)[sel[0]:sel[1]-1,*,0] ;mpt is mpt and mpt^2 concatanated
;mpt2     = (*self.data)[sel[0]:sel[1]-1,*,1] ;
cc       = info.calib_str[ sel[0]:sel[1]-1, eb_index]
effarea  = info.effarea[ eb_index]
p1s      = sz_prf.n_dimensions eq 2 ? prf[*, eb_index] : prf
if noscale then p1s = sz_prf.dimensions[0] eq (sel[1]-sel[0]) ? p1s  : p1s[ sel[0]:sel[1]-1 ] 	else $
                p1s = ( sz_prf.dimensions[0] eq (sel[1]-sel[0]) ? p1s  : p1s[ sel[0]:sel[1]-1] ) * cc.scale^(1+isq)
; This way is slower
;bproj    = ~isq ? total( p1s*cc.gridtran) + (p1s*cc.gridmod)#mpt : $
;                  total( p1s*cc.gridtran^2) + (p1s*2*cc.gridmod*cc.gridtran)#mpt $
;                   + (p1s*cc.gridmod^2)#mpt2
bproj    = ~isq ? total( p1s*cc.gridtran) + (p1s*cc.gridmod)#(*self.data)[sel[0]:sel[1]-1,*,0] : $
                  total( p1s*cc.gridtran^2) + (p1s*2*cc.gridmod*cc.gridtran)#(*self.data)[sel[0]:sel[1]-1,*,0] $
                   + (p1s*cc.gridmod^2)#(*self.data)[sel[0]:sel[1]-1,*,1]


if exist(shape) then begin
	shapesize = n_elements(shape) le 2 ? shape : size(/dimension, shape)
	bproj = reform(/over, bproj, shapesize)
	endif
;profile = (  cc.gridtran * total(img) + ( mpt # img[*] ) * cc.gridmod ) * cc.scale

return, bproj
end

PRO Hsi_Mpat__Define

dummy = {Hsi_mpat, $

         INHERITS Framework }
END

;---------------------------------------------------------------------------
; End of 'hsi_modul_pattern__define.pro'.
;---------------------------------------------------------------------------