tacw.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include "tpcift.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpctac.h"
/*****************************************************************************/
 
/*****************************************************************************/
int tacIsWeighted(
  TAC *tac
) {
  if(tac==NULL) return(0);
  if(tac->weighting==WEIGHTING_UNKNOWN || tac->weighting==WEIGHTING_OFF) return(0);
  if(tac->weighting<WEIGHTING_LAST) return(1);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacWCopy(
  TAC *tac1,
  TAC *tac2,
  int i1,
  int i2
) {
  if(tac1==NULL || tac2==NULL) return(TPCERROR_FAIL);
  if(i1<0 || i1>i2) return(TPCERROR_FAIL);
  if(i2>=tac1->_sampleNr || i2>=tac2->_sampleNr) return(TPCERROR_FAIL);
  
  for(int i=i1; i<=i2; i++) {
    tac2->w[i]=tac1->w[i];
  }
  tac2->weighting=tac1->weighting;
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacWMove(
  TAC *tac,
  int ow,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(tac==NULL || tac->tacNr<1 || tac->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  int ret;
 
  if(tac->weighting==WEIGHTING_UNKNOWN) tac->weighting=WEIGHTING_OFF;
 
  /* Search for the weight column */
  int wc=-1;
  for(int i=0; i<tac->tacNr; i++) {
    if(strncasecmp(tac->c[i].name, "WEIGHT", 6)!=0) continue;
    if(wc<0) wc=i;
    else if(ow==0) {
      /* Error if found more than one */
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_DUPLICATE_DATA);
      return TPCERROR_DUPLICATE_DATA;
    }
  }
  /* Error if not found */
  if(wc<0) {
    if(ow) tac->weighting=WEIGHTING_OFF;
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_WEIGHTS);
    return TPCERROR_NO_WEIGHTS;
  }
  /* Error if overwriting not allowed */
  if(ow==0 && tacIsWeighted(tac)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_DUPLICATE_DATA);
    return TPCERROR_DUPLICATE_DATA;
  }
  /* Copy weights to their right place */
  for(int i=0; i<tac->sampleNr; i++) tac->w[i]=tac->c[wc].y[i];
  tac->weighting=WEIGHTING_ON_GENERAL;
  /* Then delete the original column */
  ret=tacDeleteTACC(tac, wc);
  statusSet(status, __func__, __FILE__, __LINE__, ret);
  return(ret);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacWByFreq(
  TAC *tac,
  isotope isot,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(%s)\n", __func__, isotopeName(isot));
  if(tac==NULL || tac->tacNr<1 || tac->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  
  tac->weighting=WEIGHTING_ON_GENERAL;
  if(tac->sampleNr==1) {tac->w[0]=1.0; return(TPCERROR_OK);}
 
  if(tac->isframe!=0) { /* weights based on frame lengths */
    tac->weighting=WEIGHTING_ON_F;
 
    for(int i=0; i<tac->sampleNr; i++) {
      tac->w[i]=tac->x2[i]-tac->x1[i]; 
    }
 
  } else { /* weights based on sample distance */
    tac->weighting=WEIGHTING_ON_F;
 
    int i, i1, i2;
    double t, t1, t2, f;
    for(i=0; i<tac->sampleNr; i++) {
      t=t1=t2=tac->x[i];
      /* Find the closest sample time before this one */
      for(i1=i; i1>=0; i1--) {t1=tac->x[i1]; if(t1<t) break;} 
      /* Find the closest sample time after this one */ 
      for(i2=i; i2<tac->sampleNr; i2++) {t2=tac->x[i2]; if(t2>t) break;} 
      /* Mean sample distance */
      f=0.0;
      if(t1<t) f+=t-t1; else f+=t2-t;
      if(t2>t) f+=t2-t; else f+=t-t1;
      f*=0.5; if(f<=0.0) f=1.0;
      tac->w[i]=f; 
    }
 
  }
 
  /* Account for decay correction, if requested */
  if(isot!=ISOTOPE_UNKNOWN) {
    tac->weighting=WEIGHTING_ON_FD;
    double lambda=lambdaFromIsotope(isot);
    if(tac->tunit==UNIT_SEC) lambda/=60.; else if(tac->tunit==UNIT_MSEC) lambda/=60000.;
    if(tac->isframe!=0) {
      for(int i=0; i<tac->sampleNr; i++)
        tac->w[i]*=decayCorrectionFactorFromLambda(-lambda, tac->x1[i], (tac->x2[i]-tac->x1[i])); 
    } else {
      for(int i=0; i<tac->sampleNr; i++)
        tac->w[i]*=decayCorrectionFactorFromLambda(-lambda, tac->x[i], 0.0); 
    }
  }
 
  /* Scale weights so that sum of weights equals sample number */
  double sumw=0.0;
  for(int i=0; i<tac->sampleNr; i++) if(isfinite(tac->w[i])) sumw+=tac->w[i];
  if(!(sumw>0.0)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
  sumw/=(double)tac->sampleNr;
  for(int i=0; i<tac->sampleNr; i++) tac->w[i]/=sumw;
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
unsigned int tacWSampleNr(
  TAC *tac
) {
  if(tac==NULL || tac->tacNr<1 || tac->sampleNr<1) return((unsigned int)0);
  if(tac->weighting==WEIGHTING_OFF) return((unsigned int)tac->sampleNr);
  unsigned int n=0;
  for(int i=0; i<tac->sampleNr; i++) if(tac->w[i]>0.0) n++;
  if(n==0 && tac->weighting==WEIGHTING_UNKNOWN) n=(unsigned int)tac->sampleNr;
  return(n);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacWeightNorm(
  TAC *tac,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(tac==NULL || tac->tacNr<1 || tac->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(!tacIsWeighted(tac)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_WEIGHTS);
    return TPCERROR_NO_WEIGHTS;
  }
 
  double wSum=0.0; int wNr=0;
  for(int i=0; i<tac->sampleNr; i++)
    if(isfinite(tac->w[i]) && tac->w[i]>0.0) {wSum+=tac->w[i]; wNr++;}
  if(wNr==0 || !(wSum>0.0)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
  double wMean=wSum/(double)wNr;
  for(int i=0; i<tac->sampleNr; i++) 
    if(isfinite(tac->w[i]) && tac->w[i]>0.0) tac->w[i]/=wMean;
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacWeightModerate(
  TAC *tac,
  const double minprop,
  const int doZeroes,
  const int doNaNs,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(%g, %d, %d)\n", __func__, minprop, doZeroes, doNaNs);
  if(tac==NULL || tac->tacNr<1 || tac->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(!isfinite(minprop) || minprop>=1.0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
  if(!tacIsWeighted(tac)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_WEIGHTS);
    return TPCERROR_NO_WEIGHTS;
  }
 
  double wMax=nan("");
  for(int i=0; i<tac->sampleNr; i++) if(isfinite(tac->w[i])) {
    if(!isfinite(wMax)) {wMax=tac->w[i]; continue;}
    if(tac->w[i]>wMax) wMax=tac->w[i];
  }
  if(verbose>2) {printf("weight_maximum := %g\n", wMax); fflush(stdout);}
  if(!isfinite(wMax) || wMax<=0.0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_WEIGHTS);
    return TPCERROR_NO_WEIGHTS;
  }
 
  double wLimit=minprop*wMax;
  if(verbose>2) {printf("weight_limit := %g\n", wLimit); fflush(stdout);}
 
  int n=0;
  if(wLimit>0.0) {
    for(int i=0; i<tac->sampleNr; i++) {
      if(tac->w[i]==0.0) {
        if(doZeroes) {tac->w[i]=wLimit; n++;}
        continue;
      }
      if(isnan(tac->w[i])) {
        if(doNaNs) {tac->w[i]=wLimit; n++;}
        continue;
      }
      if(tac->w[i]<wLimit) {tac->w[i]=wLimit; n++;}
    }
  }
  if(verbose>1 && minprop>0.0) {printf("  %d weight(s) set to %g.\n", n, wLimit); fflush(stdout);}
 
  return(tacWeightNorm(tac, status));
}
/*****************************************************************************/
 
/*****************************************************************************/
int sifWeight(
  TAC *sif,
  isotope isot,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(%s)\n", __func__, isotopeName(isot));
  if(sif==NULL || sif->tacNr<1 || sif->sampleNr<1) {
    if(sif!=NULL && verbose>2) {
      printf("  sif.tacNr := %d\n  sif.sampleNr := %d\n", sif->tacNr, sif->sampleNr);
    }
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  sif->weighting=WEIGHTING_OFF;
  if(sif->tacNr<2 || !sif->isframe) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
    return TPCERROR_UNSUPPORTED;
  }
  /* Check times */
  double xmin, xmax;
  if(tacXRange(sif, &xmin, &xmax) || !(xmax>xmin)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return TPCERROR_INVALID_XRANGE;
  }
 
  /* Calculate trues */
  double trues[sif->sampleNr];
  for(int i=0; i<sif->sampleNr; i++) trues[i]=sif->c[0].y[i]-sif->c[1].y[i];
  if(verbose>2) {
    printf("\nTrues:\n");
    for(int i=0; i<sif->sampleNr; i++) printf("\t%g\n", trues[i]);
    fflush(stdout);
  }
  /* Correct for decay, if requested */
  if(isot!=ISOTOPE_UNKNOWN) {
    for(int i=0; i<sif->sampleNr; i++)
      trues[i]*=decayCorrectionFactorFromIsotope(isot, sif->x1[i]/60., (sif->x2[i]-sif->x1[i])/60.);
    if(verbose>2) {
      printf("\nTrues after decay correction:\n");
      for(int i=0; i<sif->sampleNr; i++) printf("\t%g\n", trues[i]);
      fflush(stdout);
    }
  }
  /* Calculate weights */
  double wMax=0.0;
  for(int i=0; i<sif->sampleNr; i++) {
    double fdur=sif->x2[i]-sif->x1[i]; if(fdur<0.0) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
      return TPCERROR_INVALID_X;
    }
    if(trues[i]>0.0) {
      sif->w[i]=fdur*fdur/trues[i];
      if(sif->w[i]>wMax) wMax=sif->w[i];
    } else {
      sif->w[i]=nan("");
    }
  }
  if(!(wMax>0.0)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
    return TPCERROR_INVALID_VALUE;
  }
  if(verbose>1) {printf("\twMax := %g\n", wMax); fflush(stdout);}
  /* Set maximum weight to frames which had no true counts */
  for(int i=0; i<sif->sampleNr; i++) if(!(trues[i]>0.0) && isnan(sif->w[i])) sif->w[i]=wMax;
  if(verbose>2) {
    printf("\nWeights before normalization:\n");
    for(int i=0; i<sif->sampleNr; i++) printf("\t%g\n", sif->w[i]);
    fflush(stdout);
  }
 
  /* Normalize weights */
  sif->weighting=WEIGHTING_ON_COUNTS;
  if(tacWeightNorm(sif, status)!=TPCERROR_OK) return(status->error);
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacSetWeights(
  TAC *tac,
  weights weightMethod,
  int weightNr,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s(tac, %d, %d)\n", __func__, weightMethod, weightNr); fflush(stdout);}
 
  if(tac==NULL || tac->tacNr<1 || tac->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(weightNr<1) weightMethod=WEIGHTING_OFF;
  if(weightNr>tac->sampleNr) weightNr=tac->sampleNr;
  int origSampleNr=tac->sampleNr;
 
  /* Get isotope from TAC data */
  isotope fisot=tacGetIsotope(tac);
  if(verbose>3) {printf("  isotope :=  %s\n", isotopeName(fisot)); fflush(stdout);}
  if(fisot==ISOTOPE_UNKNOWN) {
    if(weightMethod==WEIGHTING_ON_COUNTS || weightMethod==WEIGHTING_ON_FD) {
      tac->weighting=WEIGHTING_OFF;
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNKNOWN_ISOTOPE);
      return TPCERROR_UNKNOWN_ISOTOPE;
    }
  }
 
  if(tac->format==TAC_FORMAT_SIF && weightMethod==WEIGHTING_ON_GENERAL)
    weightMethod=WEIGHTING_ON_COUNTS;
 
  if(weightMethod==WEIGHTING_OFF) {
 
    if(verbose>2 && !tacIsWeighted(tac)) {
      printf("  Note: data did not contain weights.\n"); fflush(stdout);}
    tac->weighting=WEIGHTING_OFF;
    for(int i=0; i<tac->sampleNr; i++) tac->w[i]=0.0;
    for(int i=0; i<weightNr; i++) tac->w[i]=1.0;
    if(verbose>1) {printf("  weights removed.\n"); fflush(stdout);}
 
  } else if(weightMethod==WEIGHTING_ON_GENERAL) {
 
    if(verbose>0) {fprintf(stderr, "Error: cannot apply general weighting.\n"); fflush(stderr);}
    for(int i=0; i<tac->sampleNr; i++) tac->w[i]=0.0;
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
 
  } else if(weightMethod==WEIGHTING_UNKNOWN) {
 
    for(int i=weightNr; i<tac->sampleNr; i++) tac->w[i]=0.0;
 
    if(!tacIsWeighted(tac)) {
      if(verbose>2) {printf("  Note: data does not contain weights.\n"); fflush(stdout);}
      for(int i=0; i<weightNr; i++) tac->w[i]=1.0;
    } else {
      if(verbose>2) {printf("  Note: data contains weights.\n"); fflush(stdout);}
      tac->sampleNr=weightNr;
      int ret=tacWeightNorm(tac, status);
      tac->sampleNr=origSampleNr;
      if(ret!=TPCERROR_OK) return(ret);
    }
 
  } else if(weightMethod==WEIGHTING_ON_COUNTS) {
 
    if(verbose>2 && tacIsWeighted(tac)) {
      printf("  Note: data already contained weights.\n"); fflush(stdout);}
    tac->weighting=WEIGHTING_OFF;
    for(int i=0; i<tac->sampleNr; i++) tac->w[i]=0.0;
 
    if(tac->format==TAC_FORMAT_SIF) {
 
      /* TAC contains SIF data, just add weights */
      tac->sampleNr=weightNr;
      int ret=sifWeight(tac, fisot, status);
      tac->sampleNr=origSampleNr;
      if(ret!=TPCERROR_OK) return(ret);
      tac->weighting=WEIGHTING_ON_COUNTS;
      if(verbose>1) {printf("  weights added.\n"); fflush(stdout);}
 
    } else {
 
      double weight_moderate=100.0;
 
      /* Create temp SIF */
      if(!unitIsTime(tac->tunit) || !tacIsX(tac)) {
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
        return TPCERROR_INVALID_X;
      }
      if(!unitIsRadioactivity(tac->cunit) && !unitIsRAConc(tac->cunit)) {
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INCOMPATIBLE_UNIT);
        return TPCERROR_INCOMPATIBLE_UNIT;
      }
      TAC sif; tacInit(&sif);
      if(tacAllocate(&sif, weightNr, 3)!=TPCERROR_OK) {
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INCOMPATIBLE_DATA);
        return TPCERROR_INCOMPATIBLE_DATA;
      }
      sif.sampleNr=weightNr;
      sif.tacNr=2; // 3rd contains trues that is not saved
      sif.format=TAC_FORMAT_SIF;
      sif.isframe=1; // SIF always contains frame start and end times
      sif.tunit=tac->tunit;
      sif.cunit=tac->cunit;
      tacSetIsotope(&sif, fisot);
      tacXCopy(tac, &sif, 0, sif.sampleNr-1);
      int ret=tacXUnitConvert(&sif, UNIT_SEC, status);
      if(ret!=TPCERROR_OK) {tacFree(&sif); return(ret);}
      for(int j=0; j<sif.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
      ret=tacDecayCorrection(&sif, fisot, 0, status);
      if(ret!=TPCERROR_OK) {tacFree(&sif); return(ret);}
      if(unitIsRadioactivity(sif.cunit)) {
        ret=tacYUnitConvert(&sif, UNIT_BQ, status);
      } else {
        ret=tacYUnitConvert(&sif, UNIT_BQ_PER_ML, status);
        for(int i=0; i<sif.sampleNr; i++) sif.c[0].y[i]*=1000.; // Counts from one litre
        sif.cunit=UNIT_BQ;
      }
      if(ret!=TPCERROR_OK) {tacFree(&sif); return(ret);}
      /* 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;
      }
      /* Calculate weights based on approximated SIF data */
      ret=sifWeight(&sif, fisot, status);
      if(ret==TPCERROR_OK) ret=tacWeightNorm(&sif, status);
      if(ret!=TPCERROR_OK) {tacFree(&sif); return(ret);}
 
      double prop=0.0; if(weight_moderate>1.0) prop=1.0/weight_moderate;
      ret=tacWeightModerate(&sif, prop, 0, 0, status);
      if(ret!=TPCERROR_OK) {tacFree(&sif); return(ret);}
 
      /* Copy weights into TAC */
      {
        int j;
        for(j=0; j<sif.sampleNr; j++) tac->w[j]=sif.w[j];
        for( ; j<tac->sampleNr; j++) tac->w[j]=0.0;
      }
      tacFree(&sif);
      tac->weighting=WEIGHTING_ON_COUNTS;
      if(verbose>1) {printf("  weights added.\n"); fflush(stdout);}
 
    }
 
  } else if(weightMethod==WEIGHTING_ON_F) {
 
    if(verbose>2 && tacIsWeighted(tac)) {
      printf("  Note: data already contained weights.\n"); fflush(stdout);}
    tac->weighting=WEIGHTING_OFF;
    for(int i=0; i<tac->sampleNr; i++) tac->w[i]=0.0;
    tac->sampleNr=weightNr;
    int ret=tacWByFreq(tac, ISOTOPE_UNKNOWN, status);
    if(ret==TPCERROR_OK) ret=tacWeightNorm(tac, status);
    tac->sampleNr=origSampleNr;
    if(ret!=TPCERROR_OK) return(ret);
    if(verbose>1) {printf("  weights added.\n"); fflush(stdout);}
 
  } else if(weightMethod==WEIGHTING_ON_FD) {
 
    if(verbose>2 && tacIsWeighted(tac)) {
      printf("  Note: data already contained weights.\n"); fflush(stdout);}
    tac->weighting=WEIGHTING_OFF;
    for(int i=0; i<tac->sampleNr; i++) tac->w[i]=0.0;
    tac->sampleNr=weightNr;
    int ret=tacWByFreq(tac, fisot, status);
    if(ret==TPCERROR_OK) ret=tacWeightNorm(tac, status);
    tac->sampleNr=origSampleNr;
    if(ret!=TPCERROR_OK) return(ret);
    if(verbose>1) {printf("  weights added.\n"); fflush(stdout);}
 
  } else {
    if(verbose>0) {fprintf(stderr, "Error: invalid weighting scheme.\n"); fflush(stderr);}
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/