LALDetCharHveto.c 17.6 KB
Newer Older
1 2
#include <lal/LALDetCharHveto.h>

3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
/*
 * Okay, now we're playing a bit fast and loose. These are defined for use with
 * the SnglBurst's next pointer. Since we don't use the SnglBurst as a linked
 * list --- leaving that instead to the glib sequence --- the next pointer is
 * hijacked for use with "masking" trigs via the Prune function. Since several
 * rounds are usually required, when we "prune" a trig via its SNR, we don't
 * want to remove and then reload it later --- that takes a long time. Instead,
 * we "mask" it by assigning its next pointer one of the below values. Both
 * are near guaranteed to crash the program if someone attempts to use them
 * in the normal fashion, and are informative if looking at them in gdb.
 *
 * In short, don't use the SnglBurst next pointer for something in the hveto 
 * program. I hereby disclaim responsibility if by some magic, the program
 * actually allocates to 0xDEADBEEFDEADBEFF.
 */
#define TRIG_MASKED (SnglBurst*)0xDEADBEEFDEADBEEF
#define TRIG_NOT_MASKED (SnglBurst*)0x0

21 22 23 24 25 26 27
/*
 * Scan through a list of triggers and count the instances of each trigger type
 * and count its coincidences with the target channel. The hash tables for the
 * channel count and coincidence count must already be initialized. The
 * trig_sequence parameter must contain a full list of the triggers to be
 * scanned. The chan parameter is the 'target' channel -- usually h(t). The
 * twind parameter controls the length of the coincidence window.
28 29 30 31 32
 *
 * coinc_type = 0: Allow all coincidences
 * coinc_type = 1: Allow only one coincidence between target channel and others
 * (e.g. if two triggers from the same channel are found in coincidence, only
 * one is recorded).
33
 */
