Документ взят из кэша поисковой машины. Адрес оригинального документа : http://hea-www.harvard.edu/ChaMPlane/quantile/pre-release/grid.pl Дата изменения: Sat Jan 14 03:10:43 2006 Дата индексирования: Mon Feb 25 17:26:32 2013 Кодировка: Поисковые слова: http astrokuban.info astrokuban

#!/usr/bin/perl -w
#
# JaeSub Hong, 2003-2005, version 1.5
# Please report any problem or suggestion at jaesub@head.cfa.harvard.edu
#
# quantile grid generator
# generate grid patterns for spectra modeled by any two parameters
#
# Refer to
# J. Hong, E.M. Schlegel & J.E. Grindlay, 2004 ApJ 614 p508
# and references therein
# usage
# grid.pl
# -model modelname
# -range x.x:y.y
# -par1 name:value1,value2,value3,..:points
# -par2 name:value1,value2,value3,..:points
# -arf arf_file
# -rmf rmf_file
# [-etcp nameA:valueA,nameB:valueB,...]
# [-frac x.x,y.y,...]
# [-sherpa dir]
#
# requires
# sherpa (http://cxc.harvard.edu/sherpa/)
# grid.sl (provided along with this file)
#
# examples
#
# grid.pl \
# -model "xsphabs[A]*xspowerlaw[B]" \
# -range 0.3:8.0 \
# -par1 A.1:0.1,1,10:100 \
# -par2 B.1:1,2,3:100 \
# -etcp A.2:0.0,B.2:1 \
# -rmf test.rmf \
# -arf test.arf
#
# grid.pl \
# -model "xsphabs[A]*xspowerlaw[B]" \
# -range 0.3:8.0 \
# -par1 A.nH:0.01,0.1,1,10:200 \
# -par2 B.PhoIndx:0,1,2,3,4:200 \
# -etcp A.Redshift:0.0,B.norm:1.0 \
# -rmf test.rmf \
# -arf test.arf \
# > test_powerlaw.rdb
#
# grid.pl \
# -model "xsphabs[A]*xsbremss[B]" \
# -range 0.3:8.0 \
# -par1 A.nH:0.01,0.1,1,10:100 \
# -par2 B.kT:0.5,1.0,2.0,4.0,10.0:400 \
# -etcp A.Redshift:0.0,B.norm:1.0 \
# -rmf test.rmf \
# -arf test.arf \
# > test_bremss.rdb
#
#
# input
# Run the program without option for basic options or see below.
# For the model names and their parameters, refer Sherpa or
# xspec manual. As for the parameter name, one can simply use numbers
# (e.g. A.1 or B.1 instead of A.nh or B.PhoIndx), and the parameter
# names can be found in sherpa (load a model and type "show").
#
# output
# The result will be dumped into the screen in RDB format.
# One can save it by piping.
# e.g.) grid.pl .... > result.rdb
# For RDB format, refer
# http://hea-www.harvard.edu/MST/simul/software/docs/rdb.html
# Basically, it is an ASCII text table file where
# each field is separated by .
# The file will contains 2+n columns where n is the number of "frac"
# you request.
# First two columns are the two parameter values, and the rest will
# be the quantiles (Ex% but not Qx).
#
# limitation
# potentially many, but the accuracy of grid is limited by
# the energy scale of rmf and arf files. Need to be re-visited at
# some point.
#
#----------------------------------------------------------------------
# get the parameter files from command line
use Getopt::Long;
GetOptions(
"model=s" => \$model,
"range=s" => \$range,
"frac=s" => \$frac,
"par1=s" => \$par1,
"par2=s" => \$par2,
"etcp=s" => \$etcp,
"arf=s" => \$arf,
"rmf=s" => \$rmf,
"full" => \$full,
"sherpa=s" => \$sherpa,
"gridsl=s" => \$gridsl,
"batch=s" => \$batch,
"keep=s" => \$keep,
"raw" => \$raw,
"time=s" => \$time,
"help" => \$help,
"debug:i" => \$debug);

$debug = 0 unless defined $debug;

if (defined $help
|| !defined $model
|| !defined $par1
|| !defined $par2
|| !defined $arf
|| !defined $rmf
) {
($usage = <<" EOFHELP" ) =~ s/^\t\t//gm;
Find quantiles for given fractions
usage: $0 [-help] [-debug [X]]
-model modelname
-range X.X:Y.Y
-par1 name:value1,value2,value3,..:points
-par2 name:value1,value2,value3,..:points
-arf arf_file
-rmf rmf_file
[-frac X.X,Y.Y,...]
[-etcp nameA:valueA,nameB:valueB,...]
[-sherpa dir]
[-gridsl dir]
[-batch batch_file]
[-keep]
options
-help print this message
-debug set the level of debug
-model modelname
-range energy range
-arf arf file
-rmf rmf file
-frac fraction for quantile values
-par1 parameter properties for grid lines
-par2 parameter properties for grid lines
-etcp the rest of the parameters
-sherpa sherpa directory
-gridsl location (directory) of grid.sl
-batch temp batch file name for sherpa input
-keep if you want to keep the batch file

EOFHELP
print $usage;
exit;
}

#-----------------------------------------------------------------------
# here we go
# there's got to be a way to do all these in sherpa-slang, but
# my kung-fu is not there yet.

unless (defined $sherpa)
{ $sherpa = "sherpa"; }
else { $sherpa .= "/sherpa"; }

unless (defined $gridsl)
{ $gridsl = "grid.sl"; }
else { $gridsl .= "/grid.sl"; }

$batch = "sherpa.in" unless defined $batch;

