Download

package OMFastAnalyse;
use strict;
use English;
use Carp;
use vars
  qw(@ISA $VERSION $name $author $date $version @instList $numberOfStreams);
@ISA = qw(Module);
use ModuleResources;
use Astro::FITS::CFITSIO qw(TRUE FALSE);
#use ModuleUtil;

# Declare identity, version, author, date, etc.
$name    = __PACKAGE__;
$VERSION = '1.33';
$version = $VERSION;
$author="Duncan John Fyfe,Ian Stewart,Duncan Law-Green, Jose V Perea, Laura Tomas";
$date="2021-02-17";

#
# ChangeLog
# =========
#
# Version 1.33 - 2021-02-17 (JVP)
# ------------
#
# + Time binning for the timeseries set to the gretest value between
#   10 sec and 2*SAMPTIME/1024 from 'OM FAST MODE OSW IMAGE'
#
# Version 1.32 - 2020-12-14 (JVP)
# ------------
#
# + Ignoring omdetect, omatt and per-source processing for Grism filters
#
# Version 1.31 - 2020-02-21 (JVP)
# ------------
#
# + Typo bug fixes: readASCIIfile( $file ) -> readASCIIFile( name=> $file ) [ be aware of f -> F]
#                   Remove duplicated variable declaration ( remove duplicated my)
#
# Version 1.30 - 2020-01-14 (JVP)
# ------------
#
# + Do not exit from the module if no sources found in the first Fast exposure. Next in the foreach loop instead
#   Eventually, source level products could be created for the rest of windows if there are more windows.
#
# Version 1.29 - 2015-06-23 (LTC)
# ------------
#
# + implot (OM Fast Sky Images):
#   Create PNG product
#
#
# Version 1.28 - 2014-10-22 (JVP)
# ------------
#
# + omlcbuild (OM Light Curves):
#   Including an auxiliary imaging window data used for determining the background level around the source 
#
# Version 1.27 - 2009-05-18 (DLG)
# ------------
#
# + set omdetect value of nsigma to 3 to avoid spurious detections in fast mode
#   VY email 18 May 2009
#
# Version 1.26 - 2006-09-20 (DJF)
# ------------
#
# + Record use of exposures for source specific product creation in the database.
#
# Version 1.25 - 2006-09-05 (DJF)
# ------------
#
# + change Fast mode source list content from 'OM OSW SOURCE LIST' to 'OM OSW FAST SOURCE LIST'
# + Change the omdetect region file to a product.
#
# Version 1.24 - 2005-10-27 (IMS)
# ------------
#
# + Parameter --tests='chi' (default was 'both') added to lcplot call. This at advice of ACS 21 Oct.
#
# Version 1.23 - 2004-03-17 (DH)
# ------------
#
# + removed xcolumn and ycolumn parameters from evselects 
#   which were not creating imagesets
#
# Version 1.22 - 2002-10-25 (DH)
# ------------
#
# + omprep, change modeset=3 to modeset=1
#   got it wrong, 3 is for tracking
#
# Version 1.21 - 2002-10-25 (DH)
# ------------
#
# + omprep, change modeset=1 to modeset=3
# + omdetect, remove parameters: boxscale ,  outputregionfile  
#
# Version 1.20 - 2002-04-18 (DH)
# ------------
#
# + Change omregion parameters as per Simon Rosen email dated 18/4/02.
#
# Version 1.19 - 2001-12-05 (DFJ)
# ------------
#
# + Corrected content string for Fast mode images
#
# Version 1.18 - 2001-12-04 (DFJ)
# ------------
#
# + Added 'srcradius => -12' parameter to omregion.
#   The '-' causes a fixed radius to be assumed.  This overcomes a 
#   bug in calculating fluxes with adaptive radii
#
# Version 1.17 - 2001-11-30 (DFJ)
# ------------
#
# + Removed developers comments 
#
# Version 1.16 - 2001-11-27 (DFJ)
# ------------
#
# + Corrected condition on find of FAE files.  Was required=>true, now required=> false.
#
# Version 1.15 - 2001-11-26 (DFJ)
# ------------
#
# + Removed dependence on OMExpAnalyse generated flag.  OMFastAnalyse
#   does it's own selection of Fast Mode event lists. 
#
# Version 1.14 - 2001-11-07 (DFJ)
# ------------
#
# + Change evaluate rules to ignore if OMgetFlat is ignored.
#
# Version 1.13 - 2001-11-06 (DFJ)
# ------------
#
# + Avoid failing when there are no sources to generate source products from
#
# Version 1.12 - 2001-11-05 (DFJ)
# ------------
#
# + Fix for lcplot value VPS, was confusion over whether it should be /VPS or VPS
# + Added exception for when pstopdf conversion does not work
#
# Version 1.11 - 2001-11-05 (DFJ)
# ------------
#
# + Correct spelling mistake in parameter of lcbuild
#
# Version 1.10 - 2001-11-05 (DFJ)
# ------------
#
# + Fix for lcplot value VPS, was confusion over whether it should be /VPS or VPS
#
# Version 1.09 - 2001-11-05 (DFJ)
# ------------
#
# + Change lcplot device from default (PS) to VPS.  
#   This should change the orientation of the plot.
#
# Version 1.08 - 2001-11-05 (DFJ)
# ------------
#
# + Fix selection of osw.  Was failing to recognize 0 as a valid osw.
#
# Version 1.07 - 2001-10-02 (DFJ)
# ------------
#
# + Set a default osw of 9 when a valid value cannot be found.
#
# Version 1.06 - 2001-10-31 (DFJ)
# ------------
#
# + Change determination of osw to account for different sources (osw or osw_id)
#
# Version 1.05 - 2001-10-26 (DFJ)
# ------------
#
# + Fix second call to evselect.  Input is srcRegionFile not RegionFile
#
# Version 1.04 - 2001-10-26 (DFJ)
# ------------
#
# + Change region file format to ascii
#
# Version 1.03 - 2001-10-25 (DFJ)
# ------------
#
# + Fix omfastflat parameter fastimageset => fastimgset
#
# Version 1.02 - 2001-10-25 (DFJ)
# ------------
#
# + Fix evselect parameter outset => imageset
#
# Version 1.01 - 2001-10-25 (DFJ)
# ------------
#
# + Fix selection of Fast mode event files
#
# Version 1.00 - 2001-10-19 (DFJ)
# ------------
#
# + Initial Version (Using SAS 5.2.1)
#

