Re: displaying pha file regions

In ascds_help, you wrote:
>Hello,

>	I'm investigating a user question about how to display the binning regions 
>attached to pha2 files.  (I have limited experience, but am drawing
>some from dph's http://space.mit.edu/ASC/analysis/AGfCHRS/AGfCHRS.html)

>	This particular pha2 (found in
>	/data/ciao_demo/threads/acishetg/1451/primary)
> has 36 regions in the attached "REGION" block.  Can ds9 be used to display 
>them?  Or can I dmcopy a subset of the columns to another file and display that?
>	Here's a list of the columns in the region block:
>

Here are two versions of a slang function to plot the region.  One is
for chips, the other isis.  (There are currently some parser problems
in 2.2.1, and some isis-ciao FITS reading incompatibilities.).  

The file comments give example usage.
  
You have to import isis, because I used isis/pgplot functions to
display, but these could be changed by someone motivated to make it a
pure-chips function.  Since there is as yet no version-variable
variable to test, I simply commented out mutually exclusive sections
of the code in each file, and put in appropriate comments.

I also don't know if I can do conditional plots in chips:

  if (condition) curve x x y y

whereas in isis/pgplot I can: 

  if (condition) oplot(x,y); 

which I used to conditionally overplot background regions (they might
not have been extracted).

I did some tests on acis/hetg, acis/hetg/cc, hrcs/letg (w/ bowtie), w/
and w/out background regions.  I tested in the released ciao 2.2.1 on
debian linux. 

The following should be put into two files:

plt_tgphareg_ciao.sl
plt_tgphareg_isis.sl

Then you
chips
chips> import("isis");
chips> evalfile("plt_tgphareg_ciao.sl");
chips> plt_tgphareg( pha_file, evt_file, part, order);


I didn't do any argument or error checking, so you can break this by
giving it incompatible part and orders (e.g., part=3 only has order -1
or +1).  And I'm sure there are other ways to make the code more
elegant (like the if blocks on shape might do better w/ a case/switch
block).

Note that since this is pgplot output, subsequent chips commands,
like undo or print, won't work.


--Dave


%; Time-stamp: <2002-01-31 14:40:29 dph> 
%; MIT Directory: ~dph/libisis
%; File: plt_tg_phareg_ciao.sl
%; Author: D. Huenemoerder
%; Original version: 2002.01.31
%;====================================================================
% version: 1.0
%
% purpose: quick and dirty plot of events in diffraction coords, w/
% binning region overlayed.
%
%

define plt_tgphareg( f_pha, f_evt, part, order )
{

  variable tg_lam, tg_d, tg_part, tg_m, rowid, shape, r ;
  variable s_xverts, s_yverts;         % source region vertices
  variable mpbu, mpbd, bu_xverts, bu_yverts, bd_xverts, bd_yverts;
  variable ymin, ymax;

%+
% EXAMPLE:
%
% chips> import("isis");
% chips> evalfile("plt_tgphareg.sl");
% chips> f_pha="Pha/acisf00007_005N001_pha2.fits";
% chips> f_evt="Evt/acisf00007_005N001_evt1a.fits";
% chips> order = -1;
% chips> part = 2;
% chips> plt_tgphareg( f_pha, f_evt, part, order );
%
%-
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+ test cases:
% evalfile("/nfs/cxc/h1/dph/libisis/plt_tgphareg.sl");

% f_pha="Pha/acisf00007_005N001_pha2.fits";
% f_evt="Evt/acisf00007_005N001_evt1a.fits";
% order = -1;
% part = 2;
% plt_tgphareg( f_pha, f_evt, part, order);

% order=-1;
% part=3;
% f_pha="../../CXC/Test_ciao_2/LETGS_2.2/Pha/hrcf01715_002N002_pha2.fits";
% f_evt="../../CXC/Test_ciao_2/LETGS_2.2/Evt/hrcf01715_002N002_evt1a.fits";
% plt_tgphareg( f_pha, f_evt, part, order);


% Generic chips/slang:

% read the pha file region extension
%
%  variable f_reg = sprintf("%s[REGION]", f_pha);

%%%% incompatibilities:
% in isis, next line:
%  variable s = fits_read_table( strcat(f_pha, "[REGION]" ) );

% but in ciao 2.2.1 chips:
  variable s = readfile(strcat(f_pha,"[REGION]"));
%%%%

%%% worse incompatibility: 
%%% isis lowercases columns.  ciao varmm doesn't.
%%% put cols into variable until a version check can be done...

%% for isis:
%   tg_lam = s.tg_lam;
%   tg_d = s.tg_d;
%   tg_part = s.tg_part;
%   tg_m = s.tg_m;
%   rowid = s.rowid;
%   shape = s.shape;
%   r = s.r;
%%
%% for chips:
   tg_lam = s.TG_LAM;
   tg_d = s.TG_D;
   tg_part = s.TG_PART;
   tg_m = s.TG_M;
   rowid = s.ROWID;
   shape = s.SHAPE;
   r = s.R;
   vmessage("Region tags parsed...\n");
%%%%

  variable mps = (where( (tg_part == part) and (tg_m == order) and (rowid == "SOURCE") ))[0];
  mpbu = where( (tg_part == part) and (tg_m == order) and (rowid == "BACKGROUND_UP") );
  mpbd = where( (tg_part == part) and (tg_m == order) and (rowid == "BACKGROUND_DOWN") );
  vmessage("Found indices for part and order for source and bkg.\n");


  if (shape[mps] == "box")
  {
    s_xverts = [-0.5, 0.5, 0.5, -0.5, -0.5] * (r[mps,0])[0] + (tg_lam[mps])[0] ;
    s_yverts = [-0.5, -0.5, 0.5, 0.5, -0.5] * (r[mps,1])[0] + (tg_d[mps])[0] ;
  }
  if (shape[mps] == "polygon")
  {
    s_xverts = tg_lam[mps,*];
    s_yverts = tg_d[mps,*];
  }
  
  if (length(mpbu) > 0)   % we have a background_up region
  {
    if (shape[mpbu][0] == "box")
    {
      bu_xverts = [-0.5, 0.5, 0.5, -0.5, -0.5] * (r[mpbu,0])[0] + (tg_lam[mpbu])[0] ;
      bu_yverts = [-0.5, -0.5, 0.5, 0.5, -0.5] * (r[mpbu,1])[0] + (tg_d[mpbu])[0] ;
    }
    if (shape[mpbu][0] == "polygon")
    {
      bu_xverts = tg_lam[mpbu[0],*];
      bu_yverts = tg_d[mpbu[0],*];
    }
  }
  
  vmessage("Parsed background_up region.\n");

  if (length(mpbd) > 0)   % we have a background_down region
  {
    if (shape[mpbu[0]] == "box")
    {
      bd_xverts = [-0.5, 0.5, 0.5, -0.5, -0.5] * (r[mpbd,0])[0] + (tg_lam[mpbd])[0] ;
      bd_yverts = [-0.5, -0.5, 0.5, 0.5, -0.5] * (r[mpbd,1])[0] + (tg_d[mpbd])[0] ;
    }
    if (shape[mpbd[0]] == "polygon")
    {
      bd_xverts = tg_lam[mpbd[0],*];
      bd_yverts = tg_d[mpbd[0],*];
    }
  }  

  vmessage("Parsed background_down region.\n");

  % read the events:

%%%%% incompatibilities:
%%% in isis:
%  (tg_lam,tg_d,tg_part,tg_m) = fits_read_col(f_evt, "tg_lam", "tg_d", "tg_part", "tg_m");

%% in ciao 2.2.1 chips:
  variable e = readbintab(strcat(f_evt,"[cols tg_lam,tg_d,tg_part,tg_m]"));
  tg_lam = e.tg_lam;
  tg_d = e.tg_d;
  tg_part = e.tg_part;
  tg_m = e.tg_m;
%%%%%

  % find the events for the given order and part.
  %  part mapping:  1 => HEG;  2 => MEG;  3 => LEG.
  %  Note: order for LETG/HRC is -1 or +1.

  variable l_evt = where ( (tg_part == part ) and (tg_m == order) );
  
  vmessage("Read events; set variables. Filtered on part and order.\n");

  
% isis plot commands;

%  (could be converted to chips equivalents, if desired;
%   I plotting indexed arrays, or doing conditional plots.  Does this
%   work in chips-slang?)

  vmessage("Plotting w/ isis/pgplot functions...\n");

  variable  p1=open_plot("/xwin"); resize(18,0.6);

  xrange(0, max( tg_lam[l_evt] )*1.2);

  ymax = max( tg_d[l_evt] );
  ymin = min( tg_d[l_evt] );

  if (length(mpbu) > 0)
  {
    ymax = max( bu_yverts );
  }

  if (length(mpbd) > 0)
  {
    ymin = min( bd_yverts );
  }

  yrange( ymin*1.1, ymax*1.1);
  
  pointstyle(-1);connect_points(0);
  plot( tg_lam[l_evt], tg_d[l_evt], 1);

  connect_points(1);
  if (length(mpbu) > 0)  oplot(bu_xverts,bu_yverts,5);
  if (length(mpbd) > 0)  oplot(bd_xverts,bd_yverts,6);

  oplot(s_xverts,s_yverts,2);

}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%; Time-stamp: <2002-01-31 14:39:39 dph> 
%; MIT Directory: ~dph/libisis
%; File: plt_tg_phareg_isis.sl
%; Author: D. Huenemoerder
%; Original version: 2002.01.31
%;====================================================================
% version: 1.0
%
% purpose: quick and dirty plot of events in diffraction coords, w/
% binning region overlayed.
%
%

define plt_tgphareg( f_pha, f_evt, part, order )
{

  variable tg_lam, tg_d, tg_part, tg_m, rowid, shape, r ;
  variable s_xverts, s_yverts;         % source region vertices
  variable mpbu, mpbd, bu_xverts, bu_yverts, bd_xverts, bd_yverts;
  variable ymin, ymax;

%+
% EXAMPLE:
%
% chips> import("isis");
% chips> evalfile("plt_tgphareg_ciao.sl");   % must be in path, or change to full path.
% chips> f_pha="Pha/acisf00007_005N001_pha2.fits";   % pha Type II file
% chips> f_evt="Evt/acisf00007_005N001_evt1a.fits";  % event file
% chips> order = -1;   % diffraction order
% chips> part = 2;     % MEG
% chips> plt_tgphareg( f_pha, f_evt, part, order );
%
%-
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+ test cases:
% evalfile("/nfs/cxc/h1/dph/libisis/plt_tgphareg.sl");

% f_pha="Pha/acisf00007_005N001_pha2.fits";
% f_evt="Evt/acisf00007_005N001_evt1a.fits";
% order = -1;
% part = 2;
% plt_tgphareg( f_pha, f_evt, part, order);

% order=-1;
% part=3;
% f_pha="../../CXC/Test_ciao_2/LETGS_2.2/Pha/hrcf01715_002N002_pha2.fits";
% f_evt="../../CXC/Test_ciao_2/LETGS_2.2/Evt/hrcf01715_002N002_evt1a.fits";
% plt_tgphareg( f_pha, f_evt, part, order);


% Generic chips/slang:

% read the pha file region extension
%
%  variable f_reg = sprintf("%s[REGION]", f_pha);

%%%% incompatibilities:
% in isis, next line:
  variable s = fits_read_table( strcat(f_pha, "[REGION]" ) );

% but in ciao 2.2.1 chips:
%  variable s = readfile(strcat(f_pha,"[REGION]"));
%%%%

%%% worse incompatibility: 
%%% isis lowercases columns.  ciao varmm doesn't.
%%% put cols into variable until a version check can be done...

%% for isis:
   tg_lam = s.tg_lam;
   tg_d = s.tg_d;
   tg_part = s.tg_part;
   tg_m = s.tg_m;
   rowid = s.rowid;
   shape = s.shape;
   r = s.r;
%%
%% for chips:
%   tg_lam = s.TG_LAM;
%   tg_d = s.TG_D;
%   tg_part = s.TG_PART;
%   tg_m = s.TG_M;
%   rowid = s.ROWID;
%   shape = s.SHAPE;
%   r = s.R;
   vmessage("Region tags parsed...\n");
%%%%

  variable mps = (where( (tg_part == part) and (tg_m == order) and (rowid == "SOURCE") ))[0];
  mpbu = where( (tg_part == part) and (tg_m == order) and (rowid == "BACKGROUND_UP") );
  mpbd = where( (tg_part == part) and (tg_m == order) and (rowid == "BACKGROUND_DOWN") );
  vmessage("Found indices for part and order for source and bkg.\n");

  if (shape[mps] == "box")
  {
    s_xverts = [-0.5, 0.5, 0.5, -0.5, -0.5] * (r[mps,0])[0] + (tg_lam[mps])[0] ;
    s_yverts = [-0.5, -0.5, 0.5, 0.5, -0.5] * (r[mps,1])[0] + (tg_d[mps])[0] ;
  }
  if (shape[mps] == "polygon")
  {
    s_xverts = tg_lam[mps,*];
    s_yverts = tg_d[mps,*];
  }
  
  if (length(mpbu) > 0)   % we have a background_up region
  {
    if (shape[mpbu][0] == "box")
    {
      bu_xverts = [-0.5, 0.5, 0.5, -0.5, -0.5] * (r[mpbu,0])[0] + (tg_lam[mpbu])[0] ;
      bu_yverts = [-0.5, -0.5, 0.5, 0.5, -0.5] * (r[mpbu,1])[0] + (tg_d[mpbu])[0] ;
    }
    if (shape[mpbu][0] == "polygon")
    {
      bu_xverts = tg_lam[mpbu[0],*];
      bu_yverts = tg_d[mpbu[0],*];
    }
  }
  
  vmessage("Parsed background_up region.\n");

  if (length(mpbd) > 0)   % we have a background_down region
  {
    if (shape[mpbu[0]] == "box")
    {
      bd_xverts = [-0.5, 0.5, 0.5, -0.5, -0.5] * (r[mpbd,0])[0] + (tg_lam[mpbd])[0] ;
      bd_yverts = [-0.5, -0.5, 0.5, 0.5, -0.5] * (r[mpbd,1])[0] + (tg_d[mpbd])[0] ;
    }
    if (shape[mpbd[0]] == "polygon")
    {
      bd_xverts = tg_lam[mpbd[0],*];
      bd_yverts = tg_d[mpbd[0],*];
    }
  }  

  vmessage("Parsed background_down region.\n");

  % read the events:

%%%%% incompatibilities:
%%% in isis:
  (tg_lam,tg_d,tg_part,tg_m) = fits_read_col(f_evt, "tg_lam", "tg_d", "tg_part", "tg_m");

%% in ciao 2.2.1 chips:
%  variable e = readbintab(strcat(f_evt,"[cols tg_lam,tg_d,tg_part,tg_m]"));
%  tg_lam = e.tg_lam;
%  tg_d = e.tg_d;
%  tg_part = e.tg_part;
%  tg_m = e.tg_m;
%%%%%

  % find the events for the given order and part.
  %  part mapping:  1 => HEG;  2 => MEG;  3 => LEG.
  %  Note: order for LETG/HRC is -1 or +1.

  variable l_evt = where ( (tg_part == part ) and (tg_m == order) );
  
  vmessage("Read events; set variables. Filtered on part and order.\n");

  
% isis plot commands;
%  (could be converted to chips equivalents, if desired;
%   I plotting indexed arrays, or doing conditional plots.  Does this
%   work in chips-slang?)

  vmessage("Plotting w/ isis/pgplot functions...\n");

  variable  p1=open_plot("/xwin"); resize(18,0.6);

  xrange(0, max( tg_lam[l_evt] )*1.2);

  ymax = max( tg_d[l_evt] );
  ymin = min( tg_d[l_evt] );

  if (length(mpbu) > 0)
  {
    ymax = max( bu_yverts );
  }

  if (length(mpbd) > 0)
  {
    ymin = min( bd_yverts );
  }

  yrange( ymin*1.1, ymax*1.1);
  
  pointstyle(-1);connect_points(0);
  plot( tg_lam[l_evt], tg_d[l_evt], 1);

  connect_points(1);
  if (length(mpbu) > 0)  oplot(bu_xverts,bu_yverts,5);
  if (length(mpbd) > 0)  oplot(bd_xverts,bd_yverts,6);

  oplot(s_xverts,s_yverts,2);

}
----
You received this message because you are
subscribed to the isis-users list.
To unsubscribe, send a message to
isis-users-request_at_email.domain.hiddenwith the first line of the message as:
unsubscribe