% % file: do_streak_fit.sl in $SCI/TGCat/Streak_fwhm % % 10/17/08-dd Get and fit the streak events from an obsid. % 9/14+/12-dd Updated to use the evt2.gz file directly from TGCat. % Modified findzo_plus to use eneries < 3 keV always. % % Usage: % - edit this file to select the desired obsid, then: % - isis> .source do_streak_fit % % To update the streak-vs-time plot: % - edit focus_trends_v3.dat and enter the FWHM etc values from % the do_streak_fit output plot. % - isis> .source make_streak_plot % % - - - Creating findzo_plus.sl - - - % % Add get_streak_events to findzo by doing: % ! cat findzo-0.9.5.sl get_streak_events.sl > findzo_plus.sl % where get_streak_events.sl is just the "define findzo ..." routine % but end it prematurely with: % ... % % And draw a first region box to isolate transfer streak events % % without zero order and with energy filtering. % % % l=where(xp > (gx-boxw_strc/2) and xp < (gx+boxw_strc/2) and % yp > (gy-boxl_strc/2) and yp < (gy+boxl_strc/2) and % x > 0.0 and y > 0.0 and % hypot(xp-gx,yp-gy)>rad_strc and (en > n_eng_lo and en < n_eng_hi)); % % % offset w.r.t. zo for the returned event locations: % return xp[l]-gx, yp[l]-gy; % } % % Also: % - Keep the same size regions even if array width is smaller, % to keep scatter hallow out of streak histogram. % Comment out these lines in findzo_plus.sl : % %%% rad_strc *= 0.4; % %%% boxw_strc *= 0.5; % % - Use only energies less than 3 keV by changing: % %%variable n_eng_hi = 4000; % energy filtering (high threshold) % variable n_eng_hi = 3000; % energy filtering (high threshold) % % % - - - % Just do this once if calling do_streak_fit from a loop #iftrue % Things we'll use besides isis: require("histogram"); % findzo extended to include "get_streak_events" .load findzo_plus ; % Set the location of the tgcat dir: %%variable tgcat_dir = "/data/tgcat"; %%variable tgcat_dir = "/nfs/spectra/d2/dd/Streak_data"; % Directory where TGCat evt2 files are available: % as tgcat_dir+"/tgcid__evt2.gz" variable tgcat_dir = "/nfs/notpirx/d1/tgcat/arcA/evt2"; variable xc, yc, xs, ys; variable obsid_int, obsid_str, obs_to_do, iobs=0; variable tgid_int, tgid_str, grat_str = "m"; % read in the TGCat info - mapping tgid to obsid: variable tgc_obsids, tgc_tgids; (tgc_tgids, tgc_obsids) = readcol("tgcat_findzo_HETG.txt", [1,2]); % ------------------------------------ % Set the obisd to process: obsid_int = 15484; % ------------------------------------ #endif % make the string too obsid_str=string(obsid_int); % Lookup the tgid for this obsid: tgid_int = where(obsid_int == tgc_obsids); % Check if the obsid is found (if not, go to bottom/else section) if (length(tgid_int) > 0) { tgid_int = tgc_tgids[ tgid_int[-1] ]; tgid_str=string(tgid_int); % Comment out all the work to just check which obsids are in the lookup %%#iffalse message(""); message(" - - - - - - - "); message(" ObsId "+obsid_str+" aka TGCatId "+tgid_str); message(""); % clear out previous (findzo) plots (if re-running it) close_plot; close_plot; close_plot; % find the center with findzo using MEG (or LEG): (xc,yc) = findzo(tgcat_dir+"/tgcid_"+tgid_str+"_evt2.gz",grat_str); % now get the rotated streak events: (xs, ys) = get_streak_events(tgcat_dir+"/tgcid_"+tgid_str+"_evt2.gz",grat_str,4,xc,yc); % plot them window(1); xrange;yrange; connect_points(0); plot(xs,ys); connect_points(1); % Load a real MEG-1 spectrum (to setup a dataset): delete_data(all_data); ()=load_data("obs_14239_pha2",9); %% plot_data_counts(1); list_data; % Get the data counts strucuture: (bin_lo, bin_hi, value, err) variable dstruct; dstruct = get_data_counts(1); % create a fake counts array based on the streak values converted % to wavelength, centered at 15.0 A. % The units are sky pixels (0.492" each) variable streak_pha; streak_pha = hist1d (xs+15.0, dstruct.bin_lo); put_data_counts (1, dstruct.bin_lo, dstruct.bin_hi, streak_pha, sqrt(streak_pha)); % Now group, notice, and fit... % MEG pha has binning of 0.005 A % so group x40 for 0.2 sky-pixel bins group_data(1,40); % limit the fit to +/- 2.0 pixels: xnotice(1,13.0,17.0); xrange(10.0,20.0); xlin; ylog; %%()=open_plot("/XWIN"); ()=open_plot("/XSERV"); plot_data_counts(1); fit_fun("gauss(1)+poly(1)"); % guess the number of event per exposure time: set_par("gauss(1).area",0.5*length(xs)/get_data_exposure(1), 0); set_par("gauss(1).center",15.0, 0, 14.0,16.0); set_par("gauss(1).sigma",(38.0/24.0)/2.35, 0); % background ~ constant with exposure set_par("poly(1).a0",0.0005, 0); set_par("poly(1).a1",0.0, 1); set_par("poly(1).a2",0.0, 1); set_fit_method("lmdif"); %% ()=eval_counts; %% rplot_counts(1); ()=fit_counts; message(" Parameters for the Streak fit:"); message(""); list_par; % Get the 90% confidence limits: variable conflo, confhi; (conflo,confhi)=conf("gauss(1).sigma",1,0.001); message(" Obsid "+obsid_str+" :"); message(""); message(" Streak area (counts): "+string(get_par(1)*get_data_exposure(1))); message(""); message(" FWHM of the Streak fit = " + string(24.0*(2.35*get_par(3))) + " um" ); message(" ( "+string(24.0*(2.35*conflo)) + " to "+string(24.0*(2.35*confhi))+", 90% confidence; +/- "+ string(0.5*24.0*(2.35*(confhi-conflo)))+" )"); message(""); % Save it to .gif %%variable GIFwin=open_plot("streak_"+obsid_str+"_plot.gif/GIF"); % or if get GIF errors like: % PGPLOT, Writing new GIF image as: streak_12089_plot.gif P % Abort trap % [Dan-50] % then use .ps: variable GIFwin=open_plot("streak_"+obsid_str+"_plot.ps/CPS"); variable titstr="Obsid "+obsid_str+" : FWHM of Streak = " + string(24.0*(2.35*get_par(3))) + " um" + " +/- "+string(0.5*24.0*(2.35*(confhi-conflo)))+" (90%)"; title(titstr); errorbars(1); rplot_counts(1); close_plot(GIFwin); % Show the final plot on screen too: title(titstr); errorbars(1); rplot_counts(1); % End of "Comment out all work" above %%#endif message(" - - - - - - - "); message(""); } % if tgid_int is valid else { message(""); message(" - - - - - - - "); message(" obsid "+obsid_str+" does not have a TGCat findzo extraction."); message(" - - - - - - - "); message(""); }