34
void XLALDetCharScanTrigs( LALGHashTable *chancount, LALGHashTable *chanhist, LALGSequence* trig_sequence, const char* chan, double twind, int coinctype ){
35 36 37 38 39

	/*
	 * pointer to the position in the list for the current auxiliary trigger
	 * being examined
	 */
40
	GSequenceIter *cur_trig = XLALGSequenceBeginRaw( trig_sequence, LALGTYPE_SNGL_BURST );
41 42 43 44 45 46 47 48 49 50 51

	LIGOTimeGPS start;
	LIGOTimeGPS stop;
	XLALGPSSetREAL8( &start, 0 );
	XLALGPSSetREAL8( &stop, 1 );

	// coincidence window
	LALSeg *wind;
	wind = XLALSegCreate( &start, &stop, 0 );

	// trigger pointer
52 53
	SnglBurst *sb_target;
	GHashTable *prevcoinc;
54 55 56 57 58 59

	// Iterate through trigger list -- it should be sorted by GPS time
	// TODO: Use foreach here instead?
	for( ; !g_sequence_iter_is_end(cur_trig) ; cur_trig = g_sequence_iter_next(cur_trig) ){

		// Current trigger
60
		sb_target = (SnglBurst*)g_sequence_get( cur_trig );
61
		if( sb_target->next == TRIG_MASKED ) continue;
62 63 64 65 66

		/*
		 * Increment the trigger count for this channel.
		 * If it is not in the hash table, add it and initialize properly
		 */
67 68
		int *value;
		value = XLALGetGHashTblIntPtr( chancount, sb_target->channel );
69 70
		if( value ){
			(*value)++;
71
			XLALPrintInfo( "Count: Incrementing, %s value: %i\n", sb_target->channel, *value );
72
			//g_hash_table_insert( chancount, XLALStringDuplicate(sb_target->channel), value );
73
		} else {
Chris Pankow's avatar
Chris Pankow committed
74
			XLALPrintInfo( "Count: Adding %s with time %d.%d\n", sb_target->channel, sb_target->peak_time.gpsSeconds, sb_target->peak_time.gpsNanoSeconds );
75 76
                        XLALSetGHashTblInt( chancount, sb_target->channel, 1 );
                        value = XLALGetGHashTblIntPtr( chancount, sb_target->channel );
77 78 79 80 81 82
		}

		// Is it the channel we're looking at?
		// TODO: Consider doing this from the perspective of h(t), rather than
		// each aux channel, since we're checking coincidences of N to 1 rather
		// than 1 to N.
83 84 85 86 87
		// FIXME: check to make sure we don't have channels which are identical
		// up to a point in the string
		if( strstr( sb_target->channel, chan ) ){
			// Yes, create window segment
			start = stop = sb_target->peak_time;
88 89
			start = *XLALGPSAdd( &start, -twind/2.0 );
			stop = *XLALGPSAdd( &stop, twind/2.0 );
Chris Pankow's avatar
Chris Pankow committed
90
			XLALPrintInfo( "creating segment from %s %d.%d %d.%d\n", sb_target->channel, start.gpsSeconds, start.gpsNanoSeconds, stop.gpsSeconds, stop.gpsNanoSeconds  );
91
			XLALSegSet( wind, &start, &stop, sb_target->event_id );
92 93 94 95
		} else { // No, go to the next
			continue;
		}

96 97 98
		// FIXME: Free memory
		prevcoinc = g_hash_table_new( g_str_hash, g_str_equal );

Chris Pankow's avatar
Chris Pankow committed
99
		XLALPrintInfo( "Checking for event %d within %d %d\n", sb_target->peak_time.gpsSeconds, wind->start.gpsSeconds, wind->end.gpsSeconds );
100 101 102 103

		// This is our secondary pointer which counts out from the current
		// trigger
		GSequenceIter *trigp = cur_trig;
104
		SnglBurst* sb_aux;
105 106 107 108
		gboolean begin = FALSE;
		if( g_sequence_iter_is_begin(trigp) ){
			begin = TRUE;
		} else {
109
			trigp = g_sequence_iter_prev(trigp);
110
			sb_aux = (SnglBurst*)g_sequence_get(trigp);
111
		}
112

113 114 115
		// Sweep backward, accumulate triggers in the window until we're outside
		// of it
		while( !begin && XLALGPSInSeg( &sb_aux->peak_time, wind ) == 0 ){
116 117 118 119

			/*
			 * If we want unique coincidences, check to see if this channel has
			 * already been added.
120 121 122
			 *
			 * For now, don't use the target channel
			 * FIXME: We may want this information in the future
123
			 */
124
			if( (coinctype == 1 && g_hash_table_lookup( prevcoinc, sb_aux->channel )) || strstr( sb_aux->channel, chan ) || sb_aux->next == TRIG_MASKED ){
125 126 127 128
				if( g_sequence_iter_is_begin(trigp) ) break;
				trigp = g_sequence_iter_prev(trigp);
				sb_aux = (SnglBurst*)g_sequence_get(trigp);
				continue;
129 130 131 132
			} else {
				// FIXME: Use g_hash_table_add when compatible
				g_hash_table_insert( prevcoinc, &sb_aux->channel, &sb_aux->channel );
			}
133 134

			// TODO: Macroize?
135
			value = XLALGetGHashTblIntPtr( chanhist, sb_aux->channel );
136 137 138 139
			// If we have an entry for this channel, use it, otherwise create a
			// new one
			if( value != NULL ){
				(*value)++;
140
				XLALPrintInfo( "Left Coincidence: Incrementing, %s->%ld, time %d.%d value: %i\n", sb_aux->channel, sb_target->event_id, sb_aux->peak_time.gpsSeconds, sb_aux->peak_time.gpsNanoSeconds, *value );
141
				//g_hash_table_insert( chanhist, &sb_aux->channel, value );
142
			} else {
Chris Pankow's avatar
Chris Pankow committed
143
				XLALPrintInfo( "Left Coincidence: Adding %s->%ld with time %d.%d\n", sb_aux->channel, sb_target->event_id, sb_aux->peak_time.gpsSeconds, sb_aux->peak_time.gpsNanoSeconds );
144 145
                                XLALSetGHashTblInt( chanhist, sb_aux->channel, 1 );
                                value = XLALGetGHashTblIntPtr( chanhist, sb_aux->channel );
146
			}
147 148 149 150 151 152 153 154 155
			if( g_sequence_iter_is_begin(trigp) ) break;
			trigp = g_sequence_iter_prev(trigp);
			sb_aux = (SnglBurst*)g_sequence_get(trigp);
		}

		trigp = cur_trig;
		trigp = g_sequence_iter_next(trigp);
		if( g_sequence_iter_is_end(trigp) ) break;
		sb_aux = (SnglBurst*)g_sequence_get(trigp);
156 157 158

		// Sweep forward, accumulate triggers in the window until we're outside
		// of it
159 160 161 162 163 164
		//do {
		while( XLALGPSInSeg( &sb_aux->peak_time, wind ) == 0 ){
			//trigp = g_sequence_iter_next(trigp);
			//if( g_sequence_iter_is_end(trigp) ) break;
			//sb_aux = (SnglBurst*)g_sequence_get(trigp);
			
165
			// Not the target channel?
166
			if( strstr( sb_aux->channel, chan ) || sb_aux->next == TRIG_MASKED ){ 
167
				trigp = g_sequence_iter_next(trigp);
168
				if( g_sequence_iter_is_end(trigp) ) break;
169 170 171
				sb_aux = (SnglBurst*)g_sequence_get(trigp);
				continue;
			}
172 173

			if( coinctype == 1 && g_hash_table_lookup( prevcoinc, sb_aux->channel ) ){
174 175 176 177
				trigp = g_sequence_iter_next(trigp);
				if( g_sequence_iter_is_end(trigp) ) break;
				sb_aux = (SnglBurst*)g_sequence_get(trigp);
				continue;
178 179 180 181
			} else {
				// FIXME: Use g_hash_table_add when compatible
				g_hash_table_insert( prevcoinc, &sb_aux->channel, &sb_aux->channel );
			}
182 183

			// TODO: Macroize?
184
			value = XLALGetGHashTblIntPtr( chanhist, sb_aux->channel );
185 186 187 188
			// If we have an entry for this channel, use it, otherwise create a
			// new one
			if( value != NULL ){
				(*value)++;
189
				XLALPrintInfo( "Right Coincidence: Incrementing, %s->%ld, time %d.%d value: %i\n", sb_aux->channel, sb_target->event_id, sb_aux->peak_time.gpsSeconds, sb_aux->peak_time.gpsNanoSeconds, *value );
190
				//g_hash_table_insert( chanhist, &sb_aux->channel, value );
191
			} else {
Chris Pankow's avatar
Chris Pankow committed
192
				XLALPrintInfo( "Right Coincidence: Adding %s->%ld with time %d.%d\n", sb_aux->channel, sb_target->event_id, sb_aux->peak_time.gpsSeconds, sb_aux->peak_time.gpsNanoSeconds );
193 194
                                XLALSetGHashTblInt( chanhist, sb_aux->channel, 1 );
                                value = XLALGetGHashTblIntPtr( chanhist, sb_aux->channel );
195
			}
196 197 198 199 200
			trigp = g_sequence_iter_next(trigp);
			if( g_sequence_iter_is_end(trigp) ) break;
			sb_aux = (SnglBurst*)g_sequence_get(trigp);
		}
		//} while( XLALGPSInSeg( &sb_aux->peak_time, wind ) == 0 );
201
		g_hash_table_destroy(prevcoinc);
202 203 204 205 206 207 208 209 210 211 212
	}
	XLALFree( wind );
}

