%%%%%%%
%%%% only_afterglow function, given status, return only glow-flagged
%%%% events' indices
%%%%%%%


define only_afterglow (status)
{
   variable afterglow_bits = (0xF shl 16);   %  bits 16-19
   variable not_afterglow_bits = ~afterglow_bits;
   return (( (status & not_afterglow_bits) == 0) and (status & afterglow_bits));
}


%%%%%%%
%%%%%%%  read columns from the file
%%%%%%%

pn="/nfs/wiwaxia/d1/dph/CXO_Data/Capella/1318";
fn= pn + "/acisf01318_000N002_evt1a.fits";
(s,x,y,p,r,d,m,ml)=fits_read_col(fn,"status","x","y","tg_part","tg_r","tg_d","tg_m","tg_mlam");

connect_points(0);    % don't connect points
p1=open_plot("/xwin"); resize(18);


%%%%%%%
%%%%%%%  % "good" meg +-1st order, broad region of -1st order
%%%%%%%

wmin=-20.;wmax=-6.;
ymin=-2e-3; ymax=2e-3;

lmgood=where( (s == 0) and (p == 2) and (abs(m) == 1)
 and (ml > wmin) and (ml < wmax)
); 

		                         % afterglow in MEG +-1st order
lmglow=where( only_afterglow(s) and (p == 2) and (abs(m) == 1)
  and (ml > wmin) and (ml < wmax)
); 

frac_glow = typecast(length(lmglow), Double_Type)
                  / ( length(lmglow) + length(lmgood) );

xrange(wmin,wmax);
yrange(ymin, ymax);
label("Wavelength [A]", "tg_d [deg]", "Capella/HETGS, ObsID 1318 (repro)");
pointstyle(-1);plot(ml[lmgood],d[lmgood],1);
pointstyle(-5);oplot(ml[lmglow],d[lmglow],2);


% for output:

pgif=dup_plot("Meg_m1.gif/gif",p1);
pointstyle(-1);plot(ml[lmgood],d[lmgood],1);
pointstyle(-5);oplot(ml[lmglow],d[lmglow],2);
close_plot(pgif);
window(p1);



%%%%%%%%%%%%
%%%%%%%%%%%% select Fe XVII 15A -1st order region, "good" events
%%%%%%%%%%%%

wmin=-15.05;
wmax=-14.98;
l15=where( (p == 2) and (ml > wmin) and (ml < wmax) and (m ==-1) and (s == 0));

            % select Fe XVII 15A -1st order region, afterglow-flagged events
l15del=where( (p == 2) and (ml > wmin) and (ml < wmax) and (m ==-1) and only_afterglow(s));


xrange(wmin-0.01,wmax+0.01);
yrange(ymin, ymax);
pointstyle(-1);plot(ml[l15],d[l15],1);
pointstyle(-5);oplot(ml[l15del],d[l15del],2);


% for output:

pgif=dup_plot("Meg_m1_15A.gif/gif",p1);
pointstyle(-1);plot(ml[l15],d[l15],1);
pointstyle(-5);oplot(ml[l15del],d[l15del],2);
close_plot(pgif);
window(p1);

frac_glow_15 = typecast(length(l15del), Double_Type)
                  / ( length(l15del) + length(l15) );



%%%%%%%%%%%%
%%%%%%%%%%%% select Fe XVII 17A -1st order region, "good" events
%%%%%%%%%%%%

wmin=-17.15;
wmax=-17.0;
l17=where( (p == 2) and (ml > wmin) and (ml < wmax) and (m ==-1) and (s == 0));

            % select Fe XVII 17A -1st order region, afterglow-flagged events
l17del=where( (p == 2) and (ml > wmin) and (ml < wmax) and (m ==-1) and only_afterglow(s));

xrange(wmin-0.01,wmax+0.01);
yrange(ymin,ymax);
pointstyle(-1);plot(ml[l17],d[l17],1);
pointstyle(-5);oplot(ml[l17del],d[l17del],2);

frac_glow_17 = typecast(length(l17del), Double_Type)
                  / ( length(l17del) + length(l17) );


% for output:

pgif=dup_plot("Meg_m1_17A.gif/gif",p1);
pointstyle(-1);plot(ml[l17],d[l17],1);
pointstyle(-5);oplot(ml[l17del],d[l17del],2);
close_plot(pgif);
window(p1);


%%%%%%%%%%%%
%%%%%%%%%%%% print glow fractions:
%%%%%%%%%%%%

