1#!/usr/local/bin/perl 2# 3# Implements plotting of McStas resolution function data using PGPLOT 4# 5# 6# This file is part of the McStas neutron ray-trace simulation package 7# Copyright (C) 1997-2004, All rights reserved 8# Risoe National Laborartory, Roskilde, Denmark 9# Institut Laue Langevin, Grenoble, France 10# 11# This program is free software; you can redistribute it and/or modify 12# it under the terms of the GNU General Public License as published by 13# the Free Software Foundation; version 2 of the License. 14# 15# This program is distributed in the hope that it will be useful, 16# but WITHOUT ANY WARRANTY; without even the implied warranty of 17# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18# GNU General Public License for more details. 19# 20# You should have received a copy of the GNU General Public License 21# along with this program; if not, write to the Free Software 22# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 23 24use PDL; 25use PDL::Core; 26use PDL::Math; 27use PDL::Slatec; 28use PDL::IO::FastRaw; 29use PGPLOT; 30use PDL::Graphics::PGPLOT; 31 32# Determine the path to the McStas system directory. This must be done 33# in the BEGIN block so that it can be used in a "use lib" statement 34# afterwards. 35BEGIN { 36 ENV_HEADER 37} 38use lib $MCSTAS::perl_dir; 39 40require "mcfrontlib2D.pl"; 41 42$PI = 3.14159265358979323846; 43 44sub read_mcstas_info { 45 my ($file) = @_; 46 my $basedir; 47 $basedir = $1 if $file && $file =~ m|^(.*)/[^/]*$|; 48 my $handle = new FileHandle; 49 open $handle, $file or die "Could not open file '$file'"; 50 $info = read_simulation_info($handle); 51 close($handle); 52 return ($info); 53} 54 55sub read_mcstas_res { 56 my ($filename) = @_; 57 my ($data,$kix,$kiy,$kiz,$kfx,$kfy,$kfz,$x,$y,$z,$pi,$pf); 58 my ($size,$ki,$kf,$q,$qx,$qy,$qz,$p,$Ei,$Ef,$w); 59 my ($r,$qx_mc,$qy_mc,$qz_mc,$w_mc, $npts,$cntr,$gaus); 60 my ($ave_q,$unit_q,$unit_n,$unit_z,$tmat,$q_t); 61 my ($A,$ave_A,$mid_A,$C,$umat,$C_t,$res_mat); 62 my ($pos); 63 64 # Read data from file (either raw or ascii). 65 if($filename =~ /\.raw$/) { 66 $data = readfraw($filename); 67 ($kix,$kiy,$kiz,$kfx,$kfy,$kfz,$x,$y,$z,$pi,$pf) = dog $data; 68 } else { 69 ($kix,$kiy,$kiz,$kfx,$kfy,$kfz,$x,$y,$z,$pi,$pf) = rcols($filename); 70 $data = cat ($kix,$kiy,$kiz,$kfx,$kfy,$kfz,$x,$y,$z,$pi,$pf); 71 } 72 # Compute some basic entities 73 ($size) = $kix->dims; 74 $ki = cat($kix, $kiy, $kiz); 75 $kf = cat($kfx, $kfy, $kfz); 76 $q = $ki - $kf; 77 $Ei = 2.072*($kix*$kix+$kiy*$kiy+$kiz*$kiz); 78 $Ef = 2.072*($kfx*$kfx+$kfy*$kfy+$kfz*$kfz); 79 $w = $Ei-$Ef; 80 $p = $pi*$pf; 81 # Compute coordinate change: X along average Q vector projected 82 # into plane, Y perpendicular to X in plane, Z upwards. 83 $ave_q = sumover($q*$p->dummy(1,3)) / (sum($p)); 84 $unit_q = $ave_q->copy; 85 $unit_q->set(1,0); # Force into scattering plane. 86 $unit_q /= sqrt(inner($unit_q,$unit_q)); 87 $unit_n = pdl($unit_q->at(2), 0, -$unit_q->at(0)); 88 $unit_z = pdl(0,1,0); 89 # Build orthogonal transformation matrix, and change coordinates of Q. 90 $tmat = cat ($unit_q, $unit_n, $unit_z); 91 $q_t = xchg(PDL::Primitive::matmult($tmat,$q->dummy(2)),1,2); 92 $q_t = $q_t->clump(2); 93 ($qx,$qy,$qz) = dog $q_t; 94 95 # Now compute resolution matrix. 96 $A = append($q->transpose, $w->dummy(0)); 97 $ave_A = sumover($A->transpose*$p->dummy(1,4)) / sum($p); 98 $mid_A = $A - $ave_A->dummy(1); 99 # Get the covariance matrix in original coordinates. 100 $C = PDL::Primitive::matmult 101 ($mid_A->transpose, $mid_A*$p->dummy(0,4)) / sum($p); 102 # Change coordinates, and compute the resolution matrix. 103 $umat = transpose(append(transpose(append($tmat,pdl [0])), 104 pdl [[0],[0],[0],[1]])); 105 $C_t = inner2t($umat->transpose,$C,$umat); 106 $res_mat = $C_t->matinv; 107 print "The covariance matrix is\n"; 108 print $C_t; 109 print "and the resolution matrix is\n"; 110 print $res_mat; 111 print "Gaussian half width [Qx Qy Qz En] in Angs-1 and meV are\n"; 112 $gqx = int(2.3548/sqrt($res_mat->at(0,0))*1e4)/1e4; 113 $gqy = int(2.3548/sqrt($res_mat->at(1,1))*1e4)/1e4; 114 $gqz = int(2.3548/sqrt($res_mat->at(2,2))*1e4)/1e4; 115 $gen = int(2.3548/sqrt($res_mat->at(3,3))*1e4)/1e4; 116 print "[$gqx $gqy $gqz $gen]\n"; 117 118 return($qx,$qy,$qz,$w,$p,$C_t,$res_mat,$size); 119} 120 121sub plot_mcstas_res { 122 my ($filename,$device,$qx,$qy,$qz,$w,$p,$C_t,$res_mat,$size,$interactive,$si) = @_; 123 # Plot histograms for the four 1-d projections. 124 if (defined(&dev)) { $dev = dev "$device",4,2; } 125 else { $dev = pgopen("$device"); pgsubp(4,2); } 126 die "DEV/PGOPEN $device failed!" unless $dev > 0; 127 128 pgsch(2.1); 129 pgsci(1); 130 131 # Make a 3d visualization of the resolution elipsoid. Use MC 132 # choice to eliminate the weights. 133 $r = random $size; 134 $qx_mc = $qx->where($p > $r*max($p)); 135 $qy_mc = $qy->where($p > $r*max($p)); 136 $qz_mc = $qz->where($p > $r*max($p)); 137 $w_mc = $w->where($p > $r*max($p)); 138 $npts = $w_mc->nelem; 139 $R0 = 1; 140 $NP = $res_mat; 141 142 # plot 2D histograms, and add the gaussian ellipsoid on top of each 143 pgpanl(1,1); 144 q_hist2($qx_mc, $w_mc, "Q\\dx\\u [\\A\\u-1\\d]","\\gw [meV]",50,0); 145 mcs_proj($R0,$NP,1, $qx_mc->sum/$npts, $w_mc->sum/$npts, pdl([0,1,3]),pdl([0,3])); 146 147 pgpanl(2,1); 148 q_hist2($qy_mc, $w_mc, "Q\\dy\\u [\\A\\u-1\\d]","\\gw [meV]",50,0); 149 mcs_proj($R0,$NP,0, $qy_mc->sum/$npts, $w_mc->sum/$npts, pdl([0,1,3]),pdl([1,3])); 150 151 pgpanl(3,1); 152 q_hist2($qz_mc, $w_mc, "Q\\dz\\u [\\A\\u-1\\d]","\\gw [meV]",50,0); 153 mcs_proj($R0,$NP,0, $qz_mc->sum/$npts, $w_mc->sum/$npts, pdl([0,2,3]),pdl([2,3])); 154 155 pgpanl(4,1); 156 q_hist2($qx_mc, $qy_mc, "Q\\dx\\u [\\A\\u-1\\d]","Q\\dy\\u [\\A\\u-1\\d]",50,1); 157 mcs_proj($R0,$NP,2, $qx_mc->sum/$npts, $qy_mc->sum/$npts,pdl([0,1,3]),pdl([0,1])); 158 159 160 pgpanl(1,2); 161 my $offset=-1; 162 pgmtxt("t",$offset-0*1.2,.0,0.0,"Bragg (Gaussian) Half Widths"); 163 pgmtxt("t",$offset-1*1.2,0.2,0.0,"\\gDQ\\dx\\u = " . 164 int(2.3548/sqrt($res_mat->at(0,0))*1e4)/1e4 . " \\A\\u-1\\d"); 165 pgmtxt("t",$offset-2*1.2,0.2,0.0,"\\gDQ\\dy\\u = " . 166 int(2.3548/sqrt($res_mat->at(1,1))*1e4)/1e4 . " \\A\\u-1\\d"); 167 pgmtxt("t",$offset-3*1.2,0.2,0.0,"\\gDQ\\dz\\u = " . 168 int(2.3548/sqrt($res_mat->at(2,2))*1e4)/1e4 . " \\A\\u-1\\d"); 169 pgmtxt("t",$offset-4*1.2,0.2,0.0,"\\gD\\gw = " . 170 int(2.3548/sqrt($res_mat->at(3,3))*1e4)/1e4 . " meV"); 171 172 pgmtxt("t",$offset-6*1.2,0.0,0.0,"Resolution matrix [Q\\dx\\u Q\\dy\\u Q\\dz\\u \\gw]:"); 173 $pos = matout($offset-7*1.2, $res_mat); 174 pgmtxt("t",$pos+$offset-0*1.2,.0,0.0,"Covariance matrix [Q\\dx\\u Q\\dy\\u Q\\dz\\u \\gw]:"); 175 $pos = matout($pos+$offset-1*1.2, $C_t); 176 177 pgpanl(2,2); 178 pgmtxt("t",$offset-0*1.2,0.0,0.0,"File: $filename"); 179 my $time=gmtime; 180 pgmtxt("t",$offset-1*1.2,0.0,0.0,"Date: $time"); 181 pgmtxt("t",$offset-2*1.2,0.0,0.0,"X along <Q> in plane"); 182 pgmtxt("t",$offset-3*1.2,0.0,0.0,"Y perp. to X in plane, Z upwards"); 183 184 pgpanl(3,2); 185 186 my $i; 187 my $j=0; 188 my $shift=0.0; 189 pgmtxt("t",$offset,0.0,0.0,"Instrument simulation parameters:"); 190 foreach $i (keys %{$si->{'Params'}}) { 191 $j = $j+1; 192 pgmtxt("t",$offset-$j*1.2,$shift,0.0,$i . " = " . $si->{'Params'}{$i}); 193 if ($j > 20) { $shift = $shift+0.5; $j=0; } 194 } 195 if ($j == 0 && $shift == 0) { pgmtxt("t",$offset-2*1.2,0.0,0.0,"None"); } 196} 197 198sub chol { 199 my ($A) = @_; 200 my ($i,$j,$k,$L,$n,$n2,$li,$lj,$v); 201 202 ($n,$n2) = $A->dims; 203 die "Must be square matrix" unless $n==$n2; 204 $L = zeroes $n,$n; 205 $li = $lj = pdl []; # Handle special case for i=0 206 for($i=0; $i<$n; $i++) { 207 $li = $L->mslice([0,$i-1],[$i]) if $i; 208 $v = $A->at($i,$i) - sum($li*$li); 209 die "Not positive definite" unless $v >= 0; 210 $L->set($i,$i, sqrt($v)); 211 for($j=$i+1; $j<$n; $j++) { 212 $lj = $L->mslice([0,$i-1],[$j]) if $i; 213 $L->set($i,$j, ($A->at($i,$j) - sum($li*$lj))/$L->at($i,$i)); 214 } 215 } 216 return $L; 217} 218 219sub q_hist2 { 220 my ($x,$y,$xl,$yl,$npts,$plot_wedge) = @_; 221 ($xmin,$xmax) = minmax($x); 222 ($ymin,$ymax) = minmax($y); 223 $dx=($xmax-$xmin)/$npts; 224 $dy=($ymax-$ymin)/$npts; 225 my $tr; 226 if (defined(&label_axes)) 227 { $tr = pdl [$xmin + $dx/2, $dx, 0, $ymin + $dy/2, 0, $dy]; } 228 else 229 { $tr = cat $xmin + $dx/2, $dx, pdl(0), $ymin + $dy/2, pdl(0), $dy; } 230 $hxy = histogram2d($x, $y, $dx, $xmin, $npts, $dy, $ymin, $npts); 231 my ($min, $max) = (min($hxy), max($hxy)); 232 if ($min == $max) { 233 if($min == 0) { 234 $max = 1; 235 } else { 236 $min = 0.9*$min; 237 $max = 0.9*$max; 238 } 239 } 240 my $ramp = pdl [[ 0, 1/8, 3/8, 5/8, 7/8, 8/8], 241 [ 0, 0, 0, 1, 1, .5], 242 [ 0, 0, 1, 1, 0, 0], 243 [.5, 1, 1, 0, 0, 0]]; 244 my $numcol = 64; 245 # now do the plottings 246 pgswin($xmin,$xmax,$ymin,$ymax); 247 pgscir(16,16+$numcol-1); 248 ctab $ramp; 249 # If using the black&white postscript driver, swap foreground and 250 # background when doing the image to get more printer-friendly 251 # output. 252 my ($buf, $len); 253 my ($r0, $g0, $b0, $r1, $g1, $b1); 254 pgqinf("TYPE", $buf, $len); 255 if($buf =~ /^V?PS$/i) { 256 pgqcr(0, $r0, $g0, $b0); 257 pgqcr(1, $r1, $g1, $b1); 258 pgscr(0, $r1, $g1, $b1); 259 pgscr(1, $r0, $g0, $b0); 260 } 261 imag $hxy, $min, $max, $tr; 262 if ($plot_wedge) { pgwedg("RI", 0.5, 3.0, $min, $max, ' '); } 263 if($buf =~ /^V?PS$/i) { 264 pgscr(0, $r0, $g0, $b0); 265 pgscr(1, $r1, $g1, $b1); 266 } 267 pglab($xl, $yl,""); 268} 269 270sub matout { 271 my ($pos,$x) = @_; 272 my @lines = split("\n","$x"); 273# shift(@lines);shift(@lines);pop(@lines); 274 for(@lines) { 275 if(m'\[([^]]*)\]') { 276 pgmtxt("t",$pos,0.0,0.0,$1); 277 $pos-= 1.2; 278 } 279 } 280 return $pos; 281} 282 283# The rest of this file is converted from rescal5 matlab code. 284sub rot_elip { 285 my ($a,$b,$phi) = @_; 286 my($n,$x,$y,$s,$c,$th); 287 288 $n = 100; 289 $th = sequence($n+1)/$n*2*$PI; 290 $x = $a*cos($th); 291 $y = $b*sin($th); 292 $c = cos($phi); 293 $s = sin($phi); 294 $th = $x*$c - $y*$s; 295 $y = $x*$s + $y*$c; 296 $x = $th; 297 return ($x,$y); 298} 299 300sub rc_int { 301 my ($i,$r0,$m) = @_; 302 my ($n1,$n2,$r,$sel,$b,$mp,$new); 303 304 ($n1,$n2) = $m->dims; 305 die "Must have square input matrix" unless $n1==$n2; 306 $r = sqrt(2*$PI/$m->at($i,$i))*$r0; 307 $sel = pdl [0..$i-1,$i+1..$n1-1]; 308 $b = $m->slice(",($i)") + $m->slice("($i),"); 309 $b = $b->dice($sel); 310 311 $mp = zeroes $n1-1,$n2-1; 312 if($i > 0) { 313 $mp = $mp->ins($m->mslice([0,$i-1],[0,$i-1]),0,0); 314 } 315 if($i < $n1 - 1) { 316 $mp = $mp->ins($m->mslice([$i+1,$n1-1],[$i+1,$n2-1]),$i,$i); 317 } 318 if($i > 0 && $i < $n1 - 1) { 319 $mp = $mp->ins($m->mslice([0,$i-1],[$i+1,$n2-1]),0,$i); 320 $mp = $mp->ins($m->mslice([$i+1,$n1-1],[0,$i-1]),$i,0); 321 } 322 $new = $mp - 1/(4*$m->at($i,$i))* 323 PDL::Primitive::matmult($b->dummy(0),$b->dummy(1)); 324 return ($r, $new); 325} 326 327 328 329sub mcs_proj { 330 my ($R0,$A,$index,$x0,$y0,$sel1,$sel2) = @_; 331 my($B,$R0P,$MP,$x,$y); 332 333 $B = $A->dice($sel1,$sel1); 334 ($R0P,$MP) = rc_int($index,$R0,$B); 335 ($x,$y) = proj_elip($MP); 336 #poly($x,$y, {COLOUR => RED}); 337 hold; 338 line($x+$x0,$y+$y0,{COLOUR => BLACK}); 339 ($x,$y) = proj_elip($A->dice($sel2,$sel2)); 340 #poly($x,$y, {COLOUR => GREEN}); 341 line($x+$x0,$y+$y0,{COLOUR => BLACK, LINESTYLE => 'DOT-DASH'}); 342 rel; 343} 344 345sub proj_elip { 346 my ($MP) = @_; 347 my ($const,$theta,$S,$MP2,$hwhm_xp,$hwhm_yp,$x,$y); 348 349 $const = 1.17741; 350 $theta = 0.5*atan(2*$MP->at(0,1)/($MP->at(0,0)-$MP->at(1,1))); 351 $S = cat(cat(cos($theta), sin($theta)), 352 cat(-sin($theta), cos($theta))); 353 $MP2 = inner2t($S->transpose,$MP,$S); 354 $hwhm_xp=$const/sqrt($MP2->at(0,0)); 355 $hwhm_yp=$const/sqrt($MP2->at(1,1)); 356 ($x,$y) = rot_elip($hwhm_xp,$hwhm_yp,$theta); 357 return ($x,$y); 358} 359 360# Start of mcresplot program ===================================================== 361 362my $cc; 363my $filename; 364my $interactive=1; 365for($i = 0; $i < @ARGV; $i++) { 366 $_ = $ARGV[$i]; 367 # Options specific to mcplot. 368 if(/^-psc$/ || /^-c$/) { 369 $cc = "c"; $interactive=0; 370 } elsif(/^-ps$/ || /^-p$/) { 371 $cc = "p"; $interactive=0; 372 } elsif(/^-gif$/ || /^-g$/) { 373 $cc = "g"; $interactive=0; 374 } elsif(/^--help$/ || /^-h$/ || /^-v$/) { 375 print "mcresplot [-ps|-psc|-gif|-v] <FILE from Res_monitor>\n"; 376 print " The FILE to be used by mcresplot is generated when using Res_sample\n"; 377 print " at the sample position, and Res_monitor afterwards.\n"; 378 print " Plots the instrument resolution function (projections).\n"; 379 print " When using -ps -psc -gif, the program writes the hardcopy file\n"; 380 print " and then exits.\n"; 381 print "SEE ALSO: mcstas, mcdoc, mcplot, mcrun, mcgui, mcresplot, mcstas2vitess\n"; 382 print "DOC: Please visit http://www.mcstas.org/\n"; 383 exit; 384 } else { 385 $filename = $ARGV[$i]; 386 } 387} 388die "mcresplot <file name from Res_monitor>\n" unless @ARGV; 389 390#read resolution data 391my ($qx,$qy,$qz,$w,$p,$C_t,$res_mat,$size) = read_mcstas_res($filename); 392# now get parameter list (if any) 393my $simulation_info; 394$simulation_info =read_mcstas_info($filename); 395# now either output direct and exit, or plot xwin and wait for exit 396if ($interactive) { 397 print "Type 'P' 'C' or 'G' (in graphics window) for hardcopy, 'Q' to quit.\n"; 398} 399for (;;) { 400 if($cc =~ /[pcg]/i) { # Hardcopy? 401 my $ext="ps"; 402 my $dev = ($cc =~ /c/i) ? "cps" : "ps"; 403 if($cc =~ /g/i) { $dev = "gif"; $ext="gif"; } 404 my $fileout = "$filename.$ext"; 405 plot_mcstas_res($filename, "$fileout/$dev", $qx,$qy,$qz,$w,$p,$C_t,$res_mat,$size,0,$simulation_info); 406 print "Wrote file '$fileout' ($dev)\n"; 407 } 408 if ($interactive == 0) { $cc = "q"; } 409 else { 410 my ($ax,$ay,$cx,$cy) = (0,0,0,0); 411 plot_mcstas_res($filename, "/xwin", $qx,$qy,$qz,$w,$p,$C_t,$res_mat,$size,1,$simulation_info); 412 pgband(0, 0, $ax, $ay, $cx, $cy, $cc); 413 } 414 last if $cc =~ /[xq]/i; 415} 416if (defined(&close_window)) { close_window(); } 417 else { pgclos(); } 418 4191; 420