/*
 * Do a round of vetoes. In short, check the significance of each channel and
 * remove the highest signifiance channel, unless it falls below the
 * significance threshold. The winner is returned as the first parameter. The
 * chan parameter is the target channel. The t_ratio parameter is the ratio of
 * the veto window duration to the total examined livetime.
 */
213
double XLALDetCharVetoRound( char** winner, LALGHashTable* chancount, LALGHashTable* chanhist, const char* chan, double t_ratio ){
214
	double mu, sig, max_sig=-1;
215
	int *k;
216

217 218
	// No reference channel present in sequence?
	// This is mostly just here to prevent a segfault later.
219
	if( !XLALGHashTableKeyExists(chancount, chan) ){
220 221 222 223
		XLALPrintInfo( "Reference channel not present in channel count.\n" );
		return -1.0;
	}

224 225 226
	GHashTableIter iter;
	gpointer key, val;
	// Number of triggers in target channel
227
	int n_h = XLALGetGHashTblInt(chancount, chan);
228

229
	XLALGHashTableBeginRaw( chanhist, LALGTYPE_INT, &iter );
230 231 232 233
	// Iterate over each channel and check the significance of the
	// accumulated trigger number
	while( g_hash_table_iter_next( &iter, &key, &val ) ){
		if( strstr( key, chan ) ){
Chris Pankow's avatar
Chris Pankow committed
234
			//XLALPrintInfo( stderr, "Skipping channel %s\n", (char *)key );
235 236 237
			continue;
		}
		// Number of triggers in auxillary channel
238
		int n_aux = XLALGetGHashTblInt(chancount, key);
239
		mu = n_h*t_ratio*n_aux;
240

241
		k = (int*)val;
242

243
		XLALPrintInfo( "Total coincidences for channel %s: %i\n", (char *)key, *k );
244
		XLALPrintInfo( "Mu for channel %s: %g\n", (char *)key, mu );
245
		sig = XLALDetCharHvetoSignificance( mu, *k );
246
		XLALPrintInfo( "Significance for this channel: %g\n", sig );
247
		if( sig > max_sig && !strstr(chan, (char*)key) ){
248 249 250
			max_sig = sig;
			*winner = XLALRealloc(*winner, (strlen((char *)key) + 1) * sizeof(char));
			strcpy( *winner, (char *)key );
251 252
		}
	}
253
	XLALPrintInfo( "winner: %s\n", *winner );
254 255 256
	return max_sig;
}

