% -*- mode: SLang; mode: fold -*-
%; Time-stamp: <2013-04-25 22:19:24 dph>
%; Directory:  ~dph/h3/CXC/TG/Scripts_devel/sl_add_tg/Sl/
%; File:       add_pha_counts.sl
%; Author:     David P. Huenemoerder <dph@space.mit.edu>
%; Orig. version: 2006.06.12
%;========================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This file is part of add_pha (sum Chandra PHA spectra)
% Copyright (C) 2007 Massachusetts Institute of Technology 
%
% This software was developed by the MIT Kavli Institute for
% Astrophysics and Space Research under contract SV3-73016 from the
% Smithsonian Institution.
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either  version 2 of the
% License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% versioning...
%
private variable _version = [0, 5, 0 ]   ; % major, minor, patch
private variable add_pha_counts_version = sum( _version * [ 10000, 100, 1 ] ) ; 
private variable add_pha_counts_version_string = sprintf("%d.%d.%d", 
                                     _version[0],
                                     _version[1],
                                     _version[2] ) ; 

% version 0.5.0 - allow summed exposure or mean via top-level qualifier.
%
%
% version 0.4.2 - change copyright from MIT-integrated to MIT-standalone (GPL)
%
% version 0.4.1 - avoid cfitsio bug on Char_Type vs UChar_Type in null prim.
%
% version 0.4.0 - fix bug when responses defined in header.
%
% version:  0.3.0
%           2006.07.03 removed a few slang-2 features in order to be
%                      slang-1/ciao-3.3 compatible.  Features used
%                      were  'foreach name ( ... )'
%                            'i=1; "string$i"$'
%                            'wherefirst()'
%                      Instead of #if conditionals, just re-wrote.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% purpose:  add counts from list of pha spectra, write output to 
%           Type I pha file.  Compute mean of EXPOSURES and store in
%           header.  Compute exposure-weighted BACKSCAL sum.
%
%  Note:  Summed COUNTS and mean EXPOSURE are consistent with dmarfadd
%         for summing orders or summing observations.   Weighted sum
%         of BACKSCAL is likewise consistent and can be used w/ summed
%         background done analogously to ARFS by dmarfadd.
%
%  Output files will have minimal OGIP pha header.
%  Output will have minimal CXC pha structure, and include bin_lo,
%  bin_hi columns, if grating obs.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Ignore_PHA_Response_Keywords = 1 ; %prevent loading of responses.


private variable USE_MEAN_EXPOSURE = 1 ; 



%%%%%%%%%%%%%%%%%%%%%
%
% variable keywords, to be determined by data or params.
% This needs to be in program scope so it can be known to the header
% callback function.
%
private variable dyn_keys = struct {
%
    object,   % from 1st infiles header 
    mission,  % from 1st
    telescop, % from 1st 
    instrume, % from 1st 
    detnam,   % from 1st 
    grating,  % from 1st infiles header 
    tg_part,  % from 1st hist info 
    tg_m,     % from 1st hist info 
    exposure, % computed 
    backscal, % computed 
    detchans, % computed 
    comments  % argument 
};


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define the data-dependent keywords;
% fref is a file reference for copying static stuff from the 
%  first input file's header (assuming it's relevant to all files)
%
% oops - need data and comments...
%
private define set_data_keywords( fref,  h, comments ) 
{
%
% Compute mean exposure,  exposure-weighted backscal.
%
    variable i ;
    variable var_keys = @dyn_keys ; 

    var_keys.backscal = 0.0 ;
    var_keys.exposure = 0.0 ; 

    for (i=0; i<length(h); i++ )
    {
	var_keys.backscal += get_data_backscale( h[i] )[0] * get_data_exposure( h[i] ) ;
	var_keys.exposure += get_data_exposure( h[i] ) ;
    }

    if ( USE_MEAN_EXPOSURE ) var_keys.exposure /= length( h ) ;  % mean exposure
    var_keys.backscal /= var_keys.exposure ;   % exposure weighted backscale sum

%
% misc info, taken only from the first infile:
%
    var_keys.object  = get_data_info( h[0] ).target ;
    var_keys.tg_part = get_data_info( h[0] ).part ;
%    var_keys.tg_m    = abs( get_data_info( h[0] ).order ) ; 
    var_keys.tg_m    =  get_data_info( h[0] ).order ; 

%
% get miscellaneous key values from input file header:
%
    var_keys.mission   =  fits_read_key( fref, "MISSION" ) ; 
    var_keys.telescop  =  fits_read_key( fref, "TELESCOP" ) ; 
    var_keys.instrume  =  fits_read_key( fref, "INSTRUME" ) ; 
    var_keys.detnam    =  fits_read_key( fref, "DETNAM" ) ; 
    var_keys.grating   =  fits_read_key( fref, "GRATING" ) ;

    var_keys.detchans  =  length( get_data_counts(h[0]).value ) ; 

%
% copy the comments:
%
    var_keys.comments  = [@comments, sprintf( "Summed %d spectra",  length(h) )] ;

    return ( var_keys ) ; 
}


