#!/usr1/local/bin/perl5 -w ################################### # Robin Corbet (UMBC), 2010-08-09 # ################################### use strict 'subs'; use strict 'refs'; $VERSION = "1.2"; $help = " like_lc.pl Typing like_lc.pl -help will display this help information. Example perl script that generates Fermi LAT light curves using likelihood analysis. This is not an official part of the Fermi analysis software. WARNING! The models and other parameters used here are for DEMONSTRATION purposes only! You should carefully evaluate the choice of parameters for your own particular needs! To use this script you must be set up to run both the Fermi science tools and HEASOFT. This script checks the environment variables LHEASOFT, PFILES, and FERMI_DIR as a check on this. In addition, if you are using the auto-generation option to automatically make a model file from the 1FGL catalog, the environment variable MY_FERMI_DIR needs to be set. This is where the 1FGL catalog and the make1FGLxml.py script must be located if used. The model file generated by make1FGLxml.py uses the parameters in the 1FGL catalog. In this example only the parameters of the source of interest are left free, others are held fixed at their catalog values. Spectra are power laws. This may very well not be optimum! Among the inputs prompted for are two files: slist.dat This should contain a list of source names and their coordinates in the form: 1FGLJ0240.5+6113 40.1317, 61.2281 i.e. the 1FGL catalog name must be given because if you will use make1FGLxml.py to generate the model file. Note that there is _no_ space in the source name. i.e. \"1FGLJ0240.5+6113\" rather than \"1FGL J0240.5+6113\" Source names should not be duplicated because the output files are constructed from the source names. Optionally, an entry in the slist.dat file my have the form: A_Source 40.1317, 61.2281 super_model.xml where the 4th entry is the name of the model file to use for that source. The source name (e.g. \"A_Source\") must be contained in the model file specified. plist.dat This contains a list of photon files, one file per line. It is also necessary to provide a single spacecraft (FT2) file. Several temporary files are created while runing this script. At least most are deleted at the end. The output file for each source is an ASCII file called XXX_like_lc.dat, where \"XXX\" is the source name as given in slist.dat This file contains: time (MJD) flux (cts/cm^2/s), flux error, time bin half width (days), TS value, spectral index, index error CHANGES: 1.2 Fix for bug causing mangled output if logParabola model was used 1.1 Support for selecting either P6_V3 or P6_V11 IRFs "; $TRUE = 1; $FALSE = 0; # minimal output from tools $CHATTER = 0; # Conversion from MET (s) to MJD (days, Terrestrial Time) $MJDREF = 51910.0 + 7.428703703703703E-4; $DEFAULT_BIN_SIZE = 5; $DEFAULT_ZENITH_LIMIT = 105; $DEFAULT_BORE = 180; $DEFAULT_EVENT_CLASS_MIN = 3; $DEFAULT_MODEL_FILE = "outmodel.xml"; $DEFAULT_SOURCE_LIST = "slist.dat"; $DEFAULT_FT1_LIST = "plist.dat"; $DEFAULT_IRF_CODE = 2; # P6_V11 $DEFAULT_CATALOG_FILE = "gll_psc_v02.fit"; $DEFAULT_TS_THRESHOLD = 2; # Region size $DEFAULT_ROI = 10; # Default, "hidden" parameters $DEFAULT_EMIN = 100; $DEFAULT_EMAX = 200000; $catalog_file = $DEFAULT_CATALOG_FILE; print "like_lc.pl: V$VERSION\n"; # Check command line parameters $numArgs = $#ARGV + 1; if($numArgs > 0){print "$numArgs command-line arguments:\n"}; foreach $argnum (0 .. $#ARGV) { print "$ARGV[$argnum]\n"; if($ARGV[$argnum] =~ /-18m/) { print "Using 18 month catalog instead of 1FGL\n"; $catalog_file = "gll_psc18month_uw11b.fits"; } elsif($ARGV[$argnum] =~ /-help|-h|--help/) { print $help; exit; } } $irf_code = get_value("Which IRFs? (1 = P6_V3, 2 = P6_V11)", $DEFAULT_IRF_CODE, 1, 2); if($irf_code == 1){ $irf_prefix = "P6_V3_"; $gll_file = "gll_iem_v02.fit"; $isotropic_file = "isotropic_iem_v02.txt"; } elsif($irf_code == 2){ $irf_prefix = "P6_V11_"; $gll_file = "gll_iem_v02_P6_V11_DIFFUSE.fit"; $isotropic_file = "isotropic_iem_v02_P6_V11_DIFFUSE.txt"; } else{ print "Got an invalid IRF code ($irf_code) somehow!\n"; } $auto_model = get_yn("Auto-generate model files from 1FGL catalog?", $TRUE); ########################################################## # Check if FTOOLS are available ################################# if($ENV{'LHEASOFT'} !~/\S/ || $ENV{'PFILES'} !~/\S/ ) { print "\n You need to set up HEASOFT to use this script.\n\n"; exit(0); } if($ENV{'FERMI_DIR'} !~/\S/ ) { print "\n You need to set up the Fermi Science Tools to use this script.\n\n"; exit(0); } if($auto_model){ if($ENV{'MY_FERMI_DIR'} !~/\S/ ){ print "\n You need to set the environment variable MY_FERMI_DIR where 1FGL catalog and python script are located.\n\n"; exit(0); } } # Official Fermi directory $fermi_dir = $ENV{'FERMI_DIR'}; # My own Fermi directory. 1FGL catalog and make1FGLxml.py must be here. $my_fermi_dir = $ENV{'MY_FERMI_DIR'}; ######################################################## # Prompt for inputs ####################### if($auto_model){ print "auto model\n"; } else{ print "NOT auto model\n"; $default_model_file = get_file("Give XML model file", $DEFAULT_MODEL_FILE); } $ts_threshold = get_value("Give minimum TS value for including fluxes in output", $DEFAULT_TS_THRESHOLD, "INDEF", "INDEF"); # File containing source names and parameters and open the file $source_list = get_file_open("Give source parameter file", *SOURCE_LIST, $DEFAULT_SOURCE_LIST); # = get_file_open2("Give source parameter file", $DEFAULT_SOURCE_LIST); # Get names of files which contain lists of the # photon files and a single spacecraft file # that must cover the time range of all the photon files $ft1_file_list = get_file("Give photon file name list", $DEFAULT_FT1_LIST); # Look for spacecraft files in the current directory. opendir(DIR, "."); @files = grep(/_SC00\.fits$/,readdir(DIR)); closedir(DIR); $files = @files; #print "no. of SC files found = $files\n"; # If there's exactly one spacecraft file in the current directory suggest that as # default if($files == 1){ $ft2_file = get_file("Give spacecraft (FT2) file name", $files[0]); } else{ $ft2_file = get_file("Give spacecraft (FT2) file name"); } $roi = get_value("Give aperture radius (degrees)", $DEFAULT_ROI, 0, 180); # source radius is 10 degrees bigger than ROI $srcrad = $roi + 10; $emin_to_use = get_value("Give Emin (MeV)", $DEFAULT_EMIN, 0, 1E10); $emax_to_use = get_value("Give Emax (MeV)", $DEFAULT_EMAX, $emin_to_use, 1E10); $zenith_limit = get_value("Give Zenith limit (degrees)", $DEFAULT_ZENITH_LIMIT, 0, "INDEF"); $event_class_min = get_value("Give EVENT_CLASS minimum (use 0 for simulated data)", $DEFAULT_EVENT_CLASS_MIN, 0, 3); if($event_class_min == 3 || $event_class_min == 0){ $irfs = $irf_prefix."DIFFUSE"; } elsif($event_class_min == 2){ $irfs = $irf_prefix."SOURCE"; } elsif($event_class_min == 1){ $irfs = $irf_prefix."TRANSIENT"; } else{ print "event class minimum not allowed\n"; } $bore_limit = get_value("Give Bore limit (degrees)", $DEFAULT_BORE, 0, 360); $bin_size = get_value("Give bin size (days)", $DEFAULT_BIN_SIZE, 0, 1e10); # Convert to seconds for use in script $bin_size = $bin_size * 86400; ############################################ # Start of source loop ############### print "Getting sources\n"; while(chomp($line = )){ print "line = $line\n"; # remove any commas $line =~ s/\,/ /g; # Need to have at least source name, RA and dec # in input file @fields = split(/\s+/, $line); $length = @fields; $source = $fields[0]; $ra = $fields[1]; $dec = $fields[2]; # Check we have at least 3 columns in the source list file if($length < 3){ print "*** Error in input file!\n"; print "The source list file ($source_list) contains a line with less than 3 entries:\n"; print "$line\n"; print "Each line in this file must contain (at least) source name, RA, dec.\n"; exit(0); } # We have $default_model_file, $model_file, $model_file_in, $model_file_out # which is probably too many and we should be able to simplify! # Check if a model file was specified in the entry $model_given = $FALSE; if($length > 3) { $model_file = $fields[3]; open(TEMP, "$model_file") || die "Can't open the specified model file, $model_file\n"; close(TEMP); $model_given = $TRUE; } else { $model_file = $default_model_file; } # change any plus sign in source name to "p" # to avoid FTOOLS problems (may be better way to work # around this! # Use source name as is so can be used with xml file directly $src = $source; $src =~ s/\+/p/g; # need to escape plus signs to do matching $source_find = $source; $source_find =~ s/\+/\\\+/g; # Create output file to contain light curve $outfile = $src."_like_lc.dat"; open(OUTFILE, ">$outfile") || die "Can't open output results file: $outfile\n"; # Make names of temporary files $eventfile0 = make_temp("eventfile0.fits"); $eventfile1 = make_temp("eventfile1.fits"); $eventfile2 = make_temp("eventfile2.fits"); $expcube = make_temp("expcube.fits"); $expmap = make_temp("expmap.fits"); $results_file = make_temp("results.dat"); $specfile = make_temp("counts_spectra.fits"); @tidy_file_list = ($eventfile1, $eventfile2, $expcube, $expmap, $results_file, $specfile); # Combine all photon files together $command = "gtselect chatter=$CHATTER zmax=180 emin=$emin_to_use emax=$emax_to_use infile=\"\@$ft1_file_list\" outfile=\'$eventfile0\' ra=$ra dec=$dec rad=$roi evclsmin=$event_class_min evclsmax=10 tmin=0 tmax=0"; doit($command); if($auto_model && !$model_given){ $py_file = make_temp("pyfile.py"); # make model file (using Tyrel Johnson's script) # Make temporary perl script to run make1FGLxml $model_file_in = make_temp("myLATxmlmodel.xml"); $model_file_out = "edit_myLATxmlmodel.xml"; open(PY_FILE, ">$py_file"); print PY_FILE "import sys\n"; print PY_FILE "sys.path.append(\"$my_fermi_dir\")\n"; # print PY_FILE "from rcmake1FGLxml import *\n"; print PY_FILE "from make1FGLxml import *\n"; print PY_FILE "s1=srcList(\'$my_fermi_dir/$catalog_file\',\'$eventfile0\',\'$model_file_in\')\n"; print PY_FILE "s1.makeModel(\'$fermi_dir/refdata/fermi/galdiffuse/$gll_file\',\'gal_v02\',\'$fermi_dir/refdata/fermi/galdiffuse/$isotropic_file\',\'eg_v02\')\n"; close(PY_FILE); $command = "python $py_file"; doit($command); # Edit output of make1FGLxml to make everything constant except source of interest open(MODEL_FILE_IN, $model_file_in) || die "Can't open input model file: $model_file_in\n"; open(MODEL_FILE_OUT, ">$model_file_out") || die "Can't open output model file: $model_file_out\n"; $MY_SOURCE = $FALSE; print "looking for source: $source\n"; while(chomp($line = )){ # Check if it's a source name if($line =~ /source name/){ # print "$line\n"; # Check to see if line matches source name if($line =~ /$source_find/){ print "found the source!\n"; $MY_SOURCE = $TRUE; } else { $MY_SOURCE = $FALSE; } } # Freeze parameters for all sources at cataloged values except # for # This is for demonstration only, it is not necessarily recommended. # Time taken to perform likelihood analysis depends on how # many parameters are free. if(!$MY_SOURCE){$line =~ s/parameter free=\"1\"/parameter free=\"0\"/g;} print MODEL_FILE_OUT "$line\n"; } print MODEL_FILE_OUT "$line\n"; close(MODEL_FILE_IN); close(MODEL_FILE_OUT); $command = "rm $py_file"; doit($command); $command = "rm $model_file_in"; doit($command); } else { # Use a file which is either the one entered at a prompt, or one read from the source list file $model_file_out = $model_file; } # The file just created is the "master" combined file. # Rerun gtselect on it later to just use portions of the data file. # get TSTART and TSTOP of master photon file $command = "fkeypar \"$eventfile0\" TSTART"; doit($command); chomp($ph_tstart = `pget fkeypar value`); print "photon start time = $ph_tstart\n"; $command = "fkeypar \"$eventfile0\" TSTOP"; doit($command); chomp($ph_tstop = `pget fkeypar value`); print "photon stop time = $ph_tstop\n"; ######################################################## # Start of time loop ########################### $time1 = $ph_tstart; while ($time1 < $ph_tstop) { $time2 = $time1 + $bin_size; # gtselect - reselect time range $command = "gtselect chatter=$CHATTER zmax=180 emin=$emin_to_use emax=$emax_to_use infile=\"$eventfile0\" outfile=\'$eventfile1\' ra=$ra dec=$dec rad=$roi evclsmin=$event_class_min evclsmax=10 tmin=$time1 tmax=$time2"; doit($command); # gtmktime - good time intervals $command = "gtmktime chatter=$CHATTER scfile=\"$ft2_file\" filter=\"(DATA_QUAL==1) && (LAT_CONFIG==1) && (angsep(RA_ZENITH,DEC_ZENITH,$ra,$dec)\+$roi<$zenith_limit) && (angsep($ra,$dec,RA_SCZ,DEC_SCZ)<$bore_limit)\" roicut=n evfile=\"$eventfile1\" outfile=\"$eventfile2\""; doit($command); # Check for existence of output file # Can sometimes get: # Caught St13runtime_error at the top level: Zero rows returned from FT2 file for this filter: if(!check_file_exists("$eventfile2")){ tidy_files(@tidy_file_list); $time1 = $time2; next; } # gtltcube - live time cube $command = "gtltcube chatter=$CHATTER evfile=$eventfile2 scfile=$ft2_file outfile=$expcube dcostheta=0.025 binsz=1"; doit($command); # gtexpmap - exposure map $command = "gtexpmap chatter=$CHATTER evfile=$eventfile2 scfile=$ft2_file expcube=$expcube outfile=$expmap irfs=$irfs srcrad=$srcrad nlong=120 nlat=120 nenergies=20"; doit($command); # gtdiffrsp - calculate diffuse response columns. (If using other than DIFFUSE, which is pre-computed, this will take a long time!) $command = "gtdiffrsp chatter=$CHATTER evfile=\'$eventfile2\' scfile=$ft2_file srcmdl=$model_file_out irfs=$irfs"; doit($command); # gtlike - do the actual likelihood analysis $command = "gtlike results=$results_file specfile=$specfile chatter=$CHATTER irfs=$irfs expcube=$expcube srcmdl=$model_file_out statistic=UNBINNED optimizer=MINUIT evfile=$eventfile2 scfile=$ft2_file expmap=$expmap cmap=none bexpmap=none"; doit($command); # Extract flux from results of gtlike #open(RESULTS_FILE, $results_file) || # die "Can't open gtlike output file: $results_file\n"; # Extract flux from results of gtlike # skip to next time interval if no results file available unless(open(RESULTS_FILE, $results_file)){ print "Can't open gtlike output file: $results_file\n"; tidy_files(@tidy_file_list); $time1 = $time2; next; } # Find the line containing the source flux and write to the output light curve file $MY_SOURCE = $FALSE; $FOUND_SOURCE = $FALSE; print "Looking for: $source_find\n"; #while(chomp($line = )){ while(){ chomp($line = $_); # if($line =~ /1FGL/){ # Allow for more general model files # This is a bit flakey, relies on source name being followed by "Integral", "Prefactor", "norm", or "Normalization" if($line =~ /Integral|Prefactor|norm|Normalization/){ #print "Found Integral: $line\n"; if($line =~ /$source_find/){ print "found the source!\n"; $MY_SOURCE = $TRUE; $FOUND_SOURCE = $TRUE; #print "MY_SOURCE = TRUE\n"; } else { $MY_SOURCE = $FALSE; #print "MY_SOURCE = FALSE\n"; } } # Extract TS value if($MY_SOURCE && ($line =~ /TS value/)){ # parse out value. # This match for TS value would only work with a decimal point present... # $line =~ /(\d+\.\d+)/; $line =~ /TS value\': \'(.*)\',/; $ts_value = $1; } # Extract Index and error value if($MY_SOURCE && ($line =~ /Index/)){ # parse index and error (yes, this could be rewritten more compactly!) # but perhaps has more robustiness?? $spec_index = $line; $spec_index =~ s/\'Index\'://; $spec_index =~ s/\+\/\-//; $spec_index =~ s/\'//g; $spec_index =~ s/,//; } if($MY_SOURCE && ($line =~ /Flux/)){ # parse flux and error (yes, this could be rewritten more compactly!) # but perhaps has more robustiness?? $out = $line; $out =~ s/\'Flux\'://; $out =~ s/\+\/\-//; $out =~ s/\'//g; $out =~ s/,//; # write output. # Convert MET to MJD $time = ($time1 + $time2)/2; $time = ($time/86400) + $MJDREF; $terr = ($time2 - $time1)/(86400 * 2); # Exclude measurements with low TS values if($ts_value > $ts_threshold){ print OUTFILE "$time $out $terr $ts_value $spec_index\n"; } else{ print "TS value ($ts_value) less than threshold ($ts_threshold)\n"; } } } # tidy up files tidy_files(@tidy_file_list); #$command = "rm $eventfile1"; # doit($command); #$command = "rm $eventfile2"; # doit($command); #$command = "rm $expcube"; # doit($command); #$command = "rm $expmap"; # doit($command); #$command = "rm $results_file"; # doit($command); #$command = "rm $specfile"; # doit($command); if(!$FOUND_SOURCE){ $command = "rm $eventfile0"; doit($command); die "ERROR: Couldn't find source $src in the results from gtlike. The source specified in slist.dat is (probably) not present in the model file. Check for typos in slist.dat and/or the model file.\n";} ########################### # End of time loop for one source ########################### $time1 = $time2; } close(OUTFILE); # End of source loop } $command = "rm $eventfile0"; doit($command); ############################ # Subroutines ########## ################################################# # System call, and write to screen what it's doing # Usage: doit("string to execute"); sub doit{ print "\n$_[0]\n"; system($_[0]); } #################################################### # make a temporary filename to use that includes the # process ID to (hopefully) make it unique # Usage: $file_name = make_temp("file_name_base") sub make_temp{ # The line below created output files on /tmp # However, many systems only have limited space there. #return "/tmp/".$$.$_[00]; return "tmp_".$$.$_[00]; } ####################################### # Prompt for and get a file name. # Check that file exists before accepting its name. # Usage: $file_name = get_file("String to prompt for name"[, "default name"]); # The default file name is not required. sub get_file{ my $prompt = $_[0]; my $default = $_[1]; my $waiting = 1; while($waiting){ if ($default){ print "$prompt [default = $default]\: "; } else{ print "$prompt\: "; } chomp($file_name = ); if($file_name eq "" && $default) {$file_name = $default; print "(Using default: \"$default\")\n";} if (open(CHECK, $file_name)){ # print "opened $file_name\n"; $waiting = 0; close(CHECK); } else { print "Can't open file: $file_name\n"; } } return $file_name; } ####################################### # Prompt for and get a file name, then open it. # This should probably be combined with get_file() # Check that file exists before accepting its name. # Usage: $file_name = get_file("String to prompt for name", FILE_HANDLE, [, "default name"]); # The default file name is not required. sub get_file_open{ my $prompt = $_[0]; my $check = $_[1]; my $default = $_[2]; my $waiting = 1; while($waiting){ if ($default){ print "$prompt [default = $default]\: "; } else{ print "$prompt\: "; } chomp($file_name = ); if($file_name eq "" && $default) {$file_name = $default; print "(Using default: \"$default\")\n";} if (open($check, $file_name)){ # print "opened $file_name\n"; $waiting = 0; # close(CHECK); } else { print "Can't open file: $file_name\n"; } } return $file_name; } ############################################################ # Prompt for a Y/N response and return true or false (1 or 0) # default value is determined by argument # Usage: $response = get_yn("Prompt", $default); # $default should be 1, 0, or "INDEF" sub get_yn{ my $prompt = $_[0]; my $default = $_[1]; if($default eq 1){$default_yn = "Y";} elsif($default eq 0){$default_yn = "N";} my $INDEF = "INDEF"; my $waiting = 1; while($waiting){ if($default eq $INDEF){ print "$prompt [y/n]\: "; } else{ if ($default){ print "$prompt [y/n; default = y]\: "; } else{ print "$prompt [y/n; default = n]\: "; } } chomp($response = substr(uc(),0,1)); #print "response = $response\n"; #if($response eq "" && ($default != $INDEF)) if($response eq "") {$value = $default; print "(Using default: \"$default_yn\")\n"; $waiting = 0;} # Check response format OK elsif($response eq "Y"){ # print "got y\n"; $value = 1; $waiting = 0; } elsif($response eq "N"){ # print "got n\n"; $value = 0; $waiting = 0; } else { print "Response ($response) is not \"y\" or \"n\"\n"; } } return $value; } ############################################################ # Prompt for and get a numeric value. # If specified, check value is within specified limits. # Usage: $value = get_value("String to prompt for value", default, minimum, maximum) # If a default, min, or max value is not to be specified, call with value set to # the string "INDEF" sub get_value{ my $prompt = $_[0]; my $default = $_[1]; my $minimum = $_[2]; my $maximum = $_[3]; #print "default = $default, min = $minimum, max = $maximum\n"; my $INDEF = "INDEF"; my $waiting = 1; while($waiting){ if ($default eq $INDEF){ print "$prompt\: "; $use_default = 0; } else{ print "$prompt [default = $default]\: "; $use_default = 1; } chomp($value = ); if($value eq "" && ($use_default)) #if($value eq "") {$value = $default; print "(Using default: \"$default\")\n";} # Check data range is OK if(check_value($value, $minimum, $maximum)){ $waiting = 0; } #if (($value <= $maximum || $maximum eq $INDEF) && # ( $value > $minimum ) || $minimum eq $INDEF ){ # $waiting = 0; # } else { print "Value ($value) is not in the range $minimum to $maximum\n"; } } return $value; } ################################################### # # check for the existence of a specified file # return TRUE if it exists, otherwise FALSE sub check_file_exists{ my $file_name = shift; if (open(CHECK, $file_name)){ # print "opened $file_name\n"; close(CHECK); return $TRUE; } else{ print "File $file_name not found\n"; return $FALSE; } } ################################################ # tidy_files # delete a bunch if intermediate files sub tidy_files{ while(@_){ my $file_name = shift(@_); my $command = "rm $file_name"; doit ($command); } } ########################################## # sub check_value{ my $value = shift; my $minimum = shift; my $maximum = shift; my $result = $TRUE; if($value !~ /[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?/){ print "$value is not a number!\n"; $result = $FALSE; return $result; } #[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)? if($maximum ne "INDEF"){ if($value > $maximum){$result = $FALSE;} } if($minimum ne "INDEF"){ if($value < $minimum){$result = $FALSE;} } return $result; }