1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
|
/*****************************************************************************
* Johns Hopkins University
* Center For Astrophysical Sciences
* FUSE
*****************************************************************************
*
* Synopsis: cf_apply_dead_time(fitsfile *header, long nevents, float *time,
* float *weight, long nseconds, float *timeline,
* float *ids_dtc, float *elec_dtc)
*
* Description: Applies dead-time correction to all photons. For HIST data,
* uses mean of dead-time arrays. Writes mean dead-time
* values to file header. Issue warning if DET_DEAD > 1.5.
*
* Arguments: fitsfile *header Pointer to FITS file containing the
* header of the intermediate data file
* long nevents The number of events
* float *time An array of event times
* float *weight An array of event weights
* long nseconds The number of timeline values
* float *timeline An array of timeline times
* float *aic_rate Active Image Counter rate
* float *ids_dtc Output of cf_ids_dead_time
* float *elec_dtc Output of cf_elec_dead_time
*
* Calls:
*
* Return: 0 on success
*
* History: 08/01/03 1.1 wvd Begin work
* 08/04/03 1.2 wvd Convert aic_rate to type short.
* 11/25/03 1.3 wvd For HIST data, require that
* mean_tot_dtc be >= 1.
* 11/26/03 1.4 wvd Change aic_rate to type float.
* 02/09/04 1.5 wvd Re-write with new logic.
* 04/09/04 1.6 bjg Include stdlib.h
* 11/07/05 1.7 wvd Add a test to catch bad DTC values.
* 03/10/06 1.8 wvd Clean up i/o.
* 12/29/06 1.9 wvd Check DET_DEAD rather than TOT_DEAD.
* 04/07/07 1.10 wvd Clean up compiler errors.
*
****************************************************************************/
#include <stdlib.h>
#include <string.h>
#include "calfuse.h"
int
cf_apply_dead_time(fitsfile *header, long nevents, float *time, float *weight,
long nseconds, float *timeline, float *ids_dtc, float *elec_dtc)
{
char CF_PRGM_ID[] = "cf_apply_dead_time";
char CF_VER_NUM[] = "1.10";
char comment[FLEN_CARD], datestr[FLEN_CARD], instmode[FLEN_VALUE];
char hkexists[FLEN_VALUE];
double ctimeb, ctimee, eng_time;
int errflg=0, status=0, timeref;
int det_flag=FALSE, ids_flag=FALSE;
long j, k;
float *tot_dtc, det_dead, ids_dead, mean_tot_dtc=0.;
cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing");
if ((errflg = cf_proc_check(header, CF_PRGM_ID))) return errflg;
/* Read header keywords. Set flag if DET_DEAD > 1.5 */
FITS_read_key(header, TSTRING, "INSTMODE", instmode, NULL, &status);
FITS_read_key(header, TSTRING, "HKEXISTS", hkexists, NULL, &status);
FITS_read_key(header, TFLOAT, "DET_DEAD", &det_dead, NULL, &status);
if (det_dead > 1.5) det_flag = TRUE;
/*
* If there is no housekeeping file, AND the engineering snapshot times
* are bad, AND the mean value of any DTC array is > 1.5, set the
* combined DTC array to unity and issue a warning.
*/
if(!strncasecmp(hkexists, "N", 1)) {
FITS_read_key(header, TDOUBLE, "CTIME_B", &ctimeb, NULL, &status);
FITS_read_key(header, TDOUBLE, "CTIME_E", &ctimee, NULL, &status);
if ((eng_time = (ctimee-ctimeb) * 86400.) < 1.) {
FITS_read_key(header, TFLOAT, "IDS_DEAD", &ids_dead, NULL, &status);
if (ids_dead > 1.5) ids_flag = TRUE;
}
}
/*
* Compute total dead-time correction and its mean.
*/
tot_dtc = (float *) cf_malloc(sizeof(float) * nseconds);
if (det_flag || ids_flag) {
for (k = 0; k < nseconds; k++) tot_dtc[k] = 1.0;
mean_tot_dtc = 1.0;
}
else {
for (k = 0; k < nseconds; k++) {
tot_dtc[k] = elec_dtc[k] * ids_dtc[k];
mean_tot_dtc += tot_dtc[k];
}
mean_tot_dtc /= nseconds;
}
/* In TTAG mode, scale weights by tot_dtc as a function of time. */
if (!strncmp(instmode, "TTAG", 4)) for (j=k=0; j<nevents; j++) {
while(timeline[k+1]-FRAME_TOLERANCE < time[j] && k+1 < nseconds)
k++;
weight[j] *= tot_dtc[k];
}
/* In HIST mode, scale all weights by mean dead-time correction. */
else if (mean_tot_dtc > 1.)
for (j = 0; j < nevents; j++)
weight[j] *= mean_tot_dtc;
else mean_tot_dtc = 1.;
cf_verbose(2, "Mean total dead-time correction: %f", mean_tot_dtc);
FITS_update_key(header, TFLOAT, "TOT_DEAD", &mean_tot_dtc, NULL, &status);
if (det_flag || ids_flag) {
FITS_write_comment(header, " ", &status);
if (ids_flag) {
sprintf(comment,
"IDS_DEAD > 1.5. Setting dead-time correction to unity.");
cf_if_warning(comment);
FITS_write_comment(header, comment, &status);
}
if (det_flag) {
sprintf(comment,
"DET_DEAD > 1.5. Setting dead-time correction to unity.");
cf_if_warning(comment);
FITS_write_comment(header, comment, &status);
}
fits_get_system_time(datestr, &timeref, &status);
sprintf(comment, "CalFUSE v%s %.10s", CALFUSE_VERSION, datestr);
FITS_write_comment(header, comment, &status);
FITS_write_comment(header, " ", &status);
}
free (tot_dtc);
cf_proc_update(header, CF_PRGM_ID, "COMPLETE");
cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done processing");
return 0;
}
|
