#! /bin/sh # This is the LHEA perl script: /cvmfs/extras-fp7.egi.eu/extras/heasoft/hitomi/x86_64-unknown-linux-gnu-libc2.19-0/bin/sxsregext.pl # The purpose of this special block is to make this script work with # the user's local perl, regardless of where that perl is installed. # The variable LHEAPERL is set by the initialization script to # point to the local perl installation. #------------------------------------------------------------------------------- eval ' if [ "x$LHEAPERL" = x ]; then echo "Please run standard LHEA initialization before attempting to run /cvmfs/extras-fp7.egi.eu/extras/heasoft/hitomi/x86_64-unknown-linux-gnu-libc2.19-0/bin/sxsregext.pl." exit 3 elif [ "$LHEAPERL" = noperl ]; then echo "During LHEA initialization, no acceptable version of Perl was found." echo "Cannot execute script /cvmfs/extras-fp7.egi.eu/extras/heasoft/hitomi/x86_64-unknown-linux-gnu-libc2.19-0/bin/sxsregext.pl." exit 3 elif [ `$LHEAPERL -v < /dev/null 2> /dev/null | grep -ic "perl"` -eq 0 ]; then echo "LHEAPERL variable does not point to a usable perl." exit 3 else # Force Perl into 32-bit mode (to match the binaries) if necessary: if [ "x$HD_BUILD_ARCH_32_BIT" = xyes ]; then if [ `$LHEAPERL -V 2> /dev/null | grep -ic "USE_64_BIT"` -ne 0 ]; then VERSIONER_PERL_PREFER_32_BIT=yes export VERSIONER_PERL_PREFER_32_BIT fi fi exec $LHEAPERL -x $0 ${1+"$@"} fi ' if(0); # Do not delete anything above this comment from an installed LHEA script! #------------------------------------------------------------------------------- #!/usr/bin/perl # # File name: sxsregext.pl # Author: A. J. Sargent NASA GSFC # $Date: 2016/03/24 13:26:07 $ # Version: 0 # # Create a region file and run the extractor tool to create an SXS # image, lightcurve and spectrum file # # Tool Dependencies: # # Library Dependencies: # gen/lib/perl/ahlog # gen/lib/perl/ahapp # gen/lib/perl/ahgen # # Modification History: # # Set up use strict; use warnings; use ahlog ; use ahapp ; use ahgen qw (:ALL) ; use ahfilterlib ; # turn on AUTOFLUSH $| = 1; # For a 2D array, sort the second dimension, then the first sub sortxy { sort { $a->[1] <=> $b->[1] || $a->[0] <=> $b->[0] } @_ } ######################### # Variable Definitions ######################### my $infile =""; # Input event file my $regmode =""; # Region file type ([SKY],DET) my $region =""; # Region file name my $extract =0; # Extract image/lightcurve/spectra ([yes],no) my $outroot =""; # Root of output file names my $resolist =""; # List of ITYPE or ALL my $teldeffile =""; # Input teldef file (or CALDB) my $outexp =""; # Output exposure map my $ehkfile =""; # Input EHK file my $delta =0.0;# Range of off-axis angle (theta) in one map [arcmin] my $numphi =0; # Number of azimuth (phi) bins" my $badimgfile =""; # Input bad pixel image file (or NONE) my $pixgtifile =""; # Flickering pixel list (or NONE) my $instmap =""; # Instrument map (or CALDB) my $qefile =""; # QE file (or CALDB) my $contamifile=""; # Contamination file (or CALDB) my $obffile =""; # Optical blocking filter file for SXS (or CALDB) my $fwfile =""; # Filter wheel file for SXS (or CALDB) my $gvfile =""; # Gatevalve file for SXS (or CALDB) my $maskcalsrc =""; # Mask out calibration sources? ([yes]/no) my $fwtype =""; # Filter wheel type for SXS ([OPEN], FE55, BE, NE, or POLY) my $specmode =""; # Type of input energy ([MONO] or SPEC) my $specfile =""; # Input spectrum file my $specform =""; # Format of spectrum file ([FITS] or ASCII) my $telescop ="HITOMI"; my $instrume ="SXS"; my $outreg =""; # Output region file my $img =""; # Output image file my $lc =""; # Output lightcurve file my $pha =""; # Output spectrum file my @pixels = ( 0 ) x 35; my $sxsregexterror = 0; # sxsregext exit status ######################### # Main Code Block ######################### ahapp::startup () ; ahapp::begin_processing () ; $sxsregexterror = get_parameters () ; unless ( $sxsregexterror == 0 ) { ah_info "HIGH", ahapp::write_parameters () ; ahlog::ah_err "get_parameters" ; ahapp::end_processing($sxsregexterror); } $sxsregexterror = initialize () ; unless ( $sxsregexterror == 0 ) { ah_info "HIGH", ahapp::write_parameters () ; ahlog::ah_err "initialize" ; ahapp::end_processing($sxsregexterror); } $sxsregexterror = do_work () ; unless ( $sxsregexterror == 0 ) { ahlog::ah_err "do_work" ; ahapp::end_processing($sxsregexterror); } $sxsregexterror = finalize () ; unless ( $sxsregexterror == 0 ) { ahlog::ah_err "finalize" ; ahapp::end_processing($sxsregexterror); } # We're done. ahapp::end_processing($sxsregexterror); ######################### # Subroutines ######################### sub get_parameters { $infile = ahapp::query_parameter("infile"); $regmode = ahapp::query_parameter("regmode"); $region = ahapp::query_parameter("region"); $extract = ahapp::query_parameter("extract",1); $outroot = ahapp::query_parameter("outroot"); $resolist = ahapp::query_parameter("resolist"); $teldeffile = ahapp::query_parameter("teldeffile"); $outexp = ahapp::query_parameter("outexp"); $ehkfile = ahapp::query_parameter("ehkfile"); $delta = ahapp::query_parameter("delta"); $numphi = ahapp::query_parameter("numphi"); $badimgfile = ahapp::query_parameter("badimgfile"); $pixgtifile = ahapp::query_parameter("pixgtifile"); $instmap = ahapp::query_parameter("instmap"); $qefile = ahapp::query_parameter("qefile"); $contamifile= ahapp::query_parameter("contamifile"); $obffile = ahapp::query_parameter("obffile"); $fwfile = ahapp::query_parameter("fwfile"); $gvfile = ahapp::query_parameter("gvfile"); $maskcalsrc = ahapp::query_parameter("maskcalsrc"); $fwtype = ahapp::query_parameter("fwtype"); $specmode = ahapp::query_parameter("specmode"); $specfile = ahapp::query_parameter("specfile"); $specform = ahapp::query_parameter("specform"); return 0; } # ------------------------------------------------------------------------------ sub initialize { my $status = 0; # Input file checking if (isRequiredFileNotFound($infile)) { return 1; } if (isRequiredFileNotFound($region)) { return 1; } if ( uc $regmode eq "SKY" ) { if (isRequiredFileNotFound($ehkfile)) { return 1; } # Setup output region file name $outreg = $outroot . ".reg"; if (removeOutputFileIfClobbering($outreg,$ahapp::clobber) ) { return 1; } if (removeOutputFileIfClobbering($outexp,$ahapp::clobber)) { return 1; } } else { if ( ! $extract ) { ahlog::ah_err "Nothing to do."; return 1; } # Use the input region file for the extractor input $outreg = $region; } # Make sure we have a selection file if (isBadCALDBFileParameterValue($teldeffile,"teldeffile")) { return 1; } # Check for valid resolist parameter if ( uc $resolist eq "ALL" ) { $resolist = "0,1,2,3,4"; } # Verify that we have valid values for ITYPE my @splitresolist = split ",",$resolist; foreach my $itype ( @splitresolist ) { if ( $itype !~ /[0-4]/ ) { ahlog::ah_err "Not a valid ITYPE(=$itype), must be between 0 and 4"; return 1; } } # Set up output file names if ( $extract ) { $img = $outroot . ".img"; $lc = $outroot . ".lc"; $pha = $outroot . ".pha"; if (removeOutputFileIfClobbering($img,$ahapp::clobber) ) { return 1; } if (removeOutputFileIfClobbering($lc ,$ahapp::clobber) ) { return 1; } if (removeOutputFileIfClobbering($pha,$ahapp::clobber) ) { return 1; } } # Write all parameters to this script to the log file. ah_info "HIGH", ahapp::write_parameters () ; return 0; } # ------------------------------------------------------------------------------ sub do_work { my $status = 0; # If 'regmode' is SKY: # a. run ahexpmap # b. determine pixels to extract # c. run ahmkregion to create DET region # d. condense region file if ( uc $regmode eq "SKY" ) { # Run ahexpmap $status = ahgen::run_ftool("ahexpmap", "ehkfile=$ehkfile", "gtifile=$infile", "instrume=$instrume", "badimgfile=$badimgfile", "pixgtifile=$pixgtifile", "outfile=$outexp", "outmaptype=EXPOSURE", "delta=$delta", "numphi=$numphi", "stopsys=SKY", "instmap=$instmap", "qefile=$qefile", "contamifile=$contamifile", "vigfile=NONE", "obffile=$obffile", "fwfile=$fwfile", "gvfile=$gvfile", "maskcalsrc=$maskcalsrc", "fwtype=$fwtype", "specmode=$specmode", "specfile=$specfile", "specform=$specform", #"energy=$energy", #"evperchan=$evperchan", #"abund=$abund", #"cols=$cols", #"covfac=$covfac" ); if ( $status ) { ahlog::ah_err "ahexpmap failed."; return $status; } # Determine the number of rows (bins) in the output file my $naxis2 = ahgen::get_keyword($outexp,"OFFAXISHIST","NAXIS2"); if( ahgen::get_error_flag ) { ahlog::ah_err "NAXIS2 keyword not defined in $outexp"; return ahgen::get_error_flag; } my $ra_nom = ahgen::get_keyword($outexp,"OFFAXISHIST","RA_NOM"); if( ahgen::get_error_flag ) { ahlog::ah_err "RA_NOM keyword not defined in $outexp"; return ahgen::get_error_flag; } my $dec_nom = ahgen::get_keyword($outexp,"OFFAXISHIST","DEC_NOM"); if( ahgen::get_error_flag ) { ahlog::ah_err "DEC_NOM keyword not defined in $outexp"; return ahgen::get_error_flag; } # Process the output exposure map my @ranomxp = ahgen::read_column($outexp,"OFFAXISHIST","RANOMXP"); if($naxis2 != @ranomxp ) { ahlog::ah_err "Error reading RANOMXP column in $outexp"; return 1; } my @decnomxp = ahgen::read_column($outexp,"OFFAXISHIST","DECNOMXP"); if($naxis2 != @decnomxp ) { ahlog::ah_err "Error reading DECNOMXP column in $outexp"; return 1; } my @pa_nomxp= ahgen::read_column($outexp,"OFFAXISHIST","PA_NOMXP"); if($naxis2 != @pa_nomxp ) { ahlog::ah_err "Error reading PA_NOMXP column in $outexp"; return 1; } # Loop through each row in the output exposure map foreach my $bin ( 0 .. $naxis2 - 1 ) { # Get row element for each column my $ranomxp = $ranomxp[$bin]; my $decnomxp = $decnomxp[$bin]; my $pa_nomxp = $pa_nomxp[$bin]; # Create a temporary coordinate bin file my $binfile = "bin" . $bin . ".pnt"; ahapp::add_temp_file($binfile); # Run coordpnt to convert the input sky region to pixel region $status = ahgen::run_ftool("coordpnt", "input=$region", "outfile=$binfile", "telescop=$telescop", "instrume=$instrume", "ra=$ranomxp", "dec=$decnomxp", "roll=$pa_nomxp", "ranom=$ra_nom", "decnom=$dec_nom", "teldeffile=$teldeffile", "startsys=RADEC", "stopsys=LOWEST", "pixeltest=partial", "clobber=yes", ); if ( $status ) { ahlog::ah_err "coordpnt failed."; return $status; } # Verify that a region file was created if (isRequiredFileNotFound($binfile)) { return 1; } # Open the output region file and parse each pixel region open REGFILE, $binfile; foreach my $row ( ) { # Skip rows that are not pixel( ... ) unless ( $row =~ /pixel\(.*\)/ ) { next; } # Read the string of pixels between parentheses $row =~ /pixel\((.*)\)/; # Split the string on commas my @pixel_list = split (",",$1); # Store the pixel in the pixels array foreach my $ipix ( @pixel_list ) { $pixels[$ipix] = 1; } } close REGFILE; } # Condense the list of pixels my @pixlist = (); foreach my $ipix ( 0 .. 35 ) { push @pixlist, $ipix if $pixels[$ipix]; } # Merge the pixels array into a comma-delimited string my $pixlist = join (",",@pixlist); # Run ahmkregion to convert pixel list to DET region $status = ahgen::run_ftool("ahmkregion", "instrume=$instrume", "ra=45.0", # Not relevant "dec=30.0", # Not relevant "roll=0.0", # Not relevant "teldeffile=$teldeffile", "inregion=NONE", "outroot=$outroot", "outtextregion=no", "outxco=no", "pixlist=$pixlist", "cleanup=yes"); if ( $status ) { ahlog::ah_err "ahmkregion failed."; return $status; } # Verify that a DET BOX region file was created my $detreg = $outroot . ".DET.box.reg"; if (isRequiredFileNotFound($detreg)) { return 1; } # Create the simplified region file my @allinfo = (); my $numboxes = 0; # Open the output region file and parse each pixel region open DETREG, $detreg; foreach my $row ( ) { # Skip rows that are not pixel( ... ) unless ( $row =~ /box\(.*\)/ ) { next; } # Read the string of pixels between parentheses $row =~ /box\((.*)\)/; # Split the string on commas my ($xcen,$ycen,$xsize,$ysize,$angle) = split (",",$1); # round the center points my $xcenter = int($xcen + 0.5); my $ycenter = int($ycen + 0.5); # push the elements into a 2d array push @allinfo, [$xcenter,$ycenter, $xsize, $ysize]; $numboxes += 1; } close DETREG; unless ( $numboxes ) { ahlog::ah_err "No regions found after running ahmkregion"; return 1; } # Sort by ycenter, then xcenter my @sorted = sortxy (@allinfo); # Count the number of unique y-values my @nx = ( 0 ) x 6; # Initialize a 6-element array to zeroes foreach my $iy ( 0 .. $numboxes - 1 ) { # Get the ycenter value my $y1 = $sorted[$iy][1]; # add one to ycenter value $nx[$y1-1] += 1; } # Condense the boxes my $offset = 0; my $x1 = 0; my $y1 = 0; my $x2 = 0; my $y2 = 0; my $xs = 0; my $ys = 0; my $prevx = -1; open OUTREG, ">",$outreg; print OUTREG "physical\n"; foreach my $iy ( 1 .. 6 ) { my $nx = $nx[$iy-1]; # number of small boxes with that ycenter unless ( $nx ) { next; } # Skip to the next group if NX is zero my $start = 1; # initialize a new box foreach my $ix ( 0 .. $nx - 1 ) { # loop over all x boxes # Get the x and y center values from ahmkregion my $xcen = $sorted[$offset+$ix][0]; my $ycen = $sorted[$offset+$ix][1]; if ( $ix==0 ) { $prevx = $xcen; } # initialize the prevx variable if ( $start ) { # start a new box $x1 = $xcen - 0.5; # lower limit of box $x2 = $xcen + 0.5; # upper limit of box $ys = $sorted[$offset+$ix][3]; # ysize of box $start = 0; # start a new box } if ( ($xcen - $prevx) <= 1 ) { # adjacency condition $x2 = $xcen + 0.5 unless $ix == 0; # extend upper limit of box if ( $ix == $nx - 1 ) { # special case: last box at this y; close the box; $start = 1; # prepare to start a new box } } else { # close the box $start = 1; # prepare to start a new box } if ( $start ) { # last box at this y; write to the region file; my $newxcen = 0.5*($x1+$x2); # Use the midpoint between the first and last x-box my $newycen = $ycen; # ycenter is the current y group my $newxsiz = $x2-$x1; # calculate the x box size my $newysiz = $ys; # y box size is always 1 # Write to the region file print OUTREG "+box($newxcen,$newycen,$newxsiz,$newysiz)\n"; } # Save the previous xcenter for next x $prevx = $xcen; } # end loop over x # append to the offset $offset += $nx; } # end loop over y close OUTREG; } # end if regmode eq SKY # Run the extractor to create the image, lightcurve and spectra if ( $extract ) { # Create the ITYPE filtered image $status=ahgen::run_ftool("extractor", "exitnow=no", "filename=$infile", "eventsout=NONE", "imgfile=$img", "binf=1", "fullimage=yes", "phafile=NONE", "specbin=1", "wtmapb=yes", "wtmapfix=yes", "swmapx=no", "swmapy=no", "binh=1", "wmapver=2", "fitsbinlc=NONE", "qdpfile=NONE", "binlc=16", "lcthresh=0", "lcthwarn=3.0", "lctzero=yes", "unbinlc=NONE", "regionfile=NONE", "timefile=NONE", "gtinam=GTI", "xcolf=DETX", "ycolf=DETY", "zcolf=NONE", "xint=1.0", "yint=1.0", "tcol=TIME", "ecol=PI", "ccol=NONE", "gcol=ITYPE", "gstring=$resolist", "xcolh=X", "ycolh=Y", "gtitxt=NONE", "xronwn=NONE", "events=EVENTS", "gti=STDGTI", "timeorder=no", "timeref=40000.0", "eventkey=NONE", "phamax=TLMAX", "xfkey=TLMAX", "yfkey=TLMAX", "xhkey=TLMAX", "yhkey=TLMAX", "copyall=yes", "clobber=yes", ); if ( $status ) { ahlog::ah_err "extractor failed."; return $status; } # Create the grade and region filtered lightcurve and spectra $status=ahgen::run_ftool("extractor", "exitnow=no", "filename=$infile", "eventsout=NONE", "imgfile=NONE", "binf=1", "fullimage=yes", "phafile=$pha", "specbin=1", "wtmapb=yes", "wtmapfix=yes", "swmapx=no", "swmapy=no", "binh=1", "wmapver=2", "fitsbinlc=$lc", "qdpfile=NONE", "binlc=16", "lcthresh=0", "lcthwarn=3.0", "lctzero=yes", "unbinlc=NONE", "regionfile=$outreg", "timefile=NONE", "gtinam=GTI", "xcolf=DETX", "ycolf=DETY", "zcolf=NONE", "xint=1.0", "yint=1.0", "tcol=TIME", "ecol=PI", "ccol=NONE", "gcol=ITYPE", "gstring=$resolist", "xcolh=X", "ycolh=Y", "gtitxt=NONE", "xronwn=NONE", "events=EVENTS", "gti=STDGTI", "timeorder=no", "timeref=40000.0", "eventkey=NONE", "phamax=TLMAX", "xfkey=TLMAX", "yfkey=TLMAX", "xhkey=TLMAX", "yhkey=TLMAX", "copyall=yes", "clobber=yes", ); if ( $status ) { ahlog::ah_err "extractor failed."; return $status; } } # end if extract return 0; } # ------------------------------------------------------------------------------ sub finalize { my $status = 0; return 0; } # ------------------------------------------------------------------------------