tacdc.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 "tpcisotope.h"
#include "tpctac.h"
/*****************************************************************************/
 
/*****************************************************************************/
int tacGetIsotope(
  TAC *tac
) {
  if(tac==NULL) return(ISOTOPE_UNKNOWN);
  char buf[MAX_ISOTOPE_LEN+1];
  if(tacGetHeaderIsotope(&tac->h, buf, NULL)!=TPCERROR_OK) return(ISOTOPE_UNKNOWN);
  return(isotopeIdentify(buf));
}
/*****************************************************************************/
 
/*****************************************************************************/
void tacSetIsotope(
  TAC *tac,
  int isotope
) {
  if(tac!=NULL) tacSetHeaderIsotope(&tac->h, isotopeName(isotope));
}
/*****************************************************************************/
 
/*****************************************************************************/   
int tacDecayCorrection(
  TAC *tac,
  int isotope,
  int mode,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  /* check that data exists */
  if(tac==NULL || tac->sampleNr<1 || tac->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  /* check time unit */
  if(!unitIsTime(tac->tunit)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNKNOWN_UNIT);
    return TPCERROR_UNKNOWN_UNIT;
  }
  /* check isotope and get halflife */
  double f, lambda, halflife=isotopeHalflife(isotope); // halflife in min
  if(isotope==ISOTOPE_UNKNOWN || isnan(halflife)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNKNOWN_ISOTOPE);
    return TPCERROR_UNKNOWN_ISOTOPE;
  }  
  /* convert halflife to the same unit that the TAC data has */
  f=unitConversionFactor(UNIT_MIN, tac->tunit);
  if(!isnan(f)) halflife*=f;
  /* calculate the lambda */
  lambda=lambdaFromHalflife(halflife);
  if(mode==0) lambda=-lambda;
  if(verbose>2) {
    printf("halflife := %g\n", halflife);
    printf("lambda := %g\n", lambda);
  }
 
 
  /*
   *  Decay correction / removal
   */
  if(verbose>2) {
    if(mode!=0) printf("decay correction\n"); else printf("removing decay correction\n");
    fflush(stdout);
  }
  for(int i=0; i<tac->sampleNr; i++) {
    /* Calculate decay correction factor */
    if(tac->isframe) {
      if(isnan(tac->x1[i]) || isnan(tac->x2[i])) continue;
      f=decayCorrectionFactorFromLambda(lambda, tac->x1[i], tac->x2[i]-tac->x1[i]);
    } else {
      if(isnan(tac->x[i])) continue;
      f=decayCorrectionFactorFromLambda(lambda, tac->x[i], 0.0);
    }
    if(verbose>6) printf("  %10.4f  ->  %e\n", tac->x[i], f);
    /* Correct all TACs */
    for(int j=0; j<tac->tacNr; j++) tac->c[j].y[i]*=f;
    /* Opposite correction for weights, if available */
    if(tac->weighting==WEIGHTING_ON_GENERAL || tac->weighting==WEIGHTING_ON_COUNTS ||
       tac->weighting==WEIGHTING_ON_FD)
      tac->w[i]/=f;
  }
  if(tacIsWeighted(tac)) tacWeightNorm(tac, NULL);
  if(mode==0) tacSetHeaderDecayCorrection(&tac->h, DECAY_NOTCORRECTED);
  else tacSetHeaderDecayCorrection(&tac->h, DECAY_CORRECTED);
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
} 
/*****************************************************************************/
 
/*****************************************************************************/