257
void XLALDetCharPruneTrigs( LALGSequence* trig_sequence, const LALSegList* onsource, double snr_thresh, const char* refchan ){
258 259 260 261 262 263 264
	// FIXME: Actually, this should prune all triggers
	if( onsource->length == 0 ){
		return;
	}
	size_t i = 0;

	GSequenceIter *trigp, *tmp;
265
	trigp = XLALGSequenceBeginRaw( trig_sequence, LALGTYPE_SNGL_BURST );
266 267 268 269

	LALSeg onseg = onsource->segs[0];
	while( !g_sequence_iter_is_end( trigp ) ){
		SnglBurst* sb = g_sequence_get( trigp );
Chris Pankow's avatar
Chris Pankow committed
270

271 272 273 274 275 276 277 278
		int pos = XLALGPSInSeg( &sb->peak_time, &onseg );

		/*
		 * The case for this sentinel:
		 * Last trigger was within the segment previous to this
		 * This trigger is in the subsequent segment, thus we need to advance
		 * the pointer to account for it
		 */
279
		if( pos > 0 && i < onsource->length-1 ){
280 281
			i++;
			onseg = onsource->segs[i];
282 283 284 285 286 287 288 289
			continue;
		/*
		 * We're beyond the extent of the onsource list
		 */
		} else if( pos > 0 && i >= onsource->length ) {
			tmp = trigp;
			trigp = g_sequence_iter_next( trigp );
			g_sequence_remove(tmp);
290 291 292 293 294 295 296 297
			continue;
		}

		if( pos != 0 ){
			tmp = trigp;
			trigp = g_sequence_iter_next( trigp );
			g_sequence_remove(tmp);
			continue;
298 299
		}

Chris Pankow's avatar
Chris Pankow committed
300 301 302 303 304 305 306 307
		/*
		 * This is the reference channel, and we don't check it's snr to remove
		 * it. If no refchan is given, then remove anything according to 
		 * criteria.
		 */
		gboolean isrefchan = (refchan != NULL);
		if( isrefchan ) isrefchan = (strstr( refchan, sb->channel ) != NULL);

308 309 310 311 312
		/*
		 * Mark the trigger as unused, but don't delete it. This is a common
		 * where the snr threshold is raised, but we don't want to remove the
		 * trigger. See warnings and disclaimers at the top of the file.
		 */
313
		if( sb->snr < snr_thresh && !isrefchan ){
314
			sb->next = TRIG_MASKED;
315
		} else {
316
			sb->next = TRIG_NOT_MASKED;
317
		}
318
		trigp = g_sequence_iter_next( trigp );
319 320 321 322 323 324
	}
}

