;+ ;NAME: ; hxt_mcaccum ;PURPOSE: ; plots hxt data, prompts for start and end times for data ; accumulation and background, and returns an array with the ; accumulated counts, for each interval and channel. Then ; this will generate N sets of observations and sigmas ; using poisson statistics. Restricted to 1 time interval. ; Totals intervals before choosing samples, the old hxt_mcaccum ; samples before totaling. ;CALLING SEQUENCE: ; Hxt_mcaccum, index, data, index1, hxi_index, counts, sigma, bcounts, $ ; bsigma, trials=trials, img_trange=img_trange, $ ; bck_trange=bck_trange, channel=channel, $ ; accum_cnts=accum_cnts, ch_accum=ch_accum, $ ; serr_only=serr_only, norm_error=norm_error, $ ; no_serr=no_serr, st_times = st_times, $ ; int_times = int_times ;INPUT: ; index=data info from hda file ; data=counts in each subcollimator for each channel ;KEYWORDS: ; trials = no. of MC trials for this interval. Default is 100 ; img_trange = a time range for the images, any yohkoh format ; bck_range = a time_range for the background accumulation ; channel= a channel number for the plots, default is to prompt the user ; serr_only = if set, only use the systematic error term to ; obtain samples, the default is to add poisson ; variation and the systematic error term. The ; default systematic error is 0.03, and can be ; set using the sys_err keyword in HXT_MCIMG. ; norm_error = if set, vary the SC counts using a normal distribution, ; with a standard deviation equal to norm_error*counts ; no_serr = if set, do not add systematic error to the variation ; So: counts = poidev(counts0)+sys_err*counts0*(2*(randomu(64,4)-0.5)) ; is the default variation. If norm_error is set, ; then counts = counts0(1.0+norm_error(0)*randomn(64,4)+ ; sys_err*counts0*(2*(randomu(64,4)-0.5) ; If /serr_only is set, then ; counts = sys_err*counts0*(2*(randomu(64, 4)-0.5) ; And if the /no_serr keyword is set, the sys_err term is not included. ; st_times= an array of accumulation start times, in any yohkoh format ; int_times= interval times corresponding to st_times, may be a scalar, ; if this is not set, and st_times is set, the routine will ; use the differences between st_times. This must be an integer ; number of the appropriate resolution in seconds, if the ; value of int_times is less than a frame value. ;OUTPUT: ; index1= index structure for the first frame of each image interval ; hxi_index.chan= channel used ; =0 L0 (15-24.4keV) ; =1 M1 (24.4-35.2 keV) ; =2 M2 (35.2 - 56.8 keV) ; =3 HI (56.8 - 100.0 keV) ; hxi_index.actim= image accumulation time in 10ths of secs ; hxi_index.interval= actim ; hxi_index.bkg_time= start time for background accumulation ; (millisec OF day) ; hxi_index.bkg_day= day of start of background accumulation ; (days since 1-Jan-79) ; hxi_index.bkg_actim= background accumulation time in 10ths of secs ; hxi_index.bkg_interval= actim ; counts= the counts accumulated over the chosen intervals ; sigma= sqrt(counts) ; bcounts= the counts accumulated over the background interval ; bsigma= sqrt(bcounts) ;DESCRIPTION: ; Uses the function plot_lcur to select the data to be imaged, ; and the background. ;MODIFICATION HISTORY: ; Written Oct-93 by JMM, ; Changed to give full resolution, 19-jul-94, JMM ; 9/21/94, account for 4.0 second time lag in hda files, jmm ; 1/9/95, changed calculation of sigma to be the same as ; in Tom Metcalf's hxtaccumulate.pro, to be consistent ; in pixon calculations. ; 1/13/95, added MC code here, as old hxt_mcimg is not valid ; 4/20/95, made consistent with new versions of hxtimg_accum, hxtsbt_mult, ; And added channel and trange keywords, to enable batch runs ; 6/9/95, subtract background from counts before choosing intervals ; 4/15/96, Added serr_only, norm_error, ; 4/23/96, Added no_serr ; 1/6/97, Added st_times, int_times, noback ;- PRO Hxt_mcaccum, index, data, index1, hxi_index, counts, sigma, bcounts, $ bsigma, trials=trials, img_trange=img_trange, $ bck_trange=bck_trange, channel=channel, $ accum_cnts=accum_cnts, ch_accum=ch_accum, $ Serr_only=serr_only, Norm_error=norm_error, no_serr=no_serr, $ St_times=st_times, Int_times=int_times, noback=noback serr = hxi_index(0).syserr nc0 = 4 ;no. of channels nc01 = nc0-1 ndset = N_ELEMENTS(index) ;no. of data sets nds1 = ndset-1 ;For full resolution, 4 minor frames for each data set nds4 = ndset*4 & nds41 = nds4-1 ;No. of trials IF(KEYWORD_SET(trials)) THEN mct = fix(trials(0)) ELSE mct = 100 ;i hate to have to do this, but here we decompress the whole data array data_x = float(hx_decomp(data)) ;sum SC's to plot and you'll need a time array for each data point zz = hxt_chsums(index, data_x, /half_sec, tim_arr = tim_arr1, $ tim_sec = tim_arr, dt = dx4, /overwrite) time0 = tim_arr1(0) ;reform data to get into the proper form data_x = reform(data_x, nc0, 64, nds4, /overwrite) ;If i need to choose intervals by plotting no_time_selected = NOT (KEYWORD_SET(img_trange) OR KEYWORD_SET(st_times)) no_bck_selected = NOT (KEYWORD_SET(bck_trange) OR KEYWORD_SET(noback)) IF(no_time_selected OR no_bck_selected) THEN BEGIN IF(N_ELEMENTS(channel) GT 0) THEN ich = fix(channel(0)) ELSE BEGIN print, ' INPUT CHANNEL NUMBER FOR THE INTERVAL>>' print, ' 0 FOR LO CHANNEL' print, ' 1 FOR M1 CHANNEL' print, ' 2 FOR M2 CHANNEL' print, ' 3 FOR HI CHANNEL' read, ich ich = fix(ich) ENDELSE ENDIF ELSE ich = 0 y0 = reform(zz(ich, *))/(dx4*64.0) ;y0 is the curve used to find the intervals IF(KEYWORD_SET(noback)) THEN BEGIN print, 'BACKGROUND SET TO ZERO' bcounts = fltarr(64, 4) bsigma = bcounts ib1 = 0 & ib2 = 3 ; i may need this stuff, for the structures ib14 = ib1/4 & ib24 = ib2/4 bactim = total(dx4(ib1:ib2)) ENDIF ELSE BEGIN ;choose the background interval print, 'BACKGROUND SELECTION::' choose_interval, tim_arr1, ib1, ib2, data = y0, plot_title = 'HXT COUNTS/SECOND/SC', $ time_range = bck_trange, dt_resolution = -1, alt_tsel = alt_tsel ib1 = ib1(0) & ib2 = ib2(0) ; Check time interval for non-flare mode data ib14 = ib1/4 & ib24 = ib2/4 not_flare = where(gt_dp_mode(index(ib14:ib24)) NE 9) IF(not_flare(0) NE -1) THEN BEGIN message, /info, 'WARNING, NON-FLARE MODE BACKGROUND CHOSEN, IS ZERO FOR HI, M1, AND M2 CHANNELS' ENDIF ; time interval, bactim = total(dx4(ib1:ib2)) tcounts = data_x(*, *, ib1:ib2) ;Speed things up, using total commands on 1 dimension nibs = ib2-ib1+1 IF(nibs GT 1) THEN BEGIN bcounts = transpose(total(tcounts, 3)) bsigma = transpose(sqrt(total(tcounts, 3)+ $ total((serr*tcounts)^2, 3))) ENDIF ELSE BEGIN ;no need to total bcounts = transpose(tcounts) bsigma = transpose(sqrt(tcounts+(serr*tcounts)^2)) ;oops, jmm,1-aug-95 ENDELSE ENDELSE ;Ok, now you choose the image interval, subtract the background first ; you need a sum for the background IF(KEYWORD_SET(noback)) THEN bcps = fltarr(4) ELSE BEGIN IF(ib2 NE ib1) THEN bcps = total(zz(*, ib1:ib2), 2) ELSE bcps = zz(*, ib1) ENDELSE bcps = bcps/(64.0*bactim) ;average background rate per SC print, 'IMAGE INTERVAL SELECTION::' IF(KEYWORD_SET(accum_cnts)) THEN BEGIN IF(accum_cnts EQ 1) THEN acc_cnts = 200.0 ELSE acc_cnts = accum_cnts(0) IF(N_ELEMENTS(ch_accum) GT 0) THEN ich0 = ch_accum(0) ELSE ich0 = 0 y0 = reform(zz(ich0, *)/64.0) y0 = y0-bcps(ich0)*dx4 choose_interval, tim_arr1, ss_st, ss_en, data = y0, accum_cnts = acc_cnts, $ time_range = img_trange, plot_title = 'BACKGROUND SUBTRACTED HXT COUNTS/SC' ENDIF ELSE BEGIN choose_interval, tim_arr1, ss_st, ss_en, data = y0-bcps(ich), $ plot_title = 'BACKGROUND SUBTRACTED HXT COUNTS/SEC/SC', $ time_range = img_trange, dt_res = -1, st_times = st_times, int_times = int_times ENDELSE IF(N_ELEMENTS(ss_st) GT 1) THEN $ print, 'ONLY ONE INTERVAL IS ALLOWED, THE FIRST INTERVAL CHOSEN IS USED' ss_st = ss_st(0) ss_en = ss_en(0) ; Check time interval for non-flare mode data ix14 = ss_st/4 & ix24 = ss_en/4 not_flare = where(gt_dp_mode(index(ix14:ix24)) NE 9) IF(not_flare(0) NE -1) THEN BEGIN message, /info, 'WARNING, NON-FLARE MODE DATA IN CHOSEN INTERVAL IS ZERO FOR HI, M1, AND M2 CHANNELS' ENDIF ;You'll not have an index structure for each minor frame, so for ;each minor frame, find the appropriate index, and reset .time ;and .day ss_st_index = ss_st/4 & index1 = index(ss_st_index) index1.gen.time = tim_arr1(ss_st).time index1.gen.day = tim_arr1(ss_en).day ;And replicate index1 index1 = replicate(index1(0), mct+1) ; now accumulate data, dti will be the time interval length, in seconds counts = fltarr(64, nc0, mct+1) sigma = counts dti = total(dx4(ss_st:ss_en)) dt_intv = tim_arr(ss_en)-tim_arr(ss_st)+dx4(ss_en) ;go ahead and set the seed variable for the random generation, ;for some reason, poidev does not work in this routine without ;a seed being set. tcounts = data_x(*, *, ss_st:ss_en) dnn = ss_en-ss_st+1 j0 = randomu(seed) FOR j = 0, mct DO BEGIN IF(dnn GT 1) THEN counts_x0 = transpose(total(tcounts, 3)) $ ELSE counts_x0 = transpose(tcounts) IF(j GT 0) THEN BEGIN IF(KEYWORD_SET(no_serr)) THEN serr_term = 0.0 $ ELSE serr_term = serr*counts_x0*(2.0*(randomu(seed, 64, nc0)-0.5)) IF(KEYWORD_SET(norm_error)) THEN BEGIN counts(*, *, j) = counts_x0*(1.0+norm_error*randomn(seed, 64, nc0))+serr_term ENDIF ELSE IF(KEYWORD_SET(serr_only)) THEN BEGIN counts(*, *, j) = counts_x0+serr_term ENDIF ELSE BEGIN counts(*, *, j) = poidev(counts_x0, seed = seed)+serr_term ENDELSE ENDIF ELSE counts(*, *, j) = counts_x0 ENDFOR counts = counts > 0.0 ;gotta have it sigma = sqrt(counts+(serr*counts)^2) ; Replicate and fill hxi_index hxi_index.bkg_time = tim_arr1(ib1).time hxi_index.bkg_day = gt_day(tim_arr1(ib1)) hxi_index.bkg_actim = fix(10.0*bactim) hxi_index.bkg_interval = fix(10.0*(tim_arr(ib2)-tim_arr(ib1)+dx4(ib2))) hxi_index = replicate(hxi_index(0), nc0, mct+1) FOR i = 0, nc01 DO hxi_index(i, *).chan = i FOR j = 0, mct DO BEGIN hxi_index(*, j).actim = fix(10.0*dti) hxi_index(*, j).interval = fix(10.0*dt_intv) ENDFOR data_x = 0.0 ; setgra gives the subscripts for the subcollimator switch setgra, sn2gra, gra2sn bcounts = bcounts(gra2sn, *) counts = counts(gra2sn, *, *) bsigma = bsigma(gra2sn, *) sigma = sigma(gra2sn, *, *) RETURN END