tacdc.c File Reference

Decay correction for TAC data. More...

#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"

Go to the source code of this file.

Functions
int	tacGetIsotope (TAC *tac)
void	tacSetIsotope (TAC *tac, int isotope)
int	tacDecayCorrection (TAC *tac, int isotope, int mode, TPCSTATUS *status)

Detailed Description

Decay correction for TAC data.

Definition in file tacdc.c.

Function Documentation

tacDecayCorrection()

int tacDecayCorrection	(	TAC *	tac,
		int	isotope,
		int	mode,
		TPCSTATUS *	status )

Correct TAC data for physical decay, or remove the decay correction. Sample weights, if available, are modified, too.

See also

tacGetHeaderDecayCorrection, tacGetIsotope, tacGetHeaderIsotope, decayCorrectionFactorFromIsotope, tacDelay

Returns

enum tpcerror (TPCERROR_OK when successful).

Author

Vesa Oikonen

Parameters

tac	Pointer to TAC structure; status of decay correction in TAC is not verified, but set in this function. TAC must contain valid sample time unit.
isotope	Isotope code.
mode	0=Remove decay correction; 1=Correct for decay.
status	Pointer to status data; enter NULL if not needed.

Definition at line 59 of file tacdc.c.

  {
  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);
} 

Referenced by tacSetWeights().

tacGetIsotope()

int tacGetIsotope

(

TAC *

tac

)

Read isotope code from TAC header.

See also

tacGetHeaderIsotope, tacSetIsotope, isotopeIdentify, tacDecayCorrection

Author

Vesa Oikonen

Returns

Returns isotope code, or ISOTOPE_UNKNOWN, if not found.

Parameters

tac	Pointer to source TAC structure.

Definition at line 25 of file tacdc.c.

  {
  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));
}

Referenced by imgFromSIF(), tacReadModelingData(), tacReadReference(), and tacSetWeights().

tacSetIsotope()

void tacSetIsotope	(	TAC *	tac,
		int	isotope )

Write isotope code into TAC header, overwriting any previous field.

See also

tacSetHeaderIsotope, tacGetIsotope

Author

Vesa Oikonen

Parameters

tac	Pointer to target TAC structure.
isotope	Isotope code.

Definition at line 41 of file tacdc.c.

  {
  if(tac!=NULL) tacSetHeaderIsotope(&tac->h, isotopeName(isotope));
}

Referenced by imgToSIF(), tacRead4DM(), tacReadModelingData(), tacReadReference(), and tacSetWeights().