abssbkg.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[] = {
  "Subtract constant background radioactivity in on-line blood sampler data.",
  "The average of coincidence cps from 0 to specified time is calculated",
  "and subtracted from the data; in case of Scanditronics/GEMS this is",
  "a separate process for the two detector pairs, and in case of Allogg",
  "separate for the two channels.",
  "The background radioactivity is assumed to be stable during the assay",
  "(not decaying); do not use this program in case of momentary background",
  "radiation.",
  "Notice that the original data file is overwritten by default.",
  " ",
  "Usage: @P [options] abssfile time",
  " ",
  "Options:",
  " -o=<filename>",
  "     Name for corrected file; by default the original file is overwritten.",
  " --dry",
  "     Dry run: background is shown but data is not changed.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example 1: calculate background from 0-60 s and subtract it.",
  "     @P -o=us1328_corr.bld us1328.bld 60",
  "Example 2: calculate mean radioactivity [cps] from the measurement.",
  "     @P --dry us1328.bld 999999",
  " ",
  "See also: absszero, absscal, absstime, abssfch, abssexam, fitdelay",
  " ",
  "Keywords: input, blood, background, 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  bkgTime=nan("");
  int     dryRun=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(strcasecmp(cptr, "DRY")==0) {
      dryRun=1; 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], &bkgTime) || bkgTime<1.0) {
      fprintf(stderr, "Error: invalid bkg time '%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(isnan(bkgTime)) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(!outfile[0]) strcpy(outfile, abssfile);
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("abssfile := %s\n", abssfile);
    printf("bkgTime := %g\n", bkgTime);
    printf("dryRun := %d\n", dryRun);
    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);
  }
 
 
  /* 
   *  Calculate the background
   */
  if(verbose>1) {printf("calculating the background\n"); fflush(stdout);}
  int i, n1=0, n2=0;
  double fdur, a1, a2, bkg1, bkg2; 
  bkg1=bkg2=0.0;
  if(abss.format==TAC_FORMAT_ABSS_SCANDITRONICS ||
     abss.format==TAC_FORMAT_ABSS_GEMS)
  {
    for(i=0, n1=n2=0; i<abss.sampleNr; i++) {
      a1=abss.c[0].y[i]; a2=abss.c[3].y[i];
      if(verbose>5) {
        printf("%g %g    %g %g\n", abss.x1[i], abss.x2[i], a1, a2);
        fflush(stdout);
      }
      fdur=abss.x2[i]-abss.x1[i];
      if(isnan(fdur) || fdur<1.0E-06 || isnan(abss.x[i])) continue;
      if(abss.x[i]>bkgTime) break;
      a1/=fdur; a2/=fdur;
      if(!isnan(a1)) {bkg1+=a1; n1++;}
      if(!isnan(a2)) {bkg2+=a2; n2++;}
    }
    bkg1/=(double)n1;
    bkg2/=(double)n2;
    if(verbose>0 || dryRun!=0) {
      printf("background_ch1[cps] := %g\n", bkg1);
      printf("background_ch2[cps] := %g\n", bkg2);
    }
  } else if(abss.format==TAC_FORMAT_ABSS_ALLOGG_OLD ||
            abss.format==TAC_FORMAT_ABSS_ALLOGG)
  {
    for(i=0, n1=n2=0; i<abss.sampleNr; i++) {
      a1=abss.c[0].y[i]; a2=abss.c[1].y[i];
      if(verbose>5) {
        printf("%g %g    %g %g\n", abss.x1[i], abss.x2[i], a1, a2);
        fflush(stdout);
      }
      fdur=abss.x2[i]-abss.x1[i];
      if(isnan(fdur) || fdur<1.0E-06 || isnan(abss.x[i])) continue;
      if(abss.x[i]>bkgTime) break;
      a1/=fdur; a2/=fdur;
      if(!isnan(a1)) {bkg1+=a1; n1++;}
      if(!isnan(a2)) {bkg2+=a2; n2++;}
    }
    bkg1/=(double)n1;
    bkg2/=(double)n2;
    if(verbose>0 || dryRun!=0) {
      printf("background_ch1[cps] := %g\n", bkg1);
      printf("background_ch2[cps] := %g\n", bkg2);
    }
  }
  if(verbose>0 || dryRun!=0) {fflush(stdout);}
 
 
  /* If dry run, then that's it */
  if(dryRun!=0) {
    tacFree(&abss); return(0);
  }
 
 
  /*
   *  Otherways, subtract the background
   */
  if(verbose>1) {
    printf("subtracting background from samples\n"); fflush(stdout);}
  if(abss.format==TAC_FORMAT_ABSS_SCANDITRONICS ||
     abss.format==TAC_FORMAT_ABSS_GEMS)
  {
    for(i=0; i<abss.sampleNr; i++) {
      fdur=abss.x2[i]-abss.x1[i];
      if(isnan(fdur) || fdur<1.0E-06) continue;
      abss.c[0].y[i]-=bkg1*fdur;
      abss.c[3].y[i]-=bkg2*fdur;
    }
  } else if(abss.format==TAC_FORMAT_ABSS_ALLOGG_OLD ||
            abss.format==TAC_FORMAT_ABSS_ALLOGG)
  {
    for(i=0; i<abss.sampleNr; i++) {
      fdur=abss.x2[i]-abss.x1[i];
      if(isnan(fdur) || fdur<1.0E-06 || isnan(abss.x[i])) continue;
      abss.c[0].y[i]-=bkg1*fdur;
      abss.c[1].y[i]-=bkg2*fdur;
    }
  }
  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);
}
/*****************************************************************************/
 
/*****************************************************************************/