Files
Estia-McStas/bin/map_integrator.pl
T

93 lines
2.5 KiB
Perl
Executable File

#!/usr/bin/perl -w
#
use Math::Trig ;
use Getopt::Std ;
my %options=();
getopts("f:l:L:t:T:a:b:", \%options);
# otions:
# -f <input file>
# -l <lambda-range>
# -L
# -t <theta-range>
# -T
# -a <resolution>
# -b <background>
unless ($options{f}) {die "*** no input filename is given (-f <file>)\n"} ;
if ($options{l}) {$lamin = $options{l}} else {$lamin = 1} ;
if ($options{L}) {$lamax = $options{L}} else {$lamax = 20} ;
if ($options{t}) {$thmin = $options{t}} else {$thmin = 0.01} ;
if ($options{T}) {$thmax = $options{T}} else {$thmax = 100} ;
if ($options{a}) {$dqdq = $options{a}} else {$dqdq = 0.01} ;
if ($options{b}) {$bg = $options{b}} else {$bg = 0.0} ;
$q0 = 0.001 ;
print "\nmap_integrator\n" ;
$normalisedfile = $options{f}."/Rlt" ;
$reflectivityfile = $options{f}."/Rq" ;
open OUT, ">$reflectivityfile" ;
print OUT "# reflectivity R(q_z) calulated by map_integrator.pl from the output of \n" ;
unless (open IN, "<$normalisedfile") {die "*** $normalisedfile does not exist!\n\n"} ;
print " reading $normalisedfile\n" ;
$qmin = 10000 ;
$qmax = 0 ;
while (<IN>) {
if (/#/) {
print OUT "$_" ;
} else {
chomp ;
@line = split / */ ;
if ( $line[3] ) {
(undef,undef,undef,$lambda,$theta,undef,undef,$value,$sigma,undef,undef,$qzerr,undef) = @line ;
$qz = 4*pi * sin(deg2rad($theta)) / $lambda ;
#
if ($qz > 0 and $lambda >= $lamin and $lambda <= $lamax and $theta >= $thmin and $theta <= $thmax and $sigma < 1e3) {
$q = int( 0.5 + log($qz/$q0) / log(1+$dqdq) ) ;
if ($qmax<$q){$qmax=$q} ;
if ($qmin>$q){$qmin=$q} ;
$qzerr[$q] += $qzerr / $sigma ;
$sum[$q] += $value / $sigma ;
$norm[$q] += 1 / $sigma ;
$err[$q] += 1 / $sigma**2 ;
#
$nn[$q] += 1 ;
} ;
#
} ;
} ;
} ;
close IN ;
print " writing to $reflectivityfile " ;
for ($q=$qmin;$q<=$qmax;$q++) {
#if ($q*$dq>0.022) {
# $roq = 1 - ($q*$dq/0.022-1)/15 ;
#} else {
# $roq = 1 ;
#} ;
unless ($norm[$q]) {
$sum[$q] = 1e-12 ;
$norm[$q] = 1 ;
$nn[$q] = 1e-24 ;
} ;
if ( $sum[$q]/$norm[$q] > 1e-11 ){
printf OUT " %8.6f %11.5e %8.6f %11.5e\n"
, $q0*(1+$dqdq)**$q
, $sum[$q] / $norm[$q] - $bg
, 0.5 * $qzerr[$q] / $norm[$q]
, 0.5 * sqrt( 1 / $err[$q] )
;
#, 0.5*$q0*(1+$dqdq)**$q * $dqdq/sqrt(1+$dqdq) # qz bin width
#, 0.5 * sqrt( $nn[$q] ) / $norm[$q]
} ;
} ;
close OUT ;
print " ...done\n\n" ;