/*
 * Remove triggers in window centered around peak times for triggers from vchan.
 * This is generally called after a veto round.
325 326 327
 *
 * TODO: Merge vetolist creation here
 * TODO: Can we also decrement the count / coincidences efficiently here?
328
 */
329
void XLALDetCharRemoveTrigs( LALGSequence* trig_sequence, LALGSequence** tbd, const LALSeg veto, const char* vchan, const char* refchan, double snr_thresh ){
330 331

	size_t vetoed_events = 0;
332
	size_t de_vetoed_events = 0;
333
	size_t nevents = XLALGetGSequenceLength(trig_sequence);
Adam Mercer's avatar
Adam Mercer committed
334
	XLALPrintInfo( "nevents: %zu\n", nevents );
Chris Pankow's avatar
Chris Pankow committed
335
	XLALPrintInfo( "Channel to veto: %s\n", vchan );
336 337

	// Pointer to the current position in the trigger list
338
	GSequenceIter* trigp = XLALGSequenceBeginRaw( trig_sequence, LALGTYPE_SNGL_BURST );
339 340 341
	// Pointer to the trigger under examination
	SnglBurst* sb;

342
	// Store the triggers to be deleted
343 344
	*tbd = XLALCreateGSequence(LALGTYPE_SNGL_BURST);
	GSequenceIter* tbdit = XLALGSequenceBeginRaw( *tbd, LALGTYPE_SNGL_BURST );
345 346 347 348 349

	// Loop over our list, looking to remove trigs
	while( !g_sequence_iter_is_end(trigp) ){
		sb = (SnglBurst*)g_sequence_get(trigp);
		if( !sb ){
Chris Pankow's avatar
Chris Pankow committed
350
			XLALPrintError( "Invalid pointer for top level iterator!\n" );
351 352 353
		}

		// Is it the channel to veto on?
354
		if( (sb->snr < snr_thresh) | !strstr( sb->channel, vchan ) ){
355 356 357 358 359 360 361 362
			trigp = g_sequence_iter_next(trigp);
			continue; // no, move on
		}
		LALSeg *trig_veto = XLALSegCreate( &veto.start, &veto.end, (int)sb->event_id );
		// Center the window on the trigger
		XLALGPSAddGPS( &trig_veto->start, &sb->peak_time);
		XLALGPSAddGPS( &trig_veto->end, &sb->peak_time);

363
		GSequenceIter *st = trigp, *end = trigp;
364

365 366 367
		gboolean begin = g_sequence_iter_is_begin(st);
		if( !begin ){
			st = g_sequence_iter_prev(st);
368
			sb = (SnglBurst*)g_sequence_get(st);
369
		}
370

371 372 373 374
		// Backwards
		while( !begin & (XLALGPSInSeg( &sb->peak_time, trig_veto ) == 0) ){
			GSequenceIter *tmp;
			tmp = st;
375 376
			if( g_sequence_iter_is_begin(st) ){
				break;
377
			}
378 379 380
			st = g_sequence_iter_prev(st);
			g_sequence_move(tmp, tbdit);
			tbdit = g_sequence_iter_next(tbdit);
381
			vetoed_events++;
382
			XLALPrintInfo( "Backwards, deleting %s (%d.%d) id %ld\n", sb->channel, sb->peak_time.gpsSeconds, sb->peak_time.gpsNanoSeconds, sb->event_id );
383 384 385 386
			if( strstr(refchan, sb->channel) ){
				de_vetoed_events++;
			}
			sb = (SnglBurst*)g_sequence_get(st);
387 388 389 390 391 392 393 394 395
		}

		// Check to make sure that we're not at the end of the list
		if( g_sequence_iter_is_end(end) ){
			break;
		} else {
			sb = (SnglBurst*)g_sequence_get(end);
		}

396
		GSequenceIter *tmp;
397 398 399 400
		// Forwards
		while( XLALGPSInSeg( &sb->peak_time, trig_veto ) == 0 ){
			// don't invalidate the top level iterator
			tmp = end;
401 402 403 404 405 406
			end = g_sequence_iter_next(end);
			// Delete this trigger -- but don't delete anything we'd use later
			if( !strstr(sb->channel, vchan) ){
				g_sequence_move(tmp, tbdit);
				tbdit = g_sequence_iter_next(tbdit);
				vetoed_events++;
407
				XLALPrintInfo( "Forwards, deleting %s (%d.%d) id %ld\n", sb->channel, sb->peak_time.gpsSeconds, sb->peak_time.gpsNanoSeconds, sb->event_id );
408 409 410 411 412
				if( strstr(refchan, sb->channel) ){
					de_vetoed_events++;
				}
			}

413 414 415
			if( g_sequence_iter_is_end(end) ){
				break;
			}
416
			sb = (SnglBurst*)g_sequence_get(end);
417
		}
418 419 420 421 422 423
		tmp = trigp;
		sb = (SnglBurst*)g_sequence_get(tmp);
		trigp = g_sequence_iter_next(trigp);
		g_sequence_move( tmp, tbdit );
		tbdit = g_sequence_iter_next(tbdit);
		vetoed_events++;
424
		XLALPrintInfo( "Veto trig, deleting %s (%d.%d) id %ld\n", sb->channel, sb->peak_time.gpsSeconds, sb->peak_time.gpsNanoSeconds, sb->event_id );
425 426 427 428

		// FIXME: Add to veto list
		XLALFree( trig_veto );

Adam Mercer's avatar
Adam Mercer committed
429
		XLALPrintInfo( "%zu events deleted so far.\n", vetoed_events );
430
		nevents = XLALGetGSequenceLength(trig_sequence);
Adam Mercer's avatar
Adam Mercer committed
431
		XLALPrintInfo( "%zu events remain\n", nevents );
432
	}
Adam Mercer's avatar
Adam Mercer committed
433 434
	XLALPrintInfo( "Done, total events removed %zu\n", vetoed_events );
	XLALPrintInfo( "Done, ref channel total events removed %zu\n", de_vetoed_events );
435 436 437

}

