/***************************************************************** * Name of source file: j_ima_iros_work.c * Version of source file: 5.1.30 * Parent component: None * Programmer: Niels Lund & Niels J. Westergaard * Affiliation: Danish National Space Center * Contact address: Juliane Maries Vej 30 * DK-2100 Copenhagen, Denmark * Phone: +45 35325705 FAX: +45 35362475 * njw@dnsc.dk * Purpose: Initiate the backprojection and source * finding processes * Origin date: 041122 * Update history: 1.0.0 041130 First version * 1.0.4 050126 Test version for undue crash * 1.1.0 050131 Now reading FITS sky/backproj files * 1.1.1 050211 Keywords adapted to needs of j_ima_mosaic * 1.1.2 050303 Extra keyw:RLIMIDX, new cleaning * 1.1.3 050307 New DMAP: 51 radii * 1.1.4 050317 Align with PC-version as of 20050317 (NL) * 1.1.4 050328 'illum' cast to 'int' in 'if'-statement * 1.1.4 050407 New clean, etc. * 1.1.4 050407 New clean, etc. * 1.1.5 050415 New clean again, updated flux determination * 1.1.6 050418 Now uses JMXi-DBKG-MOD as default for empty field shds * 1.2.0 050425 Deals with SPR 4102,4108 and reduces RAM requirement * 1.3.0 050503 Deals with SPR 4120 and copes with up to 6 shadowgrams * 1.3.1 050517 Deals with SPR 4143 and initializes ima_sigma properly * 1.3.2 050519 Closes SPR 4157 - bad size of three arrays * 1.3.3 050602 Changes some arrays to static (SPR 4202) * 1.4.0 050620 Include user defined sources (SCREW 1740) * 1.4.1 051125 Transfer x_rms, y_rms and abs_len from jmx_id to backproj * 1.5.1 051129 Separating functions to jmx_lib_pif * 1.5.2 051202 Separating functions to jmx_lib_pif (2) * 1.5.3 051202 Give offX, x2x values in mm in IMOD-GRP * 1.5.4 060104 Bug in sh_depth corrected, detectmap now in mm * 1.5.5 060106 Introduced one more digit in logfile for src position * 1.5.6 060228 Removed parameter detMapFile and some obsolete variables * 1.5.7 060606 Shifted detCondMap to struct jmx_id and status_map * 1.5.8 060807 SPR 04551 Bad assignment of RLIMIDX corrected * 1.5.9 060807 SPR 04533 Resetting status variable has been removed * 1.5.10 060807 SPR 04523 Declaring testsum as float iso. int * 1.5.11 060807 SPR 04458 A warning log line has been added * 1.5.12 060907 SPR 04301 A relative VIGNDOL path is given * 1.6.0 060911 SCREW 1918 Image correction map and parameters introduced * 1.6.1 060912 SPR 4441 Shadowgrams with too few counts are excluded * 1.6.2 060915 SCREW 1842 A significance map can now be obtained * 1.6.3 060915 An exposure (vignetting) map is now a product * 1.6.4 060928 Included test for zero divisor in peaktransfer.c * 1.6.5 061002 New defaults for par. skyImagesOut * 1.7.0 061110 Produces images with few or zero events (SPR 4594) * 1.7.1 061114 SPR 4603 Image unit (BUNIT) is now cts/dm2/s * 1.8.0 061213 SPR 4640 Error in findpeak has been corrected * 1.8.1 070123 SPR 4590 int changed to long for 64 bit code * 1.8.2 070306 modified for new irosloop2 logic. * 2.1.0 070730 Active use of input source catalog * 2.1.0 070809 SCREWs 2021-2023 * only three fixed energy bands analyzed with IROS logic. * 2.1.1 070918 SPR 4729: Floating Point Exception * 2.1.2 070920 SPR 4731: Missing inialization fixed * 2.1.3 070920 SPR 4734: Another issing inialization fixed * 2.1.4 070921 SPR 4736: Fixing use of isnan function * 2.1.5 070925 SPR 4737: Bad array dimensioning fixed * 2.2.2 080111 Makes HK values available * 3.0.0 081111 SCREW 2080+1918 Electronic eff. corrections * 3.0.1 090114 SCREW 2098+2102+2105 and SPR 4812 * 3.0.2 090121 New table of sigma(PSF) is read * 3.0.X 090320 Modified 'peaktransfer' and error estimation * 3.0.3 090324 SPR 4837 Floating point exceptions removed * 3.1.0 090401 SPR 4737,4846 floating exceptions and negative fluxes * 3.2.0 090422 SCREW 2135 Better time correction implemented * 3.2.2 090519 SPR 4866+4869 Uninitialized variable problem * 3.2.3 090602 SPR 4875 Error in peak reintroduction * 3.2.4 090720 SPR 4886 Fix a divide by zero (BEO) * 3.2.6 091204 SCREW 2180/2181 Improve source finding/remove ghosts * 4.0.2 120321 SPR ?: Bug reported by JC, LP fixed * 4.0.6 120628 SPR ?: Reduced memory allocation demands * 4.0.7 120710 SPR ?: Flux error estimates as in version 3.2.7 * 4.0.9 120820 SPR ?: Improve FITS header content * 4.0.10 120830 SPR ?: Initialize variables as suggested by valgrind * 5.0.0 121001 Introduction of light curve generation * 5.0.1 150510 Introduction of time-averaging of burst-background estimate * 5.1.2 151106 Running version with LC generation * 5.1.9 160415 Update of light curve output and keywords according to agreement with ISDC * 5.1.10 160829 Circumvent analysis blocking caused by small amount of data in 'errorfit' * 5.1.11 160913 Update 'write_srcl_res' and 'work2d' to include 'pure burst' sources. * 5.1.12 160914 Update 'j_iir_Write_LC' to avoid floating point exc. for very short observation times * 5.1.14 160916 Update '*work2d.c' to avoid floating point exc. for very few time bins * 5.1.15 160925 Update 'j_iir_write_srclres' to include all sources checked by j_ima_iros * 5.1.16 160930 Update 'clean2', ex2trEvTslc', 'initPIF_LC' and 'work2d' to avoid crashes due to marginal user srcs * 5.1.19 161002 Update 'clean2' and 'srcidentify'o avoid potential index errors * 5.1.20 161006 Update 'work2d'. Protect against illegal values in sqrt call in line 3808 * 5.1.22 161215 Correct work2d.c. Wrong index in image transfer (line 3810). Update 'srcidentify'. * 5.1.23 161220 Correct peaktransfer & final_sky. Endless loop in peaktransfer for sources near edge of FOV * 5.1.24 161227 Correct work2d. Handling of 'edge'-bursts, Handling of bursts after 'summary_p', * 5.1.25 161228 Correct work2d. Initialization of Tslices[].threshold. 'Eliminate bad print in 'ex2trEvTslc' * 5.1.26 170102 Correct work2d. Handling of small number of time bins. * 5.1.28 170214 Replace 'srcidentify' (hard coded source list) with 'srcidentify2' reading the ISDC catalog * 5.1.29 170310 Multiply 'rateBuff, errBuff, bRateBuff and bErrBuff' by 100 before output in 'j_iir_Write_LC.c' * 5.1.30 170804 Modify max number of light curve energy bins (user def. + 3) from 7 to MAX_LC_EBIN * Type of element: Function * Part of package: j_ima_iros * Name: j_ima_iros_work * Required components: DAL, PIL, RIL, DAL3GEN *****************************************************************/ #include /* QQQQQQ */ int j_ima_iros_work( dal_element *theSWG, /* pointer to the SWG */ Auxdata auxdata, /* such as gain, GTI, and greyfilter */ struct instr_data *jmx_id, struct log_con *logger, swg_pointing *point, Events *bufevents, struct backpro *backproj, struct shadowgrams *sh, struct backpro_lists *bp, struct scwlist *SCWnow, goodsrc *GOODSRC, userbin *USERBIN, pifmap *PIFmap, timeslices *Tslices, b_list *B_list, Shd_adhoc shd_adhoc, Shd_adhoc shd_burst, /* lc_auxData *LC_AuxData, */ int chatter, int status ) { Instrument jemxNum; char ID_string[13], burst_txt[32], MPdec, Stext[12], chsign; char f_exttext[2][J_IIR_TEXTSIZE], regsrcnam[64]; char elmntName[DS_NAME_SIZE], sourcename[32], ISDCname[17]; static char SHAPE[64][10], COLOR[64][8], SourceName[32], regtext[128]; int eeb, eeuser, eelim, kMAX, i5 = 5, allsrc, reasrc, allWEpos, allWEzero; int jb, ix, iy, oldaccessum, newaccessum, jemx, addnew, old_b_listnum; int ch1, ch2, ch3, ch4, ch5, ch6, ch7, I_RA, I_dec, ibu, N_SRC2orig=0; int imin=0, jlim, jj, xix, yiy, newflag, tstsrc; int ibu2, Uibu[20], Uibu2[20], NdontUse=0, MMM=1; int preset_temp = 0, n_additional, xsrc, presrc, firstAddSource, oldBURST_NUM=0; int sky_index, sky_index0, ghoststep[2], srclist_index, prevent=0, jjee, jee; int iki, kk_offs=0, k_offset=0, kk_off2=0, srcID, doneID[20]; int preset_src_orig, catflag, ChanMax=0, chanlim=5; long Llength = 0L; static float xx[SCW_SRCMAX], yy[SCW_SRCMAX]; static float goodbet[SCW_SRCMAX]; static float backcnt[EE256+1][SCW_SRCMAX]; static float shadow[detdim]; static float peakval[2*EE256][SCW_SRCMAX]; static float shadow_r[detdim], shdg_fits[detdim]; static float burst_skyim[20][MAX_LC_EBIN][SKYDIM], skyim_mask[SKYDIM], localima[SKYDIM]; static float s_reject[EE256]; static float s_accept[EE256]; static float take_away_var[SCW_SRCMAX]; static int ubin_shd[detdim], ubin_totcnt; static int ubin_shd2[detdim], ubin_totcnt2; static int src_id[SCW_SRCMAX]; static int I_sxy[SCW_SRCMAX]; static int miniCorner[6][2]; static int burst_ima_qual[20][MAX_LC_EBIN]; static int burst_lc_output[100][MAX_LC_EBIN]; static int src_lc_output[16][MAX_LC_EBIN]; static int preset_ibu[MAX_NUM_PRESET_SRCS], ibu_preset[20], ibu_N_SRC2[20]; static short used_ee[EE256][N_SRCMAX], used2[N_SRCMAX], preset[SCW_SRCMAX]; static short AA[256]; static double solut_e[EE256][BASIC_SRCMAX]; static double solut[SCW_SRCMAX]; static double orig_solut[4][4], orig_rms[4]; static float maps[300 * sizeof(MAPS_SRC) / sizeof (float)]; MAPS_SRC *maps_src = NULL; double sh_rms[EE256+1][4]; double amx[6], sky_rms, val; double evRA, evdec, catRA, catdec; int map_n, xvalstart, yvalstart, n, s_field; int ylow, yhigh, ShNNN, skymax, miniIndex=0, ee_mask, NumUnknown; int ii, sxy_skymax; int chanmax[128], burst_num, ipr=0, bstsrc; int maxchan=-1, stepplot=0; float wbkg_pixel_area, wbkg_pixel_wgt, fitval, plotv, plsigma, contrib, sumE; float m0=0.0, m1=0.0, m2=0.0, fluxmax=0.0, meanFlux=0.0, rmsFlux=0.0, rmsFlux0=0.0; float correl, threshold[10]; float fitrms, factor, pkval, ShSUM, old_avGain, Isignif, sflux; float diff_x, diff_y, x, y, r, events, events2; float cmx, cmy, avr, aveffi, avillu; float funval, transmis, Emin, Emax, burst_sum; float sumX, sumY, sumgood, T_burst, T_total; float s_dec, s_RA, maxval; float avX0, avY0, xmm, ymm, Gmean[MAX_LC_EBIN]; float peakcount[2], xxx, yyy; float fk, sum_str, tot_area, sn_best, sn_temp; float temptot, tempgood; float sn_new, normfact, fact1, fact2, xx_dif, yy_dif; float s_sum, c_sum; float hot_reject, best_signif; float testsum = 0.0; float eff_imag_min=1000.0, eff_imag_max=-1000.0; float mean_bkg[2] = {0.0, 0.0}; float e_minim, e_maxim, meanE, shsum, wsum, sum3, sum2; float abs_len, escape, remaining, signif; float rawcount, Apix, Anoise = 0.5; float Skymax, Skymin, m01; /* float pifnum, pifmax, pifmean, pifmin, pifval, pifzero; */ float edgeXmm[9], edgeYmm[9]; float X_rms[256], Y_rms[256], flux, RMSfact, evX, evY; float RMS_correction, rms_factor, x_rms, y_rms, rms, mean, var, mean2, mean3; float mean1, rms1, rr, NBsigniflim = 4.0, av_expo; float allEvents, allEEvents, allSrcEvents, allBkgEvents, allWSrcEvents, allWBkgEvents; float maxflux=0.0, xscale=1.0; int i, j, k, ej, kkkk, num_flag1_srcs, num_flag3_srcs, ichan, kscale=1, kscale2=1; int chk_no, skysize, CHANmax, IMA_qual=0; int pixels, j_src, psrc, pd, current_start[20], current_end[20], current_ibu; int bitpix, dim1, dim2, ee_lc; int isrc, jsrc, n_src1, n_src2; int n_common, n_common2; int iee, iiee, vigncorr = 1; /* default: vignetting correction will be applied ! */ int n_sources, n_sources_ini, fitcycles; int e1, e2, ksrc, sxy; int n_src, flag1src = 0, Bburst_num=-1; int ij, ij_out, burst_number[10]; int baseErr, listlim; int e1_stop, e2_stop, ic_stop, jc_stop; int PIlop, BlistNum, ALLsrc, ntb, ilc; int ntime, nevent, ixy, edgsrc, lstsrc, iisum, iicat, reg_mask = 0; int edgeI, sourceID, last = 0, NOTIMEBINS = 0, FEWTIMEBINS = 0; double quality, diff, sum; double sh_rms2; double xrms, yrms, burstTime; double sumtot, cdelt, rad_ang; double divisor = 0.0, exposure = 0.0; double nn_f_min, nn_f_max, varim_min, varim_max; static struct fit_struc FS1, FS2; struct fit_struc *f_struc1 = {&FS1}; struct fit_struc *f_struc2 = {&FS2}; static struct candsrc CANDSRC[N_SRCMAX]; /* structure containing information about source candidates found */ int pif_status, old_psrc = 0; FILE *par_chk, *gnufil1, *gnufil2, *regfil1, *regfil2, *isdcat; char calibname[100], gnuname1[64], gnutext[256], regname1[64], textline[200], regnam2[64]; char gnuname2[64]; char *ptrA, *ptrB; maps_src = (MAPS_SRC*)maps; long nEvents; static int func_local = -1, local_calling_num=0; clock_t calling_TT, TT_now; /* function entry timing start */ if (func_local < 0) { func_local = logger->func_used_num++; sprintf(logger->func_names[func_local], "j_ima_iros_work"); if (logger->func_used_num > 199) { sprintf(logger->tmptxt, "func_used_num > 199"); logger->logstat = logprint(logger, logger->tmptxt, 6); logger->func_used_num = 199; } } logger->func_calls[func_local]++; logger->TT = TT_now = clock(); calling_TT = logger->TT; local_calling_num = logger->func_calling_num; logger->func_calling_num = func_local; /* function entry timing end */ /* 111111 */ /* QQQQQQ */ if (logger->trace) traces(func_local, 0, logger); sprintf(ID_string, "%04d%04d%04d", backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv); strcpy(logger->scw_file_name, ID_string); i5 = logger->i5; /* NL 20120122 */ threshold[0] = 3.5; /* current low threshold */ threshold[1] = logger->thrs; /* current high threshold */ threshold[2] = 3.5; /* low threshold for found sources */ threshold[3] = logger->thrs; /* high threshold for found sources */ threshold[4] = 3.5; /* low threshold for user sources */ threshold[5] = logger->thrs; /* high threshold for user sources */ threshold[6] = 3.5; /* low threshold for global sources */ threshold[7] = logger->thrg; /* high threshold for global sources */ threshold[8] = 3.5; /* low threshold for edge sources */ threshold[9] = logger->thrg; /* high threshold for edge sources */ logger->srcIDlimit = 0.05; logger->b_listnum = 0; old_b_listnum = 0; for (ibu=0; ibu<20; ibu++) { for (iee=0; iee<8; iee++) burst_ima_qual[ibu][iee] = 0; Uibu[ibu] = -1; Uibu2[ibu] = -1; doneID[ibu] = 0; ibu_preset[ibu] = 0; ibu_N_SRC2[ibu] = 0; } for (ibu=0; ibu<100; ibu++) for (iee=0; iee<8; iee++) burst_lc_output[ibu][iee] = 0; for (ibu=0; ibu<16; ibu++) for (iee=0; iee<8; iee++) src_lc_output[ibu][iee] = 0; for (i=0; i<64; i++) B_list[i].sourceID = -1; i = MAX_LC_EBIN; sprintf(logger->tmptxt, "In main: MAXLCEBIN %d, LCMASK: %d\n", i, logger->LCMAXMASK); logger->logstat = logprint(logger, logger->tmptxt, 0); for (i=0; i<10; i++) { sprintf(logger->tmptxt, "threshold[%1d]: %4.1f", i, threshold[i]); logger->logstat = logprint(logger, logger->tmptxt, 0); burst_number[i] = 0; } status = ISDC_OK; logger->new_detmap = 1; for (i=0; iNumTimeBins; i++) Tslices[i].burstflag = 0; for (i=0; isrc_illum_cnts[i][j] = 0.0; for (i=0; iorig_data_peak[j][i] = 0.0; for (i=0; i<3*SCW_SRCMAX; i++) backproj->blindspots[i] = 0; for (i=0; i<64; i++) {sprintf(COLOR[i], "white"); sprintf(SHAPE[i], "cross");} for (i=0; i<64; i++) logger->allsrc_ptr[i] = 0; sprintf(COLOR[2], "black"); sprintf(COLOR[3], "cyan"); sprintf(COLOR[4], "green"); sprintf(COLOR[5], "magenta"); sprintf(COLOR[6], "yellow"); sprintf(COLOR[7], "red"); sprintf(SHAPE[0], "circle"); sprintf(SHAPE[1], "diamond"); sprintf(SHAPE[2], "x"); sprintf(SHAPE[3], "arrow"); sprintf(SHAPE[4], "cross"); miniCorner[3][0] = 2; miniCorner[3][1] = 469; miniCorner[2][0] = 2; miniCorner[2][1] = 424; miniCorner[1][0] = 2; miniCorner[1][1] = 379; miniCorner[0][0] = 2; miniCorner[0][1] = 334; /* *************************************************************************** */ /* */ /* BEGIN BLOCK PROCESSING ONE SCIENCE WINDOW */ /* */ /* *************************************************************************** */ backproj->map_cnt = 0; for (i=0; is_detect[i] = 0; for (i=0; i<256; i++) AA[i] = 1; for (i=0; icposres[PIlop] / 2.36; Y_rms[PIlop] = jmx_id->bposres[PIlop] / 2.36; RMS_correction = 0.5; /* QQQQQ originally: 0.5; Niels Lund 20080421 */ if (jmx_id->jmx_unit == 0) { if (jmx_id->escale[PIlop] < 11.0) rms_factor = 1.0 + RMS_correction; if ((jmx_id->escale[PIlop] >= 11.0) && (jmx_id->escale[PIlop] < 35.0)) rms_factor = 1.0 + RMS_correction * (35.0 - jmx_id->escale[PIlop]) / 24.0; if (jmx_id->escale[PIlop] >= 35.0) rms_factor = 1.0; X_rms[PIlop] /= rms_factor; Y_rms[PIlop] /= rms_factor; } } for (i=0; inumBkgShds; i++) { meanE = 0.0; k = 0; for (j=0; j<256; j++) { if ((jmx_id->escale[j] >= sh->bkg_eval[i][0]) && (jmx_id->escale[j] <= sh->bkg_eval[i][1])) { k++; meanE += jmx_id->escale[j]; } } sh->E_mean_bkg[i] = meanE / (float)(k); } /* if (jmx_id->jmx_unit == 0) fine_anoder(AA); */ /* get NJW selected anodes into AA */ ch1 = ch2 = ch3 = ch4 = ch5 = ch6 = ch7 = 0; yhigh = 0; ylow = 256; for (i=0; istatus_map[i] > 0) { ch1++; if ((i%256) > yhigh) yhigh = i % 256; if ((i%256) < ylow ) ylow = i % 256; if ((jmx_id->status_map[i] <= sh->detAccLimit) && (jmx_id->status_map[i] != 32)) ch2++; } if (sh->soft_shdgr[0][i] > 0.0) ch3++; if (sh->soft_shdgr[1][i] > 0.0) ch4++; if (sh->soft_shdgr[2][i] > 0.0) ch5++; if (sh->soft_shdgr[3][i] > 0.0) ch6++; /* prepare common select map to be used for all energy bands */ ix = i / 256; sh->common_select_map[i] = 0; if ((jmx_id->status_map[i] > 0) && (AA[ix] == 1) && (jmx_id->status_map[i] <= sh->detAccLimit) && (jmx_id->status_map[i] != 32)) { if ((sh->soft_shdgr[0][i] > 0.0) && (sh->soft_shdgr[1][i] > 0.0) && (sh->soft_shdgr[2][i] > 0.0) && (sh->soft_shdgr[3][i] > 0.0)) sh->common_select_map[i] = 1; } if (sh->common_select_map[i] > 0) ch7++; } for (j=0; jempty[j][i] = 0.0; /* ************************************************************************* */ /* * * */ /* * QQQQQQ * */ /* * Make loop over search energy bins and generate corresponding * */ /* * background shadowgrams * */ /* * * */ /* * * */ /* ************************************************************************* */ backproj->additional = -1; for (eeuser=0; eeuseree_basic; eeuser++) { status = j_fitprep(eeuser, theSWG, jmx_id, logger, backproj, sh, bp, SCWnow, USERBIN, chatter, status); if (logger->trace) traces(func_local, 55, logger); if (status < 0) goto exittrace; if (logger->trace) traces(func_local, 155+eeuser*1000, logger); sprintf(logger->tmptxt, "work2d: after fitprep: E: %2d, totcnt: %7.1f", eeuser, sh->totcnt[eeuser]); logger->logstat = logprint(logger, logger->tmptxt, i5); } if (logger->trace) traces(func_local, 155, logger); /* prepare header data for fits files */ sprintf(ID_string, "%04d%04d%04d", backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv); strcpy(logger->scw_file_name, ID_string); if (logger->trace) traces(func_local, 165, logger); if ((logger->develop & 1) == 1) { sprintf(regname1, "J%1d_%s.reg", jmx_id->jmx_unit+1, ID_string); regfil1 = fopen(regname1, "wt"); fprintf(regfil1, "# Region file format: DS9 version 4.0\n"); fprintf(regfil1, "global color=black font=\"helvetica 10 normal\" select=1 highlite=1 edit=1 move=1 delete=1 include=1 fixed=0 source\n"); fprintf(regfil1, "# RA/dec/R: %7.4f %7.4f %4.0f, jmxRA/dec: %7.4f %7.4f, ObsT: %6.1lf s\n", point->j_RA / GTORAD, point->j_dec / GTORAD, point->j_posangle / GTORAD, point->jmxRA, point->jmxdec, sh->accumT); fprintf(regfil1, "IMAGE;text(382.0, 510.0) # text={JEM-X%1d_%s} font=\"helvetica 14 normal\"\n", jmx_id->jmx_unit+1, ID_string); fprintf(regfil1, "IMAGE;text(450.0, 490.0) # text={ObsT: %6.1lf s} font=\"helvetica 8 normal\" \n", sh->accumT); sprintf(logger->tmptxt, "Burst thresholds (sigma):"); for (i=4; i<10; i+=2) { sprintf(logger->tmptxt2, " %4.1f %4.1f,", threshold[i], threshold[i+1]); if ((strlen(logger->tmptxt) + strlen(logger->tmptxt2)) < MAX_STR_LEN) strcat(logger->tmptxt, logger->tmptxt2); } fprintf(regfil1, "IMAGE;text(360.0, -10.0) # text={%s} font=\"helvetica 8 normal\" \n", logger->tmptxt); } sprintf(logger->tmptxt2, " Energy bands: "); if ((strlen(logger->tmptxt) + strlen(logger->tmptxt2)) < MAX_STR_LEN) strcat(logger->tmptxt, logger->tmptxt2); eelim = sh->numShds - backproj->ee_basic; for (j=backproj->ee_basic; jnumShds; j++) { sprintf(logger->tmptxt2, " %4.1f %4.1f,", jmx_id->escale[sh->pival[j][0]], jmx_id->escale[sh->pival[j][1]]); if ((strlen(logger->tmptxt) + strlen(logger->tmptxt2)) < MAX_STR_LEN) strcat(logger->tmptxt, logger->tmptxt2); } if ((logger->develop & 1) == 1) fprintf(regfil1, "# %s\n", logger->tmptxt); if (logger->trace) traces(func_local, 255, logger); if (jmx_id->jmx_unit == 0) { sprintf(f_struc1->instru, "JMX1"); jemxNum = JMX1; strcpy(elmntName, "JMX1-EVTS-SHD"); jemx = 0; } else { sprintf(f_struc1->instru, "JMX2"); jemxNum = JMX2; strcpy(elmntName, "JMX2-EVTS-SHD"); jemx = 1; } f_struc1->naxis1 = J_IIR_DIM_STDSHDG; f_struc1->naxis2 = J_IIR_DIM_STDSHDG; f_struc1->bitpix = -32; f_struc1->jemx = jmx_id->jmx_unit; f_struc1->orbit = backproj->aux_orbit; f_struc1->pid = backproj->aux_pid; f_struc1->pidv = sh->aux_pidv; f_struc1->ctrl = 0; dim1 = dim2 = backproj->FOV_radius * 2 + 1; /* midisky image dimensions */ bitpix = -32; /* use float representation in FITS file */ if (jmx_id->jmx_unit == 0) sprintf(f_struc2->instru, "JMX1"); else sprintf(f_struc2->instru, "JMX2"); f_struc2->naxis1 = dim1; f_struc2->naxis2 = dim2; f_struc2->bitpix = -32; f_struc2->jemx = jmx_id->jmx_unit; f_struc2->orbit = backproj->aux_orbit; f_struc2->pid = backproj->aux_pid; f_struc2->pidv = sh->aux_pidv; f_struc2->crval1 = SCWnow->jmxRA; f_struc2->crval2 = SCWnow->jmxdec; f_struc2->rot_ang = SCWnow->roll_angle + jmx_id->phase2; f_struc2->crpix1 = (double)(backproj->FOV_radius+1); f_struc2->crpix2 = (double)(backproj->FOV_radius+1); f_struc2->exposure = sh->accumT; f_struc2->varia_id = 100; f_struc2->ctrl = 1; strcpy(f_struc2->bkg_model, logger->ref_file_name); for (i=0; itrace) traces(func_local, 355, logger); /* QQQQQQ check for user defined sources */ sprintf(logger->tmptxt, "##8777## In work: preset_src: %d\n", backproj->preset_src); logger->logstat = logprint(logger, logger->tmptxt, i5); j = backproj->preset_src + SCW_SRCMAX + 18; for (i=0; ipointing_RA = point->j_RA / GTORAD; SCWnow->pointing_dec = point->j_dec / GTORAD; SCWnow->roll_angle = point->j_posangle / GTORAD; SCWnow->jmxRA = point->jmxRA; SCWnow->jmxdec = point->jmxdec; /* dummy call to find catalog sources in FOV */ srcidentify2(theSWG, SCWnow->jmxRA, SCWnow->jmxdec, sourcename, ISDCname, &catflag, jmx_id, point, logger); ISDCname[16] = 0; num_flag1_srcs = 0; num_flag3_srcs = 0; n = backproj->preset_src; if (logger->trace) traces(func_local, 455, logger); /* Eliminate offending user specified sources */ if (backproj->preset_src > 0) { for (i=0; ipreset_src; i++) { if ((fabs(backproj->preset_src_RA[i] - logger->XRA) < 0.03) && (fabs(backproj->preset_src_dec[i] - logger->Xdec) < 0.03)) { if (backproj->preset_src_flag[i] == 0) backproj->preset_src_flag[i] = 3; } next_try: /* next_try */ status = getxy(backproj->preset_src_RA[i], backproj->preset_src_dec[i], jmx_id->jmx_unit, &xmm, &ymm, jmx_id, point, logger); if (status < 0) { /* remove offending source (outside FOV) */ for (j=i; j<(backproj->preset_src-1); j++) { backproj->preset_src_RA[j] = backproj->preset_src_RA[j+1]; backproj->preset_src_dec[j] = backproj->preset_src_dec[j+1]; backproj->preset_src_flag[j] = backproj->preset_src_flag[j+1]; } backproj->preset_src--; sprintf(logger->tmptxt, "Source removed: RA/dec: %f %f, flag: %d\n", backproj->preset_src_RA[i], backproj->preset_src_dec[i], backproj->preset_src_flag[i]); logger->logstat = logprint(logger, logger->tmptxt, i5); status = 0; if (backproj->preset_src > i) goto next_try; } else { backproj->preset_src_xx[i] = xmm / backproj->pix2mm + backproj->FOV_radius; backproj->preset_src_yy[i] = ymm / backproj->pix2mm + backproj->FOV_radius; backproj->preset_src_rr[i] = sqrt(xmm*xmm + ymm*ymm) / backproj->pix2mm; if ((backproj->preset_src_flag[i] & 3) == 3) num_flag3_srcs++; if ((backproj->preset_src_flag[i] & 3) == 1) num_flag1_srcs++; logger->srcIDlimit = 0.035; srcidentify2(theSWG, backproj->preset_src_RA[i], backproj->preset_src_dec[i], sourcename, ISDCname, &catflag, jmx_id, point, logger); logger->srcIDlimit = 0.05; ISDCname[16] = 0; sprintf(logger->tmptxt, "Catalog Source in FOV, X/Y/R: %6.1f %6.1f %6.1f, RAdec: %7.3f %7.3f flag: %d, ZZZ%s", backproj->preset_src_xx[i]-backproj->FOV_radius, backproj->preset_src_yy[i]-backproj->FOV_radius, backproj->preset_src_rr[i], backproj->preset_src_RA[i], backproj->preset_src_dec[i], backproj->preset_src_flag[i], sourcename); logger->logstat = logprint(logger, logger->tmptxt, 0); } } } if (logger->trace) traces(func_local, 555, logger); sprintf(logger->tmptxt, "Original number of catalog sources: %3d, inside FOV: %d, flag-1/3 sources: %2d %2d", n, backproj->preset_src, num_flag1_srcs, num_flag3_srcs); logger->logstat = logprint(logger, logger->tmptxt, 0); if (num_flag3_srcs > SCW_SRCMAX-backproj->backmodels) num_flag3_srcs = SCW_SRCMAX-backproj->backmodels; sprintf(logger->tmptxt, "Total number of accepted catalog sources: %3d, flag-1/3 preset src: %2d %2d", backproj->preset_src, num_flag1_srcs, num_flag3_srcs); logger->logstat = logprint(logger, logger->tmptxt, 0); /* ************************************************************************* */ /* * * */ /* * QQQQQQ * */ /* * check for preset source type-3 and introduce these in source list * */ /* * * */ /* * * */ /* ************************************************************************* */ n_sources = backproj->backmodels; for (i=0; isky_ydim - 294; ghoststep[1] = -318 * backproj->sky_ydim + 37; } else { ghoststep[0] = -126 * backproj->sky_ydim + 294; ghoststep[1] = 318 * backproj->sky_ydim - 37; } for (i=0; ipreset_src; i++) { if ((backproj->preset_src_flag[i] & 3) != 3) continue; /* introduce preset sources type-3 into list */ xx[n_sources] = backproj->preset_src_xx[i]; yy[n_sources] = backproj->preset_src_yy[i]; preset[n_sources] = 3; sprintf(logger->tmptxt, "##261x## work2: preset src3: %2d, xx/yy: %6.2f %6.2f", backproj->blindnum_preset, xx[n_sources], yy[n_sources]); logger->logstat = logprint(logger, logger->tmptxt, i5); /* make irosloop 'blind' to the preset sources */ sky_index = (int)(yy[n_sources] + 0.5) + (int)(xx[n_sources] + 0.5) * backproj->sky_ydim;; backproj->blindspots[backproj->blindnum_preset] = sky_index; n_sources++; if (n_sources >= SCW_SRCMAX) n_sources = SCW_SRCMAX - 1; else { backproj->blindnum_preset++; /* check for potential ghosts, eliminate these as well. Note: one more instance! NL 2009-12-03 */ if (((sky_index + ghoststep[0]) > 0) && ((sky_index + ghoststep[0]) < backproj->skydim)) { backproj->blindspots[backproj->blindnum_preset] = sky_index + ghoststep[0]; backproj->blindnum_preset++; } if (((sky_index + ghoststep[1]) > 0) && ((sky_index + ghoststep[1]) < backproj->skydim)) { backproj->blindspots[backproj->blindnum_preset] = sky_index + ghoststep[1]; backproj->blindnum_preset++; } } } if (logger->trace) traces(func_local, 655, logger); /* find common names for preset sources */ j = k = 0; /* isdcat = fopen("/r6/jemx/catalogs/XrBursters_cat.reg", "rt"); */ for (i=0; ipreset_src; i++) { sprintf(backproj->preset_src_name[i], "PresetSrc%04d_%02d", iicat, i); /* if ((isdcat != NULL) && (logger->develop)) { */ evRA = backproj->preset_src_RA[i]; evdec = backproj->preset_src_dec[i]; iicat = 9999; sprintf(sourcename, "PresetSrc%04d_%02d", iicat, i); /* status = catnib(isdcat, evRA, evdec, sourcename, &iicat, logger, status); if (iicat > 0) sprintf(sourcename, "PresetSrc%04d_%02d", iicat, i); */ logger->srcIDlimit = 0.035; srcidentify2(theSWG, evRA, evdec, sourcename, ISDCname, &catflag, jmx_id, point, logger); logger->srcIDlimit = 0.05; ISDCname[16] = 0; sprintf(backproj->preset_src_name[i], "%s", sourcename); /* } */ sprintf(logger->tmptxt, "PresetSRC%02d: RA/dec: %7.3lf %7.3lf, flag: %4d, name: %s", i, evRA, evdec, backproj->preset_src_flag[i], backproj->preset_src_name[i]); logger->logstat = logprint(logger, logger->tmptxt, 0); } /* if (isdcat != NULL) fclose(isdcat); */ isdcat = NULL; if (logger->trace) traces(func_local, 755, logger); /* ************************************************************************* */ /* ************************************************************************* */ /* * * */ /* * * */ /* * Here starts the source finding loop over the basic energy bands * */ /* * (ends approximately in line 1050) * */ /* * * */ /* * * */ /* ************************************************************************* */ /* ************************************************************************* */ /* 222222 */ /* QQQQQQ */ n_sources_ini = n_sources; for (eeb=0; eebee_basic; eeb++) { /* eeb = 0 to 2 indicates the energy band */ n_sources = n_sources_ini; if (logger->trace) traces(func_local, 90000+eeb, logger); SCWnow->cand_lst_start[eeb] = logger->N_SRC; /* save start index value of source candidate list */ for (i=0; iskydim; i++) { /* backproj->fitsimage[eeb][0][i] = USERBIN[eeb].variance[i]; */ backproj->illumx[i] = USERBIN[eeb].illumnew[i]; backproj->access[i] = 0; if (backproj->illumx[i] >= (int)(backproj->signiflim)) backproj->access[i] = 1; } for (i=0; ieffic[i] = USERBIN[eeb].effic_shdg[i]; if ((sh->common_select_map[i] > 0) && (USERBIN[eeb].effic_shdg[i] > 0.0)) backproj->soft_gr[i] = USERBIN[eeb].soft_gr[i]; else backproj->soft_gr[i] = -1.0; /* this pixel should not be used */ } /* * if( chatter > J_CHATTY_MAXIMAL ) { * RILlogMessage( NULL, Log_0, "##34## Dump of backproj->soft_gr"); * status = fits_test_output( "soft_gr_", backproj->soft_gr, 256, logger, chatter, status ); * if( status != ISDC_OK ) goto exittrace; * } */ /* transfer the raw shadowgram to be analyzed into "shadow_r[]" */ s_sum = s_accept[eeb] = s_reject[eeb] = 0.0; pixels = 0; backproj->act_pix[eeb] = 0.0; for (i=0; isoft_gr[i] > 0.0) { shadow_r[i] = USERBIN[eeb].f_shadgram[i]; shadow[i] = USERBIN[eeb].shadow[i]; s_accept[eeb] += USERBIN[eeb].f_shadgram[i]; s_sum += shadow_r[i]; backproj->act_pix[eeb] += 1.0; pixels++; } else { shadow_r[i] = 0.0; s_reject[eeb] += USERBIN[eeb].f_shadgram[i]; } } hot_reject = 0.0; /* hot_reject = hotpix_chk(shadow_r, eeb, backproj, logger); */ s_accept[eeb] -= hot_reject; s_sum -= hot_reject; s_reject[eeb] += hot_reject; sprintf(logger->tmptxt, "##35## Accepted/rejected shadowgram counts: E%2d: %8.1f %8.1f, hot_rejects: %6.1f", eeb, s_accept[eeb], s_reject[eeb], hot_reject); logger->logstat = logprint(logger, logger->tmptxt, 2); sh_rms[eeb][0] = rmscalc(shadow_r, backproj->soft_gr, logger); if( pixels == 0.0 ) { RILlogMessage( NULL, Warning_1,"##36## Attempt to divide with zero! reset to 1"); pixels = 1.0; } SCWnow->mean_pix[eeb] = s_accept[eeb] / (float)(pixels); if (logger->shdg_out_ctrl == 1) { /* output status map as FITS */ for(i=0; i<256-logger->fits_shft; i++){ for (j=0; j<256; j++) { ij = i + j*256; shdg_fits[ij] = (float)(jmx_id->status_map[ij]); } } strcpy(f_exttext[0], "status map"); strcpy(f_exttext[1], "condition map shadowgram"); status = shdg_to_fits(theSWG, eeb , shdg_fits, jemxNum, elmntName, f_exttext, f_struc1, logger, eeb, sh, chatter, status ); if (status != ISDC_OK) goto exittrace; /* output raw input shadowgram as FITS */ strcpy(f_exttext[0], "RAWSHDIN_T"); strcpy(f_exttext[1], "Raw input shadowgram - transposed"); status = shdg_to_fits(theSWG, eeb , shadow_r, jemxNum, elmntName, f_exttext, f_struc1, logger, eeb, sh, chatter, status ); if (status != ISDC_OK) goto exittrace; /* output pixelfolded shadowgram as FITS */ strcpy(f_exttext[0], "PIXFOLD_SHDGR"); strcpy(f_exttext[1], "Pixelfold shadowgram"); status = shdg_to_fits(theSWG, eeb , shadow, jemxNum, elmntName, f_exttext, f_struc1, logger, eeb, sh, chatter, status ); if (status != ISDC_OK) goto exittrace; } if (USERBIN[eeb].totcnt <= 1.0) { sh_rms[eeb][1] = 0.0; SCWnow->cand_lst_no[eeb] = 0; continue; } sh_rms[eeb][1] = rmscalc(shadow, backproj->soft_gr, logger); fk = 0.0; j = 0; for (i=0; isoft_gr[i] > 0.0) fk += shadow[i]; if (shadow[i] > 0.0) j++; backproj->shadw[i] = shadow[i]; } sprintf(logger->tmptxt, "##40## JEM-X%1d. Cnts/rms before imageprep: %8.1f/%8.2lf, after: %8.1f/%8.2lf, check: %6.0f, E%02d, nonzero: %5d", jmx_id->jmx_unit+1, s_sum, sh_rms[eeb][0], USERBIN[eeb].totcnt, sh_rms[eeb][1], fk, eeb, j); logger->logstat = logprint(logger, logger->tmptxt, 2); jmx_id->absolen = USERBIN[eeb].abs_len; /* calculate the position scatter around the X-ray absorbtion point at this energy */ jmx_id->scatt_len = sh_xy(USERBIN[eeb].x_rms, USERBIN[eeb].y_rms, logger, jmx_id); if ((logger->develop & 1) == 1) { if (eeb == 0) par_chk = fopen("work_par_chk.txt", "wt"); else par_chk = fopen("work_par_chk.txt", "at"); saveparams(par_chk, eeb, jmx_id, backproj, bp, logger); fclose(par_chk); par_chk = NULL; } /* **************************************************************************************** */ /* */ /* The "irosloop2" goes through the full "Iterative Removal Of Sources" sequence */ /* for one energy value (eeb) until the loop terminates due to limited significance or */ /* too many iteration loops */ /* */ /* **************************************************************************************** */ backproj->additional = -1; sprintf(logger->tmptxt, "##41## calling irosloop"); logger->logstat = logprint(logger, logger->tmptxt, 2); backproj->stopcode[eeb] = 0; preset_temp = backproj->preset_src; if (logger->trace) traces(func_local, 91000+eeb, logger); if (logger->trace) traces(func_local, 91000+logger->dontUseShdg[eeb], logger); if (logger->dontUseShdg[eeb] > 0) { n_sources = backproj->backmodels; } else { status = irosloop2(jmx_id->jmx_unit, eeb, xx, yy, src_id, I_sxy, solut, &fitrms, &sh_rms2, &fitcycles, &n_sources, goodbet, preset, jmx_id, logger, backproj, bp, sh, chatter, status); } if (logger->trace) traces(func_local, 92000+eeb, logger); backproj->preset_src = preset_temp; if ( status != ISDC_OK ) { logger->logstat = logprint(logger, logger->error_text, 7); status = J_ERROR_CODE_START+J_IMA_IROS_ERR+J_IIR_IROSLOOPERR; goto exittrace; } if (logger->trace) traces(func_local, 93000+eeb, logger); sprintf(logger->tmptxt, "##42## returning from irosloop, ee: %1d, n_sources: %4d, stopcode: %1x, map_cnt: %3d", eeb, n_sources, backproj->stopcode[eeb], backproj->map_cnt); logger->logstat = logprint(logger, logger->tmptxt, 2); /* status = fits_test_output( "work2_usbfi_", USERBIN[eeb].fst_imag, 511, logger, chatter, status ); * the content USERBIN[eeb].fst_image and of backproj->fitsimage[eeb][2] is unchanged here */ if( chatter > J_CHATTY_VERBOSE ) { testsum = 0.0; for( i = 0; i < backproj->skydim; i++ ) testsum += backproj->fitsimage[eeb][0][i]; RILlogMessage( NULL, Log_0,"##43## SUM fitsimage[eeb][0] = %f", testsum); testsum = 0.0; for( i = 0; i < backproj->skydim; i++ ) testsum += backproj->fitsimage[eeb][1][i]; RILlogMessage( NULL, Log_0,"##44## SUM fitsimage[eeb][1] = %f", testsum); testsum = 0.0; for( i = 0; i < backproj->skydim; i++ ) testsum += backproj->fitsimage[eeb][2][i]; RILlogMessage( NULL, Log_0,"##45## SUM fitsimage[eeb][2] = %f", testsum); testsum = 0.0; for( i = 0; i < backproj->skydim; i++ ) if (backproj->access[i] > 0) testsum += backproj->fitsimage[eeb][0][i]; RILlogMessage( NULL, Log_0,"##47## SUM fitsimage[eeb][0] = %f", testsum); testsum = 0.0; for( i = 0; i < backproj->skydim; i++ ) if (backproj->access[i] > 0) testsum += backproj->fitsimage[eeb][1][i]; RILlogMessage( NULL, Log_0,"##48## SUM fitsimage[eeb][1] = %f", testsum); testsum = 0.0; for( i = 0; i < backproj->skydim; i++ ) if (backproj->access[i] > 0) testsum += backproj->fitsimage[eeb][2][i]; RILlogMessage( NULL, Log_0,"##49## SUM fitsimage[eeb][2] = %f", testsum); } if (logger->test == 99) {n_sources = 1000; continue;} /* test witout requiring backpro8_d */ SCWnow->cand_lst_no[eeb] = n_sources - backproj->backmodels; SCWnow->fitcycles[eeb] = fitcycles; for (isrc=backproj->backmodels; isrcN_SRC].scw_ID = logger->N_SCW; CANDSRC[logger->N_SRC].ee = eeb; CANDSRC[logger->N_SRC].xx = xx[isrc]; CANDSRC[logger->N_SRC].yy = yy[isrc]; CANDSRC[logger->N_SRC].preset = preset[isrc]; CANDSRC[logger->N_SRC].src_id = src_id[isrc]; CANDSRC[logger->N_SRC].sxy = I_sxy[isrc]; CANDSRC[logger->N_SRC].src_strength = goodbet[isrc]; sprintf(logger->tmptxt, "CANDSRC[%4d]: ee: %1d, X/Y: %f %f, preset: %3d, src_id: %3d", logger->N_SRC, CANDSRC[logger->N_SRC].ee, CANDSRC[logger->N_SRC].xx, CANDSRC[logger->N_SRC].yy, CANDSRC[logger->N_SRC].preset, CANDSRC[logger->N_SRC].src_id); logger->logstat = logprint(logger, logger->tmptxt, i5); if (logger->tot_cnt_norm[eeb] == 0.0 ) { RILlogMessage( NULL, Warning_1,"##51## Attempt to divide with zero! reset to 1"); logger->tot_cnt_norm[eeb] = 1.0; } normfact = USERBIN[eeb].totcnt / logger->tot_cnt_norm[eeb]; logger->N_SRC++; if (logger->N_SRC > N_SRCMAX) { sprintf(logger->tmptxt, "##52## Too many source candidates: %5d", logger->N_SRC); logger->logstat = logprint(logger, logger->tmptxt, 7); status = J_ERROR_CODE_START+J_IMA_IROS_ERR+J_IIR_TOOMANYCAND; if (logger->trace) traces(func_local, 94000+eeb, logger); goto exittrace; } } SCWnow->backgr_str[eeb][0] = goodbet[0]; SCWnow->backgr_str[eeb][1] = goodbet[1]; SCWnow->backgr_area[eeb] = backproj->act_pix[eeb] / 100.0; SCWnow->fitrms[eeb] = fitrms; sh_rms[eeb][2] = sh_rms2; if (n_sources == 1000) continue; /* special signal for irosloop return */ /* ******************************************************** */ /* */ /* Analysis of this energy band completed. Store results */ /* for later comparison with the other energy bands. */ /* */ /* ******************************************************** */ sprintf(logger->tmptxt, " "); /* print newline */ logger->logstat = logprint(logger, logger->tmptxt, i5); sprintf(logger->tmptxt, "##53## JMX%1d E: %02d, t: %6.0f s, Const. fitpwr, %8.1lf, Backgr. fitpwr, %8.1lf", jmx_id->jmx_unit+1, eeb, sh->accumT, SCWnow->backgr_str[eeb][0], SCWnow->backgr_str[eeb][1]); logger->logstat = logprint(logger, logger->tmptxt, i5); for (isrc=backproj->backmodels; isrccand_lst_start[eeb] + isrc - backproj->backmodels; if(backproj->norm_d[isrc][eeb] == 0.0 ) { RILlogMessage( NULL, Warning_1,"##54## Attempt to divide with zero! reset to 1"); backproj->norm_d[isrc][eeb] = 1.0; } CANDSRC[n_src].src_pixel_area = 0.01 / backproj->norm_d[isrc][eeb]; if (logger->pointdata > 0) { /* get the background area corresponding to each source */ status = FUN(isrc, 0, eeb, xx, yy, jmx_id->jmx_unit, 0, &funval, jmx_id, logger, backproj, chatter, status); if (status != ISDC_OK) { logger->logstat = logprint(logger, logger->error_text, 7); status = J_ERROR_CODE_START+J_IMA_IROS_ERR+J_IIR_FUNERR; if (logger->trace) traces(func_local, 95000+eeb, logger); goto exittrace; } backcnt[eeb][isrc] = 2.0 * funval; /* determine RA and dec of the source candidates found at this energy */ xmm = (xx[isrc] - backproj->FOV_radius) * backproj->pix2mm; ymm = (yy[isrc] - backproj->FOV_radius) * backproj->pix2mm; getRAdec(xmm, ymm, jmx_id->jmx_unit, &s_RA, &s_dec, point, logger, jmx_id); CANDSRC[n_src].RA = s_RA; CANDSRC[n_src].dec = s_dec; CANDSRC[n_src].src_backgr_area = backcnt[eeb][isrc]; sprintf(logger->tmptxt, "##55## src: %1d, fitpwr: %8.1lf, XY: %+6.1f %+6.1f E: %2d RA/dec: %+7.3f %+7.3f, bkg_area: %6.1f, src_id: %3d", isrc, CANDSRC[n_src].src_strength, CANDSRC[n_src].xx, CANDSRC[n_src].yy, CANDSRC[n_src].ee, CANDSRC[n_src].RA, CANDSRC[n_src].dec, backcnt[eeb][isrc], CANDSRC[n_src].src_id); logger->logstat = logprint(logger, logger->tmptxt, i5); } else { sprintf(logger->tmptxt, "##56## src: %1d, fitpwr: %8.1lf, XY: %+6.1f %+6.1f E: %2d", isrc, CANDSRC[n_src].src_strength, CANDSRC[n_src].xx, CANDSRC[n_src].yy, CANDSRC[n_src].ee); logger->logstat = logprint(logger, logger->tmptxt, i5); } } if (logger->trace) traces(func_local, 96000+eeb, logger); } /* *************************************************************************** */ /* *************************************************************************** */ /* * * */ /* * END OF BLOCK PROCESSING 3 BASIC ENERGY BANDS FROM ONE SCIENCE WINDOW * */ /* * (Begins approximately in line 725) * */ /* * * */ /* *************************************************************************** */ /* *************************************************************************** */ /* 333333 */ /* QQQQQQ */ if (n_sources == 1000) { /* special signal for irosloop return */ status = J_ERROR_CODE_START+J_IMA_IROS_ERR+J_IIR_SPECIAL; goto exittrace; } /* software development test prints (ends approx line 960) */ for (eeb=0; eebee_basic; eeb++) { if (USERBIN[eeb].totcnt <= 1.0) { n_sources = backproj->backmodels; sh_rms[eeb][1] = 0.0; continue; } /* sprintf(logger->tmptxt, "##57## ngrp[%1d]: best srcnum: %2d, best fitcycles: %2d, best rms: %6.3f", eeb, SCWnow->cand_lst_no[eeb], SCWnow->fitcycles[eeb], SCWnow->fitrms[eeb]); logger->logstat = logprint(logger, logger->tmptxt, i5); */ sprintf(logger->tmptxt, "##58## Best fit coefficients : "); for (j=0; jcand_lst_no[eeb] + backproj->backmodels; j++) { n_src = SCWnow->cand_lst_start[eeb] + j - backproj->backmodels; if (j < backproj->backmodels) sprintf(logger->tmptxt2, " %7.1f", SCWnow->backgr_str[eeb][j]); else sprintf(logger->tmptxt2, " %7.1f", CANDSRC[n_src].src_strength); if ((strlen(logger->tmptxt) + strlen(logger->tmptxt2)) < MAX_STR_LEN) strcat(logger->tmptxt, logger->tmptxt2); } logger->logstat = logprint(logger, logger->tmptxt, i5); sprintf(logger->tmptxt, "##59## Corresp. eff. areas (cm2): "); for (j=0; jcand_lst_no[eeb] + backproj->backmodels; j++) { n_src = SCWnow->cand_lst_start[eeb] + j - backproj->backmodels; if (j < backproj->backmodels) sprintf(logger->tmptxt2, " %7.1f", SCWnow->backgr_area[eeb]); else sprintf(logger->tmptxt2, " %7.1f", CANDSRC[n_src].src_pixel_area); if ((strlen(logger->tmptxt) + strlen(logger->tmptxt2)) < MAX_STR_LEN) strcat(logger->tmptxt, logger->tmptxt2); } logger->logstat = logprint(logger, logger->tmptxt, i5); sprintf(logger->tmptxt, "##60## Integration time: %8.1lf (s), Absorbtion length: %5.2lf (mm)", sh->accumT, backproj->abs_len[eeb]); logger->logstat = logprint(logger, logger->tmptxt, i5); } for (i=0; iee_basic; e1++) used_ee[e1][i] = 0; used2[i] = 0; } sprintf(logger->tmptxt, "##61## Number of sources per E-band:"); for (eeb=0; eebee_basic; eeb++) { if (USERBIN[eeb].totcnt <= 1.0) continue; sprintf(logger->tmptxt2, " %2d", SCWnow->cand_lst_no[eeb]); if ((strlen(logger->tmptxt) + strlen(logger->tmptxt2)) < MAX_STR_LEN) strcat(logger->tmptxt, logger->tmptxt2); } logger->logstat = logprint(logger, logger->tmptxt, i5); e1_stop = e2_stop = ic_stop = jc_stop = 0; /* *************************************************************************** */ /* *************************************************************************** */ /* * * */ /* * * */ /* * START OF BLOCK IDENTIFYING COMMON SOURCES BETWEEN BASIC ENERGY BANDS * */ /* * (Ends approximately in line 1430) * */ /* * * */ /* * * */ /* *************************************************************************** */ /* *************************************************************************** */ /* 444444 */ /* QQQQQQ */ /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ kkkk = 0; logger->N_SRC2 = 0; for (e1=0; e1ee_basic; e1++) { initGOOD(e1, GOODSRC, logger, status); if (SCWnow->cand_lst_no[e1] == 0) continue; for (isrc=0; isrccand_lst_no[e1]; isrc++) { n_src1 = SCWnow->cand_lst_start[e1] + isrc; if (CANDSRC[n_src1].src_id > 0) backproj->s_detect[CANDSRC[n_src1].src_id] |= (1<<((2-e1)*4)) * (1+isrc); } } /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ /* e1-loop isrc-loop */ for (e1=0; e1ee_basic-1; e1++) { e1_stop = 1; if (SCWnow->cand_lst_no[e1] == 0) continue; e1_stop = 9; for (isrc=0; isrccand_lst_no[e1]; isrc++) { ic_stop = 1; if (used_ee[e1][isrc] > 0) continue; ic_stop = 2; n_src1 = SCWnow->cand_lst_start[e1] + isrc; if (SCWnow->backgr_area[e1] > 0.0) { divisor = sqrt((double)(USERBIN[e1].totcnt * CANDSRC[n_src1].src_pixel_area / SCWnow->backgr_area[e1])); if (divisor > 0.0) quality = (double)(CANDSRC[n_src1].src_strength) / divisor; else { quality = 0.0; sprintf(logger->tmptxt, "Normalization problems: ee: %2d, src#: %2d, Pixel area: %f, Total counts: %8.1f", e1, n_src1, CANDSRC[n_src1].src_pixel_area, USERBIN[e1].totcnt); logger->logstat = logprint(logger, logger->tmptxt, i5); RILlogMessage(NULL,Warning_1, "##63## divisor: ee: %2d, src#: %2d, totalcnt: %8.1f, pix_area: %f", e1, n_src1, USERBIN[e1].totcnt, CANDSRC[n_src1].src_pixel_area); } } else { quality = 0.0; sprintf(logger->tmptxt, "Normalization problems: SCWnow->backgr_area: %f", SCWnow->backgr_area[e1]); logger->logstat = logprint(logger, logger->tmptxt, 2); RILlogMessage(NULL,Warning_1,"##62## divisor: %f", SCWnow->backgr_area[e1]); } ic_stop = 9; sumX = CANDSRC[n_src1].xx * CANDSRC[n_src1].src_strength; sumY = CANDSRC[n_src1].yy * CANDSRC[n_src1].src_strength; sum_str = CANDSRC[n_src1].src_strength; sumgood = sum_str; tot_area = CANDSRC[n_src1].src_backgr_area; sumtot = tot_area; if( tot_area == 0.0 ) { sprintf(logger->tmptxt, "Normalization problems: CANDSRC->src_backgr_area: %f", tot_area); logger->logstat = logprint(logger, logger->tmptxt, 2); RILlogMessage( NULL, Warning_1,"##64## Attempt to divide with zero! reset to 1"); tot_area = 1.0; sumtot = 1.0; } sn_best = sumgood / sqrt(sumtot); CANDSRC[n_src1].best_sn = sn_best; GOODSRC[logger->N_SRC2].sn[e1] = sn_best; for (e2=e1+1; e2ee_basic; e2++) { e2_stop = 1; if (SCWnow->cand_lst_no[e2] == 0) continue; e2_stop = 9; /* e1-loop isrc-loop e2-loop jsrc-loop */ for (jsrc=0; jsrccand_lst_no[e2]; jsrc++) { jc_stop = 1; if (used_ee[e2][jsrc] > 0) continue; n_src2 = SCWnow->cand_lst_start[e2] + jsrc; divisor = SCWnow->backgr_area[e2]; if( divisor < 1.e-7 ) { sprintf(logger->tmptxt, "Normalization problems (e2): SCWnow->backgr_area: %f", divisor); logger->logstat = logprint(logger, logger->tmptxt, 2); RILlogMessage(NULL,Warning_1,"##65## background area < 1.e-7; STOP !"); status = -35; goto exittrace; } jc_stop = 2; if ((double)(CANDSRC[n_src2].src_strength) < sqrt((double)(USERBIN[e2].totcnt * CANDSRC[n_src2].src_pixel_area / divisor)) && (CANDSRC[n_src1].preset == 0)) { sprintf(logger->tmptxt, "##971## e2: %1d, jsrc: %1d, src.no: %2d, source strength: %f", e2, jsrc, n_src2, CANDSRC[n_src2].src_strength); logger->logstat = logprint(logger, logger->tmptxt, 2); sprintf(logger->tmptxt, "##1971## e2: %1d, jsrc: %1d, Quality problems: Pixel area: %f, Total count: %8.1f, backgr.area: %f", e2, jsrc, CANDSRC[n_src2].src_pixel_area, USERBIN[e2].totcnt , SCWnow->backgr_area[e2]); logger->logstat = logprint(logger, logger->tmptxt, 2); continue; } jc_stop = 9; xx_dif = CANDSRC[n_src1].xx - CANDSRC[n_src2].xx; yy_dif = CANDSRC[n_src1].yy - CANDSRC[n_src2].yy; /* e1-loop isrc-loop e2-loop jsrc-loop sqrt(dx2+dy2) < 1.5 */ if (sqrt((double)(xx_dif * xx_dif + yy_dif * yy_dif)) < 1.5) { tempgood = sum_str + CANDSRC[n_src2].src_strength; temptot = tot_area + CANDSRC[n_src2].src_backgr_area; if( temptot == 0.0 ) { sprintf(logger->tmptxt, "Normalization problems: total backgr area: %f", temptot); logger->logstat = logprint(logger, logger->tmptxt, 2); RILlogMessage( NULL, Warning_1,"##66## Attempt to divide with zero! reset to 1"); temptot = 1.0; } sn_temp = tempgood / sqrt((double)(temptot)); if( CANDSRC[n_src2].src_backgr_area <= 0.0 ) { sprintf(logger->tmptxt, "Normalization problems: CANDSRC->src_backgr_area: %f", CANDSRC[n_src2].src_backgr_area); logger->logstat = logprint(logger, logger->tmptxt, 2); RILlogMessage( NULL, Warning_1,"##67## Attempt to divide with zero! reset to 1"); CANDSRC[n_src2].src_backgr_area = 1.0; } sn_new = CANDSRC[n_src2].src_strength / sqrt((double)(CANDSRC[n_src2].src_backgr_area)); if ((sn_temp >= sn_best) && (sn_temp >= sn_new )) { sn_best = sn_temp; sum_str = tempgood; tot_area = temptot; } if ((sn_new >= sn_best) && (sn_new >= sn_temp)) { sn_best = sn_new; sum_str = CANDSRC[n_src2].src_strength; tot_area = USERBIN[e2].totcnt; } sumgood += CANDSRC[n_src2].src_strength; sumX += CANDSRC[n_src2].xx * CANDSRC[n_src2].src_strength; sumY += CANDSRC[n_src2].yy * CANDSRC[n_src2].src_strength; if( sumgood == 0.0 ) { sprintf(logger->tmptxt, "Normalization problems: sumgood: %f", sumgood); logger->logstat = logprint(logger, logger->tmptxt, 2); RILlogMessage( NULL, Warning_1,"##68## Attempt to divide with zero! reset to 1"); sumgood = 1.0; } avX0 = sumX / sumgood; avY0 = sumY / sumgood; logger->tmptxt[0] = 0; if (logger->pointdata > 0) { /* e1-loop isrc-loop e2-loop jsrc-loop sqrt(dx2+dy2) < 1.5 pointing data available */ xmm = (avX0 - backproj->FOV_radius) * backproj->pix2mm; ymm = (avY0 - backproj->FOV_radius) * backproj->pix2mm; getRAdec(xmm, ymm, jmx_id->jmx_unit, &s_RA, &s_dec, point, logger, jmx_id); sprintf(logger->tmptxt, "##69##"); sprintf(logger->tmptxt2, " RA/dec: %+7.3f %+7.3f E%1d fit: %7.1lf, E%1d fit: %7.1lf ", s_RA, s_dec, CANDSRC[n_src1].ee, CANDSRC[n_src1].src_strength, CANDSRC[n_src2].ee, CANDSRC[n_src2].src_strength); if ((strlen(logger->tmptxt) + strlen(logger->tmptxt2)) < MAX_STR_LEN) strcat(logger->tmptxt, logger->tmptxt2); } else { sprintf(logger->tmptxt2, "##70## XY: %+7.1f %+7.1f strengths: %7.1lf %7.1lf ", avX0-(float)(backproj->FOV_radius), avY0-(float)(backproj->FOV_radius), CANDSRC[n_src1].src_strength, CANDSRC[n_src2].src_strength); if ((strlen(logger->tmptxt) + strlen(logger->tmptxt2)) < MAX_STR_LEN) strcat(logger->tmptxt, logger->tmptxt2); } logger->logstat = logprint(logger, logger->tmptxt, 2); CANDSRC[n_src1].best_sn = sn_best; used_ee[e2][jsrc] = 1; used2[isrc] = 1; backproj->com_RA[kkkk] = s_RA; backproj->com_dec[kkkk] = s_dec; GOODSRC[logger->N_SRC2].xx = avX0; GOODSRC[logger->N_SRC2].yy = avY0; GOODSRC[logger->N_SRC2].RA = s_RA; GOODSRC[logger->N_SRC2].dec = s_dec; GOODSRC[logger->N_SRC2].preset |= CANDSRC[n_src1].preset; /* GOODSRC[logger->N_SRC2].src_id = CANDSRC[n_src1].src_id; */ /* src_id for verified sources should be identified on the basis of the average coordinates avX0, avY0 */ GOODSRC[logger->N_SRC2].srcnum = CANDSRC[n_src1].srcnum; GOODSRC[logger->N_SRC2].scw_ID = CANDSRC[n_src1].scw_ID; GOODSRC[logger->N_SRC2].sxy = CANDSRC[n_src1].sxy; GOODSRC[logger->N_SRC2].sn[e2] = sn_best; GOODSRC[logger->N_SRC2].detects |= (1<<((2-e1)*4)) * (1+isrc); GOODSRC[logger->N_SRC2].detects |= (1<<((2-e2)*4)) * (1+jsrc); GOODSRC[logger->N_SRC2].preset |= (4<N_SRC2].preset |= (4<N_SRC2].preset |= 32; sprintf(logger->tmptxt, "##71## Energy scales: %d %d, kkkk: %2d, scw_ID: %3d, N_SRC: %4d, N_SRC2: %3d, XY0: %6.1f %6.1f, XY1: %6.1f %6.1f RA/dec: %7.3lf %7.3lf", CANDSRC[n_src1].ee, CANDSRC[n_src2].ee, kkkk, GOODSRC[logger->N_SRC2].scw_ID, logger->N_SRC, logger->N_SRC2, CANDSRC[n_src1].xx, CANDSRC[n_src1].yy, CANDSRC[n_src2].xx, CANDSRC[n_src2].yy, GOODSRC[logger->N_SRC2].RA, GOODSRC[logger->N_SRC2].dec); logger->logstat = logprint(logger, logger->tmptxt, i5); jc_stop = 9; /* e1-loop isrc-loop e2-loop jsrc-loop sqrt(dx2+dy2) < 1.5 */ } jc_stop = 4; } e2_stop = 9; } if (used2[isrc] > 0) { used_ee[e1][isrc] = 1; logger->N_SRC2++; /* initGOOD(e1, GOODSRC, logger, status); */ if (logger->trace) traces(func_local, 10000+logger->N_SRC2, logger); kkkk++; if (kkkk > SCW_SRCMAX) {kkkk = SCW_SRCMAX; goto sourcefull;} } used2[isrc] = 0; } e1_stop = 9; } /* *************************************************************************** */ /* *************************************************************************** */ /* * * */ /* * * */ /* * END OF BLOCK IDENTIFYING COMMON SOURCES BETWEEN BASIC ENERGY BANDS * */ /* * (Begins approximately in line 1175) * */ /* * * */ /* * * */ /* *************************************************************************** */ /* *************************************************************************** */ sourcefull: /* sourcefull */ sprintf(logger->tmptxt, "##961## end of e1 loop: %2d, no of common sources found: %2d, stopcodes: %1d %1d", e1, logger->N_SRC2, e1_stop, e2_stop); logger->logstat = logprint(logger, logger->tmptxt, i5); n_common = kkkk; /* 555555 */ /* QQQQQQ */ /* *************************************************************************** */ /* *************************************************************************** */ /* * * */ /* * * */ /* * START OF BLOCK IDENTIFYING STRONG SOURCES IN SINGLE ENERGY BANDS * */ /* * (Ends approximately in line 1560) * */ /* * * */ /* * * */ /* *************************************************************************** */ /* *************************************************************************** */ for (e1=0; e1ee_basic; e1++) { if (SCWnow->cand_lst_no[e1] > 0) { for (isrc=0; isrccand_lst_no[e1]; isrc++) { if (used_ee[e1][isrc] > 0) continue; n_src1 = SCWnow->cand_lst_start[e1] + isrc; if(USERBIN[e1].totcnt <= 0.0 ) { RILlogMessage( NULL, Warning_1,"##72## Attempt to divide with zero totcnt!to 1"); USERBIN[e1].totcnt = 1.0; } quality = (double)(CANDSRC[n_src1].src_strength) / sqrt((double)(USERBIN[e1].totcnt)); sprintf(logger->tmptxt, "##73## Single band source: E%1d, quality: %5.1lf, limit: %5.1f", e1, quality, logger->qual_limit); if (quality >= logger->qual_limit) sprintf(logger->tmptxt2, " Accepted"); else sprintf(logger->tmptxt2, " NOT accepted"); if ((strlen(logger->tmptxt) + strlen(logger->tmptxt2)) < MAX_STR_LEN) strcat(logger->tmptxt, logger->tmptxt2); logger->logstat = logprint(logger, logger->tmptxt, 2); if (quality > logger->qual_limit) { sprintf(logger->tmptxt, "##74## Single band strong source: quality: %5.1lf ", quality); if (logger->pointdata > 0) { xmm = (CANDSRC[n_src1].xx - backproj->FOV_radius) * backproj->pix2mm; ymm = (CANDSRC[n_src1].yy - backproj->FOV_radius) * backproj->pix2mm; getRAdec(xmm, ymm, jmx_id->jmx_unit, &s_RA, &s_dec, point, logger, jmx_id); sprintf(logger->tmptxt2, "##75## RA/dec: %+7.3f %+7.3f ee: %1d fit: %7.1lf ", s_RA, s_dec, e1, CANDSRC[n_src1].src_strength); if ((strlen(logger->tmptxt) + strlen(logger->tmptxt2)) < MAX_STR_LEN) strcat(logger->tmptxt, logger->tmptxt2); backproj->com_RA[kkkk] = s_RA; backproj->com_dec[kkkk] = s_dec; } else { sprintf(logger->tmptxt2, " XY: %+7.1f %+7.1f strengths: %7.1lf ", CANDSRC[n_src1].xx - (float)(backproj->FOV_radius), CANDSRC[n_src1].yy - (float)(backproj->FOV_radius), CANDSRC[n_src1].src_strength); if ((strlen(logger->tmptxt) + strlen(logger->tmptxt2)) < MAX_STR_LEN) strcat(logger->tmptxt, logger->tmptxt2); backproj->com_RA[kkkk] = 0.0; backproj->com_dec[kkkk] = 0.0; } sprintf(logger->tmptxt2, "CR2"); if ((strlen(logger->tmptxt) + strlen(logger->tmptxt2)) < MAX_STR_LEN) strcat(logger->tmptxt, logger->tmptxt2); logger->logstat = logprint(logger, logger->tmptxt, 2); sumgood = CANDSRC[n_src1].src_strength; sumtot = CANDSRC[n_src1].src_backgr_area; if( sumtot == 0.0 ) { RILlogMessage( NULL, Warning_1,"##76## Attempt to divide with zero! reset to 1"); sumtot = 1.0; } sn_best = sumgood / sqrt(sumtot); CANDSRC[n_src1].best_sn = sn_best; GOODSRC[logger->N_SRC2].sn[e1] = sn_best; GOODSRC[logger->N_SRC2].xx = CANDSRC[n_src1].xx; GOODSRC[logger->N_SRC2].yy = CANDSRC[n_src1].yy; GOODSRC[logger->N_SRC2].RA = CANDSRC[n_src1].RA; GOODSRC[logger->N_SRC2].dec = CANDSRC[n_src1].dec; GOODSRC[logger->N_SRC2].scw_ID = CANDSRC[n_src1].scw_ID; GOODSRC[logger->N_SRC2].sxy = CANDSRC[n_src1].sxy; GOODSRC[logger->N_SRC2].preset |= CANDSRC[n_src1].preset; /* GOODSRC[logger->N_SRC2].src_id = CANDSRC[n_src1].src_id; * src_id for verified sources should be identified on the basis * of the average coordinates avX0, avY0 */ GOODSRC[logger->N_SRC2].detects |= (1<<((2-e1)*4)) * (1+isrc); GOODSRC[logger->N_SRC2].preset |= (4<tmptxt, "##59a# Energy scale: %d scw_ID: %3d, N_SRC: %4d, N_SRC2: %3d, XY0: %6.1f %6.1f ", e1, GOODSRC[logger->N_SRC2].scw_ID, n_src1, logger->N_SRC2, CANDSRC[n_src1].xx, CANDSRC[n_src1].yy); logger->logstat = logprint(logger, logger->tmptxt, 2); kkkk++; if (kkkk > SCW_SRCMAX) {kkkk = SCW_SRCMAX; goto sourcefull2;} logger->N_SRC2++; /* initGOOD(e1, GOODSRC, logger, status); */ if (logger->trace) traces(func_local, 20000+logger->N_SRC2, logger); } } } } /* ************************************************************************* */ /* *************************************************************************** */ /* * * */ /* * * */ /* * END OF BLOCK IDENTIFYING STRONG SOURCES IN SINGLE ENERGY BANDS * */ /* * (Begins approximately in line 1465) * */ /* * * */ /* * * */ /* *************************************************************************** */ /* *************************************************************************** */ /* 666666 */ /* QQQQQQ */ sourcefull2: /* sourcefull2 */ n_common2 = kkkk; sprintf(logger->tmptxt, "##98765## After sourcefind, n_common2: %2d, N_SRC2: %2d", n_common2, logger->N_SRC2); logger->logstat = logprint(logger, logger->tmptxt, i5); for (i=0; iN_SRC2; i++) { sprintf(logger->tmptxt, "##98765## i: %2d, RA/dec: %7.3f %7.3f, name: %30s", i, GOODSRC[i].RA, GOODSRC[i].dec, GOODSRC[i].name); logger->logstat = logprint(logger, logger->tmptxt, i5); } /* *********************************************************************************** */ /* *********************************************************************************** */ /* * * */ /* * * */ /* * Identify catalog sources among the sources found by ima_iros * */ /* * (for such sources use the coordinates from the catalog for the flux fitting) * */ /* * (ends approx in line 1650) * */ /* * * */ /* *********************************************************************************** */ /* *********************************************************************************** */ sprintf(logger->tmptxt, "Number of catalog sources: %3d, number of found sources: %3d", backproj->preset_src, n_common2); logger->logstat = logprint(logger, logger->tmptxt, 0); if (backproj->preset_src > 0) { for (j=0; jpreset_src; j++) { for (i=0; iN_SRC2; i++) { diff_x = backproj->preset_src_xx[j] - GOODSRC[i].xx; diff_y = backproj->preset_src_yy[j] - GOODSRC[i].yy; diff = diff_x * diff_x + diff_y *diff_y; diff = sqrt(diff); if (diff < 1.5) { backproj->preset_src_flag[j] |= 8; GOODSRC[i].catxx = backproj->preset_src_xx[j]; GOODSRC[i].catyy = backproj->preset_src_yy[j]; GOODSRC[i].catRA = backproj->preset_src_RA[j]; GOODSRC[i].catdec = backproj->preset_src_dec[j]; GOODSRC[i].src_id = backproj->preset_src_id[j]; if ((GOODSRC[i].preset & 3) == 3) { GOODSRC[i].diff_x = 0.0; GOODSRC[i].diff_y = 0.0; GOODSRC[i].preset |= 64; } else { GOODSRC[i].diff_x = diff_x; GOODSRC[i].diff_y = diff_y; GOODSRC[i].preset |= 64 + 1024 + (backproj->preset_src_flag[j] & 3); } goto next_preset_src; } else { if ((GOODSRC[i].preset & 3) != 3) GOODSRC[i].preset |= 1024; } } next_preset_src: ; /* next_preset_src */ } } if (n_common2 > n_common) { sprintf(logger->tmptxt, "##77## Analysis using %d common and %d strong source(s) between %d energy bands", n_common, n_common2-n_common, backproj->ee_basic); } else { RILlogMessage(NULL,Log_0,"##78## n_common2 = %d", n_common2 ); sprintf(logger->tmptxt, "##79## Analysis using %d common source(s) between %d energy bands", n_common, backproj->ee_basic); } logger->logstat = logprint(logger, logger->tmptxt, 2); if (logger->trace) traces(func_local, 33, logger); ksrc = backproj->backmodels + n_common2; for (n_src=0; n_srcbackmodels+n_src] = GOODSRC[n_src].catxx; yy[backproj->backmodels+n_src] = GOODSRC[n_src].catyy; } else { xx[backproj->backmodels+n_src] = GOODSRC[n_src].xx; yy[backproj->backmodels+n_src] = GOODSRC[n_src].yy; } sprintf(logger->tmptxt, "source %1d: X: %6.2f, Y: %6.2f", n_src, xx[backproj->backmodels+n_src], yy[backproj->backmodels+n_src]); logger->logstat = logprint(logger, logger->tmptxt, i5); } if (n_common2 == 0) { sprintf(logger->tmptxt, "##80## No common sources in the four energy bands"); logger->logstat = logprint(logger, logger->tmptxt, 2); } if( sh->accumT == 0.0 ) { RILlogMessage( NULL, Warning_1, "##81## Attempt to divide with accumT: (%f)! reset to 1", sh->accumT); sh->accumT = 1.0; } /* ***************************************************************************** */ /* ***************************************************************************** */ /* * * */ /* * * */ /* * FITTING LOOP FOR BASIC SOURCE SET OVER SEARCH ENERGY INTERVALS * */ /* * (Ends approximately in line 1910) * */ /* * * */ /* * * */ /* ***************************************************************************** */ /* ***************************************************************************** */ /* 777777 */ /* QQQQQQ */ for (i=0; isrclist[i] = -1; for (i=0; isrclist[i] = i; ShNNN = 0; for (i=0; isoft_gr[i] > 0.0) ShNNN++; for (eeuser=0; eeuseree_basic; eeuser++) { if (eeuser == 0) logger->new_detmap = 1; else logger->new_detmap = 0; if (logger->trace) traces(func_local, (int)(6300000+eeuser), logger); backproj->additional = -1; status = final_fit( eeuser, n_common2, solut_e, jmx_id, logger, backproj, sh, SCWnow, GOODSRC, USERBIN, chatter, status); if (logger->trace) traces(func_local, 99180, logger); kMAX = backproj->backmodels + n_common2; if (kMAX > BASIC_SRCMAX) kMAX = BASIC_SRCMAX; if (logger->trace) traces(func_local, 1000+kMAX, logger); /* output fit functions as FITS */ for (allsrc=0; allsrctmptxt2, "Src %1d E %1d FITmap * 10**4", allsrc, eeuser); strcpy(f_exttext[1], logger->tmptxt2); ShSUM = 0.0; for(i=0; i<256; i++){ for (j=0; j<256; j++) { ij = i + j*256; if (backproj->soft_gr[ij] <= 0.0) shdg_fits[ij] = NAN; else { shdg_fits[ij] = logger->fit_funcs[allsrc][ij] * solut_e[eeuser][allsrc]; ShSUM += shdg_fits[ij]; } } } if (logger->trace) traces(func_local, 181+allsrc, logger); sprintf(logger->tmptxt, "##EE: %2d, fit_function%1d: sum: %7.1f, solut_e: %7.1f, used_pixels: %5d", eeuser, allsrc, ShSUM/solut_e[eeuser][allsrc], solut_e[eeuser][allsrc], ShNNN); logger->logstat = logprint(logger, logger->tmptxt, 0); status = shdg_to_fits(theSWG, eeuser, shdg_fits, jemxNum, elmntName, f_exttext, f_struc1, logger, eeuser, sh, chatter, status ); if (status != ISDC_OK) {status= -70; goto exittrace; } } /* output data shadowgram as FITS */ strcpy(f_exttext[0], "INPUT SHDGRM "); sprintf(logger->tmptxt2, "InputShdg. E %2d", eeuser); strcpy(f_exttext[1], logger->tmptxt2); ShSUM = 0.0; for(i=0; i<256; i++){ for (j=0; j<256; j++) { ij = i + j*256; if (backproj->soft_gr[ij] <= 0.0) shdg_fits[ij] = NAN; else { shdg_fits[ij] = logger->fit_funcs[BASIC_SRCMAX+0][ij]; ShSUM += shdg_fits[ij]; } } } if (logger->trace) traces(func_local, 190, logger); sprintf(logger->tmptxt, "EE: %2d, input shadgrm: sum: %7.1f, used_pixels: %5d", eeuser, ShSUM, ShNNN); logger->logstat = logprint(logger, logger->tmptxt, i5); status = shdg_to_fits(theSWG, eeuser, shdg_fits, jemxNum, elmntName, f_exttext, f_struc1, logger, eeuser, sh, chatter, status ); if (status != ISDC_OK) {status= -71; goto exittrace; } /* output residual shadowgram as FITS */ strcpy(f_exttext[0], "RESIDUAL SHDGRM "); sprintf(logger->tmptxt2, "Residual E %2d", eeuser); strcpy(f_exttext[1], logger->tmptxt2); ShSUM = 0.0; for(i=0; i<256; i++){ for (j=0; j<256; j++) { ij = i + j*256; if (backproj->soft_gr[ij] <= 0.0) shdg_fits[ij] = NAN; else { shdg_fits[ij] = logger->fit_funcs[BASIC_SRCMAX+1][ij]; ShSUM += shdg_fits[ij]; } } } if (logger->trace) traces(func_local, 191, logger); sprintf(logger->tmptxt, "EE: %2d, resid shadgrm: sum: %7.1f, used_pixels: %5d", eeuser, ShSUM, ShNNN); logger->logstat = logprint(logger, logger->tmptxt, i5); status = shdg_to_fits(theSWG, eeuser, shdg_fits, jemxNum, elmntName, f_exttext, f_struc1, logger, eeuser, sh, chatter, status ); if (status != ISDC_OK) {status= -72; goto exittrace; } /* output efficiency shadowgram as FITS */ strcpy(f_exttext[0], "EFFICIENCY SHDGRM "); sprintf(logger->tmptxt2, "Efficiency E %2d", eeuser); strcpy(f_exttext[1], logger->tmptxt2); ShSUM = 0.0; for(i=0; i<256; i++){ for (j=0; j<256; j++) { ij = i + j*256; if (backproj->soft_gr[ij] <= 0.0) shdg_fits[ij] = NAN; else { shdg_fits[ij] = logger->fit_funcs[BASIC_SRCMAX+2][ij]; ShSUM += shdg_fits[ij]; } } } if (logger->trace) traces(func_local, 192, logger); /* sprintf(logger->tmptxt, "EE: %2d, effic shadgrm: sum: %7.1f, used_pixels: %5d", eeuser, ShSUM, ShNNN); logger->logstat = logprint(logger, logger->tmptxt, i5); */ status = shdg_to_fits(theSWG, eeuser, shdg_fits, jemxNum, elmntName, f_exttext, f_struc1, logger, eeuser, sh, chatter, status ); if (status != ISDC_OK) {status= -72; goto exittrace; } for (jsrc=0; jsrcbackmodels+n_common2; jsrc++) { sprintf(logger->tmptxt, "##3217a## E: %2d, src: %3d, solut_e: %f ", eeuser, solut_e[eeuser][jsrc]); logger->logstat = logprint(logger, logger->tmptxt, i5); } sh_rms[eeuser][3] = rmscalc(sh->resid, USERBIN[eeuser].soft_gr, logger); if (logger->shdg_out_ctrl == 1) { /* output residual shadowgram as FITS */ strcpy(f_exttext[0], "RESIDUAL_SHDGR"); strcpy(f_exttext[1], "Residual shadowgram"); status = shdg_to_fits(theSWG, eeuser, sh->resid, jemxNum, elmntName, f_exttext, f_struc1, logger, eeb, sh, chatter, status ); if (status != ISDC_OK) goto exittrace; } /* ************************************************************** */ /* */ /* Use final residual shadowgram to form sky image */ /* */ /* ************************************************************** */ /* QQQQQQ */ /* 757575 */ if (logger->trace) traces(func_local, 4000+logger->detImagesOut, logger); if (logger->trace) traces(func_local, 5000+logger->userImagesOut, logger); if (logger->trace) traces(func_local, 6000+eeuser, logger); logger->detImagesOut = 0; logger->userImagesOut = 1; for (i=0; ilocalima[i] = 0.0; status = final_sky( eeuser, n_common2, peakval, theSWG, jmx_id, logger, backproj, sh, bp, SCWnow, GOODSRC, USERBIN, chatter, status ); sprintf(logger->tmptxt, "Return from final_sky A, status: %d", status); logger->logstat = logprint(logger, logger->tmptxt, 0); if (status != ISDC_OK) goto exittrace; if (logger->trace) traces(func_local, 6000, logger); if (logger->trace) traces(func_local, 66000 + status, logger); status = fluxerror( eeuser, n_common2, peakval, logger, backproj, sh, GOODSRC, status ); sprintf(logger->tmptxt, "Return from fluxerror A, E: %2d, status: %d", eeuser, status); logger->logstat = logprint(logger, logger->tmptxt, 0); if (logger->trace) traces(func_local, 6100, logger); if (logger->trace) traces(func_local, 6100 + status, logger); kMAX = ksrc; if (kMAX > BASIC_SRCMAX) kMAX = BASIC_SRCMAX; if (logger->trace) traces(func_local, 99180, logger); if (logger->trace) traces(func_local, 2000+kMAX, logger); } /* end of ee-loop over basic energy bands */ /* *************************************************************************** */ /* *************************************************************************** */ /* * * */ /* * * */ /* * END OF BLOCK PROCESSING BASIC SOURCE SET FOR SEARCH ENERGY BANDS * */ /* * (begins approximately in line 1713 * */ /* * * */ /* * * */ /* *************************************************************************** */ /* *************************************************************************** */ /* 888888 */ /* QQQQQQ */ /* identify found sources and find 'common names' */ sprintf(logger->tmptxt, "Start source ID loop"); logger->logstat = logprint(logger, logger->tmptxt, 0); for (allsrc=backproj->backmodels; allsrcbackmodels; sprintf(GOODSRC[ALLsrc].name, "US%02d", ALLsrc+1); evRA = GOODSRC[ALLsrc].RA; evdec = GOODSRC[ALLsrc].dec; sprintf(logger->tmptxt, "Before srcidentify2, ALLsrc: %2d, RA/dec: %7.3f %7.3f, %s", ALLsrc, evRA, evdec, GOODSRC[ALLsrc].name); logger->logstat = logprint(logger, logger->tmptxt, 0); srcidentify2(theSWG, evRA, evdec, GOODSRC[ALLsrc].name, GOODSRC[ALLsrc].ISDCname, &catflag, jmx_id, point, logger); sprintf(logger->tmptxt, "Return from srcidentify2, ALLsrc: %2d, RA/dec: %7.3f %7.3f, %s", ALLsrc, evRA, evdec, GOODSRC[ALLsrc].name); logger->logstat = logprint(logger, logger->tmptxt, 0); GOODSRC[ALLsrc].ISDCname[16] = 0; sflux = GOODSRC[ALLsrc].flux[0] + GOODSRC[ALLsrc].flux[1] + GOODSRC[ALLsrc].flux[2]; /* save text for 'regfil1' */ cmx = GOODSRC[ALLsrc].xx - backproj->FOV_radius; cmy = GOODSRC[ALLsrc].yy - backproj->FOV_radius; rr = sqrt(cmx*cmx + cmy*cmy); status = getxy(evRA, evdec, jmx_id->jmx_unit, &xmm, &ymm, jmx_id, point, logger); if ((fabs(xmm - cmx * backproj->pix2mm) > 0.1) || (fabs(ymm - cmy * backproj->pix2mm) > 0.1)) { sprintf(logger->tmptxt, "Inconsistent RA/dec A"); logger->logstat = logprint(logger, logger->tmptxt, 0); } sprintf(logger->tmptxt, "Return from getxy, ALLsrc: %2d", ALLsrc); logger->logstat = logprint(logger, logger->tmptxt, 0); sprintf(GOODSRC[ALLsrc].regsrctxt, "fk5;point(%8.4lf, %8.4lf) # point=cross 20 text={%s} color=white", evRA, evdec, GOODSRC[ALLsrc].name); if ((logger->develop & 1) == 1) { fprintf(regfil1, "%s\n", GOODSRC[ALLsrc].regsrctxt); fprintf(regfil1, "#fk5x IJDstart/stop %8.2lf %8.2lf, Obstime: %7.1f, cm2: %5.1f, R: %6.2f, flux: %5.1f, XY: %5.1f %5.1f\n", logger->gtiStartIJD, logger->gtiStopIJD, sh->accumT, GOODSRC[ALLsrc].src_pixel_area[0], rr, sflux, GOODSRC[ALLsrc].xx, GOODSRC[ALLsrc].yy); } GOODSRC[ALLsrc].reg_used = 0xFFFF; sprintf(logger->tmptxt, " Found source. SCW ID#: %02d, RA/dec: %7.3f %7.3f, X/Y/R: %5.1f %5.1f %5.1f, name: %s %s", allsrc, evRA, evdec, GOODSRC[ALLsrc].xx, GOODSRC[ALLsrc].yy, rr, GOODSRC[ALLsrc].name, GOODSRC[ALLsrc].ISDCname); logger->logstat = logprint(logger, logger->tmptxt, 0); if (eeuser == backproj->ee_basic) { logger->allsrc_ptr[allsrc-backproj->backmodels] = allsrc - backproj->backmodels;; } } if ((logger->develop & 1) == 1) fclose(regfil1); /* *************************************************************************** */ /* *************************************************************************** */ /* * * */ /* * * */ /* * FITTING LOOP OVER USER ENERGY INTERVALS * */ /* * first perform source fits over basic source set * */ /* * next perform new fit including one extra preset source in addition * */ /* * to the basic set for every step of the loop * */ /* * (ends approximately in line 3360) * */ /* * * */ /* * * */ /* *************************************************************************** */ /* *************************************************************************** */ /* 999999 */ /* QQQQQQ */ /* LOOP OVER USER ENERGY BANDS -- ENDS in approx. line 3560 */ for (eeuser=backproj->ee_basic; eeusernumShds; eeuser++) { for (ntime=0; ntimeNumTimeBins; ntime++) { Tslices[ntime].dontUse &= 1; /* mask away 'used time bin' flag when starting on a new energy band */ } sprintf(regnam2, "J%1d_%04d%04d%04d_E%02d.reg", jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, eeuser); if ((logger->develop & 1) == 1) { if ((regfil2 = fopen(regnam2, "wt")) != NULL) { fprintf(regfil2, "global color=red font=\"helvetica 8 normal\" \n"); fclose(regfil2); } } prevent = 0; logger->eeuser = eeuser; status = j_fitprep(eeuser, theSWG, jmx_id, logger, backproj, sh, bp, SCWnow, USERBIN, chatter, status); if (logger->trace) traces(func_local, 45, logger); if (logger->trace) traces(func_local, 45000+eeuser, logger); if (logger->trace) traces(func_local, 45200+sh->numShds, logger); sprintf(logger->tmptxt, "work2d: after fitprep: E: %2d, totcnt: %7.1f", eeuser, sh->totcnt[eeuser]); logger->logstat = logprint(logger, logger->tmptxt, i5); backproj->additional = -1; status = final_fit( eeuser, n_common2, solut_e, jmx_id, logger, backproj, sh, SCWnow, GOODSRC, USERBIN, chatter, status); kMAX = backproj->backmodels + n_common2; if (kMAX > BASIC_SRCMAX) kMAX = BASIC_SRCMAX; if (logger->trace) traces(func_local, 99180, logger); if (logger->trace) traces(func_local, 3000+kMAX, logger); if (logger->trace) traces(func_local, (int)(GOODSRC[logger->N_SRC2-1].src_pixel_area[eeuser])+2, logger); if ((logger->LC == 0) || (logger->numEvents < logger->NumTimeBins)) { if (logger->trace) traces(func_local, 1000000+logger->numEvents, logger); if (logger->trace) traces(func_local, 1000000+logger->NumTimeBins, logger); goto NoTimeBins1; } /* ************************************************************** */ /* CALCULATE LIGHT CURVES FOR THE BASIC SOURCE SET */ /* ************************************************************** */ /* QQQQQQ */ /* calculate electronic efficiency for each timeslice */ old_avGain =0.0; for (ntime=0; ntimeNumTimeBins; ntime++) { if ((fabs(Tslices[ntime].avGain - old_avGain) < 0.1) || (Tslices[ntime].dontUse > 0) || (Tslices[ntime].avGain < 0.5)) { if (ntime == 0) for (i=0; itrace) traces(func_local, 7771, logger); status = ex2trEvTslc(0, eeuser, &(Tslices[0]), &(bufevents[0]), &(PIFmap[0]), sh, backproj, USERBIN, logger, status); backproj->avEE_events = sum2 = 0.0; j = 0; for (ntime=0; ntimeNumTimeBins; ntime++) { if (Tslices[ntime].dontUse) continue; if (Tslices[ntime].exposure < logger->timestep * 0.01) {Tslices[ntime].dontUse = 1; continue;} contrib = Tslices[ntime].Events * logger->timestep * 86400.0 / Tslices[ntime].exposure ; backproj->avEE_events += contrib; sum2 += contrib * contrib; j++; } sprintf(logger->tmptxt, "E%1d, NumTimeBinsA0: %6d, sumEvents: %8.1f, sum-square: %e, bins: %4d", eeuser, logger->NumTimeBins, backproj->avEE_events, sum2, j); logger->logstat = logprint(logger, logger->tmptxt, i5); /**/ if (j == 0) {NOTIMEBINS = 1; goto NoTimeBins1;} /**/ backproj->avEE_events /= (float)(j); if (j < 5) rmsFlux0 = sqrt(backproj->avEE_events); else rmsFlux0 = sqrt(sum2 / (float)(j) - backproj->avEE_events * backproj->avEE_events); if (j < 5) goto FewTimeBins0; sprintf(logger->tmptxt, "E%1d, NumTimeBinsA: %6d, meanEvents: %6.1lf, rms: %6.2lf, bins: %4d", eeuser, logger->NumTimeBins, backproj->avEE_events, rmsFlux0, j); logger->logstat = logprint(logger, logger->tmptxt, i5); /* re-calculate mean and rms after excluding strongly deviant time bins */ j = 0; sum = wsum = sum2 = 0.0; for (ntime=0; ntimeNumTimeBins; ntime++) { if (Tslices[ntime].dontUse) continue; contrib = Tslices[ntime].Events * logger->timestep * 86400.0 / Tslices[ntime].exposure; if (fabs(contrib - backproj->avEE_events) > 3.0 * rmsFlux0) continue; /* do not include strongly deviant time bins in calculation of mean count rate and rmsFlux0 */ sum += contrib; sumE += Tslices[ntime].Events; sum2 += contrib * contrib; j++; } if (j == 0) {NOTIMEBINS = 1; goto NoTimeBins1;} backproj->avEE_events = sum / (float)(j); if (j == 1) rmsFlux0 = sqrt(sumE); else rmsFlux0 = sqrt(sum2 / (float)(j) - backproj->avEE_events * backproj->avEE_events); FewTimeBins0: Gmean[eeuser%MAX_LC_EBIN] = backproj->avEE_events; for (ntime=0; ntimeNumTimeBins; ntime++) { if (Tslices[ntime].dontUse) continue; Tslices[ntime].Gevents[eeuser%MAX_LC_EBIN] = Tslices[ntime].Events; } sprintf(logger->tmptxt, "E%1d, NumTimeBinsB: %6d, meanEvents: %6.1lf, rms: %6.2lf, bins: %4d", eeuser, logger->NumTimeBins, backproj->avEE_events, rmsFlux0, j); logger->logstat = logprint(logger, logger->tmptxt, i5); /* QQQQQQ Begin loop over basic sources to extract light curves. Ends aprox 2360 */ for (allsrc=backproj->backmodels; allsrcbackmodels; for (i=0; i<128; i++) chanmax[i] = -1; for (i=0; iNumTimeBins; i++) Tslices[i].dontUse &= MMM; /* mask away flags from prev. source */ if (eeuser == 3) PIFmap[allsrc].BURST_NUM = 0; burst_num = oldBURST_NUM = PIFmap[allsrc].BURST_NUM; if (logger->trace) traces(func_local, 100000+allsrc*1000 + eeuser*100+ PIFmap[allsrc].BURST_NUM , logger); iisum =0.0; for (i=0; isoft_gr[i] > 0.0) iisum += backproj->pif_b_pickmap[allsrc][i]; ALLsrc = allsrc - backproj->backmodels; sprintf(logger->tmptxt, "##4343A## initPIF_LC E%2d, allsrc: %2d, allsrc: %2d, ALLsrc: %2d, RA/dec: %7.3f %7.3f, %20s %s", eeuser, allsrc, allsrc, ALLsrc, GOODSRC[ALLsrc].RA, GOODSRC[ALLsrc].dec, GOODSRC[ALLsrc].name, GOODSRC[ALLsrc].ISDCname); logger->logstat = logprint(logger, logger->tmptxt, i5); /* transfer PIF-maps to structure PIFmap */ status = initPIF_LC(allsrc, allsrc, &(solut_e[eeuser][0]), &(PIFmap[allsrc]), backproj, sh, logger, status); if (status != ISDC_OK) continue; /* goto next source */ allEvents = allEEvents = allSrcEvents = allBkgEvents = allWSrcEvents = allWBkgEvents = allWEpos = allWEzero = 0.0; /* Extract events for current source & energy from event list and place in timeslices */ if (logger->trace) traces(func_local, 7772, logger); status = ex2trEvTslc(allsrc,eeuser,&(Tslices[0]),&(bufevents[0]),&(PIFmap[0]),sh,backproj,USERBIN,logger,status); /* QQQQQQ */ /* here starts the burst search in the light curves of the found sources */ if (logger->trace) traces(func_local, 200001, logger); m0 = m1 = m2 = 0.0; for (ntime=0; ntimeNumTimeBins; ntime++) { if (Tslices[ntime].dontUse) continue; flux = Tslices[ntime].wflux[allsrc]; m0 += 1.0; m1 += flux; m2 += flux * flux; } if (logger->trace) traces(func_local, 3000000+(int)(1000.0*m0), logger); if (logger->trace) traces(func_local, 3000000+(int)(1000.0*m1), logger); if (logger->trace) traces(func_local, 3000000+(int)(m2), logger); if ((m0 > 0.0) && (m1 > 0.0)) { meanFlux = m1 / m0; if ((m2 / m0 - meanFlux * meanFlux) > 0.0) rmsFlux = sqrt(m2 / m0 - meanFlux * meanFlux); else rmsFlux = 1.0; if (rmsFlux == 0.0) {rmsFlux = m1; goto FewTimeBins1;} } else goto FewTimeBins1; if (logger->trace) traces(func_local, 200002, logger); sprintf(logger->tmptxt, "E%1d, src%1d m0: %lf, meanFlux: %f, rms: %f", eeuser, allsrc, m0, meanFlux, rmsFlux); logger->logstat = logprint(logger, logger->tmptxt, i5); if ((logger->NumTimeBins < 8) || (eeuser > MAX_LC_EBIN)) { FEWTIMEBINS = 1; goto FewTimeBins1; } m0 = m1 = m2 = 0.0; for (ntime=0; ntimeNumTimeBins; ntime++) { if (Tslices[ntime].dontUse) continue; allEvents += Tslices[ntime].Events; allEEvents += Tslices[ntime].wEvents; allWSrcEvents += Tslices[ntime].srcWCounts; allSrcEvents += Tslices[ntime].srcCounts; allWBkgEvents += Tslices[ntime].bkgWCounts; allBkgEvents += Tslices[ntime].bkgCounts; if (Tslices[ntime].av_WEff > 0.0) allWEpos++; else allWEzero++; flux = Tslices[ntime].wflux[allsrc]; if (fabs(flux - meanFlux) > 3.0*rmsFlux) continue; m0 += 1.0; m1 += flux; m2 += flux * flux; } if (logger->trace) traces(func_local, 200003, logger); if ((m0 > 0.0) && (m1 > 0.0)) { meanFlux = m1 / m0; if ((m2 / m0 - meanFlux * meanFlux) > 0.0) rmsFlux = sqrt(m2 / m0 - meanFlux * meanFlux); else rmsFlux = 1.0; } else goto FewTimeBins1; if (logger->trace) traces(func_local, 200004, logger); sprintf(logger->tmptxt, "E%1d, src%02d allev: %7.1f, allEEv: %7.1f, allWSrcEv: %7.1f, allWBkgEv: %7.1f, av_WEff>0/=0: %3d, %3d, meanFlux: %f, rms: %f", eeuser, allsrc, allEvents, allEEvents, allWSrcEvents, allWBkgEvents, allWEpos, allWEzero, meanFlux, rmsFlux); logger->logstat = logprint(logger, logger->tmptxt, i5); evRA = GOODSRC[allsrc-backproj->backmodels].RA; evdec = GOODSRC[allsrc-backproj->backmodels].dec; PIFmap[allsrc].burstRA[24] = evRA; PIFmap[allsrc].burstdec[24] = evdec; sprintf(PIFmap[allsrc].name, "%s", GOODSRC[ALLsrc].name); threshold[0] = threshold[2]; threshold[1] = threshold[3]; for (ntime=0; ntimeNumTimeBins-0; ntime++) Tslices[ntime].threshold = threshold[0]; more_bursts_1: /* more_bursts_1 */ fluxmax = 0.0; chanmax[burst_num] = -1; j = k = 0; for (ntime=0; ntimeNumTimeBins-0; ntime++) { if (Tslices[ntime].dontUse) { j += (Tslices[ntime].dontUse & 1); k += (Tslices[ntime].dontUse & 2); continue; } if (Tslices[ntime].wflux[allsrc] > fluxmax) { fluxmax = Tslices[ntime].wflux[allsrc]; chanmax[burst_num] = ntime; } } if (chanmax[burst_num] < 0) goto FewTimeBins1; if ((chanmax[burst_num] < chanlim) || (abs(logger->NumTimeBins - chanmax[burst_num]) < chanlim)) { Tslices[chanmax[burst_num]].dontUse |= 2; goto more_bursts_1; } if (backproj->avEE_events > 0.0) { RMSfact = sqrt(Tslices[chanmax[burst_num]].Events / backproj->avEE_events); if ((RMSfact > 0.0) && (rmsFlux > 0.0)) signif = (fluxmax - meanFlux) / (rmsFlux * RMSfact); else signif = 0.0; } else goto FewTimeBins1; if (logger->trace) traces(func_local, 200005, logger); /* QQQQQQ if threshold exceeded call TslBurstChk (for the basic source set) */ if ((signif > Tslices[chanmax[burst_num]].threshold) && (logger->BUsrch > 0)) { sprintf(logger->tmptxt, "FBXE%1d S%02d bst_num %2d chmax %4d flxmax %6.2lf mean %6.2lf signi %4.1lf NoUse %2d %4d RMSfact %5.1lf", eeuser, allsrc, burst_num, chanmax[burst_num], fluxmax, meanFlux, signif, j, k/2, RMSfact); status = TslBurstChk(chanmax[burst_num], meanFlux, rmsFlux, threshold, allsrc, eeuser, evRA, evdec, &(Tslices[0]), &(PIFmap[0]), &(B_list[0]), backproj, logger, status); } ChanMax = chanmax[burst_num]; NdontUse = 0; for (i=0; iNumTimeBins; i++) if (Tslices[i].dontUse > 0) NdontUse++; if (logger->trace) traces(func_local, 1000000+1000*allsrc+100*eeuser+PIFmap[allsrc].BURST_NUM, logger); if (logger->trace) traces(func_local, 2000000+1000*allsrc+100*eeuser+burst_num, logger); if (PIFmap[allsrc].BURST_NUM > burst_num) { burst_num = PIFmap[allsrc].BURST_NUM; if (logger->trace) traces(func_local, 3000000+1000*allsrc+100*eeuser+PIFmap[allsrc].BURST_NUM, logger); BlistNum = logger->b_listnum - 1; sprintf(logger->tmptxt, "Burst found (found src): %2d, RA/dec: %8.4f %8.4f, name: %s !%s!, BlNum: %d, NdntUse: %d, chnmax %4d %4d", allsrc, evRA, evdec, GOODSRC[ALLsrc].name, GOODSRC[ALLsrc].ISDCname, BlistNum, NdontUse, ChanMax, chanmax[burst_num]); logger->logstat = logprint(logger, logger->tmptxt, i5); B_list[BlistNum].RA = evRA; B_list[BlistNum].dec = evdec; status = getxy(evRA, evdec, jmx_id->jmx_unit, &xmm, &ymm, jmx_id, point, logger); B_list[BlistNum].xvalstart = xmm / backproj->pix2mm + backproj->FOV_radius; B_list[BlistNum].yvalstart = ymm / backproj->pix2mm + backproj->FOV_radius; /* save data from burst in 'found source' for 'regfil2' */ sprintf(B_list[BlistNum].regsrctxt, "fk5;point(%8.4f, %8.4f) # point=%s %d text={} color=%s", evRA, evdec, SHAPE[1], (int)(20+GOODSRC[ALLsrc].bnum*2), COLOR[eeuser]); GOODSRC[ALLsrc].reg_used = 0xFFFF; B_list[BlistNum].chanmax = chanmax[burst_num] = ChanMax; sprintf(B_list[BlistNum].regtxtlin, "time chans: %3d %3d, chmax: %4d", B_list[BlistNum].start_chan, B_list[BlistNum].end_chan, chanmax[burst_num]); sprintf(regnam2, "J%1d_%04d%04d%04d_E%02d.reg", jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, eeuser); cmx = GOODSRC[ALLsrc].xx - backproj->FOV_radius; cmy = GOODSRC[ALLsrc].yy - backproj->FOV_radius; rr = sqrt(cmx*cmx + cmy*cmy); status = getxy(evRA, evdec, jmx_id->jmx_unit, &xmm, &ymm, jmx_id, point, logger); if ((fabs(xmm - cmx * backproj->pix2mm) > 0.1) || (fabs(ymm - cmy * backproj->pix2mm) > 0.1)) { sprintf(logger->tmptxt, "Inconsistent RA/dec B"); logger->logstat = logprint(logger, logger->tmptxt, 0); } GOODSRC[ALLsrc].LC_signif = B_list[BlistNum].SLCsignif[eeuser%MAX_LC_EBIN]; if ((logger->develop & 1) == 1) { if ((regfil2 = fopen(regnam2, "at")) != NULL) { fprintf(regfil2, "%s\n", B_list[BlistNum].regsrctxt); fprintf(regfil2, "#fk5y IJDstart/stop %8.2lf %8.2lf, Obstime: %7.1f, cm2: %5.1f, R: %6.2f, signif: %5.1f, %s, XY: %5.1f %5.1f\n", logger->gtiStartIJD, logger->gtiStopIJD, sh->accumT, GOODSRC[ALLsrc].src_pixel_area[eeuser], rr, B_list[BlistNum].SLCsignif[eeuser%MAX_LC_EBIN], B_list[BlistNum].regtxtlin, GOODSRC[ALLsrc].xx, GOODSRC[ALLsrc].yy); fclose(regfil2); } } GOODSRC[ALLsrc].burstptr[GOODSRC[ALLsrc].bnum] = BlistNum; GOODSRC[ALLsrc].bnum++; if (PIFmap[allsrc].BURST_NUM < 24) goto more_bursts_1; } else chanmax[burst_num] = 0; FewTimeBins1: /* FewTimeBins1 */ /* output gnuplot lightcurves for basic source set */ txtConvert(GOODSRC[ALLsrc].name, SourceName, logger); if (logger->trace) traces(func_local, 200006, logger); if (logger->trace) traces(func_local, 200000+eeuser, logger); sprintf(gnuname1, "J%1d_lc_%04d%04d%04d_%s_E%02d.gnu", jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, SourceName, eeuser); sprintf(logger->tmptxt, "GNUNAVN: ALLsrc: %2d N_SRC2: %2d %20s J%1d_lc_%04d%04d%04d_%s_E%02d.gnu", ALLsrc, logger->N_SRC2, GOODSRC[ALLsrc].name, jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, SourceName, eeuser); logger->logstat = logprint(logger, logger->tmptxt, i5); /* save light curve for later use */ if ((gnufil1 = fopen(gnuname1, "wt")) != NULL) { ilc = allsrc; if (eeuser == backproj->ee_basic) { logger->N_SRC3 = logger->N_SRC2; /* N_SRC3 is used to control output of fits-lightcurves */ logger->allsrc_lc[allsrc-backproj->backmodels] = allsrc - backproj->backmodels;; } /* fact1 = GOODSRC[ALLsrc].src_pixel_area[eeuser]; fact2 = GOODSRC[ALLsrc].src_backgr_area[eeuser]; for (ntime=0; ntimeNumTimeBins; ntime++) { if (Tslices[ntime].dontUse) continue; Tslices[ntime].wflux[allsrc] *= fact1; Tslices[ntime].wfErr[allsrc] *= fact1; Tslices[ntime].wbflux[allsrc] *= fact2; Tslices[ntime].wbfErr[allsrc] *= fact2; } rmsFlux *= GOODSRC[ALLsrc].src_pixel_area[eeuser]; NLNLNL 20160415 */ status = gnuTslice(allsrc, eeuser, rmsFlux, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, Tslices, chanmax, logger, gnufil1, status); fclose(gnufil1); gnufil1 = NULL; /* rmsFlux /= (GOODSRC[ALLsrc].src_pixel_area[eeuser] + 0.00001); NLNLNL 20160415 */ } if (logger->trace) traces(func_local, 200007, logger); if (PIFmap[allsrc].BURST_NUM > oldBURST_NUM) { if (logger->trace) traces(func_local, 200008, logger); old_b_listnum = logger->b_listnum; } } /* QQQQQ end of loop over basic source list for light curve extraction (begins app. line 2110) */ NoTimeBins1: ; /* QQQQQQ */ /* ************************************************************** */ /* */ /* LOOP OVER ADDITIONAL SOURCES */ /* */ /* ************************************************************** */ if (logger->trace) traces(func_local, 200009, logger); flag1src = 0; for (psrc=0; psrcpreset_src; psrc++) if ((backproj->preset_src_flag[psrc] & 3) == 1) flag1src++; for (psrc=0; psrcpreset_src; psrc++) { /* beginning of loop over preset sources ('psrc' loop variable) - ends in 2945 */ for (i=0; i<128; i++) chanmax[i] = -1; for (i=0; iNumTimeBins; i++) Tslices[i].dontUse &= MMM; /* mask away flags from prev. source */ if ((backproj->preset_src_flag[psrc] & 3) != 1) continue; /* this source should not be testet in this loop */ if ((backproj->preset_src_flag[psrc] & 8) == 8) continue; /* this source should not be testet in this loop */ if (backproj->preset_src_N_SRC2[psrc] < 0) { addnew = 1; logger->new_detmap = 1; backproj->preset_src_N_SRC2[psrc] = logger->N_SRC2; logger->allsrc_ptr[logger->N_SRC2] = logger->N_SRC2; sprintf(logger->tmptxt, "##4545## n_common2: %1d, psrc: %2d, X/Y: %7.2f %7.2f, backproN_SRC2: %2d, loggerN_SRC2: %2d, name: %s", n_common2, psrc, backproj->preset_src_xx[psrc], backproj->preset_src_yy[psrc], backproj->preset_src_N_SRC2[psrc], logger->N_SRC2, backproj->preset_src_name); logger->logstat = logprint(logger, logger->tmptxt, i5); GOODSRC[backproj->preset_src_N_SRC2[psrc]].xx = backproj->preset_src_xx[psrc]; GOODSRC[backproj->preset_src_N_SRC2[psrc]].yy = backproj->preset_src_yy[psrc]; GOODSRC[backproj->preset_src_N_SRC2[psrc]].RA = backproj->preset_src_RA[psrc]; GOODSRC[backproj->preset_src_N_SRC2[psrc]].dec = backproj->preset_src_dec[psrc]; GOODSRC[backproj->preset_src_N_SRC2[psrc]].catxx = backproj->preset_src_xx[psrc]; GOODSRC[backproj->preset_src_N_SRC2[psrc]].catyy = backproj->preset_src_yy[psrc]; evRA = GOODSRC[backproj->preset_src_N_SRC2[psrc]].catRA = backproj->preset_src_RA[psrc]; evdec = GOODSRC[backproj->preset_src_N_SRC2[psrc]].catdec = backproj->preset_src_dec[psrc]; GOODSRC[backproj->preset_src_N_SRC2[psrc]].src_id = backproj->preset_src_id[psrc]; sxy = backproj->preset_src_xx[psrc] * backproj->sky_ydim + backproj->preset_src_yy[psrc]; GOODSRC[backproj->preset_src_N_SRC2[psrc]].sxy = sxy; GOODSRC[backproj->preset_src_N_SRC2[psrc]].preset |= 1; sprintf(GOODSRC[backproj->preset_src_N_SRC2[psrc]].name, "%s", backproj->preset_src_name[psrc]); sprintf(logger->tmptxt, "Preset source. SCW ID#: %02d, RA/dec: %7.3f %7.3f, X/Y/R: %5.1f %5.1f %5.1f, name: %s", logger->N_SRC2 + backproj->backmodels, backproj->preset_src_RA[psrc], backproj->preset_src_dec[psrc], backproj->preset_src_xx[psrc], backproj->preset_src_yy[psrc], backproj->preset_src_rr[psrc], backproj->preset_src_name[psrc]); logger->logstat = logprint(logger, logger->tmptxt, i5); } else { addnew = 0; logger->new_detmap = 0; } evRA = GOODSRC[backproj->preset_src_N_SRC2[psrc]].catRA = backproj->preset_src_RA[psrc]; evdec = GOODSRC[backproj->preset_src_N_SRC2[psrc]].catdec = backproj->preset_src_dec[psrc]; backproj->srclist[n_common2] = backproj->preset_src_N_SRC2[psrc]; backproj->additional = backproj->preset_src_N_SRC2[psrc]; status = final_fit( eeuser, n_common2, solut_e, jmx_id, logger, backproj, sh, SCWnow, GOODSRC, USERBIN, chatter, status); if (logger->trace) traces(func_local, 470000+status, logger); if (logger->trace) traces(func_local, 470000+eeuser, logger); if (addnew) { logger->N_SRC2++; if (logger->trace) traces(func_local, 471000+logger->N_SRC2, logger); initGOOD(eeuser, GOODSRC, logger, status); addnew = 0; } if (logger->trace) traces(func_local, 623000+logger->N_SRC2, logger); if (logger->trace) traces(func_local, 624000+n_common2, logger); for (jsrc=0; jsrcbackmodels+n_common2+1; jsrc++) { sprintf(logger->tmptxt, "##3217b## E: %2d, solut_e: %f", eeuser, solut_e[eeuser][jsrc]); logger->logstat = logprint(logger, logger->tmptxt, i5); } if (logger->trace) traces(func_local, 623100+eeuser, logger); if (NOTIMEBINS) goto NoTimeBins4; /* ************************************************************** */ /* THIS MAY BE THE RIGHT PLACE FOR CALCULATION OF LIGHT CURVES */ /* FOR THE USER DEFINED SOURCE SET */ /* ************************************************************** */ /* QQQQQQ */ presrc = backproj->backmodels + backproj->preset_src_N_SRC2[psrc]; lstsrc = backproj->backmodels + n_common2; for (i=0; i<128; i++) chanmax[i] = -1; if (eeuser == 3) PIFmap[presrc].BURST_NUM = 0; burst_num = oldBURST_NUM = PIFmap[presrc].BURST_NUM; if (logger->trace) traces(func_local, 623120+eeuser, logger); if ((logger->LC == 0) || (logger->numEvents < logger->NumTimeBins)) goto NoTimeBins2; if (logger->trace) traces(func_local, 3700000+eeuser*100+ psrc, logger); if (logger->trace) traces(func_local, 3700000+eeuser*100+ presrc, logger); if (logger->trace) traces(func_local, 3700000+eeuser*100+ lstsrc, logger); if (firstAddSource == 0) { firstAddSource = 1; if (logger->trace) traces(func_local, 623400+eeuser, logger); status = source_weight2(backproj, logger, sh, PIFmap[presrc].wgt_map, status); for (i=0; i 1.0) PIFmap[presrc].wgt_map[i] = 1.0; if (logger->trace) traces(func_local, 623500+eeuser, logger); } else { for (i=0; i 1.0) PIFmap[presrc].wgt_map[i] = 1.0; if (PIFmap[presrc-1].wgt_map[i] < 0.0) PIFmap[presrc].wgt_map[i] = 0.0; PIFmap[presrc].wgt_map[i] = PIFmap[presrc-1].wgt_map[i]; } } if (logger->trace) traces(func_local, 7720000+presrc, logger); if( status != ISDC_OK ) { if (logger->trace) traces(func_local, 7720099, logger); RILlogMessage( NULL, Warning_0, "##273## source_weight2 returned: %d, status reset!", status); status = ISDC_OK; } if (logger->trace) traces(func_local, 7730000+logger->NumTimeBins, logger); sprintf(logger->tmptxt, "##4343B## initPIF_LC E%2d, lstsrc %2d, presrc: %2d, psrc: %2d, RA/dec: %7.3f %7.3f, %20s", eeuser, lstsrc, presrc, psrc, backproj->preset_src_RA[psrc], backproj->preset_src_dec[psrc], backproj->preset_src_name[psrc]); logger->logstat = logprint(logger, logger->tmptxt, i5); /* transfer PIF-maps to structure PIFmap */ status = initPIF_LC(lstsrc, presrc, &(solut_e[eeuser][0]), &(PIFmap[presrc]), backproj, sh, logger, status); if (status != ISDC_OK) { if (logger->trace) traces(func_local, 7100000+presrc, logger); status = ISDC_OK; continue; /* goto next user-defined source */ } if (logger->trace) traces(func_local, 7700000+presrc, logger); if (PIFmap[presrc].src_pixel_area < 2.0) { if (logger->trace) traces(func_local, 7100000+(int)(1000.0*PIFmap[presrc].src_pixel_area)+presrc, logger); continue; /* do not derive light curves for very marginal sources */ } if (logger->trace) traces(func_local, 7710000+presrc, logger); /* * Renormalization removed 20160219 NLNLNL PIFmap[presrc].wsrc_pixel_area *= 0.01; PIFmap[presrc].wbkg_pixel_area *= 0.01; PIFmap[presrc].src_pixel_area *= 0.01; PIFmap[presrc].bkg_pixel_area *= 0.01; */ PIFmap[presrc].burstRA[24] = backproj->preset_src_RA[psrc]; PIFmap[presrc].burstdec[24] = backproj->preset_src_dec[psrc]; sprintf(PIFmap[presrc].name, "%s", backproj->preset_src_name[psrc]); /* if (firstAddSource == 0) { firstAddSource = 1; if (logger->trace) traces(func_local, 623400+eeuser, logger); status = source_weight2(backproj, logger, sh, PIFmap[presrc].wgt_map, status); if (logger->trace) traces(func_local, 623500+eeuser, logger); } else { for (i=0; itrace) traces(func_local, 7720000+presrc, logger); if( status != ISDC_OK ) { if (logger->trace) traces(func_local, 7720099, logger); RILlogMessage( NULL, Warning_0, "##273## source_weight2 returned: %d, status reset!", status); status = ISDC_OK; } if (logger->trace) traces(func_local, 7730000+logger->NumTimeBins, logger); */ /* Extract events for current source & energy from event list and place in timeslices */ /* Calculate fluxes, corrected for exposure, deadtime, greyfilter and efficiency */ for (ntime=0; ntimeNumTimeBins; ntime++) { Tslices[ntime].srcWCounts = 0.0; Tslices[ntime].srcCounts = 0.0; Tslices[ntime].bkgWCounts = 0.0; Tslices[ntime].bkgCounts = 0.0; } if (logger->trace) traces(func_local, 7773, logger); status = ex2trEvTslc(presrc,eeuser,&(Tslices[0]),&(bufevents[0]),&(PIFmap[0]),sh,backproj,USERBIN,logger,status); if( status != ISDC_OK ) { RILlogMessage( NULL, Warning_0, "##273## ex2trEvTslc returned: %d, status reset!", status); status = ISDC_OK; } if ((logger->NumTimeBins < 10) || (eeuser > MAX_LC_EBIN)) goto gnu_out_2; if (logger->trace) traces(func_local, 7740000+presrc, logger); m0 = m1 = m2 = 0.0; for (ntime=0; ntimeNumTimeBins; ntime++) { if (Tslices[ntime].dontUse) continue; flux = Tslices[ntime].wflux[presrc]; m0 += 1.0; if (!jirisnan(flux)) { m1 += flux; m2 += flux * flux; } else {if (logger->trace) traces(func_local, 7730000+(int)(m0), logger);} if (logger->trace) traces(func_local, 7730000+(int)(m0), logger); if (logger->trace) traces(func_local, (int)(fabs(m1*100.0))+2, logger); } if (logger->trace) traces(func_local, 7750000+presrc, logger); if (logger->trace) traces(func_local, 7750000+(int)(m0), logger); if (logger->trace) traces(func_local, (int)(fabs(m1))+2, logger); /* Calculate mean flux, rmsFlux, and find maxflux and maxchan */ if (m0 > 0.0) { /* if ((m1 <= 0.0) && (PIFmap[presrc].BURST_NUM > 0))) */ if (logger->trace) traces(func_local, 7000000+(int)(100.0*m1), logger); meanFlux = m1 / m0; if ((m2 / m0 - meanFlux * meanFlux) > 0.0) rmsFlux = sqrt(m2 / m0 - meanFlux * meanFlux); else rmsFlux = 1.0; } else continue; if (logger->trace) traces(func_local, 7760000+eeuser, logger); m0 = m1 = m2 = 0.0; for (ntime=0; ntimeNumTimeBins; ntime++) { if (Tslices[ntime].dontUse) continue; allEEvents += Tslices[ntime].wEvents; allWSrcEvents += Tslices[ntime].srcWCounts; allSrcEvents += Tslices[ntime].srcCounts; allWBkgEvents += Tslices[ntime].bkgWCounts; allBkgEvents += Tslices[ntime].bkgCounts; flux = Tslices[ntime].wflux[presrc]; if (fabs(flux-meanFlux) > 3.0*rmsFlux) continue; m0 += 1.0; m1 += flux; m2 += flux * flux; } if (logger->trace) traces(func_local, 7770000+presrc, logger); /* Calculate mean flux, rmsFlux, and find maxflux and maxchan */ if (m0 > 5) { meanFlux = m1 / m0; if ((m2 / m0 - meanFlux * meanFlux) > 0.0) rmsFlux = sqrt(m2 / m0 - meanFlux * meanFlux); else rmsFlux = 1.0; if (logger->trace) traces(func_local, 7780000+(int)(100.0*m1), logger); } else { if (logger->trace) traces(func_local, 7790000+presrc, logger); continue; } if (logger->trace) traces(func_local, 25000+presrc, logger); sprintf(logger->tmptxt, "ZBE%1d, src%02d allev: %7.1f, EEv: %7.1f, SrcEv: %7.1f, BkgEv: %7.1f, WSrcEv: %7.1f, WBkgEv: %7.1f", eeuser, presrc, allEvents, allEEvents, allSrcEvents, allBkgEvents, allWSrcEvents, allWBkgEvents); logger->logstat = logprint(logger, logger->tmptxt, i5); more_bursts_2: threshold[0] = threshold[4]; threshold[1] = threshold[5]; for (ntime=0; ntimeNumTimeBins-0; ntime++) Tslices[ntime].threshold = threshold[0]; fluxmax = -10000.0; chanmax[burst_num] = -1; j = k = 0; for (ntime=0; ntimeNumTimeBins-0; ntime++) { if (Tslices[ntime].dontUse) { j += (Tslices[ntime].dontUse & 1); k += (Tslices[ntime].dontUse & 2); continue; } flux = Tslices[ntime].wflux[presrc]; if (flux > fluxmax) { fluxmax = flux; chanmax[burst_num] = ntime; } } if (chanmax[burst_num] < 0) goto gnu_out_2; if ((chanmax[burst_num] < chanlim) || (abs(logger->NumTimeBins - chanmax[burst_num]) < chanlim)) { Tslices[chanmax[burst_num]].dontUse |= 2; goto more_bursts_2; } if (logger->trace) traces(func_local, 8111000+burst_num, logger); if (logger->trace) traces(func_local, 8112000+chanmax[burst_num], logger); if (Tslices[chanmax[burst_num]].srcWCounts < 6.0) { if (logger->trace) traces(func_local, 8120000+100*presrc+(int)(Tslices[chanmax[burst_num]].srcWCounts), logger); goto gnu_out_2; /* do not accept burst triggers based on very few counts */ } if (backproj->avEE_events > 0.0) { RMSfact = sqrt(Tslices[chanmax[burst_num]].Events / backproj->avEE_events); if ((RMSfact > 0.0) && (rmsFlux > 0.0)) signif = (fluxmax - meanFlux) / (rmsFlux * RMSfact); else signif = 0.0; } else signif = 0.0; /* if (signif > Tslices[chanmax[burst_num]].threshold) { */ if ((signif > Tslices[chanmax[burst_num]].threshold) && (logger->BUsrch > 0)) { /* QQQQQQ if threshold exceeded call TslBurstChk (for the user source set) */ sprintf(logger->tmptxt, "UBXE%1d, S%02d, brst_nm: %2d, chnmx: %4d, flxmx: %6.2lf, mean: %6.2lf, signf: %4.1lf, dntUs bns: %2d %4d, RMSfact: %5.1lf", eeuser, presrc, burst_num, chanmax[burst_num], fluxmax, meanFlux, signif, j, k/2, RMSfact); logger->logstat = logprint(logger, logger->tmptxt, i5); Bburst_num = burst_num; status = TslBurstChk(chanmax[burst_num], meanFlux, rmsFlux, threshold, presrc, eeuser, evRA, evdec, &(Tslices[0]), &(PIFmap[0]), &(B_list[0]), backproj, logger, status); } else goto gnu_out_2; if( status != ISDC_OK ) { RILlogMessage( NULL, Warning_0, "##273## TslBurstChk returned: %d, source %2d skipped!", status, presrc); if (logger->trace) traces(func_local, 8140000+100*presrc+Tslices[chanmax[burst_num]].srcWCounts, logger); goto gnu_out_2; } NdontUse = 0; for (i=0; iNumTimeBins; i++) if (Tslices[i].dontUse > 0) NdontUse++; if (logger->trace) traces(func_local, 1000*presrc+NdontUse, logger); if (logger->trace) traces(func_local, 3000000+1000*presrc+100*eeuser+PIFmap[presrc].BURST_NUM, logger); if (logger->trace) traces(func_local, 7000000+1000*psrc+100*eeuser+PIFmap[presrc].BURST_NUM, logger); reg_mask = 1< burst_num) { burst_num = PIFmap[presrc].BURST_NUM; if (logger->trace) traces(func_local, 4000000+1000*presrc+100*eeuser+burst_num, logger); ALLsrc = backproj->preset_src_N_SRC2[psrc]; BlistNum = logger->b_listnum - 1; sprintf(GOODSRC[ALLsrc].name, "%s", backproj->preset_src_name[psrc]); if ((GOODSRC[ALLsrc].reg_used & reg_mask) == 0) sprintf(regsrcnam, "%s", GOODSRC[ALLsrc].name); else sprintf(regsrcnam, ""); B_list[BlistNum].RA = evRA; B_list[BlistNum].dec = evdec; B_list[BlistNum].chanmax = chanmax[Bburst_num]; status = getxy(evRA, evdec, jmx_id->jmx_unit, &xmm, &ymm, jmx_id, point, logger); B_list[BlistNum].xvalstart = xmm / backproj->pix2mm + backproj->FOV_radius; B_list[BlistNum].yvalstart = ymm / backproj->pix2mm + backproj->FOV_radius; sprintf(logger->tmptxt, "Burst found (user src): %2d, RA/dec: %8.4f %8.4f, name: %s !%s!, BlNum: %4d, NdontUse: %2d, chanmax: %4d", presrc, evRA, evdec, GOODSRC[ALLsrc].name, GOODSRC[ALLsrc].ISDCname, BlistNum, NdontUse, chanmax[Bburst_num]); logger->logstat = logprint(logger, logger->tmptxt, i5); /* save data from burst in 'user source' for 'regfil2' */ sprintf(B_list[BlistNum].regsrctxt, "fk5;point(%8.4f, %8.4f) # point=%s %d text={%s} color=%s\n", evRA, evdec, SHAPE[1], (int)(20+GOODSRC[ALLsrc].bnum*2), regsrcnam, COLOR[eeuser]); sprintf(B_list[BlistNum].regtxtlin, "time chans: %3d %3d, chmax: %4d", B_list[BlistNum].start_chan, B_list[BlistNum].end_chan, B_list[BlistNum].chanmax); GOODSRC[ALLsrc].burstptr[GOODSRC[ALLsrc].bnum] = BlistNum; GOODSRC[ALLsrc].bnum++; /* GOODSRC[ALLsrc].src_pixel_area[eeuser] = 0.01 * PIFmap[presrc].src_pixel_area;*/ /* NL 20150602 */ GOODSRC[ALLsrc].src_pixel_area[eeuser] = 0.01 * PIFmap[presrc].wsrc_pixel_area; /* NL 20160424 */ sprintf(regnam2, "J%1d_%04d%04d%04d_E%02d.reg", jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, eeuser); if ((GOODSRC[ALLsrc].reg_used & reg_mask) == 0) sprintf(sourcename, backproj->preset_src_name[psrc]); else sprintf(sourcename, ""); GOODSRC[ALLsrc].reg_used |= reg_mask; cmx = GOODSRC[ALLsrc].xx - backproj->FOV_radius; cmy = GOODSRC[ALLsrc].yy - backproj->FOV_radius; rr = sqrt(cmx*cmx + cmy*cmy); status = getxy(evRA, evdec, jmx_id->jmx_unit, &xmm, &ymm, jmx_id, point, logger); if ((fabs(xmm - cmx * backproj->pix2mm) > 0.1) || (fabs(ymm - cmy * backproj->pix2mm) > 0.1)) { sprintf(logger->tmptxt, "Inconsistent RA/dec C"); logger->logstat = logprint(logger, logger->tmptxt, 0); } GOODSRC[ALLsrc].LC_signif = B_list[BlistNum].SLCsignif[eeuser%MAX_LC_EBIN]; if ((logger->develop & 1) == 1) { if ((regfil2 = fopen(regnam2, "at")) != NULL) { fprintf(regfil2, "fk5;point(%8.4f, %8.4f) # point=%s %d text={%s} color=%s\n", evRA, evdec, SHAPE[1], (int)(20+GOODSRC[ALLsrc].bnum*2), sourcename, COLOR[eeuser]); fprintf(regfil2, "#fk5z IJDstart/stop %8.2lf %8.2lf, Obstime: %7.1f, cm2: %5.1f, R: %6.2f, signif: %5.1f, %s, XY: %5.1f %5.1f\n", logger->gtiStartIJD, logger->gtiStopIJD, sh->accumT, GOODSRC[ALLsrc].src_pixel_area[eeuser], rr, B_list[BlistNum].SLCsignif[eeuser%MAX_LC_EBIN], B_list[BlistNum].regtxtlin, GOODSRC[ALLsrc].xx, GOODSRC[ALLsrc].yy); fclose(regfil2); } } burstTime = (PIFmap[presrc].burstTstart[burst_num] + PIFmap[presrc].burstTstop[burst_num]) / 2.0; I_RA = (int)(evRA * 100.0); I_dec = (int)(fabs(evdec * 100.0)); MPdec = 'p'; if (evdec < 0.0) MPdec = 'm'; sprintf(logger->tmptxt, "BF_E%1d, src%1d, B_FLAG: %1d, oldBF: %1d, b_Tstart/stop: %lf %lf, b_Ev: %6.1f, burstTime: %lf", eeuser, presrc, burst_num, oldBURST_NUM, PIFmap[presrc].burstTstart[burst_num], PIFmap[presrc].burstTstop[burst_num], PIFmap[presrc].burstEvents[burst_num], burstTime); logger->logstat = logprint(logger, logger->tmptxt, i5); } else chanmax[burst_num] = 0; if (PIFmap[presrc].BURST_NUM < 24) goto more_bursts_2; gnu_out_2: ; /* gnu_out_2 */ if (logger->trace) traces(func_local, 5000000, logger); if ((PIFmap[presrc].BURST_NUM > oldBURST_NUM) || (logger->USER1_lc == 1) || (PIFmap[presrc].BURST_NUM > 0)) { if (logger->trace) traces(func_local, 5000001, logger); if (PIFmap[presrc].BURST_NUM > 0) { GOODSRC[ALLsrc].preset | 512; if (GOODSRC[ALLsrc].sn[eeuser%MAX_LC_EBIN] < B_list[BlistNum].SLCsignif[eeuser%MAX_LC_EBIN]) GOODSRC[ALLsrc].sn[eeuser%MAX_LC_EBIN] = B_list[BlistNum].SLCsignif[eeuser%MAX_LC_EBIN]; } txtConvert(backproj->preset_src_name[psrc], SourceName, logger); sprintf(gnuname1, "J%1d_lc_%04d%04d%04d_%s_E%02d.gnu", jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, SourceName, eeuser); sprintf(logger->tmptxt, "GNUNAVN: presrc: %2d ALLsrc: %2d, N_SRC2: %2d %20s J%1d_lc_%04d%04d%04d_%s_E%02d.gnu, Apix: %6.1f", presrc, ALLsrc, logger->N_SRC2, backproj->preset_src_name[psrc], jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, SourceName, eeuser, PIFmap[presrc].src_pixel_area); logger->logstat = logprint(logger, logger->tmptxt, i5); if ((gnufil1 = fopen(gnuname1, "wt")) != NULL) { ilc = presrc; if ((PIFmap[presrc].BURST_NUM > oldBURST_NUM) && (eeuser == backproj->ee_basic)) { logger->N_SRC3++; /* N_SRC3 controls output of fits-lc */ logger->allsrc_lc[logger->N_SRC3-1] = backproj->preset_src_N_SRC2[psrc]; GOODSRC[ALLsrc].bnum++; } /* fact1 = GOODSRC[ALLsrc].src_pixel_area[eeuser]; fact2 = GOODSRC[ALLsrc].src_backgr_area[eeuser]; for (ntime=0; ntimeNumTimeBins; ntime++) { if (Tslices[ntime].dontUse) continue; Tslices[ntime].wflux[allsrc] *= fact1; Tslices[ntime].wfErr[allsrc] *= fact1; Tslices[ntime].wbflux[allsrc] *= fact2; Tslices[ntime].wbfErr[allsrc] *= fact2; } rmsFlux *= GOODSRC[ALLsrc].src_pixel_area[eeuser]; NLNLNL 20160415 */ status = gnuTslice(presrc, eeuser, rmsFlux, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, Tslices, chanmax, logger, gnufil1, status); fclose(gnufil1); gnufil1 = NULL; /* rmsFlux /= (GOODSRC[ALLsrc].src_pixel_area[eeuser] + 0.00001); NLNLNL 20160415 */ } } else if (logger->trace) traces(func_local, 5000007, logger); if (PIFmap[presrc].BURST_NUM > oldBURST_NUM) { if (logger->trace) traces(func_local, 5000002, logger); /* save light curve for later use */ BlistNum = logger->b_listnum - 1; for (i=old_b_listnum; ib_listnum; i++) strcpy(B_list[i].regtxtlin, textline); old_b_listnum = logger->b_listnum; } else { /* if ((backproj->preset_src_flag[psrc] & 5) == 5) { */ if ((backproj->preset_src_flag[psrc] & 3) == 3) { if (logger->trace) traces(func_local, 5000003, logger); txtConvert(backproj->preset_src_name[psrc], SourceName, logger); sprintf(gnuname1, "J%1d_lc_%04d%04d%04d_%s_E%02d.gnu", jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, SourceName, eeuser); if ((gnufil1 = fopen(gnuname1, "wt")) != NULL) { ilc = presrc; /* rmsFlux *= GOODSRC[ALLsrc].src_pixel_area[eeuser]; NLNLNL 20160415 */ status = gnuTslice(presrc, eeuser, rmsFlux, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, Tslices, chanmax, logger, gnufil1, status); fclose(gnufil1); gnufil1 = NULL; /* rmsFlux /= (GOODSRC[ALLsrc].src_pixel_area[eeuser] + 0.00001); NLNLNL 20160415 */ } goto NoTimeBins4; } } /* QQQQQQ */ NoTimeBins2: ; } /* end of loop over preset sources ('psrc' loop variable) - starts in 2430 */ /* *************************************************************************** */ /* *************************************************************************** */ /* * * */ /* * End of source fits including one extra preset source * */ /* * in addition to the basic set for every step of the loop * */ /* * (begins approximately in line 2380 - about 470 lines) * */ /* * * */ /* *************************************************************************** */ /* *************************************************************************** */ sprintf(logger->tmptxt, "##98765## After introduction of preset sources, N_SRC2: %2d, N_SRC3: %2d, eeuser: %2d", logger->N_SRC2, logger->N_SRC3, eeuser); logger->logstat = logprint(logger, logger->tmptxt, i5); for (i=0; iN_SRC2; i++) { sprintf(logger->tmptxt, "##98765## i: %2d, RA/dec: %7.3f %7.3f, name: %30s, ptr: %2d", i, GOODSRC[i].RA, GOODSRC[i].dec, GOODSRC[i].name, logger->allsrc_ptr[i]); logger->logstat = logprint(logger, logger->tmptxt, i5); } if (logger->trace) traces(func_local, 6000000, logger); if ((logger->LC == 0) || (logger->NumTimeBins < 1) || (eeuser > MAX_LC_EBIN)) goto NoTimeBins4; if (logger->trace) traces(func_local, 6000001, logger); /* generate global light curve */ /* calculate weight map for global light curve (in PIFmap[0].wgt_map) */ if (eeuser == 3) PIFmap[0].BURST_NUM = 0; burst_num = oldBURST_NUM = PIFmap[0].BURST_NUM; lstsrc = backproj->backmodels + n_common2; for (i=0; i<128; i++) chanmax[i] = -1; for (i=0; iNumTimeBins; i++) Tslices[i].dontUse &= MMM; /* mask away flags from prev. source */ status = source_weight2(backproj, logger, sh, PIFmap[0].wgt_map, status); for (i=0; i 1.0) PIFmap[0].wgt_map[i] = 1.0; PIFmap[0].wsrc_pixel_area = 0.0; PIFmap[0].src_pixel_area = 0.0; for (i=0; isoft_gr[i] <= 0.0) continue; PIFmap[0].wsrc_pixel_area += PIFmap[0].wgt_map[i]; PIFmap[0].src_pixel_area += 1.0; } allEvents = allEEvents = allSrcEvents = allBkgEvents = allWSrcEvents = allWBkgEvents = allWEpos = allWEzero = 0.0; for (ntime=0; ntimeNumTimeBins; ntime++) { /* Tslices[ntime].dontUse &= 0xFFFD; */ /* Remove old "used-bin" indicators from 'dontUse' parameters */ if (Tslices[ntime].dontUse) continue; allEvents += Tslices[ntime].Events; allEEvents += Tslices[ntime].wEvents; } more_bursts_3: /* more_bursts_3 */ /* Find next peak channel in light curve */ if (logger->trace) traces(func_local, 6000002, logger); threshold[0] = threshold[6]; threshold[1] = threshold[7]; for (ntime=0; ntimeNumTimeBins-0; ntime++) Tslices[ntime].threshold = threshold[0]; fluxmax = 0.0; chanmax[burst_num] = -1; j = k = 0; for (ntime=0; ntimeNumTimeBins-0; ntime++) { Tslices[ntime].wflux[1] = 0; Tslices[ntime].wflux[0] = 0; if (Tslices[ntime].dontUse) { j += (Tslices[ntime].dontUse & 1); k += (Tslices[ntime].dontUse & 2); if ((Tslices[ntime].dontUse&1) == 0) Tslices[ntime].wflux[1] = Tslices[ntime].Events; continue; } Tslices[ntime].wflux[1] = Tslices[ntime].Events; Tslices[ntime].wflux[0] = Tslices[ntime].Events; if ((ntime < 0) || (ntime >= logger->NumTimeBins-0)) continue; if (Tslices[ntime].Events > fluxmax) { fluxmax = Tslices[ntime].Events; chanmax[burst_num] = ntime; } } if (chanmax[burst_num] < 0) goto gnu_out_3; if ((chanmax[burst_num] < chanlim) || (abs(logger->NumTimeBins - chanmax[burst_num]) < chanlim)) { Tslices[chanmax[burst_num]].dontUse |= 2; goto more_bursts_3; } if (logger->trace) traces(func_local, 6000003, logger); kscale = logger->NumTimeBins / 511; kscale2 = (int)(511/logger->NumTimeBins); if (kscale2 < 1) kscale2 = 1; kscale++; av_expo = 0.0; k = 0; for (i=0; iNumTimeBins; i+=kscale) { for (j=0; j= logger->NumTimeBins) continue; if (Tslices[i+j].dontUse&1 == 1) continue; av_expo += Tslices[i+j].final_expo; k++; } } if (k == 0) goto gnu_out_3; av_expo /= logger->TIMEstep * (float)(k); if (logger->trace) traces(func_local, 6000004, logger); if (logger->trace) traces(func_local, 1000000+(int)(av_expo*100.0), logger); meanFlux = backproj->avEE_events; if (meanFlux <= 0.0) goto no_time_bins_3; if (logger->trace) traces(func_local, 6000005, logger); if ((chanmax[burst_num] >= 0) && (meanFlux > 0.0)) { RMSfact = sqrt(Tslices[chanmax[burst_num]].Events / meanFlux); } else RMSfact = 0.0; if (RMSfact == 0.0) {RMSfact = 1.0; goto gnu_out_3;} if (rmsFlux0 > 0.0) { signif = (fluxmax - meanFlux) / (rmsFlux0 * RMSfact); } else goto gnu_out_3; if (logger->trace) traces(func_local, 6000006, logger); sprintf(logger->tmptxt, "GBXE%1d, burst: %2d, chmax: %4d, flxmax: %6.2lf, mean: %6.2lf, signif: %4.1lf, dontUse: %2d %4d, RMSfact: %5.1lf, av_expo: %5.3lf", eeuser, burst_num, chanmax[burst_num], fluxmax, meanFlux, signif, j, k/2, RMSfact, av_expo); logger->logstat = logprint(logger, logger->tmptxt, i5); if (logger->trace) traces(func_local, 6000007, logger); /* if (signif > 1.2 * Tslices[chanmax[burst_num]].threshold) { */ /* if (signif > 1.2 * threshold[0]) { */ if ((signif > 1.2 * threshold[0]) && (logger->BUsrch > 0)) { if (logger->trace) traces(func_local, 6000008, logger); status = TslABurstChk(chanmax[burst_num], meanFlux, rmsFlux0, threshold, eeuser, &(Tslices[0]), &(PIFmap[0]), &(B_list[0]), backproj, logger, status); if (logger->trace) traces(func_local, 6000009, logger); NdontUse = 0; for (i=0; iNumTimeBins; i++) if (Tslices[i].dontUse > 0) NdontUse++; if (logger->trace) traces(func_local, 1000+NdontUse, logger); if (logger->trace) traces(func_local, 3300000+100*PIFmap[0].BURST_NUM+burst_num, logger); if (PIFmap[0].BURST_NUM > burst_num) burst_num = PIFmap[0].BURST_NUM; else chanmax[burst_num] = 0; if (PIFmap[0].BURST_NUM < 24) goto more_bursts_3; } gnu_out_3: /* gnu_out_3 */ if (logger->trace) traces(func_local, 6000010, logger); if (rmsFlux0 > 0.0) { if (logger->trace) traces(func_local, 6000011, logger); kscale = logger->NumTimeBins / 511; kscale2 = (int)(511/logger->NumTimeBins); if (kscale2 < 1) kscale2 = 1; kscale++; imin = (511 - logger->NumTimeBins / kscale) / 2; if (logger->NumTimeBins < 256) imin = (511 - logger->NumTimeBins * kscale2 - kscale2 / 2) / 2; if (kscale2 > 2) stepplot = kscale2; else stepplot = 0; sprintf(regnam2, "J%1d_%04d%04d%04d_E%02d.reg", jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, eeuser); sprintf(gnuname1, "J%1d_lc_%04d%04d%04d_Global_E%02d.gnu", jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, eeuser); if ((gnufil1 = fopen(gnuname1, "wt")) != NULL) { ilc = 0; status = gnuTslice(1, eeuser, rmsFlux0, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, Tslices, chanmax, logger, gnufil1, status); fclose(gnufil1); gnufil1 = NULL; } av_expo = 0.0; k = 0; for (i=0; iNumTimeBins; i+=kscale) { for (j=0; j= logger->NumTimeBins) continue; if (Tslices[i+j].dontUse&1 == 1) continue; av_expo += Tslices[i+j].final_expo; k++; } } if (k == 0) goto gnu_out_7; av_expo /= logger->TIMEstep * (float)(k); if ((logger->develop & 1) == 1) { if ((regfil2 = fopen(regnam2, "at")) != NULL) { /* NLNLNL fprintf(regfil2, "global color=red font=\"helvetica 8 normal\" \n"); */ j = eeuser; fprintf(regfil2, "IMAGE;text( 85.0, -10.0) # text={Energy band: %4.1f to %4.1f keV}\n", jmx_id->escale[sh->pival[j][0]], jmx_id->escale[sh->pival[j][1]]); plsigma = 20.0; /* set 1-sigma y-scale for light curve plot */ if (logger->trace) traces(func_local, 6000012, logger); fprintf(regfil2, "IMAGE;polygon %d %d ", imin-1, (int)(2.0 * plsigma +0.5)); kscale = logger->NumTimeBins / 511; kscale2 = (int)(511/logger->NumTimeBins); if (kscale2 < 1) kscale2 = 1; kscale++; imin = (511 - logger->NumTimeBins / kscale) / 2; if (logger->NumTimeBins < 256) imin = (511 - logger->NumTimeBins * kscale2 - kscale2 / 2) / 2; if (kscale2 > 2) stepplot = kscale2; else stepplot = 0; for (i=0; iNumTimeBins; i+=kscale) { plotv = 0.0; for (j=0; j= logger->NumTimeBins) continue; if (Tslices[i+j].dontUse&1 == 1) continue; if (Tslices[i+j].final_expo == 0.0) continue; plotv += Tslices[i+j].wflux[1] * logger->TIMEstep / Tslices[i+j].final_expo - meanFlux / av_expo; Tslices[i+j].GlobalLC[eeuser%MAX_LC_EBIN] = Tslices[i+j].wflux[1] * logger->TIMEstep / Tslices[i+j].final_expo - meanFlux / av_expo; } plotv /= (sqrt((double)(kscale)) * rmsFlux0); plotv *= plsigma; plotv += 0.5 + 2.0 * plsigma; if (plotv > 511.0) plotv = 511.0; if (plotv < 0.0) plotv = 0.0; fprintf(regfil2, "%d %d ", (int)(i/kscale)*kscale2+imin, (int)(plotv)); if (stepplot) fprintf(regfil2, "%d %d ", (int)(i/kscale)*kscale2+imin+stepplot, (int)(plotv)); } m0 = m1 = m2 = 0.0; for (i=0; iNumTimeBins; i++) { if (Tslices[i].dontUse&1 == 1) continue; m0 += 1.0; m1 += Tslices[i].GlobalLC[eeuser%MAX_LC_EBIN]; m2 += Tslices[i].GlobalLC[eeuser%MAX_LC_EBIN] * Tslices[i].GlobalLC[eeuser%MAX_LC_EBIN]; } if (m0 > 1.0) { m1 /= m0; if ((m2 - m1*m1) > 0.0) m2 = sqrt(m2 - m1 * m1); else m2 = 1.0; } else {m1 = 0.0; m2 = 1.0;} for (i=0; iNumTimeBins; i++) { /* normalize GlobalLC for correlation calculations */ if (Tslices[i].dontUse&1 == 1) continue; Tslices[i].GlobalLC[eeuser%MAX_LC_EBIN] = (Tslices[i].GlobalLC[eeuser%MAX_LC_EBIN] - m1) / m2; } fprintf(regfil2, "%d %d\n", (int)((float)(logger->NumTimeBins*kscale2) / (float)(kscale)+imin)+1, (int)(2.0 * plsigma +0.5)); fprintf(regfil2, "IMAGE;line(%d,%d,%d,%d) # dashlist = 5 10 dash = 1\n", imin, (int)(threshold[1]*plsigma), (int)((float)(logger->NumTimeBins*kscale2) / (float)(kscale))+imin+stepplot, (int)(threshold[1]*plsigma)); fprintf(regfil2, "IMAGE;text(%d,%d) # text={0} color = black\n", imin-15, (int)(2.0 * plsigma) - 10); fprintf(regfil2, "IMAGE;text(%d,%d) # text={Channel Width: %4.1f s} color = black\n", (int)(263), (int)(2.0 * plsigma) - 35, (float)(kscale)*logger->TIMEstep); fprintf(regfil2, "IMAGE;text(%d,%d) # text={%d} color = black\n", (int)(511 - imin + 20), (int)(2.0 * plsigma) - 10, logger->NumTimeBins); fprintf(regfil2, "IMAGE;text(%d,%d) # text={Threshold} color = black\n", 475, (int)(threshold[1]*plsigma) - 10); fclose(regfil2); regfil2 = NULL; } } } /* end of Global light curve section. Starts in 2860 */ gnu_out_7: ; if (logger->trace) traces(func_local, 6000020, logger); /* save burst data for 'regfil1' */ if (PIFmap[0].BURST_NUM > oldBURST_NUM) { /* i = burst_number[0]; i %= 5; j = (eeuser - 3) % 3; sprintf(textline, "IMAGE;point(%8.4lf, %8.4lf) # point=%s %d text={} color=%s", 256.0, 256.0, SHAPE[0], 100+i+3*j, COLOR[eeuser]); burst_number[0]++; for (i=old_b_listnum; ib_listnum; i++) strcpy(B_list[i].regsrctxt, textline); */ old_b_listnum = logger->b_listnum; } no_time_bins_3: ; /* QQQQQQ */ /* ************************************************************** */ /* */ /* GENERATE LIGHT CURVES FOR DETECTOR EDGE FIELDS */ /* */ /* ************************************************************** */ /* ************************************************************** */ /* */ /* Initialize PIFmaps for edge transient search */ /* */ /* ************************************************************** */ if (logger->trace) traces(func_local, 44880000, logger); if (eeuser == 3) { edgeXmm[0] = edgeYmm[0] = edgeYmm[1] = edgeYmm[5] = edgeXmm[3] = edgeXmm[7] = 0.0; edgeXmm[1] = 237.0; edgeYmm[7] = 233.0; edgeXmm[5] = -240.0; edgeYmm[3] = -238.0; edgeXmm[2] = edgeXmm[8] = 225.0 / sqrt(2.0); edgeXmm[6] = -223.0 / sqrt(2.0); edgeYmm[6] = 223.0 /sqrt(2.0); edgeYmm[8] = 215.0 / sqrt(2.0); edgeYmm[2] = edgeYmm[4] = edgeXmm[4] = -225.0 / sqrt(2.0); } status = source_weight2(backproj, logger, sh, PIFmap[1].wgt_map, status); for (i=0; i 1.0) PIFmap[1].wgt_map[i] = 1.0; for (ej=1; ej<9; ej++) { /* beginning of edge field loop 'ej' 1-9. Ends in line 3538 */ if (logger->trace) traces(func_local, 44880000+ej, logger); edgsrc = 89 + ej; for (i=0; iNumTimeBins; i++) Tslices[i].dontUse &= MMM; /* mask away flags from prev. source */ burst_num = oldBURST_NUM = PIFmap[edgsrc].BURST_NUM; for (i=0; isoft_gr[i] <= 0.0) continue; PIFmap[edgsrc].wsrc_pixel_area += PIFmap[edgsrc].wgt_map[i]; PIFmap[edgsrc].src_pixel_area += 1.0; } edgeI = logger->useHexpos; logger->useHexpos = 100; status = FUNedge(edgsrc, eeuser, edgeXmm[ej], edgeYmm[ej], jmx_id->jmx_unit, last, jmx_id, logger, backproj, chatter, status); logger->useHexpos = edgeI; lstsrc = backproj->backmodels + n_common2; if (logger->trace) traces(func_local, 7300000+100*edgsrc+eeuser, logger); status = ISDC_OK; status = j_pif_limits2(lstsrc, edgsrc, backproj, logger, status ); if (logger->trace) traces(func_local, 7310000, logger); if (status < 0) continue; if (logger->trace) traces(func_local, 7310000+100*lstsrc, logger); if (logger->trace) traces(func_local, 7310000+100*lstsrc+n_common2, logger); for (i=0; isoft_gr[i] <= 0.0) continue; backproj->pif_s_pickmap[lstsrc][i] = backproj->pif_s_pickmap[edgsrc][i]; backproj->pif_b_pickmap[lstsrc][i] = backproj->pif_b_pickmap[edgsrc][i]; backproj->pif_maps[lstsrc][i] = backproj->pif_maps[edgsrc][i]; } if (logger->trace) traces(func_local, 7330000, logger); backproj->pif_spix[lstsrc] = backproj->pif_spix[edgsrc]; backproj->pif_bpix[lstsrc] = backproj->pif_bpix[edgsrc]; solut_e[eeuser][lstsrc] = 0.0; if (logger->trace) traces(func_local, 7340000, logger); for (ntime=0; ntimeNumTimeBins; ntime++) { Tslices[ntime].srcWCounts = 0.0; Tslices[ntime].srcCounts = 0.0; Tslices[ntime].bkgWCounts = 0.0; Tslices[ntime].bkgCounts = 0.0; } sprintf(logger->tmptxt, "##4343C## initPIF_LC E%2d, lstsrc %2d, edgsrc: %2d, ej: %2d, XRA/dc: %7.3f %7.3f, Edge%02d", eeuser, lstsrc, edgsrc, ej, evRA, evdec, edgsrc); logger->logstat = logprint(logger, logger->tmptxt, i5); status = initPIF_LC(lstsrc, edgsrc, &(solut_e[eeuser][0]), &(PIFmap[edgsrc]), backproj, sh, logger, status); if (status < 0) continue; /* Extract events for current source & energy from event list and place in timeslices */ /* Calculate fluxes, corrected for exposure, deadtime, greyfilter and efficiency */ if (logger->trace) traces(func_local, 7774, logger); status = ex2trEvTslc(edgsrc,eeuser,&(Tslices[0]),&(bufevents[0]),&(PIFmap[0]),sh,backproj,USERBIN,logger,status); if( status != ISDC_OK ) { RILlogMessage( NULL, Warning_0, "##273## ex2trEvTslc returned: %d, status reset!", status); status = ISDC_OK; } for (ntime=0; ntimeNumTimeBins; ntime++) { /* Tslices[ntime].dontUse &= 0xFFFD; */ /* Remove old "used-bin" indicators from 'dontUse' parameters */ } m0 = m1 = m2 = 0.0; for (ntime=0; ntimeNumTimeBins; ntime++) { if (Tslices[ntime].dontUse) continue; flux = Tslices[ntime].wflux[edgsrc]; m0 += 1.0; m1 += flux; m2 += flux * flux; } if (logger->trace) traces(func_local, 7750000+edgsrc, logger); /* Calculate mean flux, rmsFlux, and find maxflux and maxchan */ if ((m0 > 0.0) && (m1 > 0.0)) { meanFlux = m1 / m0; if ((m2 / m0 - meanFlux * meanFlux) > 0.0) rmsFlux = sqrt(m2 / m0 - meanFlux * meanFlux); else rmsFlux = 1.0; } else continue; if (logger->trace) traces(func_local, 7760000+edgsrc, logger); allEEvents = allSrcEvents = allBkgEvents = allWSrcEvents = allWBkgEvents = allWEpos = allWEzero = 0.0; m0 = m1 = m2 = 0.0; for (ntime=0; ntimeNumTimeBins; ntime++) { if (Tslices[ntime].dontUse) continue; allEEvents += Tslices[ntime].wEvents; allWSrcEvents += Tslices[ntime].srcWCounts; allSrcEvents += Tslices[ntime].srcCounts; allWBkgEvents += Tslices[ntime].bkgWCounts; allBkgEvents += Tslices[ntime].bkgCounts; flux = Tslices[ntime].wflux[edgsrc]; if (fabs(flux-meanFlux) > 3.0*rmsFlux) continue; m0 += 1.0; m1 += flux; m2 += flux * flux; } if (logger->trace) traces(func_local, 7770000+edgsrc, logger); /* Calculate mean flux, rmsFlux, and find maxflux and maxchan */ if (m0 > 5) { meanFlux = m1 / m0; if ((m2 / m0 - meanFlux * meanFlux) > 0.0) rmsFlux = sqrt(m2 / m0 - meanFlux * meanFlux); else rmsFlux = 1.0; } else continue; /* proceed to next edge source */ if (logger->trace) traces(func_local, 25000+edgsrc, logger); sprintf(logger->tmptxt, "EBE%1d, src%2d allev %7.1f EEv %8.1f SrcEv %7.1f BkgEv %7.1f WSrcEv %7.1f WBkgEv %7.1f m012 %5.0f %6.1f %6.1f", eeuser, edgsrc, allEvents, allEEvents, allSrcEvents, allBkgEvents, allWSrcEvents, allWBkgEvents, m0, m1, m2); logger->logstat = logprint(logger, logger->tmptxt, i5); more_bursts_4: /* more_bursts_4 */ fluxmax = -10000.0; threshold[0] = threshold[8]; threshold[1] = threshold[9]; for (ntime=0; ntimeNumTimeBins-0; ntime++) Tslices[ntime].threshold = threshold[0]; chanmax[burst_num] = -1; j = k = 0; for (ntime=0; ntimeNumTimeBins-0; ntime++) { if (Tslices[ntime].dontUse) { j += (Tslices[ntime].dontUse & 1); k += (Tslices[ntime].dontUse & 2); continue; } flux = Tslices[ntime].wflux[edgsrc]; if (flux > fluxmax) { fluxmax = flux; chanmax[burst_num] = ntime; } } if (logger->trace) traces(func_local, 25100, logger); if ((chanmax[burst_num] < 1) || (chanmax[burst_num] > logger->NumTimeBins-1)) goto gnu_out_4; sprintf(logger->tmptxt, "EXXE%1d S%02d bstnum %2d chmax %4d flmax %6.2lf mean %6.2lf sgnf %4.1lf NoUse: %2d %4d RMSfct %3.1lf thr %3.1f, cntmax: %5.1f", eeuser, edgsrc, burst_num, chanmax[burst_num], fluxmax, meanFlux, signif, j, k/2, RMSfact, 1.2*Tslices[chanmax[burst_num]].threshold, Tslices[chanmax[burst_num]].srcWCounts); logger->logstat = logprint(logger, logger->tmptxt, i5); if (logger->trace) traces(func_local, 8110000+(int)(10.0*Tslices[chanmax[burst_num]].srcWCounts), logger); if (Tslices[chanmax[burst_num]].srcWCounts < 4.0) { if (logger->trace) traces(func_local, 8120000+(int)(10.0*edgsrc+Tslices[chanmax[burst_num]].srcWCounts),logger); goto gnu_out_4; /* do not accept burst triggers based on very few counts */ } if ((chanmax[burst_num] < chanlim) || (abs(logger->NumTimeBins - chanmax[burst_num]) < chanlim)) { Tslices[chanmax[burst_num]].dontUse |= 2; goto more_bursts_4; } if (backproj->avEE_events > 0.0) RMSfact = sqrt(Tslices[chanmax[burst_num]].Events / backproj->avEE_events); else continue; if ((RMSfact > 0.0) && (rmsFlux > 0.0)) signif = (fluxmax - meanFlux) / (rmsFlux * RMSfact); else continue; sprintf(logger->tmptxt, "EBXE%1d S%02d bst_num %2d chanmax %4d flxmax %6.2lf mean %6.2lf signf %4.1lf dontUsbns: %2d %4d RMSfact %3.1lf thresh %3.1f", eeuser, edgsrc, burst_num, chanmax[burst_num], fluxmax, meanFlux, signif, j, k/2, RMSfact, 1.2*Tslices[chanmax[burst_num]].threshold); logger->logstat = logprint(logger, logger->tmptxt, i5); /* if (signif > 1.2 * Tslices[chanmax[burst_num]].threshold) { */ /* if (signif > 1.2 * threshold[0]) { */ if ((signif > 1.2 * threshold[0]) && (logger->BUsrch > 0)) { status = TslBurstChk(chanmax[burst_num], meanFlux, rmsFlux, threshold, edgsrc, eeuser, evRA, evdec, &(Tslices[0]), &(PIFmap[0]), &(B_list[0]), backproj, logger, status); } else goto gnu_out_4; if( status != ISDC_OK ) { RILlogMessage( NULL, Warning_0, "##273## TslBurstChk returned: %d, source %2d skipped!", status, edgsrc); if (logger->trace) traces(func_local, 8140000+100*edgsrc+Tslices[chanmax[burst_num]].srcWCounts, logger); goto gnu_out_4; } ChanMax = chanmax[burst_num]; NdontUse = 0; for (i=0; iNumTimeBins; i++) if (Tslices[i].dontUse > 0) NdontUse++; if (logger->trace) traces(func_local, 1000*edgsrc+NdontUse, logger); if (logger->trace) traces(func_local, 5000000+1000*edgsrc+100*eeuser+PIFmap[edgsrc].BURST_NUM, logger); if (PIFmap[edgsrc].BURST_NUM > burst_num) { burst_num = PIFmap[edgsrc].BURST_NUM; chanmax[burst_num] = ChanMax; sprintf(logger->tmptxt, "Burst %2d found for edge source#: %3d, chanmax: %3d, NdontUse: %d", burst_num+1, edgsrc, chanmax[burst_num], NdontUse); logger->logstat = logprint(logger, logger->tmptxt, i5); burstTime = (PIFmap[edgsrc].burstTstart[burst_num] + PIFmap[edgsrc].burstTstop[burst_num]) / 2.0; sprintf(logger->tmptxt, "EBF_E%1d, src%1d, B_FLAG: %1d, oldBF: %1d, b_Tstart/stop: %lf %lf, b_Ev: %6.1f, burstTime: %lf", eeuser, edgsrc, burst_num, oldBURST_NUM, PIFmap[edgsrc].burstTstart[burst_num], PIFmap[edgsrc].burstTstop[burst_num], PIFmap[edgsrc].burstEvents[burst_num], burstTime); logger->logstat = logprint(logger, logger->tmptxt, i5); } if (logger->trace) traces(func_local, 4000060, logger); Tslices[chanmax[burst_num]].dontUse |= 2; chanmax[burst_num] = 0; if (PIFmap[edgsrc].BURST_NUM < 24) goto more_bursts_4; gnu_out_4: if (PIFmap[edgsrc].BURST_NUM > oldBURST_NUM) { sprintf(gnuname1, "J%1d_lc_%04d%04d%04d_Edge%02d_E%02d.gnu", jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, edgsrc, eeuser); if ((gnufil1 = fopen(gnuname1, "wt")) != NULL) { ilc = edgsrc; status = gnuTslice(edgsrc, eeuser, rmsFlux, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, Tslices, chanmax, logger, gnufil1, status); fclose(gnufil1); gnufil1 = NULL; } /* save burst data for 'regfil1' */ evX = edgeXmm[edgsrc - 89] / backproj->pix2mm + 256.0; evY = edgeYmm[edgsrc - 89] / backproj->pix2mm + 256.0; /* i = burst_number[edgsrc - 89]; i %= 5; j = (eeuser - 3) % 3; sprintf(textline, "IMAGE;point(%8.4lf, %8.4lf) # point=%s %d text={} color=%s", evX, evY, SHAPE[0], 100+i+j*3, COLOR[eeuser]); burst_number[edgsrc - 89]++; for (i=old_b_listnum; ib_listnum; i++) strcpy(B_list[i].regsrctxt, textline); */ old_b_listnum = logger->b_listnum; } } /* end of edge field loop 'ej' 1-9. Starts in line 3281 */ /* ************************************************************** */ /* */ /* END OF LOOP OVER DETECTOR EDGE FIELDS */ /* */ /* ************************************************************** */ NoTimeBins4: ; finalSKY: ; /* ************************************************************** */ /* */ /* Backproject final residual shadowgram to form sky image */ /* */ /* ************************************************************** */ if (logger->trace) traces(func_local, 14000+logger->detImagesOut, logger); if (logger->trace) traces(func_local, 15000+logger->userImagesOut, logger); if (logger->trace) traces(func_local, 16000+eeuser, logger); backproj->additional = -1; logger->detImagesOut = 0; logger->userImagesOut = 1; for (i=0; ilocalima[i] = 0.0; status = final_sky( eeuser, n_common2, peakval, theSWG, jmx_id, logger, backproj, sh, bp, SCWnow, GOODSRC, USERBIN, chatter, status ); sprintf(logger->tmptxt, "Return from final_sky B, status: %d", status); logger->logstat = logprint(logger, logger->tmptxt, 0); if (status != ISDC_OK) {status += 10; goto exittrace;} if (logger->trace) traces(func_local, 6000, logger); if (logger->trace) traces(func_local, 67000+status, logger); if (logger->LC > 0) { if ( strstr( logger->skyImagesOut, "NONE" ) != NULL) goto noImagesOut; Skymin = 1000.0; Skymax = -1000.0; m01 = m1 = m2 = 0.0; for (i=0; i 0.0) { if (Skymin > USERBIN[eeuser].signific[i]) Skymin = USERBIN[eeuser].signific[i]; if (Skymax < USERBIN[eeuser].signific[i]) Skymax = USERBIN[eeuser].signific[i]; m01 += 1.0; m1 += USERBIN[eeuser].signific[i]; m2 += USERBIN[eeuser].signific[i] * USERBIN[eeuser].signific[i]; } else { USERBIN[eeuser].signific[i] = 0.0; } } if (logger->trace) traces(func_local, 67100, logger); if (logger->trace) traces(func_local, 1000000+m01, logger); /* if (logger->trace) traces(func_local, 1000000+(int)(Skymin*1000.0), logger); if (logger->trace) traces(func_local, 1000000+(int)(Skymax*1000.0), logger); */ if (logger->trace) traces(func_local, 67110, logger); if (m01 > 100.0) { mean1 = m1 / m01; if ((m2 / m01 - mean1 * mean1) > 0.0) rms1 = sqrt(m2 / m01 - mean1 * mean1); else rms1 = 1.0; } else { mean1 = 0.0; rms1 = 1.0; } if (logger->trace) traces(func_local, 67200, logger); if (logger->trace) traces(func_local, 1000000+(int)(mean1 ), logger); if (logger->trace) traces(func_local, 1000000+(int)(rms1 ), logger); Skymin = 1000.0; Skymax = -1000.0; m0 = m1 = m2 = 0.0; for (i=0; i 0.0) { USERBIN[eeuser].signific[i] -= mean1; USERBIN[eeuser].signific[i] /= rms1; if (Skymin > USERBIN[eeuser].signific[i]) Skymin = USERBIN[eeuser].signific[i]; if (Skymax < USERBIN[eeuser].signific[i]) Skymax = USERBIN[eeuser].signific[i]; if (fabs(USERBIN[eeuser].signific[i]) > 4.0) continue; m0 += 1.0; m1 += USERBIN[eeuser].signific[i]; m2 += USERBIN[eeuser].signific[i] * USERBIN[eeuser].signific[i]; } } if (logger->trace) traces(func_local, 67300, logger); if (m0 > 0.0) { mean = m1 / m0; if ((m2 / m0 - mean * mean) > 0.0) rms = sqrt(m2 / m0 - mean * mean); else rms = 1.0; } else { if (logger->trace) traces(func_local, 3880001, logger); } if (logger->trace) traces(func_local, 67400, logger); sprintf(logger->tmptxt, "Signi_skyim E%2d, min/max: %6.1f %6.1f, m0: %5.0f loss: %5.0f, mean: %6.3f %6.3f, rms: %5.3f %5.3f", eeuser, Skymin, Skymax, m01, m01-m0, mean1, mean, rms1, rms); logger->logstat = logprint(logger, logger->tmptxt, i5); if (logger->trace) traces(func_local, 67500, logger); } if (logger->trace) traces(func_local, (int)(GOODSRC[logger->N_SRC2-1].src_pixel_area[eeuser])+2, logger); noImagesOut: status = fluxerror( eeuser, logger->N_SRC2, peakval, logger, backproj, sh, GOODSRC, status ); } /* end of eeuser loop - begins approx in line 2010 */ backproj->additional = -1; /* QQQQQQ */ /* AAAAAA */ /* ******************************************************************* */ /* */ /* FINAL SUMMARY PRINTOUT */ /* */ /* ******************************************************************* */ /* QQQQQQ */ sprintf(logger->tmptxt, "##98765## Before summary print, N_SRC2: %2d, N_SRC3: %2d", logger->N_SRC2, logger->N_SRC3); logger->logstat = logprint(logger, logger->tmptxt, i5); for (i=0; iN_SRC2; i++) { sprintf(logger->tmptxt, "##98765## i: %2d, RA/dec: %7.3f %7.3f, name: %30s, ptr: %2d", i, GOODSRC[i].RA, GOODSRC[i].dec, GOODSRC[i].name, logger->allsrc_ptr[i]); logger->logstat = logprint(logger, logger->tmptxt, i5); } sprintf(calibname, "c%1d_%04d%04d%04d.bin", jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv); for (jsrc=0; jsrcdevelop & 1) == 1) { calibout(jsrc, n_common2, backproj->ee_basic, vigncorr, GOODSRC[jsrc].diff_x, GOODSRC[jsrc].diff_y, solut_e[0][jsrc], solut_e[1][jsrc], solut_e[2][jsrc], backproj->resid_data_peak[0][jsrc], backproj->resid_data_peak[1][jsrc], backproj->resid_data_peak[2][jsrc], calibname, theSWG, point, jmx_id, logger, backproj, SCWnow, sh, GOODSRC, chatter, status); } } n_additional = 0; summary_p(n_common2, n_additional, solut_e, jmx_id, backproj, sh, SCWnow, GOODSRC, &(PIFmap[0]), &(Tslices[0]), &(B_list[0]), chanmax, logger, chatter, status ); /* ******************************************************************* */ /* */ /* GENERATE SKY-IMAGES FOR ACCEPTED BURST TIME PERIODS */ /* (ONLY FOR BURSTS WITHOUT IDENTIFIED SOURCE ORIGIN) */ /* */ /* ******************************************************************* */ sprintf(logger->tmptxt, "Return from summary_p\n"); logger->logstat = logprint(logger, logger->tmptxt, 0); if (logger->trace) traces(func_local, 300001, logger); if ((logger->LC == 0) || (logger->b_listnum == 0)) goto NoBurstImages; for (i=0; ifitsimage[0][0][i]; /* introduce 'blindspots' corresponding to the 'found' sources */ for (tstsrc=0; tstsrcsky_ydim] = 0; skyim_mask[GOODSRC[tstsrc].sxy - backproj->sky_ydim] = 0; skyim_mask[GOODSRC[tstsrc].sxy + backproj->sky_ydim + 1] = 0; skyim_mask[GOODSRC[tstsrc].sxy - backproj->sky_ydim - 1] = 0; skyim_mask[GOODSRC[tstsrc].sxy + backproj->sky_ydim + 1] = 0; skyim_mask[GOODSRC[tstsrc].sxy - backproj->sky_ydim - 1] = 0; skyim_mask[GOODSRC[tstsrc].sxy + 1] = 0; skyim_mask[GOODSRC[tstsrc].sxy - 1] = 0; } if (logger->trace) traces(func_local, 300002, logger); listlim = logger->b_listnum; for (ibu=0; ibu logger->NumTimeBins-2)) continue; if ((B_list[ibu].end_chan < 1) || (B_list[ibu].end_chan > logger->NumTimeBins-2)) continue; ee_lc = B_list[ibu].ee & logger->LCMAXMASK; m0 = m1 = 0.0; for (ntime = B_list[ibu].start_chan; ntime<=B_list[ibu].end_chan; ntime++) { if (Tslices[ntime].dontUse&1 == 1) continue; m0 += Gmean[ee_lc ]; m1 += (Tslices[ntime].Gevents[ee_lc ] - Gmean[ee_lc ]); } if (m0 > 0.0) B_list[ibu].GLCsignif[ee_lc ] = m1 / sqrt(m0); else B_list[ibu].GLCsignif[ee_lc ] = 1.0; sprintf(logger->tmptxt, "A%2d srcID%2d E%1d %4X st/endCh: %3d %3d ", ibu, B_list[ibu].sourceID, ee_lc, B_list[ibu].ee, B_list[ibu].start_chan, B_list[ibu].end_chan); if (logger->trace) traces(func_local, 320000+ibu, logger); sprintf(logger->tmptxt2, "SLC/GLCsig: %4.1f %4.1f, RA/dec: %7.3f %7.3f rms: %5.2f b_listnum: %2d", B_list[ibu].SLCsignif[ee_lc], B_list[ibu].GLCsignif[ee_lc], B_list[ibu].RA, B_list[ibu].dec, B_list[ibu].rmsFlux, logger->b_listnum); if (logger->trace) traces(func_local, 330000+ibu, logger); if ((strlen(logger->tmptxt) + strlen(logger->tmptxt2)) < MAX_STR_LEN) strcat(logger->tmptxt, logger->tmptxt2); logger->logstat = logprint(logger, logger->tmptxt, i5); } if (logger->trace) traces(func_local, 340000, logger); /* ***************************************************** */ /* */ /* LOOP OVER BURST PERIODS */ /* (loop ends in line 4182 (approx)) */ /* */ /* ***************************************************** */ N_SRC2orig = logger->N_SRC2; preset_src_orig = backproj->preset_src; psrc = backproj->preset_src; if (listlim > 20) listlim = 20; for (ibu=0; ibutrace) traces(func_local, 350000+ibu, logger); if ((B_list[ibu].ee & 0x4000) == 0) continue; /* only use ibu-values selected in 'summary_p' */ for (i=0; i<8; i++) B_list[ibu].src_pixel_area[i] = 0.0; if ((B_list[ibu].sourceID > 0) && (B_list[ibu].sourceID < 90)) { B_list[ibu].ee |= 0x8000; /* indicate that the sourceID is established based on light curves only */ continue; /* skip all identified sources */ } if (logger->trace) traces(func_local, 351000+ibu, logger); /* loop over (first 4 at most) user energy bins */ eelim = sh->numShds - backproj->ee_basic; if (eelim > MAX_LC_EBIN - backproj->ee_basic) eelim = MAX_LC_EBIN - backproj->ee_basic; bstsrc = logger->N_SRC2; old_psrc = psrc; IMA_qual = 0; for (eeuser=backproj->ee_basic; eeuseree_basic; eeuser++) { /* loop over energy bands (ends 4179) */ if (logger->trace) traces(func_local, 352000+eeuser, logger); if (eeuser == backproj->ee_basic) PIFmap[bstsrc].BURST_NUM = 0; burst_num = oldBURST_NUM = PIFmap[bstsrc].BURST_NUM; sxy_skymax = 0; sprintf(logger->tmptxt, "Imaging loop: ibu: %2d, eeuser: %1d, trig: %02X", ibu, eeuser, B_list[ibu].triggers); logger->logstat = logprint(logger, logger->tmptxt, i5); ee_mask = 1<trace) traces(func_local, 352000+ee_mask, logger); if ((B_list[ibu].triggers & ee_mask) != ee_mask) continue; /* skip over uninteresting energy bands */ if (logger->trace) traces(func_local, 352000+eeuser, logger); if (logger->trace) traces(func_local, 300004, logger); for (i=0; itrace) traces(func_local, 300005, logger); sprintf(burst_txt, "BurstImaJ%1d_%04d%04d%04d_E%1d_C%03d_%02d_", jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, eeuser, B_list[ibu].start_chan, (B_list[ibu].end_chan - B_list[ibu].start_chan)); sprintf(logger->tmptxt, "ibu: %1d, Tchannels: %3d %3d, TimeBins: %4d, E-range: %5d %5d, buffer events: %12ld", ibu, B_list[ibu].start_chan, B_list[ibu].end_chan, logger->NumTimeBins, sh->pival[eeuser][0], sh->pival[eeuser][1], logger->numEvents); logger->logstat = logprint(logger, logger->tmptxt, i5); if (logger->trace) traces(func_local, 300007, logger); status = ISDC_OK; current_start[ibu] = B_list[ibu].start_chan; current_end[ibu] = B_list[ibu].end_chan; current_ibu = ibu; /* m0 = m1 = 0.0; for (ntime = B_list[ibu].start_chan; ntime<=B_list[ibu].end_chan; ntime++) { m0 += Gmean[eeuser]; m1 += (Tslices[ntime].Gevents[eeuser] - Gmean[eeuser]); } if (m0 > 0.0) B_list[ibu].GLCsignif[eeuser] = m1 / sqrt(m0); else B_list[ibu].GLCsignif[eeuser] = 1.0; */ /* ***************************************************** */ /* */ /* FILL BURST SHADOWGRAM */ /* */ /* ***************************************************** */ status = fill_burst_shd(Tslices[B_list[ibu].start_chan].Tstart, Tslices[B_list[ibu].end_chan].Tstop, sh->pival[eeuser][0], sh->pival[eeuser][1], bufevents, logger->numEvents, auxdata, &(shd_burst), logger, chatter, status); events = 0.0; events2 = 0.0; m0 = m1 = m2 = 0.0; for (i=0; isoft_gr[i] <= 0.0) continue; events += shd_burst.shd[i]; m0 += 1.0; m1 += shd_burst.shd[i]; m2 += shd_burst.shd[i] * shd_burst.shd[i]; } if (m1 < 100.0) goto recover1; if (m0 > 0.0) { m1 /= m0; if ((m2 / m0 - m1 * m1) > 0.0) rms = sqrt(m2 / m0 - m1 * m1); else rms = 1.0; } else goto recover1; sprintf(logger->tmptxt, "fill_burst_shd returned %6.0f %6.0f %6ld events, rms: %6.4f, scw events: %7.0lf", events, events2, shd_burst.nEvs, rms, USERBIN[eeuser].totcnt); logger->logstat = logprint(logger, logger->tmptxt, i5); if (logger->trace) traces(func_local, 300008, logger); if( status != ISDC_OK ) { RILlogMessage( NULL, Warning_0, "##273## fill_burst_shd returned: %d, status reset!", status); sprintf(logger->tmptxt, "%s: fill_burst_shd return status: %d", burst_txt, status); logger->logstat = logprint(logger, logger->tmptxt, i5); status = ISDC_OK; } if (logger->trace) traces(func_local, 300009, logger); /* ***************************************************** */ /* */ /* SUBTRACT RESCALED SHADOGRAM FOR FULL SCW */ /* */ /* ***************************************************** */ m0 = m1 = m2 = mean2 = mean3 = 0.0; for (i=0; isoft_gr[i] > 0.0) { m0 += 1.0; m1 += shd_burst.shd[i]; mean2 += ubin_shd[i]; mean3 += fabs(ubin_shd[i]); m2 += shd_burst.shd[i] * shd_burst.shd[i]; } else shd_burst.shd[i] = 0.0; } if (m0 > 0.0) { mean = m1 / m0; mean2 /= m0; mean3 /= m0; if ((m2 / m0 - mean * mean) > 0.0) rms = sqrt(m2 / m0 - mean * mean); else rms = 1.0; } else goto recover1; T_total = (Tslices[logger->NumTimeBins-1].Tstop - Tslices[0].Tstart) * 86400.0; T_burst = (Tslices[current_end[ibu]].Tstop - Tslices[current_start[ibu]].Tstart) * 86400.0; sprintf(logger->tmptxt, "Before ubin normalization: T_burst: %5.1f, T_total: %6.1f, mean_ubin_shd: %6.4f, fabs: %6.1f", T_burst, T_total, mean2, mean3); logger->logstat = logprint(logger, logger->tmptxt, i5); sprintf(logger->tmptxt, "Before ubin normalization: m0: %6.1f, mean: %6.3f, rms: %6.4f", m0, mean, rms); logger->logstat = logprint(logger, logger->tmptxt, i5); m0 = m1 = m2 = 0.0; for (i=0; isoft_gr[i] > 0.0) { shd_burst.shd[i] -= (ubin_shd[i] * T_burst / T_total); m0 += 1.0; m1 += shd_burst.shd[i]; m2 += shd_burst.shd[i] * shd_burst.shd[i]; } else shd_burst.shd[i] = 0.0; } if (m0 > 0.0) { mean = m1 / m0; if ((m2 / m0 - mean * mean) > 0.0) rms = sqrt(m2 / m0 - mean * mean); else rms = 1.0; } else goto NoPointSrc1; sprintf(logger->tmptxt, "After ubin normalization: m0: %6.1f, mean: %6.3f, rms: %6.4f", m0, mean, rms); logger->logstat = logprint(logger, logger->tmptxt, i5); m0 = m1 = m2 = 0.0; for (i=0; isoft_gr[i] > 0.0) { sh->f_shadgram[eeuser][i] = shd_burst.shd[i] - mean; m0 += 1.0; m1 += sh->f_shadgram[eeuser][i]; m2 += sh->f_shadgram[eeuser][i] * sh->f_shadgram[eeuser][i]; } else sh->f_shadgram[eeuser][i] = 0.0; } if (m0 > 0.0) { mean = m1 / m0; if ((m2 / m0 - mean * mean) > 0.0) rms = sqrt(m2 / m0 - mean * mean); else rms = 1.0; } else goto NoPointSrc1; sprintf(logger->tmptxt, "After zero normalization: m0: %6.1f, mean: %6.3f, rms: %6.4f", m0, mean, rms); logger->logstat = logprint(logger, logger->tmptxt, i5); for (i=0; isoft_gr[i] > 0.0) { sh->f_shadgram[eeuser][i] = shd_burst.shd[i]; } else sh->f_shadgram[eeuser][i] = 0.0; } /* ***************************************************** */ /* */ /* IMAGE BACKPROJECTION */ /* */ /* ***************************************************** */ logger->burstimage = 1; status = j_fitprep(eeuser, theSWG, jmx_id, logger, backproj, sh, bp, SCWnow, USERBIN, chatter, status); logger->burstimage = 0; sprintf(logger->tmptxt, "work2d: after fitprep: E: %2d, totcnt: %7.1f", eeuser, sh->totcnt[eeuser]); logger->logstat = logprint(logger, logger->tmptxt, i5); /* check quality of sky-image - if too poor, do not output mini-images */ burst_ima_qual[ibu][eeuser] = 0; Skymin = 1000.0; Skymax = -1000.0; m01 = m1 = m2 = 0.0; for (i=0; i 0.0) && (USERBIN[eeuser].variance[i] > 0.0)) { burst_skyim[ibu][eeuser][i] = USERBIN[eeuser].fst_imag[i] / sqrt(USERBIN[eeuser].variance[i] + 5.0); if (Skymin > burst_skyim[ibu][eeuser][i]) Skymin = burst_skyim[ibu][eeuser][i]; if (Skymax < burst_skyim[ibu][eeuser][i]) Skymax = burst_skyim[ibu][eeuser][i]; m01 += 1.0; m1 += burst_skyim[ibu][eeuser][i]; m2 += burst_skyim[ibu][eeuser][i] * burst_skyim[ibu][eeuser][i]; } else { skyim_mask[i] = 0.0; burst_skyim[ibu][eeuser][i] = 0.0; } } if (m01 > 100.0) { mean1 = m1 / m01; if ((m2 / m01 - mean1 * mean1) > 0.0) rms1 = sqrt(m2 / m01 - mean1 * mean1); else rms1 = 1.0; } else { if (logger->trace) traces(func_local, 3890000, logger); burst_ima_qual[ibu][eeuser] = 0; } Skymin = 1000.0; Skymax = -1000.0; m0 = m1 = m2 = 0.0; for (i=0; i burst_skyim[ibu][eeuser][i]) Skymin = burst_skyim[ibu][eeuser][i]; if (Skymax < burst_skyim[ibu][eeuser][i]) { Skymax = burst_skyim[ibu][eeuser][i]; sxy_skymax = i; } if (fabs(burst_skyim[ibu][eeuser][i]) > 4.0) continue; m0 += 1.0; m1 += burst_skyim[ibu][eeuser][i]; m2 += burst_skyim[ibu][eeuser][i] * burst_skyim[ibu][eeuser][i]; } B_list[ibu].sxy_skymax[eeuser] = sxy_skymax; status = fits_test_output(burst_txt, burst_skyim[ibu][eeuser], 511, logger, chatter, status ); if (m0 > 100.0) { mean = m1 / m0; if ((m2 / m0 - mean * mean) > 0.0) rms = sqrt(m2 / m0 - mean * mean); else rms = 1.0; } else { if (logger->trace) traces(func_local, 3890001, logger); burst_ima_qual[ibu][eeuser] = 0; } xxx = (float)((int)(sxy_skymax / backproj->sky_ydim)); yyy = (float)(sxy_skymax % backproj->sky_ydim); B_list[ibu].ISsignif[eeuser] = (Skymax - mean) / rms; sprintf(logger->tmptxt, "Bst_im %2d E%1d, min/max: %4.1f %5.2f, x/y: %5.1f %5.1f, m0: %5.0f loss: %5.0f, av: %6.3f %6.3f, rms: %5.3f %5.3f, signif: %4.1f", ibu, eeuser, Skymin, Skymax, xxx, yyy, m01, m01-m0, mean1, mean, rms1, rms, B_list[ibu].ISsignif[eeuser]); logger->logstat = logprint(logger, logger->tmptxt, i5); /* determine RA and dec of the source candidates found at this energy */ xmm = (xxx - backproj->FOV_radius) * backproj->pix2mm; ymm = (yyy - backproj->FOV_radius) * backproj->pix2mm; getRAdec(xmm, ymm, jmx_id->jmx_unit, &s_RA, &s_dec, point, logger, jmx_id); sprintf(logger->tmptxt, "Bst_im2: ibu: %1d, xmm/ymm: %6.1f %6.1f, RA/dec: %7.3f %7.3f", ibu, xmm, ymm, s_RA, s_dec); logger->logstat = logprint(logger, logger->tmptxt, i5); if (B_list[ibu].ISsignif[eeuser] > NBsigniflim) { /* ends in 4188 */ /* ***************************************************** */ /* SATISFACTORY DOMINATING SOURCE FOUND */ /* (EXCESS > 4 * LOCAL MAP RMS) */ /* ***************************************************** */ IMA_qual = 1; B_list[ibu].ee = 0xE000 + eeuser; /* check candidate against 'found' and 'preset' sources */ /* (including already accepted 'new' burst candidates) */ newflag = 1; for (tstsrc=0; tstsrcN_SRC2; tstsrc++) { diff_x = xxx - GOODSRC[tstsrc].xx; diff_y = yyy - GOODSRC[tstsrc].yy; diff = diff_x * diff_x + diff_y *diff_y; diff = sqrt(diff); if (diff < 1.5) {j = tstsrc; newflag = 0;} } if (logger->trace) traces(func_local, 3910000, logger); if (newflag) if (logger->trace) traces(func_local, 3920000+logger->N_SRC2, logger); if (B_list[ibu].ISsignif[eeuser] > NBsigniflim) if (logger->trace) traces(func_local, 3930000+logger->N_SRC2, logger); if ((newflag) || (B_list[ibu].ISsignif[eeuser] > NBsigniflim)) { if (logger->trace) traces(func_local, 3940000, logger); NBsigniflim = B_list[ibu].ISsignif[eeuser]; if (logger->trace) traces(func_local, 3950000+(int)(NBsigniflim * 10.0), logger); GOODSRC[logger->N_SRC2].xx = backproj->preset_src_xx[psrc] = xxx; GOODSRC[logger->N_SRC2].yy = backproj->preset_src_yy[psrc] = yyy; xmm = (xxx - backproj->FOV_radius) * backproj->pix2mm; ymm = (yyy - backproj->FOV_radius) * backproj->pix2mm; getRAdec(xmm, ymm, jmx_id->jmx_unit, &s_RA, &s_dec, point, logger, jmx_id); GOODSRC[logger->N_SRC2].RA = backproj->preset_src_RA[psrc] = s_RA; GOODSRC[logger->N_SRC2].dec = backproj->preset_src_dec[psrc] = s_dec; i = srcidentify2(theSWG, s_RA, s_dec, backproj->preset_src_name[psrc], backproj->preset_src_Iname[psrc], &catflag, jmx_id, point, logger); backproj->preset_src_Iname[psrc][16] = 0; strcpy(GOODSRC[logger->N_SRC2].name, backproj->preset_src_name[psrc]); strcpy(GOODSRC[logger->N_SRC2].ISDCname, backproj->preset_src_Iname[psrc]); B_list[ibu].IDima[eeuser] = logger->N_SRC2; B_list[ibu].sourceID = logger->N_SRC2; B_list[ibu].ix = (int)(xxx+0.5); B_list[ibu].iy = (int)(yyy+0.5); B_list[ibu].RA = s_RA; B_list[ibu].dec = s_dec; B_list[ibu].xvalstart = xmm / backproj->pix2mm + backproj->FOV_radius; B_list[ibu].yvalstart = ymm / backproj->pix2mm + backproj->FOV_radius; strcpy(B_list[ibu].name[eeuser], backproj->preset_src_name[psrc]); B_list[ibu].ee = 0xE000+eeuser; sprintf(logger->tmptxt, "New burst src: X/Y: %5.1f %5.1f %s, RA/dec: %7.3f %7.3f, psrc: %2d, N_SRC2: %2d, %s, ibu: %d, i: %d", xxx, yyy, backproj->preset_src_name[psrc], backproj->preset_src_RA[psrc], backproj->preset_src_dec[psrc], psrc+1, logger->N_SRC2+1, GOODSRC[logger->N_SRC2].name, ibu, i); logger->logstat = logprint(logger, logger->tmptxt, i5); preset_ibu[psrc] = ibu; ibu_preset[ibu ] = psrc; ibu_N_SRC2[ibu ] = logger->N_SRC2; logger->new_detmap = 1; backproj->preset_src_N_SRC2[psrc] = logger->N_SRC2; sprintf(logger->tmptxt, "##4547## n_common2: %1d, psrc: %2d, X/Y: %7.2f %7.2f, backproN_SRC2: %2d, loggerN_SRC2: %2d, eeuser: %1d", n_common2, psrc, backproj->preset_src_xx[psrc], backproj->preset_src_yy[psrc], backproj->preset_src_N_SRC2[psrc], logger->N_SRC2, eeuser); logger->logstat = logprint(logger, logger->tmptxt, i5); GOODSRC[backproj->preset_src_N_SRC2[psrc]].xx = backproj->preset_src_xx[psrc]; GOODSRC[backproj->preset_src_N_SRC2[psrc]].yy = backproj->preset_src_yy[psrc]; GOODSRC[backproj->preset_src_N_SRC2[psrc]].RA = backproj->preset_src_RA[psrc]; GOODSRC[backproj->preset_src_N_SRC2[psrc]].dec = backproj->preset_src_dec[psrc]; GOODSRC[backproj->preset_src_N_SRC2[psrc]].catxx = backproj->preset_src_xx[psrc]; GOODSRC[backproj->preset_src_N_SRC2[psrc]].catyy = backproj->preset_src_yy[psrc]; GOODSRC[backproj->preset_src_N_SRC2[psrc]].catRA = backproj->preset_src_RA[psrc]; GOODSRC[backproj->preset_src_N_SRC2[psrc]].catdec = backproj->preset_src_dec[psrc]; GOODSRC[backproj->preset_src_N_SRC2[psrc]].src_id = 77; sxy = backproj->preset_src_xx[psrc] * backproj->sky_ydim + backproj->preset_src_yy[psrc]; GOODSRC[backproj->preset_src_N_SRC2[psrc]].sxy = sxy; sprintf(GOODSRC[backproj->preset_src_N_SRC2[psrc]].name, "%s", backproj->preset_src_name[psrc]); backproj->preset_src++; psrc++; if (newflag) { logger->N_SRC2++; initGOOD(eeuser, GOODSRC, logger, status); newflag = 0; } if (logger->trace) traces(func_local, 3960000+logger->N_SRC2, logger); } else { sprintf(logger->tmptxt, "Source refound: %32s, RA/dec: %7.3f %7.3f, N_SRC2: %2d", GOODSRC[j].name, GOODSRC[j].RA, GOODSRC[j].dec, j); } } NoPointSrc1: recover1: for (i=0; iLCMAXMASK; for (i=0; i<8; i++) { /* sunchronize all energy bands to strongest found source */ strcpy(B_list[ibu].name[i], B_list[ibu].name[ee_lc]); B_list[ibu].sxy_skymax[i] = B_list[ibu].sxy_skymax[ee_lc]; if (logger->trace) traces(func_local, 3971000+i, logger); } } /* end of loop over energies (eeuser) (begins 3789) */ ee_lc = B_list[ibu].ee & logger->LCMAXMASK; if (IMA_qual > 0) burst_ima_qual[ibu][ee_lc] = 1; } /* end of loop over burst detections (ibu) (begins 3766) */ sprintf(logger->tmptxt, "##98765## After analysis of 'global bursts' N_SRC2: %2d, N_SRC3: %2d", logger->N_SRC2, logger->N_SRC3); logger->logstat = logprint(logger, logger->tmptxt, i5); for (i=0; iN_SRC2; i++) { sprintf(logger->tmptxt, "##98765## i: %2d, RA/dec: %7.3f %7.3f, name: %30s, ptr: %2d", i, GOODSRC[i].RA, GOODSRC[i].dec, GOODSRC[i].name, logger->allsrc_ptr[i]); logger->logstat = logprint(logger, logger->tmptxt, i5); } /* SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS */ /* generate light curves for dominant source found during burst time intervals */ /* the generation logic is similar to the one used for 'user-sources' * (maybe better 'found sources', because refined pifmaps have not been generated NL 20150423 ?): * - an outer loop over the 'user energy bands' * - an inner loop over the new, dominant, sources (max one per E-band) * */ eelim = sh->numShds - backproj->ee_basic; if (eelim > MAX_LC_EBIN - backproj->ee_basic) eelim = MAX_LC_EBIN - backproj->ee_basic; logger->N_SRC2 = N_SRC2orig; if (logger->trace) traces(func_local, 470100+logger->N_SRC2, logger); if (preset_src_orig == backproj->preset_src) goto NoBurstImages; for (i=0; iNumTimeBins; i++) Tslices[i].dontUse &= 1; /* mask away peak count ch. flags from prev. sources */ for (eeuser=backproj->ee_basic; eeuseree_basic; eeuser++) { /* loop over (< 4) E-bands (end 4185) */ if (eeuser == backproj->ee_basic) PIFmap[bstsrc].BURST_NUM = 0; burst_num = oldBURST_NUM = PIFmap[bstsrc].BURST_NUM; status = source_weight2(backproj, logger, sh, PIFmap[1].wgt_map, status); for (i=0; i 1.0) PIFmap[1].wgt_map[i] = 1.0; /* loop over new sources, found in global-lc_sky-images (end 4178) */ for (psrc=preset_src_orig; psrcpreset_src; psrc++) { ibu = preset_ibu[psrc]; backproj->srclist[n_common2] = backproj->preset_src_N_SRC2[psrc]; backproj->additional = backproj->preset_src_N_SRC2[psrc]; status = final_fit( eeuser, n_common2, solut_e, jmx_id, logger, backproj, sh, SCWnow, GOODSRC, USERBIN, chatter, status); if (logger->trace) traces(func_local, 471100+status, logger); if (logger->trace) traces(func_local, 723000+logger->N_SRC2, logger); if (logger->trace) traces(func_local, 623100+eeuser, logger); bstsrc = backproj->backmodels + backproj->preset_src_N_SRC2[psrc]; lstsrc = backproj->backmodels + n_common2; for (i=0; i<128; i++) chanmax[i] = -1; if (eeuser == 3) PIFmap[bstsrc].BURST_NUM = 0; burst_num = oldBURST_NUM = PIFmap[bstsrc].BURST_NUM; if ((logger->LC == 0) || (logger->numEvents < logger->NumTimeBins)) goto NoTimeBins2; if (logger->trace) traces(func_local, 700000+eeuser*100+ bstsrc, logger); sprintf(logger->tmptxt, "##4343D## initPIF_LC E%2d, lstsrc %2d, bstsrc: %2d, psrc: %2d, RA/dec: %7.3f %7.3f, %20s", eeuser, lstsrc, edgsrc, psrc, s_RA, s_dec, backproj->preset_src_name[psrc]); logger->logstat = logprint(logger, logger->tmptxt, i5); /* transfer PIF-maps to structure PIFmap */ status = initPIF_LC(lstsrc, bstsrc, &(solut_e[eeuser][0]), &(PIFmap[bstsrc]), backproj, sh, logger, status); if (logger->trace) traces(func_local, 310017, logger); if (status < 0) goto nextbstsrc; /* Extract events for current source & energy from event list and place in timeslices */ /* Calculate fluxes, corrected for exposure, deadtime, greyfilter and efficiency */ for (i=0; iNumTimeBins; i++) Tslices[i].dontUse &= 1; /* mask away flags from prev. sources */ if (logger->trace) traces(func_local, 7775, logger); status = ex2trEvTslc(bstsrc, eeuser, &(Tslices[0]), &(bufevents[0]), &(PIFmap[0]), sh, backproj, USERBIN, logger, status); if( status != ISDC_OK ) { if (logger->trace) traces(func_local, 399012, logger); RILlogMessage( NULL, Warning_0, "##273## ex2trEvTslc returned: %d, status reset!", status); status = ISDC_OK; } for (ntime=0; ntimeNumTimeBins; ntime++) { /* sprintf(logger->tmptxt, "extrEvt: bstsrc: %3d, ntime: %3d, wflux: %7.3f %7.3f, dontUse: %d", bstsrc, ntime, Tslices[ntime].wflux[bstsrc], Tslices[ntime].wbflux[bstsrc], Tslices[ntime].dontUse); logger->logstat = logprint(logger, logger->tmptxt, i5); */ Tslices[ntime].srcWCounts = 0.0; Tslices[ntime].srcCounts = 0.0; Tslices[ntime].bkgWCounts = 0.0; Tslices[ntime].bkgCounts = 0.0; } if (logger->trace) traces(func_local, 600002, logger); /* Calculate mean flux, rmsFlux, and find maxflux and maxchan */ maxflux = 0.0; maxchan = -1; m0 = m1 = m2 = 0.0; for (ntime=0; ntimeNumTimeBins; ntime++) { if (Tslices[ntime].dontUse) continue; flux = Tslices[ntime].wflux[bstsrc]; m0 += 1.0; m1 += flux; m2 += flux * flux; if (maxflux < flux) {maxflux = flux; maxchan = ntime;} if ((ntime > 5) && (ntime < 10)) { sprintf(logger->tmptxt, "ExtrLCburstA: ntime: %3d, flux: %5.2f, m0/m1/m2: %4.0f %6.3f %7.3f", ntime, flux, m0, m1, m2); logger->logstat = logprint(logger, logger->tmptxt, i5); } } if ((m0 > 0.0) && (m1 > 0.0)) { meanFlux = m1 / m0; if ((m2 / m0 - meanFlux * meanFlux) > 0.0) rmsFlux = sqrt(m2 / m0 - meanFlux * meanFlux); else rmsFlux = 1.0; } else { meanFlux = 0.0; rmsFlux = 99.9; goto nextbstsrc; } if (logger->trace) traces(func_local, 610003, logger); allEEvents = allSrcEvents = allBkgEvents = allWSrcEvents = allWBkgEvents = allWEpos = allWEzero = 0.0; m0 = m1 = m2 = 0.0; for (ntime=0; ntimeNumTimeBins; ntime++) { if (Tslices[ntime].dontUse) continue; allEEvents += Tslices[ntime].wEvents; allWSrcEvents += Tslices[ntime].srcWCounts; allSrcEvents += Tslices[ntime].srcCounts; allWBkgEvents += Tslices[ntime].bkgWCounts; allBkgEvents += Tslices[ntime].bkgCounts; flux= (Tslices[ntime].wflux[bstsrc] - meanFlux) / rmsFlux; if (fabs(flux) > 3.0) continue; m0 += 1.0; m1 += flux; m2 += flux * flux; if ((ntime > 5) && (ntime < 10)) { sprintf(logger->tmptxt, "ExtrLCburstB: ntime: %3d, flux: %5.2f, m0/m1/m2: %4.0f %6.3f %7.3f", ntime, flux, m0, m1, m2); logger->logstat = logprint(logger, logger->tmptxt, i5); } } if (m0 > 5) { meanFlux = m1 / m0; if ((m2 / m0 - meanFlux * meanFlux) > 0.0) rmsFlux = sqrt(m2 / m0 - meanFlux * meanFlux); else rmsFlux = 1.0; } else { meanFlux = 0.0; rmsFlux = 99.9; goto nextbstsrc; } if (logger->trace) traces(func_local, 610004, logger); sprintf(logger->tmptxt, "BBE%1d, src%2d allev: %7.1f, allEEv: %7.1f, allSrcEv: %7.1f, allBkgEv: %7.1f, allWSrcEv: %7.1f, allWBkgEv: %7.1f, m012: %5.0f %6.1f %6.1f", eeuser, bstsrc, allEvents, allEEvents, allSrcEvents, allBkgEvents, allWSrcEvents, allWBkgEvents, m0, m1, m2); logger->logstat = logprint(logger, logger->tmptxt, i5); if (logger->trace) traces(func_local, 600006, logger); B_list[ibu].src_pixel_area[eeuser] = 0.01 * PIFmap->wsrc_pixel_area; /* NLNLNL 20160224 */ strcpy(B_list[ibu].name[eeuser], backproj->preset_src_name[psrc]); txtConvert(backproj->preset_src_name[psrc], SourceName, logger); sprintf(gnuname1, "J%1d_lc_%04d%04d%04d_%s_E%02d.gnu", jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, SourceName, eeuser); if ((gnufil1 = fopen(gnuname1, "wt")) != NULL) { ilc = bstsrc; if (eeuser == (B_list[ibu].ee & logger->LCMAXMASK)) logger->N_SRC3++; /* N_SRC3 controls output of fits-lc */ logger->N_SRCT = logger->N_SRC2; if (logger->N_SRCT < logger->N_SRC3) logger->N_SRCT = logger->N_SRC3; logger->allsrc_lc[logger->N_SRC3-1] = backproj->preset_src_N_SRC2[psrc]; logger->allsrc_ptr[logger->N_SRCT-1] = backproj->preset_src_N_SRC2[psrc]; GOODSRC[logger->N_SRC3-1].bnum++; /* for (ntime=1; ntimeNumTimeBins-1; ntime++) { if (Tslices[ntime].dontUse) continue; Tslices[ntime].flux *= GOODSRC[ALLsrc].src_pixel_area[eeuser]; Tslices[ntime].fl_err *= GOODSRC[ALLsrc].src_pixel_area[eeuser]; Tslices[ntime].bflux *= GOODSRC[ALLsrc].backgr_pixel_area[eeuser]; Tslices[ntime].bfl_err *= GOODSRC[ALLsrc].backgr_pixel_area[eeuser]; } rmsFlux *= GOODSRC[ALLsrc].src_pixel_area[eeuser]; */ status = gnuTslice(bstsrc, eeuser, rmsFlux, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, Tslices, chanmax, logger, gnufil1, status); fclose(gnufil1); gnufil1 = NULL; } nextbstsrc: ; } /* end of loop over new, dominant sources (psrc) begins in line 4196 (approx) */ /* ***************************************************** */ /* */ /* END OF LOOP OVER NEW, DOMINANT, SOURCES */ /* */ /* ***************************************************** */ } /* end of loop over energies (eeuser) begins in line 4185 (approx) */ /* ******************************************************** */ /* */ /* END OF EEUSER LOOP */ /* */ /* ******************************************************** */ sprintf(logger->tmptxt, "##98765## After final loop over dominant sources, N_SRC2: %2d, N_SRC3: %2d", logger->N_SRC2, logger->N_SRC3); logger->logstat = logprint(logger, logger->tmptxt, i5); for (i=0; iN_SRC2; i++) { sprintf(logger->tmptxt, "##98765## i: %2d, RA/dec: %7.3f %7.3f, name: %30s, ptr: %2d", i, GOODSRC[i].RA, GOODSRC[i].dec, GOODSRC[i].name, logger->allsrc_ptr[i]); logger->logstat = logprint(logger, logger->tmptxt, i5); } if (logger->trace) traces(func_local, 783000, logger); /* print full list of burst detections */ /* and combine identical source candidates into short list of candidates */ /* eliminate candidates whose light curves have already been searched for bursts (such cases should not exist) */ for (ibu=0; ibub_listnum; ibu++) { if (B_list[ibu].sourceID < 0) continue; ee_lc = B_list[ibu].ee & logger->LCMAXMASK; sprintf(logger->tmptxt, "B%2d srcID%2d ID%2d IDima%2d Etrig %1d %4X st/end/maxChan: %3d %3d %3d, SGsig: %4.1f %4.1f, Ix/y/sign:", ibu, B_list[ibu].sourceID, B_list[ibu].ID[ee_lc], B_list[ibu].IDima[ee_lc], ee_lc, B_list[ibu].ee, B_list[ibu].start_chan, B_list[ibu].end_chan, B_list[ibu].chanmax, B_list[ibu].SLCsignif[ee_lc], B_list[ibu].GLCsignif[ee_lc]); for (j=0; jee_basic; sprintf(logger->tmptxt2, " %3d %3d %4.1f", (int)(B_list[ibu].sxy_skymax[jj] / backproj->sky_ydim), B_list[ibu].sxy_skymax[jj] % backproj->sky_ydim, B_list[ibu].ISsignif[jj]); if ((strlen(logger->tmptxt) + strlen(logger->tmptxt2)) < MAX_STR_LEN) strcat(logger->tmptxt, logger->tmptxt2); } logger->logstat = logprint(logger, logger->tmptxt, i5); } NoBurstImages: ; /* NoBurstImages: */ status = ISDC_OK; if (logger->LC == 0) goto slut; k = 0; jlim = sh->numShds; if (jlim > MAX_LC_EBIN) jlim = MAX_LC_EBIN; for (ibu=0; ibub_listnum; ibu++) { ee_lc = B_list[ibu].ee & logger->LCMAXMASK; if (logger->trace) traces(func_local, 1000000 + logger->b_listnum*10000 + ee_lc*100 + jlim, logger); if ((B_list[ibu].ee & 0x4000) == 0) continue; /* use only primary bursts from summary_p and from imaging */ if (B_list[ibu].end_chan == 0) continue; /* use only sensible bursts */ sprintf(logger->tmptxt, "J%1d_%04d%04d%04d E%d b# %2d, Id/S/G/I Id/I:", jmx_id->jmx_unit+1 ,backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, ee_lc, ibu); for (j=backproj->ee_basic; jtrace) traces(func_local, 2000000+j, logger); if (B_list[ibu].IDima[j] < 0) sprintf(Stext, "--"); else { if (B_list[ibu].IDima[j] < 100) sprintf(Stext, "%02d", B_list[ibu].IDima[j]); else sprintf(Stext, "U%1d", B_list[ibu].IDima[j]-100); } sprintf(logger->tmptxt2, " %02d %4.1f %4.1f %4.1f %2s %4.1f |", B_list[ibu].ID[j], B_list[ibu].SLCsignif[j], B_list[ibu].GLCsignif[j], B_list[ibu].ISsignif[j], Stext, B_list[ibu].Isignif[j]); if ((strlen(logger->tmptxt) + strlen(logger->tmptxt2)) < MAX_STR_LEN) strcat(logger->tmptxt, logger->tmptxt2); if (logger->trace) traces(func_local, 3000000+j, logger); } logger->logstat = logprint(logger, logger->tmptxt, i5); ee_lc = B_list[ibu].ee & logger->LCMAXMASK; sprintf(logger->tmptxt, "b# %02d, b_listnum: %2d, triggers: %02X, %s, EE: %04X, eeuser: %d", ibu, logger->b_listnum, B_list[ibu].triggers, B_list[ibu].name, B_list[ibu].ee, ee_lc); logger->logstat = logprint(logger, logger->tmptxt, i5); if (logger->trace) traces(func_local, 1000000 + logger->b_listnum*10000 + ee_lc*100 + jlim, logger); } sprintf(logger->tmptxt, " "); logger->logstat = logprint(logger, logger->tmptxt, i5); /* ***************************************************************** */ /* */ /* FINAL OUTPUT OF SELECTED BURST LIGHT CURVES TO REGFILE */ /* */ /* ***************************************************************** */ /* 4444444444444444444444444444444444444444444444444444444444444444444444 */ k_offset = 1; for (ibu=0; ibub_listnum; ibu++) { /* loop over burst periods [ibu] (ends in 4712) */ ee_lc = B_list[ibu].ee & logger->LCMAXMASK; sprintf(logger->tmptxt, "%2d B_list.ee: %04X %1d\n", ibu, B_list[ibu].ee, ee_lc); logger->logstat = logprint(logger, logger->tmptxt, i5); if (logger->trace) traces(func_local, 9780001, logger); if ((B_list[ibu].ee&0xF000) != 0xE000) continue; if (logger->trace) traces(func_local, 9780002, logger); if ((logger->burst_lc_out == 0) && (burst_ima_qual[ibu][ee_lc] == 0)) continue; if (logger->trace) traces(func_local, 9780003, logger); srcID = B_list[ibu].sourceID; if ((srcID <= 0) || (srcID >= 90)) continue; if (logger->trace) traces(func_local, 9780004, logger); xxx = (float)((int)(B_list[ibu].sxy_skymax[ee_lc] / backproj->sky_ydim)); yyy = (float)(B_list[ibu].sxy_skymax[ee_lc] % backproj->sky_ydim); if (logger->trace) traces(func_local, 9780014, logger); sprintf(logger->tmptxt, "%2d B_list.srcID: %2d, ix.iy.name: %5.1f %5.1f %s, x/yvalstart: %5.1f %5.1f", ibu, srcID, xxx, yyy, B_list[ibu].name[ee_lc], B_list[ibu].xvalstart, B_list[ibu].yvalstart); logger->logstat = logprint(logger, logger->tmptxt, i5); if (logger->trace) traces(func_local, 9780214, logger); eelim = sh->numShds - backproj->ee_basic; if (eelim > MAX_LC_EBIN - backproj->ee_basic) eelim = MAX_LC_EBIN - backproj->ee_basic; for (eeuser=backproj->ee_basic; eeuseree_basic; eeuser++) { /* loop over eeuser (ends in 4490) */ if (burst_lc_output[srcID][eeuser] > 0) { kk_offs = 2 * burst_lc_output[srcID][eeuser]; if (kk_offs > 63) kk_offs = 63; kk_off2 = kk_offs + 2; } else { burst_lc_output[srcID][eeuser] = k_offset; kk_offs = 2 * k_offset; if (kk_offs > 63) kk_offs = 63; kk_off2 = kk_offs + 2; } if (logger->trace) traces(func_local, 9780005, logger); if (logger->trace) traces(func_local, 9780005+eeuser*10, logger); sprintf(regnam2, "J%1d_%04d%04d%04d_E%02d.reg", jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, eeuser); printf("qxq %s\n", regnam2); GOODSRC[srcID].LC_signif = B_list[ibu].ISsignif[ee_lc]; if ((logger->develop & 1) == 1) { if ((regfil2 = fopen(regnam2, "at")) != NULL) { /* regfil2 write section (ends in 4487 */ /* if (eeuser == ee_lc) { */ fprintf(regfil2, "IMAGE;line(%d,%d,%d,%d) # dashlist = 8 6 dash = 1 color=%s\n", B_list[ibu].iy+1, B_list[ibu].ix+1, (int)(B_list[ibu].start_chan*kscale2/kscale+imin), (int)(kk_off2 * plsigma), COLOR[k_offset+2]); fprintf(regfil2, "IMAGE;line(%d,%d,%d,%d) # dashlist = 8 6 dash = 1 color=%s\n", B_list[ibu].iy+1, B_list[ibu].ix+1, (int)(B_list[ibu].end_chan*kscale2/kscale+imin), (int)(kk_off2 * plsigma), COLOR[k_offset+2]); if (logger->trace) traces(func_local, 9780006, logger); /* } */ if (kk_offs < 2 * k_offset) goto end_lightcurve; if (B_list[ibu].ID[ee_lc] > 0) goto skip_fk5_print; reg_mask = 1<trace) traces(func_local, 9780007, logger); cmx = GOODSRC[srcID ].xx - backproj->FOV_radius; cmy = GOODSRC[srcID ].yy - backproj->FOV_radius; rr = sqrt(cmx*cmx + cmy*cmy); if (logger->trace) traces(func_local, 9780017, logger); status = getxy(B_list[ibu].RA, B_list[ibu].dec, jmx_id->jmx_unit, &xmm, &ymm, jmx_id, point, logger); if ((fabs(xmm - cmx * backproj->pix2mm) > 0.1) || (fabs(ymm - cmy * backproj->pix2mm) > 0.1)) { if (logger->trace) traces(func_local, 9780027, logger); sprintf(logger->tmptxt, "Inconsistent D B_list.name: %20s, GOODSRC.name: %20s", B_list[ibu].name[ee_lc], GOODSRC[srcID].name); logger->logstat = logprint(logger, logger->tmptxt, i5); } if (logger->trace) traces(func_local, 9780037, logger); /* B_list[ibu].ISsignif[ee_lc] = 99.0; GOODSRC[srcID].src_pixel_area[0] = 999.0; */ sprintf(regtext, " time chans: %3d %3d, chmax: %4d, E: %1d, ibu: %1d, %14s", B_list[ibu].start_chan, B_list[ibu].end_chan, B_list[ibu].chanmax, (B_list[ibu].ee & logger->LCMAXMASK), ibu, B_list[ibu].name[ee_lc]); if (logger->trace) traces(func_local, 9780047, logger); fprintf(regfil2, "#fk5u IJDstart/stop %8.2lf %8.2lf, Obstime: %7.1f, cm2: %5.1f, R: %6.2f, signif: %5.1f, %s, XY: %5.1f %5.1f\n", logger->gtiStartIJD, logger->gtiStopIJD, sh->accumT, GOODSRC[srcID].src_pixel_area[ee_lc], rr, B_list[ibu].ISsignif[ee_lc], regtext, GOODSRC[srcID ].xx, GOODSRC[srcID ].yy); /* logger->gtiStartIJD, logger->gtiStopIJD, sh->accumT, B_list[ibu].src_pixel_area[ee_lc], rr, B_list[ibu].ISsignif[ee_lc], regtext, GOODSRC[srcID ].xx, GOODSRC[srcID ].yy); */ if (logger->trace) traces(func_local, 9780057, logger); skip_fk5_print: /* recover lightcurves from gnu-files */ if (logger->trace) traces(func_local, 9780010+eeuser, logger); txtConvert(B_list[ibu].name[ee_lc], SourceName, logger); sprintf(gnuname2, "J%1d_lc_%04d%04d%04d_%s_E%02d.gnu", jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, SourceName, eeuser); if (logger->trace) traces(func_local, 9780020+eeuser, logger); m0 = m1 = m2 = 0.0; if ((gnufil2 = fopen(gnuname2, "rt")) != NULL) { if (logger->trace) traces(func_local, 9780030+eeuser, logger); sprintf(logger->tmptxt,"Successfully opened %s for reading", gnuname2); logger->logstat = logprint(logger, logger->tmptxt, i5); for (ntime=0; ntimeNumTimeBins; ntime++) { if ((Tslices[ntime].dontUse & 1) > 0) continue; fgets(textline, 200, gnufil2); if (feof(gnufil2)) goto udoneAN; sscanf(textline, "%d %f", &j, &(Tslices[ntime].wflux[1])); m0 += 1.0; m1 += Tslices[ntime].wflux[1]; m2 += Tslices[ntime].wflux[1] * Tslices[ntime].wflux[1]; k = strlen(textline); textline[k-1] = 0; } if (logger->trace) traces(func_local, 9780040+eeuser, logger); udoneAN: fclose(gnufil2); gnufil2 = NULL; if (m0 > 4) { mean = m1 / m0; if ((m2 / m0 - mean * mean) > 0.0) rmsFlux = sqrt(m2 / m0 - mean * mean); else rmsFlux = 1.0; } else { mean = 0.0; rmsFlux = 99.9; goto end_lightcurve; } if (logger->trace) traces(func_local, 9780050+eeuser, logger); } else { if (logger->trace) traces(func_local, 9780060+eeuser, logger); sprintf(logger->tmptxt,"Could not open %s for reading", gnuname2); logger->logstat = logprint(logger, logger->tmptxt, 0); goto end_lightcurve; } if (logger->trace) traces(func_local, 9780070+eeuser, logger); m0 = m1 = m2 = 0.0; for (ntime=0; ntimeNumTimeBins; ntime++) { if (Tslices[ntime].dontUse & 1) continue; if (fabs(Tslices[ntime].wflux[1] - mean) / rmsFlux > 3.0) continue; m0 += 1.0; m1 += Tslices[ntime].wflux[1]; m2 += Tslices[ntime].wflux[1] * Tslices[ntime].wflux[1]; } if (logger->trace) traces(func_local, 9780080+eeuser, logger); if (m0 > 4) { mean = m1 / m0; if ((m2 / m0 - mean * mean) > 0.0) rmsFlux = sqrt(m2 / m0 - mean * mean); else rmsFlux = 1.0; for (ntime=0; ntimeNumTimeBins; ntime++) { if (Tslices[ntime].dontUse & 1) continue; Tslices[ntime].wflux[1] -= mean; } } else { mean = 0.0; rmsFlux = 99.9; goto end_lightcurve; } if (logger->trace) traces(func_local, 9780090+eeuser, logger); fprintf(regfil2, "IMAGE;polygon %d %d ", imin-1, (int)(kk_off2 * plsigma + 0.5)); if (logger->trace) traces(func_local, 9780100+eeuser, logger); kscale = logger->NumTimeBins / 511; kscale2 = (int)(511/logger->NumTimeBins); if (kscale2 < 1) kscale2 = 1; kscale++; if (logger->NumTimeBins < 256) imin = (511 - logger->NumTimeBins * kscale2 - kscale2 / 2) / 2; if (kscale2 > 2) stepplot = kscale2; else stepplot = 0; for (i=0; iNumTimeBins; i+=kscale) { plotv = 0.0; for (j=0; j= logger->NumTimeBins) continue; if (Tslices[i+j].dontUse&1 == 1) continue; plotv += Tslices[i+j].wflux[1]; } plotv /= (sqrt((double)(kscale)) * rmsFlux); plotv *= plsigma; plotv += 0.5 + kk_off2 * plsigma; if (plotv > 511.0) plotv = 511.0; if (plotv < 0.0) plotv = 0.0; fprintf(regfil2, "%d %d ", (int)(i/kscale)*kscale2+imin, (int)(plotv)); if (stepplot) fprintf(regfil2, "%d %d ", (int)(i/kscale)*kscale2+imin+stepplot, (int)(plotv)); } if (logger->trace) traces(func_local, 9780110+eeuser, logger); fprintf(regfil2, "%d %d # color=%s\n", (int)((float)(logger->NumTimeBins*kscale2) / (float)(kscale))+imin+1, (int)(kk_off2 * plsigma +0.5), COLOR[k_offset+2]); status = LCcorrelation(eeuser, Tslices, logger, B_list[ibu].start_chan, B_list[ibu].end_chan, &correl, status); fprintf(regfil2, "#fk5c Correlation: %5.3f\n", correl); B_list[ibu].correl = correl; end_lightcurve: fclose(regfil2); regfil2 = NULL; if (logger->trace) traces(func_local, 9780140+eeuser, logger); } /* end of regfile2 write sequence (start in line 4316) */ } if (logger->trace) traces(func_local, 9780150+eeuser, logger); } /* end of eeuser-loop (start in line 4295) */ if (logger->trace) traces(func_local, 9780200+k_offset, logger); if (kk_offs == 2 * k_offset) k_offset++; if (logger->trace) traces(func_local, 9780200+k_offset, logger); if (logger->trace) traces(func_local, 9780200+ibu, logger); next_ibu: ; } /* end of ibu-loop (start in line 4467) */ if (logger->trace) traces(func_local, 9780300+k_offset, logger); /* 66666666666666666666666666666666666666666666666666666666666666666666666666666 */ /* QQQQQQ */ if (logger->trace) traces(func_local, 9880000+n_common2, logger); for (isrc=0; isrctrace) traces(func_local, 9881000+isrc, logger); for (eeuser=backproj->ee_basic; eeuseree_basic; eeuser++) { /* loop over eeuser (ends in 4548) */ if (src_lc_output[isrc][eeuser] > 0) { kk_offs = 2 * src_lc_output[isrc][eeuser]; if (kk_offs > 63) kk_offs = 63; kk_off2 = kk_offs + 2; } else { src_lc_output[isrc][eeuser] = k_offset; kk_offs = 2 * k_offset; if (kk_offs > 63) kk_offs = 63; kk_off2 = kk_offs + 2; } if (logger->trace) traces(func_local, 9882000+eeuser, logger); if (logger->trace) traces(func_local, 9882000+imin , logger); sprintf(regnam2, "J%1d_%04d%04d%04d_S%02d.reg", jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, eeuser); printf("qzq %s\n", regnam2); if ((logger->develop & 1) == 1) { if ((regfil2 = fopen(regnam2, "at")) != NULL) { /* regfil2 write section (ends in 4546 */ fprintf(regfil2, "IMAGE;line(%d,%d,%d,%d) # dashlist = 8 6 dash = 1 color=%s\n", (int)(GOODSRC[isrc].yy+1), (int)(GOODSRC[isrc].xx+1), (int)(imin), (int)(kk_off2 * plsigma), COLOR[k_offset+2]); fprintf(regfil2, "IMAGE;line(%d,%d,%d,%d) # dashlist = 8 6 dash = 1 color=%s\n", (int)(GOODSRC[isrc].yy+1), (int)(GOODSRC[isrc].xx+1), (int)(imin + logger->NumTimeBins*kscale2/kscale), (int)(kk_off2 * plsigma), COLOR[k_offset+2]); if (logger->trace) traces(func_local, 9880006, logger); if (kk_offs < 2 * k_offset) goto end_lightcurve2; /* recover lightcurves from gnu-files */ if (logger->trace) traces(func_local, 9880010+eeuser, logger); txtConvert(GOODSRC[isrc].name, SourceName, logger); sprintf(gnuname2, "J%1d_lc_%04d%04d%04d_%s_E%02d.gnu", jmx_id->jmx_unit+1, backproj->aux_orbit, backproj->aux_pid, sh->aux_pidv, SourceName, eeuser); if (logger->trace) traces(func_local, 9880020, logger); if (logger->NumTimeBins < 5) goto end_lightcurve2; m0 = m1 = m2 = 0.0; if ((gnufil2 = fopen(gnuname2, "rt")) != NULL) { if (logger->trace) traces(func_local, 9880030, logger); sprintf(logger->tmptxt,"Successfully opened %s for reading", gnuname2); logger->logstat = logprint(logger, logger->tmptxt, i5); for (ntime=0; ntimeNumTimeBins; ntime++) { if ((Tslices[ntime].dontUse & 1) > 0) continue; fgets(textline, 200, gnufil2); if (feof(gnufil2)) goto udoneAN2; sscanf(textline, "%d %f", &j, &(Tslices[ntime].wflux[1])); m0 += 1.0; m1 += Tslices[ntime].wflux[1]; m2 += Tslices[ntime].wflux[1] * Tslices[ntime].wflux[1]; k = strlen(textline); textline[k-1] = 0; } if (logger->trace) traces(func_local, 9880040, logger); udoneAN2: fclose(gnufil2); gnufil2 = NULL; if (m0 > 4) { mean = m1 / m0; if ((m2 / m0 - mean * mean) > 0.0) rmsFlux = sqrt(m2 / m0 - mean * mean); else rmsFlux = 1.0; if (logger->trace) traces(func_local, 9880052, logger); } else { mean = 0.0; rmsFlux = 99.9; if (logger->trace) traces(func_local, 9880053, logger); } } else { if (logger->trace) traces(func_local, 9880060, logger); sprintf(logger->tmptxt,"Could not open %s for reading", gnuname2); logger->logstat = logprint(logger, logger->tmptxt, i5); goto next_sru; } if (logger->trace) traces(func_local, 9880070, logger); m0 = m1 = m2 = 0.0; for (ntime=0; ntimeNumTimeBins; ntime++) { if (Tslices[ntime].dontUse & 1) continue; if (fabs(Tslices[ntime].wflux[1] - mean) / (rmsFlux+0.0001) > 3.0) continue; m0 += 1.0; m1 += Tslices[ntime].wflux[1]; m2 += Tslices[ntime].wflux[1] * Tslices[ntime].wflux[1]; } if (logger->trace) traces(func_local, 9880080, logger); if (m0 > 4) { mean = m1 / m0; if ((m2 / m0 - mean * mean) > 0.0) rmsFlux = sqrt(m2 / m0 - mean * mean); else rmsFlux = 1.0; for (ntime=0; ntimeNumTimeBins; ntime++) { if (Tslices[ntime].dontUse & 1) continue; Tslices[ntime].wflux[1] -= mean; } } else { mean = 0.0; rmsFlux = 99.9; goto end_lightcurve2; goto next_sru; } if (logger->trace) traces(func_local, 9880090, logger); fprintf(regfil2, "IMAGE;polygon %d %d ", imin-1, (int)(kk_off2 * plsigma + 0.5)); kscale = logger->NumTimeBins / 511; kscale2 = (int)(511/logger->NumTimeBins); if (kscale2 < 1) kscale2 = 1; kscale++; if (logger->trace) traces(func_local, 9880100, logger); if (logger->trace) traces(func_local, 9880100+kscale, logger); for (i=0; iNumTimeBins; i+=kscale) { plotv = 0.0; for (j=0; j= logger->NumTimeBins) continue; if (Tslices[i+j].dontUse&1 == 1) continue; plotv += Tslices[i+j].wflux[1]; } plotv /= (sqrt((double)(kscale)) * (rmsFlux+0.00001)); plotv *= plsigma; plotv += 0.5 + kk_off2 * plsigma; if (plotv > 511.0) plotv = 511.0; if (plotv < 0.0) plotv = 0.0; fprintf(regfil2, "%d %d ", (int)(i/kscale)*kscale2+imin, (int)(plotv)); if (stepplot) fprintf(regfil2, "%d %d ", (int)(i/kscale)*kscale2+imin+stepplot, (int)(plotv)); } if (logger->trace) traces(func_local, 9880110, logger); fprintf(regfil2, "%d %d # color=%s\n", (int)((float)(logger->NumTimeBins*kscale2) / (float)(kscale))+imin+1, (int)(kk_off2 * plsigma +0.5), COLOR[k_offset+2]); end_lightcurve2: fclose(regfil2); regfil2 = NULL; if (logger->trace) traces(func_local, 9880140, logger); } /* end of regfile2 write sequence (start in line 4529) */ } if (logger->trace) traces(func_local, 9880150, logger); } /* end of eeuser-loop (start in line 4508) */ if (logger->trace) traces(func_local, 9880200+k_offset, logger); if (kk_offs == 2 * k_offset) k_offset++; if (logger->trace) traces(func_local, 9880200+k_offset, logger); if (logger->trace) traces(func_local, 9880200+ibu, logger); next_sru: ; } /* end of isrc-loop (start in line 4504) */ if (logger->trace) traces(func_local, 9880300+k_offset, logger); /* 77777777777777777777777777777777777777777777777777777777777777777777777777777 */ slut: ; /* slut: */ if (logger->trace) traces(func_local, 100000+logger->N_SRC2, logger); if (logger->trace) traces(func_local, 99, logger); /* function exit timing start */ TT_now = clock(); logger->func_times[func_local] += difftime(TT_now, logger->TT); logger->TT = TT_now; logger->func_calling_num = local_calling_num; /* function exit timing end */ /* Argument changed from ISDC_OK to 'status' with version 3.0.1 */ return( status ); exittrace: /* exittrace */ if (logger->trace) traces(func_local, 1234567, logger); if (logger->trace) traces(func_local, 1000-status, logger); goto slut; } /* *********************************************************************************** */ void txtConvert(char srcname[], char SrcName[], struct log_con *logger) { int i, i5; i5 = 0; for (i=0; i<31; i++) { if (srcname[i] == 0) { SrcName[i] = 0; return; } if (srcname[i] == ' ') { SrcName[i] = '_'; continue; } SrcName[i] = srcname[i]; } SrcName[i] = 0; sprintf(logger->tmptxt, "Convert %s to %s", srcname, SrcName); logger->logstat = logprint(logger, logger->tmptxt, i5); return; }