# Declare list of instruments this module is interested in
@instList=qw(om);

# Number of streams
sub numberOfStreams {
	return numberOfExposures();
}

sub evaluateRules {

	# If upstream module has been nored, so if this one
	return ignore() if ignored(module => 'OMgetFlat'
		,instrument => thisInstrument
	    ,stream => 1
	);   

	return ignore() 
		if ignored(module => 'OMExpAnalyse'
			,instrument => thisInstrument
			,stream => thisStream
		);

	# Start conditions
	start()
		if complete(module => 'OMExpAnalyse'
			,instrument => thisInstrument
			,stream => thisStream
		);

}

sub performAction {
	info("Module version number: $version");
	# Misc. variables
	my (@extraopts);

	# Find exposure ID
	my $exp_id=exposureID(instrument => thisInstrument 
		,stream => thisStream
	);

	my $grism = 0;
	my $filter=getExposureProperty(instrument => 'om',
                                       exp_id => $exp_id,
                                       name => 'filter');
	info("DEBUG --- exp_id: $exp_id, filter: $filter");
	if ($filter =~ /grism/i){ $grism = 1}

	# Get FAE File
	my @FAEList=findFile(class => 'odf'
		,instrument => thisInstrument
		,exp_id => $exp_id
		,content => 'Fast mode events'
		,required => 'false'
	);

	if (scalar(@FAEList)) {
		my $tmp = scalar(@FAEList) ;
		info("$tmp Fast Mode Windows to process");
	}
	else {
		info("No Fast Mode Windows to process");
		return success();
	}

	# Find necessary files
	my $windowFile=findFile(class => 'odf'
		,instrument => thisInstrument
		,exp_id =>	$exp_id
		,content => 'Priority window data'
		,required => 'true'
	);

	# Find periodic housekeeping file
 	my $phkFile=findFile(class => 'odf'
		,instrument => thisInstrument
		,content => 'Periodic housekeeping'
		,required => 'true'
	);

	# Find non-periodic housekeeping file
	my $nphkFile=findFile(class => 'odf'
		,instrument => thisInstrument
		,content => 'Non-periodic housekeeping'
		,required => 'true'
	);
	
	# Find omprep tracking history
	my $prepThFile=findFile(class => 'intermediate'
		,instrument => thisInstrument
		,exp_id => $exp_id
		,content => 'omprep tracking history'
	);

		# Flat created by omflatgen
	my $flatGenFile=findFile(class => 'intermediate'
		,instrument => thisInstrument,
		,content => 'OM Flatfield'
	);


	# Loop over each image (one per OSW)
	foreach my $FAEFile (@FAEList) {

		my $imInfo; 
		
		$imInfo=fileInfo(class => 'odf'
			,name => $FAEFile
			);					
		my $osw = $$imInfo{'osw'};
		$osw = $$imInfo{'osw_id'} unless (defined($osw));
		$osw = 9 unless (defined($osw));
		info("Processing OSW $osw");

		# Create a raw image
		# omprep , evselect , omfastshift , omfastflat , omdetect , ommatt (source info)
		# (foreach source: omregion , evselect ,evselect , omlcbuild , lcplot , ps2pdf

		# Run omprep on the Event List 

		my $rawImFile=newFile(class => 'intermediate'
			,instrument => thisInstrument
			,exp_id => $exp_id
			,osw_id => $osw
			,content => 'Fast Mode Raw Image'
		);

		my $EvListFile=newFile(class => 'intermediate'
			,instrument => thisInstrument
			,exp_id => $exp_id
			,osw_id => $osw
			,content => 'PPS Event List File'
		);

		info("Start omprep");
		doCommand('omprep'
			,set => $FAEFile
			,nphset => $nphkFile
			,pehset => $phkFile 
			,wdxset => $windowFile
			,outset => $EvListFile
			,modeset => '1'
		) or return exception();

		# evselect an image
		
		info("Start evselect(1)");
		doCommand('evselect'
			,table => $EvListFile
			,xcolumn => 'RAWX'
			,ycolumn => 'RAWY'
			,withimageset => 'true'
			,imageset => $rawImFile
		) or return exception();

		# omfastshift
		info("Start omfastshift");
		doCommand('omfastshift'
			,nphset => $nphkFile
			,thxset => $prepThFile
			,set => $EvListFile
		) or return exception();

		# omfastflat

		my $FastImgFile=newFile(class => 'product'
			,instrument => thisInstrument
			,exp_id => $exp_id
			,osw_id => $osw
			,content => 'OM FAST MODE OSW IMAGE'
		);

		my $oswFlatFile=newFile(class => 'intermediate'
			,instrument => thisInstrument
			,exp_id => $exp_id
			,osw_id => $osw
			,content => 'OM OSW Flat Field'
		);
		
		info("Start omfastflat");
		doCommand('omfastflat'
			,slewflatset => $flatGenFile
			,oswflatset => $oswFlatFile
			,fastimgset => $FastImgFile
			,set => $EvListFile
		) or return exception();

		# omdetect

		my ($RegionFile, $sourceListFile, $skyImage);
		if ( $grism ){
		    info("Ignoring omdetect for Grism exposures, exposure: $exp_id, filter: $filter");
		} else {

		$RegionFile=newFile(class => 'product'
			,instrument => thisInstrument
			,exp_id => $exp_id
			,osw_id => $osw
			,content => 'OM OSW FAST REGION FILE'
			,format => 'ASCII'
		);

		$sourceListFile=newFile(class => 'product'
			,instrument => thisInstrument
			,exp_id => $exp_id
			,osw_id => $osw
			,content => 'OM OSW FAST SOURCE LIST'
		);

		info("Start omdetect");
		doCommand('omdetect'
			,set => $FastImgFile
			,nsigma => 3
			,regionfile => $RegionFile
			,outset => $sourceListFile
		) or return exception();
		}

		if ( !fileExists( file => $sourceListFile)){
		    info("Ignoring omatt becasue there is no source list."); }
		else {
		# Use omatt to make the sky image.
	
		$skyImage=newFile(class => 'product'
			,instrument => thisInstrument
			,exp_id => $exp_id
			,osw_id => $osw
			,content => 'OM FAST MODE OSW SKY IMAGE'
		);

		info("Start omatt");
		doCommand('omatt'
			,set => $FastImgFile
			,sourcelistset => $sourceListFile
			,ppsoswset => $skyImage
			,usecat => 'F'
		) or return exception();
		}

               #LTC start 
               #Call to implot to produce graphic images for OM Fast mode windows

               if (fileExists( file => $skyImage))
               {


	          info("Sky image found ");

                  my $pngskyImage = newFile(
			   class => 'product'
			   , instrument => thisInstrument
			   , exp_id => $exp_id
			   , osw_id => $osw
			   , content => 'OM FAST MODE OSW SKY IMAGE'
			   , format => 'PNG'
		   );

                  if (fileExists( file => $sourceListFile))

                  {

                     doCommand(
			   'implot', set => $skyImage
			   , device => "$pngskyImage/PNG"
			   , withsrclisttab => 'yes'
                           , srclisttab => $sourceListFile.':SRCLIST'
                           , zscaletype=>'sqrt'
                           , colourmap=>'7'
		     ) or return exception();

                                 
	            info("implot has created PNG with src list");


		 }

                 else


                 {

                     doCommand(
			   'implot', set => $skyImage
			   , device => "$pngskyImage/PNG"
			   , withsrclisttab => 'no'
                           , zscaletype=>'sqrt'
                           , colourmap=>'7'
		     ) or return exception();

                                 
	            info("implot has created PNG wothout src list");

                 

		 }

              }
            else
            {
                info("Sky image not found");

               

            }

            
            #if (fileExists( file => $pngskyImage))

            #{
            #info("DEBUG: PNG file $pngskyImage created");
            #}
            #else
            #{
            #info("PNG file does not exist"); 
            # return exception();        
	    #}

            #LTC end


	if ( $grism ){
	     info("Ignoring per-source processing for Grism exposures, exposure: $exp_id, filter: $filter");
	     return success();
	}

		# per source processing


		my @RegionFile= readASCIIFile(name => $RegionFile);
		my $regionSourceCount = scalar(@RegionFile);
		$regionSourceCount = $regionSourceCount - 1;  # Exclude format line in region file

		my $sourceListCount = numberFITSRows(file => $sourceListFile, extension => 'SRCLIST');
		
		info("Per-Source Processing");
		info("Number of sources from region file $RegionFile: $regionSourceCount");
		info("Sources in sourcelist $sourceListFile         : $sourceListCount");

		if (!$regionSourceCount || !$sourceListCount) {

			info("One of the source counts is zero.  No source level products will be produced for $exp_id.");
			#return success();
			next;

		}
		
		for(my $source = 1;$source < ($regionSourceCount+1);$source++) {

			info("Source Number : $source");
			setExposureProperty( instrument => thisInstrument , exp_id => $exp_id , name => 'ssp' , value => 0+TRUE );

			my $bkgRegionFile=newFile(class => 'intermediate'
				,instrument => thisInstrument
				,exp_id => $exp_id
				,osw_id => $osw
				,src_num => $source
				,content => "OM Background Region File"
			);

			my $bkgRatesFile=newFile(class => 'intermediate'
				,instrument => thisInstrument
				,exp_id => $exp_id
				,osw_id => $osw
				,src_num => $source
				,content => "OM Background Rates File"
			);

			my $srcRegionFile=newFile(class => 'intermediate'
				,instrument => thisInstrument
				,exp_id => $exp_id
				,osw_id => $osw
				,src_num => $source
				,content => "OM Source Region File"
			);

			my $srcRatesFile=newFile(class => 'intermediate'
				,instrument => thisInstrument
				,exp_id => $exp_id
				,osw_id => $osw
				,src_num => $source
				,content => "OM Source Rates File"
			);

			#omregion
			info("Start omregion");
			doCommand('omregion'
				,set => $sourceListFile
				,srcnumber => "$source"
				,nfwhm => '3'
				,bkginner => '1.2'
				,bkgouter => '2.5'
				,bkgfile => $bkgRegionFile
				,srcfile => $srcRegionFile
				,srcradius => -6
			) or return exception();

			# evselect (source rate)
			my $FastImgFile=findFile(class => 'product'
                                ,instrument => thisInstrument
                                ,exp_id => $exp_id
                                ,osw_id => $osw
                                ,content => 'OM FAST MODE OSW IMAGE'
                                );
                        my $tbin = 10;
                        if ( fileExists( file => $FastImgFile )){
                             if ( hasFITSKeyword(file => $FastImgFile, extension => 'PRIMARY', keyword => 'SAMPTIME') ){
                                  my $samptime = readFITSKeyword(file => $FastImgFile, extension => 'PRIMARY', keyword => 'SAMPTIME');
                                  info("Sampling timebinsize of $FastImgFile = $samptime/1024 sec");
                                  if (2*$samptime/1024 > 10) { $tbin = 2*$samptime/1024; }
                             }
                        }

			info("Start evselect (2)");
			info("timebinsize for evselect = $tbin");
			doCommand('evselect'
				,table => $EvListFile
				,expression=> "((WIN_FLAG .eq. 0) .and. (region($srcRegionFile, CORR_X, CORR_Y)))"
				#,xcolumn => 'CORR_X'
				#,ycolumn => 'CORR_Y'
				,maketimecolumn => 'true' 
				,withrateset => 'T'
				,timebinsize => $tbin
				,rateset =>	$srcRatesFile
			) or return exception();

			# evselect (background rate)

			info("Start evselect (3)");
			doCommand('evselect'
				,table => $EvListFile
				,expression=> "((WIN_FLAG .eq. 0) .and. (region($bkgRegionFile, CORR_X, CORR_Y)))"
				#,xcolumn => 'CORR_X'
				#,ycolumn => 'CORR_Y'
				,maketimecolumn => 'true' 
				,withrateset => 'T'
				,timebinsize => $tbin
				,rateset =>	$bkgRatesFile
			) or return exception();

			# omlcbuild

			my $lightCurveFile=newFile(class => 'product'
				,instrument => thisInstrument
				,exp_id => $exp_id
				,osw_id => $osw
				,src_num => $source
				,content => "OM OSW SOURCE TIMESERIES"
			);
		
			info("Start omlcbuild");
						
# Searching for the auxiliary imaging window data file used for determining the background level around the source
			# Get IMI Files : 00IMI.FIT, 01IMI.FIT or 02IMI.FIT
			my @IMIList=findFile(class => 'odf'
					     ,instrument => thisInstrument
					     ,exp_id => $exp_id
					     ,content => 'Imaging mode image'
					     ,required => 'false'
					     );
			# No bkg from image if 'Imaging mode image' not found
			my $bkgFromImage = 'no';
			my $bkgImageFile = '';
				
			foreach my $IMIFile ( sort @IMIList ){		# Be careful with this 'sort' 
			
				if (fileExists(file => $IMIFile)){
				info("Imaging mode image for determining background $exp_id:$IMIFile");
				$bkgImageFile = $IMIFile;
				$bkgFromImage = 'yes';
				last; # halt the foreach loop if 'Imaging mode image' is found
				}
			}
			if ( $bkgFromImage =~ /no/){		# No bkg from image if 'Imaging mode image' not found
				info("No Imaging mode image (IMI) found for this Fast exposure: $FAEFile");
			}


			doCommand('omlcbuild'
				,srcregionset => $srcRegionFile
				,bkgregionset => $bkgRegionFile
				,srcrateset => $srcRatesFile
				,bkgrateset => $bkgRatesFile
				,sourcelistset => $sourceListFile
				,wdxset => $windowFile
				,outset => $lightCurveFile
				,bkgfromimage => $bkgFromImage
				,imageset => $bkgImageFile
			) or return exception();
##DEBUG
			info("DEBUG ++ omlcbuild complete, generated TIMESR file: $lightCurveFile");
			if (-e $lightCurveFile) {
			    info("DEBUG ++ TIMESR file $lightCurveFile exists!");
			} else {
			    info("DEBUG ++ WARNING: TIMESR file $lightCurveFile not found!");
			}
##DEBUG

			#lcplot

			my $lcPlotFile=newFile(class => 'intermediate'
				,instrument => thisInstrument
				,exp_id => $exp_id
				,osw_id => $osw
				,src_num => $source
				,format => 'PS'
				,content => "OM OSW SOURCE TIMESERIES"
			);

			info("Start lcplot");
			doCommand('lcplot'
				,set => $lightCurveFile
				,binsize => '1'
				,plotdevice=>'/VPS'
				,plotfile => $lcPlotFile
                                ,tests => 'chi'
			) or return exception();

			my $pdfLCplot=newFile(class => 'product'
				,instrument => thisInstrument
				,exp_id => $exp_id
				,osw_id => $osw
				,src_num => $source
				,'format' => 'PDF'
				,content => "OM OSW SOURCE TIMESERIES"
			) or return exception();

			info("Start PStoPDF");
			PStoPDF(source => $lcPlotFile
				,destination => $pdfLCplot
			) or return exception();

		}

	}
			
	# Everything was OK
	return success();
}

1;