#-----------------------------------------------------------------------
# write the sherpa batch input file

$range =~ s/:/,/;
$frac = "0.25,0.5,0.75" unless defined $frac;
$rmf = "ideal.rmf" unless defined $rmf;
$arf = "ideal.arf" unless defined $arf;

if (defined $raw) {
$x_axe = "channels";
$time = 1000000 unless defined $time;
$timeset = "fakeit time = $time";
$suffix = "_raw";
$rerun = "$timeset\nfakeit 1\n";
}
else {
$x_axe = "energy";
$timeset = "";
$suffix = "";
$rerun = "";
}


($input = <<"EOFINPUT" ) =~ s/^\t//gm;
paramprompt off
INSTRUMENT 1 = RSP[instrumentA]
instrumentA.rmf = $rmf
instrumentA.arf = $arf
ignore $x_axe :${range}:
source 1 = $model
$timeset
fakeit 1
frac=[$frac]
EOFINPUT

if (defined $etcp) {
@etcp = split/,/, $etcp;
foreach (@etcp) {
($name, $value) = split/:/, $_;
$input .= "$name = $value\n";
}
}

($p1n,$p1v,$p1p) = split/:/, $par1;
($p2n,$p2v,$p2p) = split/:/, $par2;

@p1v = split/,/, $p1v;
@p2v = split/,/, $p2v;

#($p1_min, $p1_max) = ($p1v[0], $p1v[$#p1v]);
#($p2_min, $p2_max) = ($p2v[0], $p2v[$#p2v]);

# yes, switching (1<->2) is right
#$p1s = ($p1_max-$p1_min)/$p2p;
#$p2s = ($p2_max-$p2_min)/$p1p;

# go for the par1 first

if (defined $full) {
$gid = 0;
foreach $j (1 .. $#p1v) {
($p1_min, $p1_max) = ($p1v[$j-1], $p1v[$j]);
# yes, switching (1<->2) is right
$p1s = ($p1_max-$p1_min)/$p2p * ($#p1v);
foreach ($p1=$p1_min;$p1<=$p1_max;$p1 +=$p1s) {
$input .="$p1n = $p1\n";

foreach $i (1 .. $#p2v) {
($p2_min, $p2_max) = ($p2v[$i-1], $p2v[$i]);
# yes, switching (1<->2) is right
$p2s = ($p2_max-$p2_min)/$p1p * ($#p2v);
foreach ($p2=$p2_min;$p2<=$p2_max;$p2 +=$p2s) {
$input .="$p2n = $p2\n$rerun"
."quant = mod_quantile${suffix}(1,frac)\n"
."quant_ = pri_quantile(quant)\n"
."print(\"==gx$gid\t\"+string($p1)+\"\t\"+string($p2)+\"\t\"+quant_)\n";
}
}
$gid++;
}
}
} else {
$gid = 0;
foreach (@p1v) {
$p1 = $_;
$input .="$p1n = $p1\n";

foreach $i (1 .. $#p2v) {
($p2_min, $p2_max) = ($p2v[$i-1], $p2v[$i]);
# yes, switching (1<->2) is right
$p2s = ($p2_max-$p2_min)/$p1p * ($#p2v);
foreach ($p2=$p2_min;$p2<=$p2_max;$p2 +=$p2s) {
$input .="$p2n = $p2\n$rerun"
."quant = mod_quantile${suffix}(1,frac)\n"
."quant_ = pri_quantile(quant)\n"
."print(\"==gx$gid\t\"+string($p1)+\"\t\"+string($p2)+\"\t\"+quant_)\n";
}
}
$gid++;
}

$gid = 0;
foreach (@p2v) {
$p2 = $_;
$input .="$p2n = $p2\n";

foreach $i (1 .. $#p1v) {
($p1_min, $p1_max) = ($p1v[$i-1], $p1v[$i]);
# yes, switching (1<->2) is right
$p1s = ($p1_max-$p1_min)/$p2p * ($#p1v);
foreach ($p1=$p1_min;$p1<=$p1_max;$p1 +=$p1s) {
$input .="$p1n = $p1\n$rerun"
."quant = mod_quantile${suffix}(1,frac)\n"
."quant_ = pri_quantile(quant)\n"
."print(\"==gy$gid\t\"+string($p1)+\"\t\"+string($p2)+\"\t\"+quant_)\n";
}
}
$gid++;
}
}

open(BAT,"> $batch") or die "can't find $batch\n";
print BAT $input;
close(BAT);

#-----------------------------------------------------------------------
# now go, really.

die "can't find $gridsl\n" unless -e $gridsl;
@output = `$sherpa --slscript $gridsl --batch $batch`;
`rm -rf $batch` unless defined $keep;

#-----------------------------------------------------------------------
# let's try to understand the results

@frac = split/,/,$frac;

$tit ="gridID\t$p1n\t$p2n";
$unit ="S\tN\tN";
for ($i=0;$i<=$#frac;$i++){
$ifrac = sprintf("%d",$frac[$i]*100.);
$tit .= "\tE$ifrac";
$unit .= "\tN";
}

$tit =~ s/\./_/g;

print "#\n";
print "# $0 \\ \n";
print "# -model \"$model\" \\ \n";
print "# -range $range \\ \n";
print "# -frac $frac \\ \n";
print "# -par1 $par1 \\ \n";
print "# -par2 $par2 \\ \n";
print "# -etcp $etcp \\ \n" if defined $etcp;
print "# -rmf $rmf \\ \n";
print "# -arf $arf \\ \n";
print "#\n";
print "$tit\n";
print "$unit\n";
foreach (@output){
next unless /^==(.*)/;
print "$1\n";
}

exit;