438 439 440
/*
 * Turn all the peak times for channel vchan into a segment list of vetoes.
 */
441
void XLALDetCharTrigsToVetoList( LALSegList* vetoes, LALGSequence* trig_sequence, const LALSeg veto, const char* vchan ){
442 443 444 445

	float wind = XLALGPSDiff(&veto.end, &veto.start);

	// Pointer to the current position in the trigger list
446
	GSequenceIter* trigp = XLALGSequenceBeginRaw( trig_sequence, LALGTYPE_SNGL_BURST );
447 448 449 450 451 452
	// Pointer to the trigger under examination
	SnglBurst* sb;

	while( !g_sequence_iter_is_end(trigp) ){
		sb = (SnglBurst*)g_sequence_get(trigp);
		if( !sb ){
Chris Pankow's avatar
Chris Pankow committed
453
			XLALPrintError( "Invalid pointer for top level iterator!\n" );
454 455 456 457 458 459 460 461 462 463 464 465 466 467 468
		}

		if( strstr( sb->channel, vchan ) ){
            LALSeg vetotmp;
            LIGOTimeGPS start = sb->peak_time;
            LIGOTimeGPS stop = sb->peak_time;
            XLALGPSSetREAL8( &start, -wind/2.0 );
            XLALGPSSetREAL8( &stop, wind/2.0 );
            XLALSegSet( &vetotmp, &start, &stop, sb->event_id );
			XLALSegListAppend( vetoes, &vetotmp );
		}
		trigp = g_sequence_iter_next(trigp);
	}
}

469
/*
470
 * Calculate the significance of a set of triggers from the Poisson survival
471 472
 * function given the expected number of triggers.
 */
473
double XLALDetCharHvetoSignificance( double mu, int k ){
474 475 476 477 478 479 480
	if( mu < 0 ){
		XLALPrintError( "%s(): attempt to calculate significance with a negative mu (%f)", __func__, mu );
	}
	if( k < 0 ){
		XLALPrintWarning( "%s(): attempt to calculate significance with a negative k (%d)", __func__, k );
	}

481 482 483 484 485
	double sig = gsl_sf_gamma_inc_P(k, mu);
	if( sig != 0.0 ){
		return -log10(sig);
	} else {
		return -k*log10(mu) + mu*log10(exp(1)) + gsl_sf_lngamma(k+1)/log(10);
486
	}
487
}