/*=============================================================================
find_coinc - A simple utility to find coincidences between event lists
Written Aug 2001 by Peter Shawhan
Uses the "Metaio" parsing code by Philip Charlton
=============================================================================*/

#include <stdio.h>
#include <string.h>
#include "metaio.h"

#define MAXSTACK 100

/*===========================================================================*/
int main( int argc, char **argv )
{
    char file1[256], file2[256];
    double window;
    int status, status2, file2eof=0, i, ncoinc=0;
    int isec1, ins1, isec2, ins2;
    int sec1, ns1;
    double time1, time2, delta;

    struct MetaioParseEnvironment parseEnv1, parseEnv2;
    const MetaioParseEnv env1 = &parseEnv1;
    const MetaioParseEnv env2 = &parseEnv2;

    /* Stack of rows from the second file */
    int nstack = 0;
    int ssec[MAXSTACK], sns[MAXSTACK];
    int newlow;

    /*------ Beginning of code ------*/

    /*-- Check number of command-line arguments --*/
    if ( argc != 4 ) {
      printf( "Look for coincidences in two event lists\n" );
      printf( "Usage:  find_coinc <file1> <file2> <window>\n" );
      return 3;
    }

    /*-- Parse command-line arguments --*/
    strncpy( file1, argv[1], sizeof(file1) ); file1[sizeof(file1)-1] = 0;
    strncpy( file2, argv[2], sizeof(file2) ); file2[sizeof(file2)-1] = 0;
    if ( sscanf( argv[3], "%lf", &window ) < 1 ) {
      printf( "Error parsing window size %s\n", argv[3] );
      return 1;
    }

    /*-- Open the files and read up to the beginning of each Stream --*/
    status = MetaioOpen( env1, file1 );
    if ( status != 0 ) {
      printf( "Error opening file %s\n", file1 );
      MetaioClose( env1 );
      return 2;
    }

    status = MetaioOpen( env2, file2 );
    if ( status != 0 ) {
      printf( "Error opening file %s\n", file2 );
      MetaioClose( env2 );
      return 2;
    }

    /*-- Locate the "start_time" and "start_time_ns" columns in each file" --*/
    isec1 = MetaioFindColumn( env1, "start_time" );
    ins1 = MetaioFindColumn( env1, "start_time_ns" );
    if ( isec1 < 0 || ins1 < 0 ) {
      printf( "File %s does not contain start_time and/or start_time_ns\n",
	      file1 );
      MetaioClose( env1 ); MetaioClose( env2 ); return 4;
    }

    isec2 = MetaioFindColumn( env2, "start_time" );
    ins2 = MetaioFindColumn( env2, "start_time_ns" );
    if ( isec2 < 0 || ins2 < 0 ) {
      printf( "File %s does not contain start_time and/or start_time_ns\n",
	      file2 );
      MetaioClose( env1 ); MetaioClose( env2 ); return 4;
    }

    /*--------------------------------------*/
    /*-- Loop over rows in the first file --*/
    while ( (status=MetaioGetRow(env1)) == 1 ) {

      /*-- Get the time of this event --*/
      sec1 = env1->ligo_lw.table.elt[isec1].data.int_4s;
      ns1 = env1->ligo_lw.table.elt[ins1].data.int_4s;

      /*-- Check for coincidences (within the window) from the second file --*/
      /*-- We check the stack, enlarging it as needed --*/
      newlow = 0;
      for ( i=0; i<nstack+1; i++ ) {

	/*-- If i==nstack, then we need to enlarge the stack --*/
	if ( i == nstack ) {
	  /*-- If we hit the end of file2, break --*/
	  if ( file2eof ) { break; }

	  /*-- Make sure we have more room in the stack --*/
	  if ( nstack == MAXSTACK ) {
	    printf( "Ran out of space in the event stack (size=%d)\n",
		    MAXSTACK );
	    MetaioClose( env1 ); MetaioClose( env2 ); return 5;
	  }

	  /*-- Read a row --*/
	  status2 = MetaioGetRow(env2);
	  if ( status2 == -1 ) {
	    printf( "Error while getting row from file %s\n", file2 );
	    MetaioClose( env1 ); MetaioClose( env2 ); return 6;
	  } else if ( status2 == 0 ) {
	    /*-- Reached end of file --*/
	    file2eof = 1;
	    break;
	  }

	  /*-- Add this event time to the stack --*/
	  ssec[nstack] = env2->ligo_lw.table.elt[isec2].data.int_4s;
	  sns[nstack] = env2->ligo_lw.table.elt[ins2].data.int_4s;
	  nstack++;
	}

	delta = (double) (ssec[i]-sec1) + 0.000000001 * (double) (sns[i]-ns1);

	/*-- If stack row is earlier by more than the size of the window,
	  get ready to delete it --*/
	if ( delta < -window ) {
	  newlow++;
	} else if ( delta <= window ) {
	  /*-- Found a coincidence! --*/
	  ncoinc++;
	  printf( "Coinc: time1=%d.%09d time2=%d.%09d delta=%.9f\n",
		  sec1, ns1, ssec[i], sns[i], delta );
	} else {
	  /*-- Stack row is too far in the future --*/
	  break;
	}
      }
      /*-- If necessry, compress the stack to get rid of far-past rows --*/
      if ( newlow > 0 ) {
	for ( i=newlow; i<nstack; i++ ) {
	  ssec[i-newlow] = ssec[i];
	  sns[i-newlow] = sns[i];
	}
	nstack -= newlow;
      }

    }

    /*-- Close files --*/
    MetaioClose(env1);
    MetaioClose(env2);

    printf( "Total coincidences found: %d\n", ncoinc );

    return 0;
}