LAL  7.5.0.1-8083555
ComputeDataQualityVector.c
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2008, 2009 Jordi Burguet-Castell, Xavier Siemens
3  *
4  * This program is free software; you can redistribute it and/or modify
5  * it under the terms of the GNU General Public License as published by
6  * the Free Software Foundation; either version 2 of the License, or
7  * (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with with program; see the file COPYING. If not, write to the
16  * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
17  * MA 02110-1301 USA
18  */
19 
20 #include <complex.h>
21 #include <lal/ComputeDataQualityVector.h>
22 
23 #include <math.h> /* to use isnan() and isinf() */
24 
25 
26 /**
27  * Compute the Data Quality Vector as defined in
28  * https://www.lsc-group.phys.uwm.edu/daswg/wiki/S6OnlineGroup/CalibratedData
29  *
30  * Copied for reference:
31  *
32  * State vector channel: IFO:DMT-STATE_VECTOR
33  *
34  * The data type is float(!) so they have to be first converted to int.
35  * Then, bits 0-4 define the states. Bits 5-15 are always set to 1,
36  * such that the word 0xffff means science mode:
37  *
38  * bit0=SCI (operator set to go to science mode)
39  * bit1=CON conlog unsets this bit is non-harmless epics changes
40  * bit2=UP (set by locking scripts)
41  * bit3=!INJ Injections unset this bit
42  * bit4=EXC Unauthorized excitations cause this bit to be unset
43  *
44  * Channel Name: IFO:DMT-DATA_QUALITY_VECTOR where IFO is one of (H1, H2, L1)
45  * Sample Rate: 1 Hz, for any state to be true, it must be true every
46  * sample of h(t) within the second
47  * Type: INT4
48  *
49  * The quality channel will use the following bitmask definition
50  *
51  * SCIENCE 1 // SV_SCIENCE & LIGHT
52  * INJECTION 2 // Injection: same as statevector
53  * UP 4 // SV_UP & LIGHT
54  * CALIBRATED 8 // SV_UP & LIGHT & (not TRANSIENT)
55  * BADGAMMA 16 // Calibration is bad (outside 0.8 < gamma < 1.2)
56  * LIGHT 32 // Light in the arms ok
57  * MISSING 64 // Indication that data was dropped in DMT
58  *
59  * All the arrays must have been previously allocated.
60  * r_x is the number of x_value(s) that are in a single DQ sample (=
61  * x.length/n_dq). So it is the "rate" of x with respect to DQ samples.
62  *
63  * Other special meanings:
64  * t_bad_left: time (in s) of the last NOT-UP event in the Data Quality
65  * t_bad_right: time (in s) of the next NOT-UP event in the Data Quality
66  * wings: duration (in s) of the wings used for calibration
67  *
68  * If t_bad_left < wings then it doesn't matter how much less it
69  * is. Same thing for t_bad_right > wings.
70  */
71 int XLALComputeDQ(REAL4* sv_data, int r_sv,
72  REAL4* lax_data, REAL4* lay_data, int r_light,
73  COMPLEX16* gamma_data, int r_gamma,
74  int t_bad_left, int t_bad_right, int wings,
75  int missing,
76  int* dq_data, int n_dq) /* output */
77 {
78  int i, j; /* counters */
79  float sum_x, sum_y;
80  int light; /* is there light in the arms? */
81  int sv_science, sv_injection, sv_up; /* state vector info */
82  int science, injection, up; /* DQ (not identical to the sv_*) */
83  int calibrated;
84  int badgamma;
85  int dq_value;
86 
87  /* Fill one by one the contents of the DQ vector */
88  for (i = 0; i < n_dq; i++) {
89  /* light */
90  sum_x = 0;
91  sum_y = 0;
92  for (j = 0; j < r_light; j++) {
93  sum_x += lax_data[i*r_light + j];
94  sum_y += lay_data[i*r_light + j];
95  }
96 
97  light = (sum_x/r_light > 100 && sum_y/r_light > 100);
98  /* "is the mean higher than 100 for both arms?" */
99 
100  /* science, injection, up (stuff coming from the state vector) */
101  sv_science = 1; /* in science mode */
102  sv_injection = 0; /* with no injection going on */
103  sv_up = 1; /* and IFO is up */
104  for (j = 0; j < r_sv; j++) { /* go over SV samples in a DQ sample */
105  int s = (int) sv_data[i*r_sv + j]; /* convert from float to int */
106  if ((s & (1 << 0)) == 0) sv_science = 0;
107  if ((s & (1 << 3)) == 0) sv_injection = 1;
108  if ((s & (1 << 2)) == 0) sv_up = 0;
109  }
110 
111  science = sv_science && light;
112  injection = sv_injection;
113  up = sv_up && light; /* these are the definitions for the DQ vector */
114 
115 
116  /* calibrated */
117  /* Because we will have to compute UP for the Data Quality
118  * everywhere before anything, to know if something funny
119  * happens within a "wings" distance from this data, we will
120  * leave the computation of the calibrated flag for next
121  * loop.
122  */
123  /* calibrated = up && (! transient); */
124 
125  /* badgamma */
126  badgamma = 0;
127 
128  for (j = 0; j < r_gamma; j++) {
129  REAL8 re = creal(gamma_data[i*r_gamma + j]);
130  if (re < 0.8 || re > 1.2 || isnan(re) || isinf(re))
131  badgamma = 1;
132  }
133 
134  /* data quality */
135  dq_value = 0;
136  if (science) dq_value += (1 << 0);
137  if (injection) dq_value += (1 << 1);
138  if (up) dq_value += (1 << 2);
139  /* if (calibrated) dq_value += (1 << 3); */ /* we'll do that later */
140  if (badgamma) dq_value += (1 << 4);
141  if (light) dq_value += (1 << 5);
142  if (missing) dq_value += (1 << 6); /* directly from the argument */
143 
144  dq_data[i] = dq_value;
145  }
146 
147  /* Now look for the transients and fill the "calibrated" bit. */
148  for (i = 0; i < n_dq; i++) {
149  calibrated = 1;
150 
151  if (i - wings < t_bad_left || i + wings > n_dq + t_bad_right)
152  calibrated = 0;
153 
154  for (j = 0; j < wings; j++) {
155  int pos = i - j;
156  if (pos > 0) {
157  if ((dq_data[pos] & (1 << 2)) == 0) /* if not up */
158  calibrated = 0;
159  }
160  pos = i + j; /* note that this includes dq_data[i] having UP=1 */
161  if (pos < n_dq) {
162  if ((dq_data[pos] & (1 << 2)) == 0) /* if not up */
163  calibrated = 0;
164  }
165  }
166 
167  if (calibrated) dq_data[i] += (1 << 3);
168  /* we checked that we were DQ up all the time in +/- wings
169  * seconds, so that also includes sv_up && light */
170  }
171 
172  return 0;
173 }
int XLALComputeDQ(REAL4 *sv_data, int r_sv, REAL4 *lax_data, REAL4 *lay_data, int r_light, COMPLEX16 *gamma_data, int r_gamma, int t_bad_left, int t_bad_right, int wings, int missing, int *dq_data, int n_dq)
Compute the Data Quality Vector as defined in https://www.lsc-group.phys.uwm.edu/daswg/wiki/S6OnlineG...
double complex COMPLEX16
Double-precision floating-point complex number (16 bytes total)
double REAL8
Double precision real floating-point number (8 bytes).
float REAL4
Single precision real floating-point number (4 bytes).