%
%  process the files, and fill in the var_keys structure:
%
private define _add_pha_counts( infiles, outfile, comments )
{
    variable i ;
    variable h = Integer_Type[ 0 ] ;
    variable data, var_keys ; 
	
% Load the files.

    for( i=0; i<length(infiles); i++)
	h =[h, load_data( infiles[i] ) ] ;

% combine datasets and obtain summed counts.
% If any grids mis-match, there will be an error.
%
    variable g = combine_datasets( h ) ;
    if ( -1 != g )
    {
	data = get_combined( g, &get_data_counts ) ;
    }
    else
    {
	delete_data( h ) ; 
	error("%% List of infiles must specify commensurate spectra.");
    }

%
% set header values:
%    
    var_keys = set_data_keywords( infiles[0], h, comments ) ; 

%
% and clean-up
%
    delete_data( h ) ; % clean up
%
    return( data, var_keys ) ; 
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% write Type I PHA file, given an appropriate structure.
% Need a cxc header
%

private variable pha_record_tg = struct { channel, bin_lo, bin_hi, counts, } ;
%private variable pha_record_im = struct { channel, pi, counts } ;
private variable pha_record_im = struct { channel, counts } ;

%
% define static values:
%
private define write_static_keywords( fp )
{
    fits_write_date (fp);
    fits_update_key( fp, "HDUNAME",  "SPECTRUM" );
    fits_update_key( fp, "HDUCLASS", "OGIP" ); 
    fits_update_key( fp, "HDUCLAS1", "SPECTRUM" ) ; 
    fits_update_key( fp, "HDUCLAS2", "TOTAL", "TOTAL | NET | BDG" ) ; 
    fits_update_key( fp, "HDUCLAS3", "TYPE:I", ":I->single,col; :II->multiple,rows" ) ; 
    fits_update_key( fp, "HDUCLAS4", "COUNT", "COUNT | RATE" ) ; 
    fits_update_key( fp, "HDUVERS",  "1.2.1" ) ; 
    fits_update_key( fp, "CREATOR",  "add_pha_counts.sl v" + add_pha_counts_version_string ) ; 
    fits_update_logical( fp, "POISSERR", 1, NULL ) ;
    fits_update_key( fp, "CHANTYPE",  "PI" ) ; 
    fits_update_key( fp, "BACKFILE",  "NONE" ) ; 
    fits_update_key( fp, "CORRFILE",   "NONE" ) ; 
    fits_update_key( fp, "AREASCAL",   1.0 ) ; 
    fits_update_key( fp, "CORRSCAL",   1.0 ) ; 
    fits_update_key( fp, "RESPFILE",  "NONE" ) ; 
    fits_update_key( fp, "ANCRFILE",  "NONE" ); 
    fits_update_key( fp, "GROUPING",  0, "0 => ungrouped"  ); 
    fits_update_key( fp, "QUALITY",   0, "0 => all good" ) ; 
}


% write the NULL PRIMARY keywords which are known a priori:
%
private define write_static_null_hdr_keywords( fp )  
{
    fits_update_key( fp, "ORIGIN", "ASC", NULL ) ;
    fits_update_key( fp, "CREATOR",  "add_pha_counts.sl v" + add_pha_counts_version_string ) ; 
    fits_write_comment( fp, "" ); 
    fits_update_key( fp, "TITLE", "WARNING! Summed counts spectrum; use with caution." ) ; 
    fits_write_comment( fp, "" ); 
    fits_write_date (fp);
}

% write the NULL PRIMARY keywords which are also in the interesting
% extension, and which have been copied from the first input pha file:
%
% revise for slang-1 compatibility:
private define write_var_null_hdr_keywords( fp, var_keys )
{
    foreach ( ["mission", "telescop", "instrume"] )
    {
	variable field = () ; 
	fits_update_key( fp, field, get_struct_field( var_keys, field )	) ;
    }
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
private define write_keywords_callback( fp, var_keys )
{
%
% write the a priori known keywords:
%
    write_static_keywords( fp ) ;

% we now need to write the computed keywords in the var_keys struct.
% we can loop over all struct keys.  But some should also have custom
% comment fields for some keys.  Store in an assoc. array.

    variable cmts = Assoc_Type[ String_Type ] ;

    cmts[ "tg_part"] = "1=>HEG, 2=>MEG, 3=>LEG; 0=>zero order or direct" ; 
%    cmts[ "tg_m" ] = "absolute value of grating order" ;
    cmts[ "tg_m" ] = "grating order" ;
    if ( USE_MEAN_EXPOSURE ) cmts[ "exposure" ] = "Mean exposure;  T = sum(t_i) / N" ;
    else  cmts[ "exposure" ] = "Summed exposure;  T = sum(t_i)" ;
    cmts[ "backscal" ] = "Exposure-weighted backscal; sum( t_i Back_i )/T";

% now update each keyword in the structure:
%
    foreach ( get_struct_field_names( var_keys ) )
    {
	variable n = () ;
	variable v = get_struct_field( var_keys, n ) ;
	variable c = NULL ; 
	if ( assoc_key_exists( cmts, n ) ) c = cmts[ n ] ;

%
% "history" or "comment" fields are special cases.  They may be
% arrays, and they have a different update function:
%
	variable vv ; 

	if ( n == "history")
	{
	    fits_write_comment( fp, "" ) ; 
	    foreach ( v )
	    {
		vv = () ;
		fits_write_history( fp, vv ) ;
	    }
	}
	else if ( n == "comments")
	{
	    fits_write_comment( fp, "" ) ; 
	    foreach ( v )
	    {
		vv = () ;
		fits_write_comment( fp, vv ) ;
	    }
	}
	else  % update the named keyword, given the value and comment:
	{
	    fits_update_key( fp, n, v, c ) ;
	}
    }
}

%
% Write the binary table.  We conditionally reverse the arrays because
% ISIS uses wavelength order for gratings, but the OGIP spec is in
% ascending energy.
%
private define write_pha( file, data, var_keys )
{
    variable p;

    if ( var_keys.tg_part == 0 )
    {
	p = @pha_record_im ;
	p.counts = @data.value ;   % no reverse for imaging mode (no grid)
    }
    else
    {
	p = @pha_record_tg ;
	p.bin_lo = reverse( data.bin_lo ) ;
	p.bin_hi = reverse( data.bin_hi ) ;
	p.counts = reverse( data.value ) ; 
    }

    p.channel = [ 1:length(data.value) ] ;
%    p.pi = double( p.channel ) ;

    variable fp = fits_open_file( file, "c" ) ; 

% v0.4.1: cfitsio bug; change Char_Type to UChar_Type  ***   BUG ***
%    fits_create_image_hdu (fp, NULL, Char_Type, Int_Type[0]);
% ***
    fits_create_image_hdu (fp, NULL, UChar_Type, Int_Type[0]);
% ***
    write_static_null_hdr_keywords( fp )  ;
    write_var_null_hdr_keywords( fp, var_keys ) ; 
    fits_write_chksum (fp);

    fits_write_binary_table( fp, "SPECTRUM", p, &write_keywords_callback, var_keys ) ;

    if ( var_keys.tg_part > 0 )
    {
        variable fields = get_struct_field_names (p);
% revise for slang-1 compatibility
%   "foreach field ("  ->  "foreach ( ... { variable field = () ..."
%   wherefirst() -> where()[0]
%   "TUNIT$i"$ -> sprintf()
%
	foreach ( ["bin_lo", "bin_hi"] )
	  {
	      variable field = () ; 
	      variable i = (where( fields == field ))[0] + 1 ;
	      fits_update_key( fp, sprintf("TUNIT%d",i), "angstrom") ;
	  }
    }
    fits_write_chksum (fp);
    fits_close_file (fp);
}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
define add_pha_counts()    %  infiles, outfile[, comments] )
{
%        infiles = String_Type[]
%        outfile = String_Type,
%        comments = String_Type[]
%
%%%%%%%
%   Input files can have cfitsio virtual file syntax.  It is assumed
%   that all input files reference the same type of data (grid,
%   absolute order, grating type, resolution).  The output grid will
%   be determined by the first infile.  Other header attributes
%   (TG_PART, TG_M, OBJECT, OBS_ID), if present, will be determined
%   from the first input file.
%
%   No binning or grouping is supported; no ARF or RMF associations
%   are made; no background is processied.  Those can be done in an
%   analysis session, or with other tools.
%
%   Since this function sums counts, no error column will be written.
%   In accordance with this, the header will contain POISSERR = T
%
%%%%%%%%%%%%%%%%%%%%%
%
%% Parameters:
%
%   "infile" is a list of input file names, which may be cfitsio
%            virtual file specifications to select particular rows and
%            columns.
%
%   "outfile" is an arbitrary file name.
%
%   "comments" is optional array of strings to write as FITS COMMENT
%              records. 
%
%   Example:
%
%     infiles = [ "pha-1.fits[tg_part==2 && (tg_m==-1 || tg_m==1)]",
%                 "pha-2.fits[#row==9 || #row==10]",
% 	          "pha-3.fits" ];
%
%     add_pha_counts( infiles, "pha-123.fits", ["Capella", "MEG -1,+1"] );
%
%%%%%%%%
%   Notes on procedure:
%
%        Data will be read with the ISIS load_data() function, which
%        performs minimal checking and stores metadata.  Further
%        checking is done when data are collected via
%        combine_datasets() and get_combined(), which requires grids
%        to match.  No grid-matching is attempted.  This can be done
%        in an analysis session, such as with ISIS'
%        match_dataset_grids(), which will ensure that counts, arfs,
%        and rmfs match and are used consistently.
%     
%
    variable infiles, outfile, comments ;

    switch( _NARGS )
    {
	case 2:
	( infiles, outfile ) = () ;
	comments = [""]  ; 
    }
    {
	case 3:
	( infiles, outfile, comments ) = () ;
	comments = [ comments ];
    }
    {
	message("\n%% USAGE: add_pha_counts( infiles, outfile[,comments]	);");
	message("\n%% PURPOSE: add counts from 2 or more input	spectra.");
	message("\n%% EXAMPLE:\n" 
         +  "infiles = [ \"pha-1.fits[tg_part==2 && (tg_m==-1 || tg_m==1)]\",\n"
         +  "            \"pha-2.fits[#row==9 || #row==10]\",\n"
 	 +  "            \"pha-3.fits\" ];\n"
	 +  "add_pha_counts( infiles, \"pha-123.fits\", [\"Capella\", \"MEG -1,+1\"] );\n" 
         ) ;
	return ; 
    }

    USE_MEAN_EXPOSURE = 1 ;
    if ( qualifier_exists( "sum" ) )  USE_MEAN_EXPOSURE = 0 ;
    if ( qualifier_exists( "mean" ) ) USE_MEAN_EXPOSURE = 1 ;

    if ( length( infiles ) == 1 ) infiles = [ infiles ] ; %force to array.

    variable dat, kwd ; 

    (dat, kwd) = _add_pha_counts( infiles, outfile, comments ) ;

    write_pha( outfile, dat, kwd ) ; 
}

provide( "add_pha_counts" ) ;