abssexam.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpctac.h"
#include "tpcstatist.h"
#include "tpcabss.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Scanditronics and GEMS on-line blood sampler have sometimes failed to",
  "correctly record the count data from either of the two channels, seen as",
  "substantially lower or zero counts.",
  "Use this program to examine the channel coincidence count ratio from",
  "multiple files. Program examines individual ABSS data files whether",
  "the channel ratio exceeds a predermined limit, and if it does, then prints",
  "the following information:",
  "  column #1: measurement date,",
  "  column #2: Lower-to-higher channel ratio (the closer to 1 the better),",
  "  column #3: channel 1 mean cps from ten highest overall measurements,",
  "  column #4: channel 2 mean cps from ten highest overall measurements,",
  "  column #5: file name.",
  " ",
  "Usage: @P [options] abssfile(s)",
  " ",
  "Options:",
  " -limit=<limit for channel ratio>",
  "     Limit for the lower-to-higher channel ratio; by default, 0.7.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "     @P *.bld",
  " ",
  "See also: abssfch, absscal, abssbkg, absszero, absstime",
  " ",
  "Keywords: input, blood, calibration, ABSS",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* Turn on the globbing of the command line, since it is disabled by default in
   mingw-w64 (_dowildcard=0); in MinGW32 define _CRT_glob instead, if necessary;
   In Unix&Linux wildcard command line processing is enabled by default. */
