tacweigh.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpccsv.h"
#include "tpcift.h"
#include "tpctac.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Add or remove the sample (time frame) weighting information to TAC file for",
  "parameter estimations and curve fitting using formula (Mazoyer et al, 1986):",
  "  weight=(frame duration)^2 / (decay corrected trues in a frame)",
  "or optionally:",
  "  weight=(frame duration),",
  "or:",
  "  weight=(frame duration)*exp(-t*ln(2)/halflife),",
  " ",
  "TACs are assumed to be corrected for decay, and SIFs not corrected.",
  "The relative weights are adjusted using a scan information file (SIF),",
  "or TACs in the file (volume weighted average of all TACs or given TAC);",
  "in the latter case the units in TAC file must be set correctly.", 
  " ",
  "Usage: @P [Options] tacfile [sif | tacname]",
  " ",
  "Options:",
  " -rm",
  "     Existing weights are removed.",
  " -list",
  "     Weights are not calculated, but existing weights are printed on screen.",
  "     Return code is non-zero if data does not contain weights.",
  " -wf | -wfd",
  "     With -wf Weights are based only on frame length or sampling interval.",
  "     With -wfd the late frames are given less weight by using formula:",
  "     weight=(frame duration)*exp(-t*ln(2)/halflife) (Thiele et al, 2008).",
  " -i=<Isotope code>",
  "     Isotope, for example C-11, in case it is not found inside SIF or TAC",
  "     file. Isotope is only needed with SIF, and with option -wfd.",
  " -moderate=<ratio>",
  "     Weights are moderated by limiting the ratio of maximum and minimum",
  "     weights to the given ratio (100 by default). Weights that are lower",
  "     than maximum/ratio are set to maximum/ratio. Enter zero to not",
  "     moderate the weights. Otherwise ratio must be larger than one.",
  " -sif=<Filename>",
  "     SIF data based on TAC file is written in given file; cannot be used",
  "     if SIF is given as argument.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Note that absolute weights cannot be calculated. Relative weights are",
  "scaled so that average weight is 1.0 for frames that have weight above 0.",
  " ",
  "Reference:",
  "1. Mazoyer BM, Huesman RH, Budinger TF, Knittel BL. Dynamic PET data",
  "   analysis. J Comput Assist Tomogr 1986; 10:645-653.",
  "2. Thiele F, Buchert R. Evaluation of non-uniform weighting in non-linear", 
  "   regression for pharmacokinetic neuroreceptor modelling.", 
  "   Nucl Med Commun. 2008; 29:179-188.",
  " ",
  "See also: sifcat, sifisot, taclist, imgweigh, tacframe, imghead, tacdecay",
  " ",
  "Keywords: TAC, SIF, modelling, weighting, fitting",
  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     siffile[FILENAME_MAX], tacfile[FILENAME_MAX], newsif[FILENAME_MAX];
  double   weight_moderate=100.0;
  int      weight_method=WEIGHTING_ON_COUNTS;
  int      remove_weights=0;
  int      print_weights=0;
  isotope  isot=ISOTOPE_UNKNOWN;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  siffile[0]=tacfile[0]=newsif[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    if(strncasecmp(cptr, "LIST", 1)==0) {
      print_weights=1; continue;
    } else if(strcasecmp(cptr, "RM")==0 || strcasecmp(cptr, "DEL")==0) {
      remove_weights=1; continue;
    } else if(strncasecmp(cptr, "I=", 2)==0) {
      if((isot=isotopeIdentify(cptr+2))==ISOTOPE_UNKNOWN) {
        fprintf(stderr, "Error: invalid isotope '%s'.\n", cptr+2); return(1);}
      continue;
    } else if(strncasecmp(cptr, "SIF=", 4)==0) {
      cptr+=4; strlcpy(newsif, cptr, FILENAME_MAX); if(strlen(newsif)>0) continue;
    } else if(strncasecmp(cptr, "moderate=", 9)==0) {
      cptr+=9; weight_moderate=atofVerified(cptr); 
      if(weight_moderate>1.0) continue;
      if(weight_moderate<=0.0) continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weight_method=WEIGHTING_ON_F; continue;
    } else if(strcasecmp(cptr, "WFD")==0) {
      weight_method=WEIGHTING_ON_FD; continue;
    }
    fprintf(stderr, "Error: invalid option '%s'\n", argv[ai]);
    return(1);
  } else break; // tac name argument may start with '-'
  
  /* 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);}
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-1;
 
  /* Process other arguments, starting from the first non-option */
  if(ai<argc) {strlcpy(tacfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {strlcpy(siffile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  if(siffile[0] && (weight_method==WEIGHTING_ON_F || weight_method==WEIGHTING_ON_FD)) {
    fprintf(stderr, "Warning: using options -wf and -wfd with SIF is not recommended.\n");
  }
 
  /* Is something missing? */
  if(!tacfile[0]) {
    fprintf(stderr, "Error: missing file name; use option --help\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile := %s\n", tacfile);
    //printf("siffile := %s\n", siffile); // may actually be TAC name
    if(newsif[0]) printf("newsif := %s\n", newsif);
    printf("print_weights := %d\n", print_weights);
    printf("remove_weights := %d\n", remove_weights);
    printf("weight_method := %d\n", weight_method);
    if(isot!=ISOTOPE_UNKNOWN) printf("isotope := %s\n", isotopeName(isot));
    printf("moderate := %g\n", weight_moderate);
    fflush(stdout);
  }
 
 
 
  /*
   *  Read TAC data
   */
  if(verbose>1) {fprintf(stdout, "reading %s\n", tacfile); fflush(stdout);}
  TAC tac; tacInit(&tac);
  if(tacRead(&tac, tacfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), tacfile);
    tacFree(&tac); return(2);
  }
  if(verbose>2) {
    printf("weighting := %d\n", tac.weighting);
    printf("fileformat := %s\n", tacFormattxt(tac.format));
    printf("tacNr := %d\n", tac.tacNr);
    printf("sampleNr := %d\n", tac.sampleNr);
    printf("xunit := %s\n", unitName(tac.tunit));
    printf("yunit := %s\n", unitName(tac.cunit));
    if(tac.isframe) printf("frames := yes\n"); else printf("frames := no\n");
    fflush(stdout);
  }
  if(verbose>2) {
    printf("contains_weights := ");
    if(tacIsWeighted(&tac)) printf("yes\n"); else printf("no\n");
    fflush(stdout);
  }
 
 
  /*
   *  Print existing weights, if required, and quit
   */
  if(print_weights!=0) {
    if(verbose>2) printf("printing existing weights\n");
    if(!tacIsWeighted(&tac)) {
      printf("contains_weights := no\n");
      //printf("%s does not contain weights.\n", tacfile);
      tacFree(&tac); return(100);
    } else if(tac.isframe) {
      printf("start[]\tend[]\tweight\n");
      for(int i=0; i<tac.sampleNr; i++) printf("%g\t%g\t%.4e\n", tac.x1[i], tac.x2[i], tac.w[i]);
      tacFree(&tac); return(0);
    } else {
      printf("time[]\tweight\n");
      for(int i=0; i<tac.sampleNr; i++) printf("%g\t%.4e\n", tac.x[i], tac.w[i]);
      tacFree(&tac); return(0);
    }
  }
 
  /* Check if datafile contains information on decay correction */
  decaycorrection fdc=tacGetHeaderDecayCorrection(&tac.h);
  if(fdc==DECAY_UNKNOWN) {
    if(verbose>1) printf("Note: status of current decay correction is not known.\n");
  } else if(fdc==DECAY_CORRECTED) {
    if(verbose>1) fprintf(stderr, "Note: physical decay is corrected in %s.\n", tacfile);
  } else if(fdc==DECAY_NOTCORRECTED) {
    fprintf(stderr, "Note: physical decay is not corrected in %s.\n", tacfile);
  }
 
  /* If datafile contains valid isotope, then check that it is the same as given by user. */
  /* Not all methods need isotope. */
  /* No problem if isotope is not yet known, because it may be in SIF */
  {
    isotope fisot=tacGetIsotope(&tac);
    if(verbose>3) {printf("tac.isotope :=  %s\n", isotopeName(fisot)); fflush(stdout);}
    if(isot==ISOTOPE_UNKNOWN && fisot==ISOTOPE_UNKNOWN) {
      if(!siffile[0] && (weight_method==WEIGHTING_ON_COUNTS || weight_method==WEIGHTING_ON_FD)) {
        fprintf(stderr, "Error: valid isotope not specified.\n");
        tacFree(&tac); return(1);
      }
    } else if(isot==ISOTOPE_UNKNOWN && fisot!=ISOTOPE_UNKNOWN) {
      isot=fisot;
      if(verbose>1) printf("isotope := %s\n", isotopeName(isot));
    } else if(isot!=ISOTOPE_UNKNOWN && fisot==ISOTOPE_UNKNOWN) {
      fisot=isot;
      tacSetIsotope(&tac, fisot);
      //if(verbose>1) printf("isotope := %s\n", isotopeName(isot));
    }
    if(isot!=fisot && (weight_method==WEIGHTING_ON_COUNTS || weight_method==WEIGHTING_ON_FD)) {
      fprintf(stderr, "Error: different isotope in %s\n", tacfile);
      tacFree(&tac); return(3);
    }
  }
 
 
 
  /*
   *  Remove the stored weights, if requested
   */
  if(remove_weights!=0) {
    if(verbose>1) printf("removing weights\n");
    if(!tacIsWeighted(&tac)) {
      printf("%s does not contain weights.\n", tacfile);
      tacFree(&tac); return(0);
    }
    tac.weighting=WEIGHTING_OFF;
    if(verbose>1) printf("saving modified %s\n", tacfile);
    FILE *fp; fp=fopen(tacfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing\n");
      tacFree(&tac); return(11);
    }
    int ret=tacWrite(&tac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); return(12);
    }
    /* Tell user what we did */
    if(verbose>0) printf("Weights removed in %s\n", tacfile);
    tacFree(&tac); return(0);
  }
 
 
  /*
   *  Try to read SIF or identify specified 'head' TAC
   */
  TAC sif; tacInit(&sif);
  int headIndex=-1;
  if(siffile[0]) { /* User gave either SIF or TAC name */
    /* Check if this is SIF file */
    if(verbose>1) {fprintf(stdout, "trying to read %s\n", siffile); fflush(stdout);}
    int ret=tacRead(&sif, siffile, &status);
    if(ret==TPCERROR_OK) {
      if(verbose>2) {
        printf("sif.weighting := %d\n", sif.weighting);
        printf("sif.fileformat := %s\n", tacFormattxt(sif.format));
        printf("sif.tacNr := %d\n", sif.tacNr);
        printf("sif.sampleNr := %d\n", sif.sampleNr);
        printf("sif.xunit := %s\n", unitName(sif.tunit));
        printf("sif.yunit := %s\n", unitName(sif.cunit));
        if(sif.isframe) printf("sif.frames := yes\n"); else printf("sif.frames := no\n");
        fflush(stdout);
      }
      if(sif.format!=TAC_FORMAT_SIF) {
        fprintf(stderr, "Error: file is not in SIF format.\n");
        tacFree(&tac); tacFree(&sif); return(4);
      }
      if(verbose>1) printf("%s is SIF\n", siffile);
      /* tacReadSIF() adds column for trues but does not include it in tacNr */
      if(sif.tacNr!=2 || sif._tacNr!=3) {
        fprintf(stderr, "Error: invalid SIF.\n");
        tacFree(&tac); tacFree(&sif); return(4);
      }
      double promptMax, randomMax;
      if(tacYRange(&sif, 0, NULL, &promptMax, NULL, NULL, NULL, NULL) ||
         tacYRange(&sif, 1, NULL, &randomMax, NULL, NULL, NULL, NULL))
      {
        fprintf(stderr, "Error: invalid SIF.\n");
        tacFree(&tac); tacFree(&sif); return(4);
      }
      if(verbose>2) printf("  promptMax := %g\n  randomMax := %g\n", promptMax, randomMax);
      if(!(promptMax>0.0)) {
        if(weight_method==WEIGHTING_ON_COUNTS) {
          fprintf(stderr, "Error: SIF does not contain prompts.\n");
          tacFree(&tac); tacFree(&sif); return(4);
        }
        fprintf(stderr, "Warning: SIF does not contain prompts.\n"); fflush(stderr);
      }
      if(!(randomMax>0.0)) {
        fprintf(stderr, "Warning: SIF does not contain randoms.\n"); fflush(stderr);
      }
    } else { // siffile is not an existing file; check if it is a TAC name
      if(verbose>1) fprintf(stdout, "trying to find TAC %s\n", siffile);
      int n=tacSelectTACs(&tac, siffile, 1, &status);
      if(n<1) {
        fprintf(stderr, "Error: no TAC matching '%s' was found.\n", siffile);
        tacFree(&tac); return(4);
      }
      if(n>1) {
        if(verbose>1) printf(" %d TACs match '%s'.\n", n, siffile);
        headIndex=tacSelectBestReference(&tac);
      } else {
        for(int i=0; i<tac.tacNr; i++) if(tac.c[i].sw) {headIndex=i; break;}
      }
      if(headIndex>=0) {
        if(verbose>0) printf("headTAC := %s\n", tac.c[headIndex].name);
      } else {
        fprintf(stderr, "Error: no single TAC matching '%s' was found.\n", siffile);
        tacFree(&tac); return(4);
      }
      strcpy(siffile, "");
    }
  }
 
 
  /*
   *  If we have a SIF, then try to find isotope inside it,
   *  and check that time frames match the TAC data.
   */
  if(sif.sampleNr>0) {
    isotope fisot=tacGetIsotope(&sif);
    if(verbose>3) {printf("sif.isotope :=  %s\n", isotopeName(fisot)); fflush(stdout);}
    if(isot==ISOTOPE_UNKNOWN && fisot==ISOTOPE_UNKNOWN) {
      fprintf(stderr, "Error: valid isotope not specified.\n");
      tacFree(&tac); tacFree(&sif); return(3);
    }
    if(isot==ISOTOPE_UNKNOWN && fisot!=ISOTOPE_UNKNOWN) {
      isot=fisot;
      tacSetIsotope(&tac, fisot);
      if(verbose>1) printf("isotope := %s\n", isotopeName(isot));
    } else if(isot!=ISOTOPE_UNKNOWN && fisot==ISOTOPE_UNKNOWN) {
      fisot=isot;
      tacSetIsotope(&sif, fisot);
    }
    if(isot!=fisot) {
      fprintf(stderr, "Error: different isotope in %s\n", siffile);
      tacFree(&tac); tacFree(&sif); return(4);
    }
    if(isot==ISOTOPE_UNKNOWN && (weight_method==WEIGHTING_ON_COUNTS || weight_method==WEIGHTING_ON_FD))
    {
      fprintf(stderr, "Error: isotope is required with SIF and selected weighting method.\n");
      tacFree(&tac); tacFree(&sif); return(4);
    }
    if(sif.sampleNr!=tac.sampleNr || !tacXMatch(&tac, &sif, verbose-5)) {
      fprintf(stderr, "Error: frame times in TAC and SIF do not match.\n");
      tacFree(&tac); tacFree(&sif); return(4);
    }
  }
 
  /* If user wants to save SIF, it must be saved with no decay decay correction, and therefore
     we must know the isotope */
  if(newsif[0] && isot==ISOTOPE_UNKNOWN) {
    fprintf(stderr, "Error: isotope is required for making SIF.\n");
    tacFree(&tac); tacFree(&sif); return(1);
  }
 
 
 
  /*
   *  If SIF was not given, then create SIF data from regional data.
   *  Count-related data is computed here for the SIF even if user just
   *  wants the weights to be based on frame lengths.
   */
  if(!siffile[0]) {
    if(verbose>1) printf("creating SIF data from regional data\n");
    /* Make sure that TAC time units are known */
    if(!unitIsTime(tac.tunit)) {
      /* Try to guess the unit */
      double xmin, xmax;
      if(tacXRange(&tac, &xmin, &xmax)) {
        fprintf(stderr, "Error: invalid sample times.\n");
        tacFree(&tac); return(2);
      }
      if(xmax<=60.0) {
        tac.tunit=UNIT_MIN;
        fprintf(stderr, "Warning: unknown time unit; assumed min.\n");
      } else if(xmax>360.0) {
        tac.tunit=UNIT_SEC;
        fprintf(stderr, "Warning: unknown time unit; assumed sec.\n");
      } else {
        fprintf(stderr, "Error: unknown time unit.\n");
        tacFree(&tac); return(2);
      }
    }
    /* Make sure that TAC calibration units are known, if we need them */
    if(weight_method==WEIGHTING_ON_COUNTS) {
      if(tac.cunit==UNIT_UNKNOWN) {
        fprintf(stderr, "Error: unknown calibration unit.\n");
        tacFree(&tac); return(2);
      }
      if(!unitIsRadioactivity(tac.cunit) && !unitIsRAConc(tac.cunit)) {
        fprintf(stderr, "Error: invalid unit %s.\n", unitName(tac.cunit));
        tacFree(&tac); return(2);
      }
    } else if(newsif[0]) { // needed for making SIF
      if(tac.cunit==UNIT_UNKNOWN) {
        fprintf(stderr, "Error: unknown calibration unit.\n");
        tacFree(&tac); return(2);
      }
      if(!unitIsRadioactivity(tac.cunit) && !unitIsRAConc(tac.cunit)) {
        fprintf(stderr, "Error: invalid unit %s.\n", unitName(tac.cunit));
        tacFree(&tac); return(2);
      }
    }
    /* Allocate memory for SIF data */
    if(tacAllocate(&sif, tac.sampleNr, 3)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot prepare data for SIF.\n");
      tacFree(&tac); tacFree(&sif); return(6);
    }
    sif.sampleNr=tac.sampleNr;
    sif.tacNr=2; // 3rd contains trues that is not saved
    /* Set SIF compatible titles into TAC struct */
    strcpy(sif.c[0].name, "Prompts");
    strcpy(sif.c[1].name, "Randoms");
    strcpy(sif.c[2].name, "Trues"); // computed, not saved
    /* Set basic header information */
    sif.format=TAC_FORMAT_SIF;
    sif.isframe=tac.isframe; // For now, but SIF always contains frame start and end times
    sif.tunit=tac.tunit;
    sif.cunit=tac.cunit;
    tacSetIsotope(&sif, isot);
    /* Study number is needed when SIF is saved */
    {
      char studynr[MAX_STUDYNR_LEN+1];
      if(tacGetHeaderStudynr(&tac.h, studynr, &status)==TPCERROR_OK)
        tacSetHeaderStudynr(&sif.h, studynr);
    }
    /* Copy times, and convert to seconds */
    tacXCopy(&tac, &sif, 0, tac.sampleNr-1);
    if(tacXUnitConvert(&sif, UNIT_SEC, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&sif); return(6);
    }
    /* Guess frame start and end times, if not present in data */
    if(tacSetX(&sif, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&sif); return(6);
    }
    sif.isframe=1; // SIF always contains frame start and end times
    /* Estimate SIF counts */
    if(headIndex>=0) {
      /* If 'head' region name was given, then copy it to sif */
      for(int i=0; i<tac.sampleNr; i++) sif.c[0].y[i]=tac.c[headIndex].y[i];
    } else {
      /* otherwise we must calculate the size weighted average of all TACs */
      for(int j=0; j<tac.sampleNr; j++) {
        double sum=0.0, wsum=0.0;
        for(int i=0; i<tac.tacNr; i++) {
          if(!isfinite(tac.c[i].y[j])) continue;
          double w=tac.c[i].size; if(!(w>0.0)) w=1.0;
          sum+=w*tac.c[i].y[j]; wsum+=w;
        }
        if(wsum>0.0) sif.c[0].y[j]=sum/wsum; else sif.c[0].y[j]=0.0;
      } // next time frame
    }
    /* Remove decay correction, if isotope is known. If isotope is known, 'trues' are decay
       corrected for weight calculation in sifWeight(), although values in SIF will not 
       be changed. */
    if(isot!=ISOTOPE_UNKNOWN) {
      if(tacDecayCorrection(&sif, isot, 0, &status)!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s\n", errorMsg(status.error));
        tacFree(&tac); tacFree(&sif); return(6);
      }
    }
    if(unitIsRadioactivity(sif.cunit)) {
      if(tacYUnitConvert(&sif, UNIT_BQ, &status)!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s\n", errorMsg(status.error));
        tacFree(&tac); tacFree(&sif); return(6);
      }
    } else if(unitIsRAConc(sif.cunit)) {
      if(tacYUnitConvert(&sif, UNIT_BQ_PER_ML, &status)!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s\n", errorMsg(status.error));
        tacFree(&tac); tacFree(&sif); return(6);
      }
      for(int i=0; i<sif.sampleNr; i++) sif.c[0].y[i]*=1000.; // Counts from one litre
      sif.cunit=UNIT_BQ;
    }
    /* Multiply with frame duration (Bq/s -> Bq/frame) */
    for(int i=0; i<sif.sampleNr; i++) sif.c[0].y[i]*=(sif.x2[i]-sif.x1[i]);
    sif.cunit=UNIT_COUNTS;
    /* Set prompts, randoms, and trues */
    for(int i=0; i<sif.sampleNr; i++) {
      sif.c[2].y[i]=sif.c[0].y[i];
      sif.c[1].y[i]=0.0;
    }
    /* Save the SIF, if requested */
    if(newsif[0]) {
      if(verbose>1) printf("writing %s\n", newsif);
      FILE *fp; fp=fopen(newsif, "w");
      if(fp==NULL) {
        fprintf(stderr, "Error: cannot open file for writing (%s)\n", newsif);
        tacFree(&tac); tacFree(&sif); return(11);
      }
      if(tacWriteSIF(&sif, fp, 0, &status)!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s\n", errorMsg(status.error));
        fclose(fp); tacFree(&tac); tacFree(&sif); return(12);
      }
      fclose(fp);
      if(verbose>=0) printf("%s saved.\n", newsif);
    }
  }
 
 
  /*
   *  Calculate weights into SIF
   */
  if(verbose>1) {printf("calculating weights from SIF data\n"); fflush(stdout);}
  if(weight_method==WEIGHTING_ON_COUNTS) {
    if(sifWeight(&sif, isot, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot calculate weights.\n");
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&sif); return(4);
    }
  } else if(weight_method==WEIGHTING_ON_F) {
    if(tacWByFreq(&sif, ISOTOPE_UNKNOWN, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot calculate weights.\n");
      fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), tacfile);
      tacFree(&tac); tacFree(&sif); return(4);
    }
  } else if(weight_method==WEIGHTING_ON_FD) {
    if(tacWByFreq(&sif, isot, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot calculate weights.\n");
      fprintf(stderr, "Error: %s (%s)\n", errorMsg(status.error), tacfile);
      tacFree(&tac); tacFree(&sif); return(4);
    }
  } else {
    fprintf(stderr, "Error: invalid settings.\n");
    tacFree(&tac); tacFree(&sif); return(4);
  }
  /* Moderate weights, if required */
  if(weight_moderate>0.0) {
    if(verbose>1) {printf("moderate and normalize weights\n"); fflush(stdout);}
    double prop=0.0; if(weight_moderate>1.0) prop=1.0/weight_moderate;
    if(tacWeightModerate(&sif, prop, 0, 0, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&sif); return(8);
    }
  } else {
    /* Just normalize weights */
    if(verbose>1) {printf("normalize weights\n"); fflush(stdout);}
    if(tacWeightNorm(&sif, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&sif); return(9);
    }
  }
  /* Copying weights into TAC */
  tacWCopy(&sif, &tac, 0, tac.sampleNr-1);
  tacFree(&sif);
 
  /*
   *  Save TAC file with weights
   */
  {
    /* Some file formats cannot save weights */
    if(tac.format==TAC_FORMAT_SIMPLE) tac.format=TAC_FORMAT_PMOD;
    if(verbose>1) printf("saving modified %s\n", tacfile);
    FILE *fp; fp=fopen(tacfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing\n");
      tacFree(&tac); return(21);
    }
    int ret=tacWrite(&tac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); return(22);
    }
    /* Tell user what we did */
    if(verbose>0) printf("Weights added in %s\n", tacfile);
  }
 
  tacFree(&tac);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/