% % file: do_streak_fit.sl % % Get and fit the streak events from an obsid. 10/17/08 - dd % % Usage: % - edit this file to select the desired obsid, then: % - hydra> .source do_streak_fit % % To update the streak plot: % - edit values into focus_trends_v2.dat % - isis> .source make_streak_plot % % Note that the default findzo energy range of 0.5 to 4 keV is used. % % - - - % Add get_streak_events to findzo by doing: % ! cat findzo.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; % % - - - .load findzo_plus % Set the location of the tgcat dir: variable tgcat_dir = "/data/tgcat"; variable xc, yc, xs, ys; variable obsid_str, obs_to_do, iobs=0; variable grat_str = "m"; % ----- Loop over a bunch... %% obs_to_do=[702,1921,2741,2739,2708]; %%obs_to_do=[3814,3823,4572,4584,5040,6601,7449,7451]; %%obs_to_do=[1010,1103,1318,2583,3674,5955,6471]; %%obs_to_do=[337,1014,1705,3167,3706,3708,5173,7291]; %%obs_to_do=[105,106,1252,1451]; %%obs_to_do=[706,609,4552,7280]; %%obs_to_do=[1928,4574,15,425,3504,]; %%_for iobs (0,length(obs_to_do)-1,1) %%{ %%obsid_str = string(obs_to_do[iobs]); % - - - - - Or just do one % Select the obsid from available ones: % 1235, 62538 %obsid_str="1252"; %obsid_str="106"; %obsid_str="1103"; %obsid_str="1318"; % 16, 1450 %obsid_str="1451"; %obsid_str="105"; %obsid_str="337"; % 1456, 457, 1481 %obsid_str="609"; % 1456b, 459, 1511, 605, 639, 373, 837, 92 %obsid_str="706"; % 57, 705, 629, 640, 604 %obsid_str="425"; %obsid_str="660"; % --- Large errors, not used, Not on Mac --- %obsid_str="1706"; % --- Large errors, not used, Not on Mac --- %obsid_str="1714"; % --- Large errors, not used, Not on Mac --- %obsid_str="1705"; %obsid_str="704"; % --- ZO windowed out, not used, Not on Mac --- %obsid_str="1712"; % --- 3C 273, No HETG -- not used. Not on Mac--- %obsid_str="15"; %obsid_str="646"; % --- Large errors, not used, Not on Mac --- %obsid_str="702"; %obsid_str="17"; % --- Large errors, not used, Not on Mac --- %obsid_str="632"; % --- Large errors, not used, Not on Mac --- %obsid_str="1943"; % --- Large errors, not used, Not on Mac --- %obsid_str="1928"; %obsid_str="1010"; %obsid_str="1019"; % --- Large errors, not used, Not on Mac --- %obsid_str="2429"; % --- Large errors, not used, Not on Mac --- %obsid_str="1014"; %obsid_str="2463"; % --- Large errors, not used, Not on Mac --- %obsid_str="1921"; %obsid_str="3167"; %obsid_str="2741"; %obsid_str="2583"; %obsid_str="3662"; % --- Her X-1 No HETG -- not used. Not on Mac--- %obsid_str="2739"; %obsid_str="2708"; %obsid_str="3708"; %obsid_str="3706"; %obsid_str="4420"; % very high counts %obsid_str="3814"; %obsid_str="3504"; %obsid_str="4439"; % --- Large errors, not used, Not on Mac --- %obsid_str="4430"; % --- Large errors, not used, Not on Mac --- %obsid_str="3567"; % --- Large errors, not used, Not on Mac --- %obsid_str="3824"; % --- Large errors, not used, Not on Mac --- %obsid_str="3568"; % --- Large errors, not used, Not on Mac --- %obsid_str="3674"; %obsid_str="3823"; %obsid_str="4574"; %%%obsid_str="4572"; % skip: double source ~ aligned by roll %obsid_str="4760"; % --- Large errors, not used, Not on Mac --- %obsid_str-4U1624; % --- Large errors, not used, Not on Mac --- %obsid_str="5169"; % --- Large errors, not used, Not on Mac --- %obsid_str="4584"; %obsid_str="4552"; %obsid_str="5040"; %obsid_str="5173"; %obsid_str="5955"; %obsid_str="6088"; % --- Large errors, not used, Not on Mac --- %obsid_str="5514"; % --- Large errors, not used, Not on Mac --- %obsid_str="6187"; % --- Large errors, not used, Not on Mac --- %obsid_str="6601"; %obsid_str="7280"; %obsid_str="6471"; %obsid_str="6926"; % --- Large errors, not used, Not on Mac --- %obsid_str="7291"; %obsid_str="6613"; % --- Large errors, not used, Not on Mac --- %obsid_str="7449"; %obsid_str="7451"; %obsid_str="8380"; % --- Large errors, not used, Not on Mac --- %obsid_str="8426"; % --- Large errors, not used, Not on Mac --- %obsid_str="7485"; %obsid_str="7511"; %obsid_str="8525"; %obsid_str="9638"; %obsid_str="9703"; %obsid_str="9076"; %obsid_str="9858"; %obsid_str="9705"; %obsid_str="9712"; %obsid_str="10659"; % 4U 1957+11 Y Offset = 0.0 arc min %obsid_str="10661"; % " %obsid_str="10695"; % 4U1556-60 Y Offset = 0.0 arc min %obsid_str="10857"; % 4U 0614+091 Y Offset = 0.0 arc min %obsid_str="10679"; % EV Lac --- very low counts in streak --- %obsid_str="10599"; % Capella: Y Offset = 0.33 arc min %%%% 9947 141.3 HD 110432 Jose T. 2009-11-13 %obsid_str="11931"; % Capella: Y Offset = 0.33 arc min % %obsid_str="11074"; % LMC X-1, evt2 E < 3keV YOffset = 0.0 %obsid_str="12072"; % LMC X-1, evt2 E < 3keV YOffset = 0.0 %obsid_str="12069"; % LMC X-1, evt2 E < 3keV YOffset = 0.0 %obsid_str="12070"; % LMC X-1, evt2 E < 3keV YOffset = 0.0 %obsid_str="11044"; % very bright Cyg X-1, evt2 E < 3keV YOffset = 0.0 %obsid_str="11058"; % 4U 1626-67 (Schulz is Obs) [2010-01-14] YOffset = 0.0 %obsid_str="11987"; % LMC X-1, evt2 E < 3keV YOffset = 0.0 %obsid_str="12089"; % LMC X-1, evt2 E < 3keV YOffset = 0.0 %%%%obsid_str="12090"; % LMC X-1, evt2 E < 3keV 13.2 2010-02-26 YOffset = 0.0 % %obsid_str="10663"; % Mkn 421 YOffset = 0.167 %obsid_str="10670"; % Mkn 421 YOffset = 0.167 % %%%% HD 42054 (11021) Testa [126.2, 2010-06-21] %%%% 4U 1702-429 (11045) GO [2010-06-29] YOffset = 1.0* %%%% GX 9+9 (11072) GO 2010-07-13 YOffset = 0.0 %%%% Cyg X-1 (12313,'14) GO 2010-07-22 YOffset = 0.0 %%%% GX 17+2 (11088,'888) GO 2010-07-25 YOffset = 0.0 % %obsid_str="11815"; % GX 13+1, evt2 E < 3keV YOffset = 0.0 %obsid_str="11816"; % GX 13+1, evt2 E < 3keV YOffset = 0.0 %obsid_str="11814"; % GX 13+1, evt2 E < 3keV YOffset = 0.0 obsid_str="11817"; % GX 13+1, evt2 E < 3keV YOffset = 0.0 % %%%% Capella 13089 2010-11-29 *** Y Offset = 0.33 arc min *** % - - - - - % LETG-ACIS for fun: %grat_str="l"; %obsid_str="8319"; % - - - - - % clear out previous plots close_plot; close_plot; close_plot; % find the center with findzo using MEG (or LEG): (xc,yc) = findzo(tgcat_dir+"/obs_"+obsid_str+"/evt2",grat_str); % now get the rotated streak events: (xs, ys) = get_streak_events(tgcat_dir+"/obs_"+obsid_str+"/evt2",grat_str,4,xc,yc); % plot them window(1); xrange;yrange; connect_points(0); plot(xs,ys); connect_points(1); % Load the real MEG-1 data: delete_data(all_data); if (grat_str == "l") { ()=load_data(tgcat_dir+"/obs_"+obsid_str+"/pha2",4); } else { ()=load_data(tgcat_dir+"/obs_"+obsid_str+"/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... if (grat_str == "l") { % LEG group_data(1,20); xnotice(1,12.0,18.0); } else { % 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); rplot_counts(1); close_plot(GIFwin); % Show the final plot on screen too: title(titstr); rplot_counts(1); % Loop over a bunch end %%}