print("MEG +-1 glow fraction = " + string(frac_glow));
print("MEG -1 15A glow fraction = " + string(frac_glow_15));
print("MEG -1 17A glow fraction = " + string(frac_glow_17));




%%%%%%%%%%%%
%%%%%%%%%%%% Look at bigger field, detector coords
%%%%%%%%%%%%

(s,dx,dy,p,m)=fits_read_col(fn,"status","detx","dety","tg_part","tg_m");


lgood=where( (s == 0)
); 

lglow=where( only_afterglow(s)
); 

limits;
pointstyle(-1);
label("DetX", "DetY", "Capella/HETGS, ObsID 1318 (repro)");

xrange(800,1600);
pointstyle(-1);plot(dx[lgood],dy[lgood],1);
pointstyle(-4);oplot(dx[lglow],dy[lglow],2);


% for output:

pgif=dup_plot("S0_detxy.gif/gif",p1);
pointstyle(-1);plot(dx[lgood],dy[lgood],1);
pointstyle(-4);oplot(dx[lglow],dy[lglow],2);
close_plot(pgif);
window(p1);


%%%%%%%%%%%%
%%%%%%%%%%%%   full field, detxy
%%%%%%%%%%%%

xrange;
%pointstyle(-1);plot(dx[lgood],dy[lgood],1);
%pointstyle(-4);oplot(dx[lglow],dy[lglow],2);

pointstyle(-1);plot(dx[lglow],dy[lglow],1);


% for output:

pgif=dup_plot("Full_Glow_detxy_.gif/gif",p1);
pointstyle(-1);plot(dx[lglow],dy[lglow],1);
close_plot(pgif);
window(p1);


%%%%%%%%%%%%
%%%%%%%%%%%%   obsid 105:
%%%%%%%%%%%%

pn = "/nfs/cxc/h1/dph/d1/CXO_Data/HETG_DIRS/obsid_105/primary";
fn = pn + "/acisf00105_001N002_evt1a.fits";

(s,x,y,p,r,d,m,ml)=fits_read_col(fn,"status","x","y","tg_part","tg_r","tg_d","tg_m","tg_mlam");

wmin=-30.;wmax=-30.;
ymin=-2e-3; ymax=2e-3;

lmgood=where( (s == 0) and (p == 2) and (abs(m) == 1)
 and (ml > wmin) and (ml < wmax)
); 

		                         % afterglow in MEG +-1st order
lmglow=where( only_afterglow(s) and (p == 2) and (abs(m) == 1)
  and (ml > wmin) and (ml < wmax)
); 

frac_glow = typecast(length(lmglow), Double_Type)
                  / ( length(lmglow) + length(lmgood) );

frac_glow;   % echo to screen (if interactive)

xrange(wmin,wmax);
yrange(ymin, ymax);
label("Wavelength [A]", "tg_d [deg]", "HETGS, ObsID 105 (repro)");
pointstyle(-1);plot(ml[lmgood],d[lmgood],1);
pointstyle(-1);oplot(ml[lmglow],d[lmglow],2);


% for output:

pgif=dup_plot("Meg_m1_105.gif/gif",p1);
pointstyle(-1);plot(ml[lmgood],d[lmgood],1);
pointstyle(-1);oplot(ml[lmglow],d[lmglow],2);
close_plot(pgif);
window(p1);



%%%%%%%%%%%%
%%%%%%%%%%%%  HEG
%%%%%%%%%%%%

wmin=-15.;wmax=15.;
ymin=-2e-3; ymax=2e-3;

lhgood=where( (s == 0) and (p == 1) and (abs(m) == 1)
 and (ml > wmin) and (ml < wmax)
); 

		                         % afterglow in HEG +-1st order
lhglow=where( only_afterglow(s) and (p == 2) and (abs(m) == 1)
  and (ml > wmin) and (ml < wmax)
); 

frac_glow = typecast(length(lhglow), Double_Type)
                  / ( length(lhglow) + length(lhgood) );

frac_glow;   % echo to screen (if interactive)

xrange(wmin,wmax);
yrange(ymin, ymax);
label("Wavelength [A]", "tg_d [deg]", "HETGS, ObsID 105 HEG (repro)");
pointstyle(-1);plot(ml[lhgood],d[lhgood],1);
pointstyle(-1);oplot(ml[lhglow],d[lhglow],2);


% for output:

pgif=dup_plot("Heg_m1_105.gif/gif",p1);
pointstyle(-1);plot(ml[lhgood],d[lhgood],1);
pointstyle(-1);oplot(ml[lhglow],d[lhglow],2);
close_plot(pgif);
window(p1);