abssfch.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 "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 calculate the channel1-to-channel2 ratio of coincidence",
  "counts, or, to correct the data from failed channel using user-provided",
  "correct channel1-to-channel2 ratio.",
  " ",
  "Usage: @P [options] abssfile [ratio]",
  " ",
  "Options:",
  " -o=<filename>",
  "     Name for corrected file; by default the original file is overwritten.",
  " -fix=<channel>",
  "     Set the channel (1 or 2) to fix; by default the channel with lower",
  "     counts is fixed.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example 1: calculate channel1-to-channel2 ratio.",
  "     @P us1328.bld",
  "Example 2: correct the failed data with given channel1-to-channel2 ratio.",
  "     @P -o=us1328_corr.bld us1328.bld 0.982",
  " ",
  "See also: abssexam, 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;
  char    abssfile[FILENAME_MAX], outfile[FILENAME_MAX];
  char   *cptr;
  double  chRatio=nan("");
  int     chFix=0;
  int     ret;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  abssfile[0]=outfile[0]=(char)0;
  /* 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, "O=", 2)==0 && strlen(cptr)>2) {
      strlcpy(outfile, cptr+2, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "FIX=", 4)==0 && strlen(cptr)>4) {
      chFix=atoi(cptr+4); if(chFix==1 || chFix==2) 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 */
  if(ai<argc) strlcpy(abssfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    if(atofCheck(argv[ai], &chRatio) || chRatio<=0.0) {
      fprintf(stderr, "Error: invalid channel ratio '%s'.\n", argv[ai]);
      return(1);
    }
    ai++;
  }
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  /* Did we get all the information that we need? */
  if(!abssfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(!outfile[0]) strcpy(outfile, abssfile);
  if(chFix!=0 && isnan(chRatio)) {
    fprintf(stderr, "Error: missing channel ratio.\n");
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("abssfile := %s\n", abssfile);
    if(!isnan(chRatio)) printf("chRatio := %g\n", chRatio);
    if(chFix!=0) printf("chFix := %d\n", chFix);
    printf("outfile := %s\n", outfile);
    fflush(stdout);
  }
 
  /*
   *  Read ABSS file
   */
  if(verbose>1) printf("reading %s\n", abssfile);
  TAC abss; tacInit(&abss);
  ret=tacRead(&abss, abssfile, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&abss); return(2);
  }
  if(verbose>1) {
    printf("sampleNr := %d\n", abss.sampleNr);
  }
  if(verbose>2) {
    printf("fileformat := %d\n", abss.format);
    printf("tacNr := %d\n", abss.tacNr);
  }
  if(abss.format!=TAC_FORMAT_ABSS_ALLOGG && 
     abss.format!=TAC_FORMAT_ABSS_ALLOGG_OLD &&
     abss.format!=TAC_FORMAT_ABSS_GEMS &&
     abss.format!=TAC_FORMAT_ABSS_SCANDITRONICS)
  {
    fprintf(stderr, "Error: not valid ABSS format.\n");
    tacFree(&abss); return(2);
  }
  if(verbose>100) abssWrite(&abss, stdout, NULL);
  if(abss.tacNr<1 || abss.sampleNr<1) {
    fprintf(stderr, "Error: not valid ABSS format.\n");
    tacFree(&abss); return(2);
  }
  /* This program is useful only for Scanditronics/GEMS sampler */
  if(abss.format!=TAC_FORMAT_ABSS_GEMS &&
     abss.format!=TAC_FORMAT_ABSS_SCANDITRONICS)
  {
    fprintf(stderr, "Error: not applicable to provided ABSS system.\n");
    tacFree(&abss); return(2);
  }
 
  /* 
   *  Check for dead channels
   */
  if(verbose>1) printf("checking for dead channels\n");
  int n1, n2, n;
  n1=n2=n=0;
  ret=abssAboveZero(&abss, &n1, &n2, &n);
  if(ret!=0) {
    fprintf(stderr, "Error: not valid ABSS format.\n");
    tacFree(&abss); return(3);
  }
  if(verbose>1) {
    printf("positive_samples[1] := %d\n", n1);
    printf("positive_samples[2] := %d\n", n2);
    printf("positive_samples := %d\n", n);
    fflush(stdout);
  }
  if(n<1) {
    fprintf(stderr, "Error: bad ABSS data.\n"); fflush(stderr);
    tacFree(&abss); return(4);
  }
  if(n<abss.sampleNr/10) {
    fprintf(stderr, "Warning: bad ABSS data.\n"); fflush(stderr);
  }
  int c1ok, c2ok;
  if(n1<=n/2) c1ok=0; else c1ok=1;
  if(n2<=n/2) c2ok=0; else c2ok=1;
  /* If both channels have problems, then we can do nothing */
  if(!c1ok && !c2ok) {
    fprintf(stderr, "Error: bad ABSS data.\n"); fflush(stderr);
    tacFree(&abss); return(4);
  }
  if(verbose>0 && !c1ok) printf("channel1 is dead\n");
  if(verbose>0 && !c2ok) printf("channel2 is dead\n");
  /* Check that user did not want to fix the only decent channel */
  if(chFix==1 && !c2ok) {
    fprintf(stderr, "Error: channel2 cannot be used to fix channel1.\n");
    tacFree(&abss); return(4);
  }
  if(chFix==2 && !c1ok) {
    fprintf(stderr, "Error: channel1 cannot be used to fix channel2.\n");
    tacFree(&abss); return(4);
  }
 
 
  /*
   *  If neither channel is dead, then check which one gave more often
   *  lower counts
   */
  if(c1ok && c2ok) {
    ret=abssHigherCounts(&abss, &n1, &n2);
    if(ret!=0) {
      fprintf(stderr, "Error: not valid ABSS format.\n");
      tacFree(&abss); return(3);
    }
    if(verbose>1) {
      printf("higher_samples[1] := %d\n", n1);
      printf("higher_samples[2] := %d\n", n2);
      fflush(stdout);
    }
  }
 
 
  /* 
   *  Calculate the channel ratio, and (usually) print it in stdout
   */
  if(verbose>1) {printf("calculating the ratio\n"); fflush(stdout);}
  /* Ratio can only be calculated if both channels are okeish */
  if(isnan(chRatio) && (c1ok+c2ok)<2) {
    fprintf(stderr, "Error: cannot compute channel ratio.\n"); fflush(stderr);
    tacFree(&abss); return(4);
  }
  double ratio=1.0;
  if((c1ok+c2ok)==2) {
    int nratio;
    ret=abssCalculateRatio(&abss, &ratio, &nratio);
    if(ret==0) {
      if(verbose>1 || nratio<10) {
        printf("channel ratio calculated from %d samples.\n", nratio);
        fflush(stdout);
      }
      if(isnan(chRatio) || verbose>0) {
        fprintf(stdout, "Ch1/Ch2 := %.6f\n", ratio);
        fflush(stdout);
      }
    } else {
      if(isnan(chRatio)) {
        fprintf(stderr, "Error: cannot compute channel ratio.\n"); 
        fflush(stderr); tacFree(&abss); return(5);
      }
      ratio=1.0;
    }
  }
 
  /* If user did not give the true Ch1/Ch2, then quit */
  if(isnan(chRatio)) {
    if(verbose>1) printf("Ch1/ch2 not given; no correction done.\n");
    tacFree(&abss); return(0);
  }
 
  /*
   *  Decide which channel to fix, if user did not tell it
   */
  if(chFix==0) {
    if((c1ok+c2ok)==1) {
      if(c1ok) chFix=2; else chFix=1;
    } else if(n1>n2 && ratio>1.0) {
      chFix=2;
    } else if(n2>n1 && ratio<1.0) {
      chFix=1;
    } else {
      fprintf(stderr, "Error: not sure which channel to fix.\n");
      fflush(stderr); tacFree(&abss); return(6);
    }
    if(verbose>1) printf("chFix := %d\n", chFix);
  }
 
  /* 
   *  Fix one of the channels
   */
  if(verbose>1) {printf("fixing channel %d\n", chFix); fflush(stdout);}
  ret=abssFixChannel(&abss, chFix, chRatio);
  if(ret!=0) {
    fprintf(stderr, "Error: cannot fix channel counts.\n"); 
    fflush(stderr); tacFree(&abss); return(7);
  }
  if(verbose>3) {
    /* Compute the channel ratio again */
    if(abssCalculateRatio(&abss, &ratio, NULL)==0) {
      printf("channel ratio after correction := %g\n", ratio);
      fflush(stdout);
    }
  }
 
  /* Write the corrected data */
  if(verbose>1) {
    printf("writing corrected data file in %s\n", outfile); fflush(stdout);}
  FILE *fp;
  fp=fopen(outfile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open file for writing.\n");
    tacFree(&abss); return(11);
  }
  ret=abssWrite(&abss, fp, &status);
  fclose(fp);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&abss); return(12);
  }
  if(verbose>0) {
    printf("corrected data written in %s\n", outfile); 
    fflush(stdout);
  }
 
  tacFree(&abss);  
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/