/*
#undef _CRT_glob
#define _CRT_glob -1
*/
int _dowildcard = -1;
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int     ai, help=0, version=0, verbose=1;
  int     ret, fileNr=0, printNr, firstfile=0;
  char   *cptr;
  double  ratioLimit=0.70;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    if(strncasecmp(cptr, "LIMIT=", 6)==0 && strlen(cptr)>6) {
      ratioLimit=atofVerified(cptr+6); if(!isnan(ratioLimit)) continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-3;
  
  /* Print help or version? */
  if(help==2) {tpcHtmlUsage(argv[0], info, ""); return(0);}
  if(help) {tpcPrintUsage(argv[0], info, stdout); return(0);}
  if(version) {tpcPrintBuild(argv[0], stdout); return(0);}
 
  /* Process other arguments, starting from the first non-option */
  for(; ai<argc; ai++) {
    if(firstfile==0) firstfile=ai;
    fileNr++;
  }
  /* Did we get all the information that we need? */
  if(fileNr==0) {tpcPrintUsage(argv[0], info, stderr); return(1);}
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("fileNr := %d\n", fileNr);
    printf("ratioLimit := %g\n", ratioLimit);
    fflush(stdout);
  }
 
 
  /*
   *  Process each file
   */
  char abssfile[FILENAME_MAX];
  TAC abss; tacInit(&abss);
  fileNr=printNr=0;
  double r=1.0;
  for(ai=firstfile; ai<argc; ai++) {
 
    /*
     *  Read file
     */
    strlcpy(abssfile, argv[ai], FILENAME_MAX); 
    if(verbose>5) {fprintf(stdout, "%s\n", abssfile); fflush(stdout);}
    ret=tacRead(&abss, abssfile, &status);
    if(ret!=TPCERROR_OK) {
      if(verbose>0) {
        fprintf(stdout, "%s\n", abssfile); fflush(stdout);
        fprintf(stderr, "Error: %s\n", errorMsg(status.error));
        fflush(stderr);
      }
      tacFree(&abss); continue;
    }
    /* This program is useful only for Scanditronics/GEMS sampler */
    if(abss.format!=TAC_FORMAT_ABSS_GEMS &&
       abss.format!=TAC_FORMAT_ABSS_SCANDITRONICS)
    {
      if(verbose>0) {
        fprintf(stderr, "bypassed %s : %s format.\n", 
                abssfile, tacFormattxt(abss.format));
        fflush(stderr);
      }
      tacFree(&abss); continue;
    }
 
    /* 
     *  Examine the data
     */
    /* Divide coincidents by frame duration */
    ret=abssCalculateCps(&abss);
    if(ret!=0) {
      if(verbose>0) {
        fprintf(stderr, "bypassed %s : invalid contents.\n", abssfile);
        fflush(stderr);
      }
      tacFree(&abss); continue;
    }
    /* Allocate extra column for the channel mean */
    if(tacAllocateMore(&abss, 1)!=TPCERROR_OK) {
      fprintf(stderr, "bypassed %s : invalid contents.\n", abssfile);
      tacFree(&abss); continue;
    }
    abss.tacNr++;
    strcpy(abss.c[abss.tacNr-1].name, "mean");
    /* Calculate the channel mean */
    ret=abssChannelMean(&abss, abss.c[abss.tacNr-1].y);
    /* Sort data by the mean column */
    if(ret==0) ret=tacSortByConc(&abss, abss.tacNr-1, &status);
 
    /* Calculate the mean from the first 10 frames from the two channels */
    //printf("x1=%g\n", abss.x1[0]);
    unsigned int n, m1=0, m2=0;
    int c1=0, c2=3;
    double avg1, avg2;
    n=10; if(n>(unsigned int)abss.sampleNr) n=(unsigned int)abss.sampleNr;
    ret=statMeanSD(abss.c[c1].y, n, &avg1, NULL, &m1);
    if(ret==0) ret=statMeanSD(abss.c[c2].y, n, &avg2, NULL, &m2);
    if(ret!=0 || m1!=m2) {
      if(verbose>0) {
        fprintf(stderr, "bypassed %s : invalid contents.\n", abssfile);
        fflush(stderr);
      }
      tacFree(&abss); continue;
    }
    /* and then the smaller-to-larger -ratio */
    r=0.0;
    if(avg1>avg2) {
      if(fabs(avg1)>1.0E-10) r=avg2/avg1;
    } else if(avg1<avg2) {
      if(fabs(avg2)>1.0E-10) r=avg1/avg2;
    } else if(avg1==avg2) {
      r=1.0;
    }
    //printf("r=%g\n", r);
 
    /* Get the date for printing */
    char abssdate[128];
    if(tacGetHeaderScanstarttime(&abss.h, abssdate, NULL)==TPCERROR_OK) {
      abssdate[10]=(char)0;
    } else {
      strcpy(abssdate, "1900-01-01");
    }
 
    /* If ratio was not lower than the limit, then continue with next file
       or print it in verbose mode
     */
    if(r>=ratioLimit) {
      if(verbose>1) {
        if(printNr==0)
          fprintf(stdout, "%-10s\t%-s\t%-s\t%-s\t%-s\n",
                  "date", "ratio", "ch1", "ch2", "file");
        fprintf(stdout, "%s\t%5.3f\t\%6.1f\t%6.1f\t%s\tOK\n", 
               abssdate, r, avg1, avg2, abssfile);
        printNr++;
      }
      fflush(stdout);
      tacFree(&abss);
      fileNr++;
      continue;
    }
 
    /* If ratio was lower than its limit, then report it */
    if(printNr==0 && verbose>0)
      fprintf(stdout, "%-10s\t%-s\t%-s\t%-s\t%-s\n",
              "date", "ratio", "ch1", "ch2", "file");
    fprintf(stdout, "%s\t%5.3f\t\%6.1f\t%6.1f\t%s", 
           abssdate, r, avg1, avg2, abssfile);
    if(verbose>1) printf("\tFAILED\n"); else printf("\n");
    fflush(stdout);
    tacFree(&abss);
    fileNr++; printNr++;
  } /* next ABSS file */
 
  /* Check that any files could be successfully processed */
  if(fileNr<1) {
    fprintf(stderr, 
      "Error: no on-line sampler datafile was successfully examined.\n");
    return(1);
  } else {
    if(verbose>0)
      fprintf(stdout, "Ready: %d datafiles were examined.\n", fileNr);
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/