delay.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpctac.h"
#include "tpcli.h"
#include "tpclinopt.h"
/*****************************************************************************/
#include "tpctacmod.h"
/*****************************************************************************/
 
/*****************************************************************************/ 
int tacDelay(
  TAC *tac,
  double dt,
  int ti,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(tac, %g, %d)\n", __func__, dt, ti);
  /* Check the function input */
  if(tac==NULL || tac->sampleNr<1 || tac->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(ti>=tac->tacNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(fabs(dt)<1.0E-12) {
    if(verbose>1) printf("  requested delay time is zero; nothing done.\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return(TPCERROR_OK);
  }
  /* Check that input does not have any missing values */
  if(tacXNaNs(tac)>0 || tacYNaNs(tac, ti)>0) {
    if(verbose>1) printf("  data contains NaN(s).\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
    return TPCERROR_MISSING_DATA;
  } 
 
  /* Get and check the x range */
  double xmin=nan(""), xmax=nan("");
  if(tacXRange(tac, &xmin, &xmax)!=0 || dt>=xmax || dt<=-xmax) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return TPCERROR_INVALID_XRANGE;
  }
  if(verbose>2) printf("  xrange: %g - %g\n", xmin, xmax);
 
 
  /* If TAC contains just frame mid times, then simple dot-to-dot interpolation */
  if(tac->isframe==0) {
    /* Make temporary x array */
    double tx[tac->sampleNr];
    for(int i=0; i<tac->sampleNr; i++) tx[i]=tac->x[i]+dt;
    /* One TAC at a time */
    for(int j=0; j<tac->tacNr; j++) if(ti<0 || ti==j) {
      /* Make temporary y array */
      double ty[tac->sampleNr];
      for(int i=0; i<tac->sampleNr; i++) ty[i]=tac->c[j].y[i];
      /* Interpolate temporary data to original times */
      if(liInterpolate(tx, ty, tac->sampleNr, tac->x, tac->c[j].y, NULL, NULL, tac->sampleNr, 0, 1, 
                       verbose-10)!=0)
      {
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
        return TPCERROR_INVALID_VALUE;
      }
    }
    /* That's it then */
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return(TPCERROR_OK);
  }
 
  /* We have frame start and end times */
  /* Convert the data into stepwise dot-to-dot data */
  TAC stac; tacInit(&stac);
  {
    int ret=tacFramesToSteps(tac, &stac, NULL);
    if(ret!=TPCERROR_OK) {
      tacFree(&stac);
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(ret);
    }
  }
  /* Change times of the stepwise data */
  for(int i=0; i<stac.sampleNr; i++) stac.x[i]+=dt;
  /* Interpolate stepwise data into original frame times, one TAC at a time */
  for(int j=0; j<tac->tacNr; j++) if(ti<0 || ti==j) {
    if(liInterpolateForPET(stac.x, stac.c[j].y, stac.sampleNr, tac->x1, tac->x2, tac->c[j].y, 
                           NULL, NULL, tac->sampleNr, 0, 1, verbose-10))
    {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return TPCERROR_INVALID_VALUE;
    }
  }
  tacFree(&stac);
 
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
void initDelayFitData(
  DELAYFITDATA *d
) {
  if(d==NULL) return;
  d->iNr=d->tNr=d->moveNr=0;
  tacInit(&d->dtac); tacInit(&d->ditac); tacInit(&d->diitac); tacInit(&d->ddtac); tacInit(&d->dddtac);
  d->mode=0; d->model=0;
  d->llsq_n=d->llsq_m=0;
  d->llsq_mat=NULL;
  d->llsq_a=NULL;
  d->llsq_b=d->llsq_x=d->llsq_wp=d->llsq_zz=NULL; 
  d->llsq_i=NULL;
  return;
}
/*****************************************************************************/
void freeDelayFitData(
  DELAYFITDATA *d
) {
  if(d==NULL) return;
  d->iNr=d->tNr=d->moveNr=0;
  tacFree(&d->dtac); tacFree(&d->ditac); tacFree(&d->diitac); tacFree(&d->ddtac); tacFree(&d->dddtac);
  d->mode=0; d->model=0;
  d->llsq_n=d->llsq_m=0;
  free(d->llsq_a); free(d->llsq_mat);
  free(d->llsq_b); free(d->llsq_x); free(d->llsq_wp); free(d->llsq_zz); 
  free(d->llsq_i);
  return;
}
/*****************************************************************************/
 
/*****************************************************************************/ 
int tacDelayFit(
  TAC *btac,
  TAC *ttac,
  int ci,
  double dtmin,
  double dtmax,
  double fitend,
  double dtstep,
  double *dt,
  int mode, 
  int model, 
  DELAYFITDATA *tdata,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>2) {
    printf("%s(inp, tis, %d, %g, %g, %g, %g, %d, ...)\n", __func__, ci, dtmin, dtmax, fitend, dtstep, mode);
    fflush(stdout);
  }
 
 
  /* Check provided data */
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
  if(ttac==NULL || (btac==NULL && tdata==NULL)) return(TPCERROR_FAIL);
  if(btac!=NULL && (btac->tacNr<1 || btac->sampleNr<4)) return(TPCERROR_FAIL);
  if(ci<0 || ci>=ttac->tacNr || ttac->sampleNr<4) return(TPCERROR_FAIL);
  if(model<0 || model>2) return(TPCERROR_FAIL);
  if(dt!=NULL) *dt=nan("");
 
  /* Initiate (possibly reused) data structure */
  DELAYFITDATA *b, ltdata; initDelayFitData(&ltdata);
  if(tdata!=NULL) {
    b=tdata;
    if(verbose>3) {
      if(b->moveNr==0) printf("  given temp data will be filled\n"); 
      else printf("  filled temp data will be used\n");
      fflush(stdout);
    }
  } else {
    b=&ltdata;
    if(verbose>4) {printf("  local temp data created\n"); fflush(stdout);}
  }
 
  /* Check and set temp data when run first time with the given TTAC */
  if(btac!=NULL || b->moveNr==0) {
    /* Delay range and steps*/
    if(!(dtmax>dtmin) || !(dtstep>0.001)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
    b->dtmin=dtmin; b->dtmax=dtmax; b->dtstep=dtstep;
    /* How many delayed curves needed */
    b->moveNr=1+(b->dtmax-b->dtmin)/b->dtstep;
    if(verbose>3) printf("  moveNr := %d\n", b->moveNr);
    if(b->moveNr<5 || b->moveNr>5000) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
    /* Get x range */
    if(tacXRange(ttac, &b->ttac_xmin, &b->ttac_xmax)!=0 || tacXRange(btac, &b->btac_xmin, &b->btac_xmax)!=0) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
    if(verbose>3) printf("  TTAC range := %g - %g\n  BTAC range := %g - %g\n",
                         b->ttac_xmin, b->btac_xmax, b->ttac_xmin, b->btac_xmax);
    /* Check that time units do match */
    if(btac->tunit!=ttac->tunit) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
      return(TPCERROR_INVALID_X);
    }
    /* Check that input curve is long enough to fit duration */
    if(!(fitend>0.0)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
    b->fitend=fitend;
    if(b->fitend>b->ttac_xmax) b->fitend=b->ttac_xmax;
    if(b->fitend>b->btac_xmax) b->fitend=b->btac_xmax;
    if(b->dtmin<0.0 && b->btac_xmax<(b->fitend+b->dtmin))
      b->fitend+=b->dtmin;
    if(verbose>3) printf("fitend=%g -> %g\n", fitend, b->fitend);
    b->iNr=0; for(unsigned int i=0; i<(unsigned int)btac->sampleNr; i++) if(btac->x[i]<=b->fitend) b->iNr++; else break;
    b->tNr=0; for(unsigned int i=0; i<(unsigned int)ttac->sampleNr; i++) if(ttac->x[i]<=b->fitend) b->tNr++; else break;
    //printf("iNr=%u tNr=%u\n", b->iNr, b->tNr);
    if(b->iNr<5 || b->tNr<5) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
 
    /* Further, when run the first time, make input TACs moved in time, interpolated and integrated 
       to tissue times */
    if(verbose>5) {printf("  preparing delayed input data\n"); fflush(stdout);}
    b->mode=mode;
    b->model=model;
    int ret;
    ret=tacAllocate(&b->dtac, b->tNr, b->moveNr+1); // additional place for tissue curve
    if(ret!=TPCERROR_OK) {
      freeDelayFitData(&ltdata);
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(TPCERROR_OUT_OF_MEMORY);
    }
    b->dtac.sampleNr=b->tNr;
    b->dtac.tacNr=b->moveNr+1;  // tissue as the last curve
    b->dtac.isframe=ttac->isframe;
    for(int i=0; i<b->dtac.sampleNr; i++) {
      b->dtac.x1[i]=ttac->x1[i]; b->dtac.x[i]=ttac->x[i]; b->dtac.x2[i]=ttac->x2[i];}
    ret=tacDuplicate(&b->dtac, &b->ditac);
    if(ret==TPCERROR_OK) ret=tacDuplicate(&b->dtac, &b->diitac); 
    if(ret==TPCERROR_OK) ret=tacDuplicate(&b->dtac, &b->ddtac); 
    if(ret==TPCERROR_OK) ret=tacDuplicate(&b->dtac, &b->dddtac); 
    if(ret!=TPCERROR_OK) {
      freeDelayFitData(&ltdata);
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(ret);
    }
    /* Integrate/derivate input data, change times, and interpolate to tissue data */
    double bd[btac->sampleNr], bdd[btac->sampleNr];
    double bi[btac->sampleNr], bii[btac->sampleNr];
    if(liDerivate3(btac->x, btac->c[0].y, btac->sampleNr, bd, bdd, 0)) {
      freeDelayFitData(&ltdata);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_VALUE);
      return(TPCERROR_INVALID_VALUE);
    }
    if(btac->isframe==0) {
      if(liIntegrate(btac->x, btac->c[0].y, btac->sampleNr, bi, 3, 0) ||
         liIntegrate(btac->x, bi, btac->sampleNr, bii, 3, 0))
        ret=TPCERROR_INVALID_VALUE;
    } else {
      if(liIntegrateFE(btac->x1, btac->x2, btac->c[0].y, btac->sampleNr, bi, bii, 0))
        ret=TPCERROR_INVALID_VALUE;
    }
    if(ret) {
      freeDelayFitData(&ltdata);
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(ret);
    }
    double dx[btac->sampleNr], dx2[btac->sampleNr], *dxp;
    if(btac->isframe==0) dxp=dx; else dxp=dx2;
    for(int di=0; di<b->dtac.tacNr-1; di++) {
      for(int i=0; i<btac->sampleNr; i++) {
        dx[i]=btac->x[i]+b->dtmin+di*b->dtstep;
        dx2[i]=btac->x2[i]+b->dtmin+di*b->dtstep;
      }
      if(b->dtac.isframe==0) {
        ret=liInterpolate(dx, btac->c[0].y, btac->sampleNr, b->dtac.x, b->dtac.c[di].y, 
                            NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
    
        if(!ret) ret=liInterpolate(dxp, bi, btac->sampleNr, b->dtac.x, b->ditac.c[di].y, 
                            NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolate(dxp, bii, btac->sampleNr, b->dtac.x, b->diitac.c[di].y, 
                            NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolate(dx, bd, btac->sampleNr, b->dtac.x, b->ddtac.c[di].y, 
                            NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolate(dx, bdd, btac->sampleNr, b->dtac.x, b->dddtac.c[di].y, 
                            NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
      } else {
        ret=liInterpolateForPET(dx, btac->c[0].y, btac->sampleNr, b->dtac.x1, b->dtac.x2, 
                            b->dtac.c[di].y, NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolateForPET(dxp, bi, btac->sampleNr, b->dtac.x1, b->dtac.x2, 
                            b->ditac.c[di].y, NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolateForPET(dxp, bii, btac->sampleNr, b->dtac.x1, b->dtac.x2, 
                            b->diitac.c[di].y, NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolateForPET(dx, bd, btac->sampleNr, b->dtac.x1, b->dtac.x2, 
                            b->ddtac.c[di].y, NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
        if(!ret) ret=liInterpolateForPET(dx, bdd, btac->sampleNr, b->dtac.x1, b->dtac.x2, 
                            b->dddtac.c[di].y, NULL, NULL, b->dtac.sampleNr, 3, 1, 0);
      }
      if(ret) {
        freeDelayFitData(&ltdata);
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
        return(TPCERROR_INVALID_XRANGE);
      }
    }
    if(verbose>6) {
      printf("integrated input curve 1:\n");
      for(int i=0; i<b->dtac.sampleNr; i++) 
        printf("%g  %g %g %g\n", b->dtac.x[i], b->dtac.c[0].y[i], b->ditac.c[0].y[i], b->diitac.c[0].y[i]);
      fflush(stdout);
    }
 
 
    /* Allocate memory for linear LSQ */
    if(verbose>5) {printf("  preparing data for LLSQ\n"); fflush(stdout);}
    b->llsq_n=3+b->model; 
    b->llsq_m=b->dtac.sampleNr;
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    b->llsq_mat=(double*)malloc((b->llsq_n*b->llsq_m)*sizeof(double));
    if(b->llsq_mat==NULL) {freeDelayFitData(&ltdata); return(TPCERROR_OUT_OF_MEMORY);}
    b->llsq_a=(double**)malloc(b->llsq_n*sizeof(double*));
    if(b->llsq_a==NULL) {freeDelayFitData(&ltdata); return(TPCERROR_OUT_OF_MEMORY);}
    for(int ni=0; ni<b->llsq_n; ni++) b->llsq_a[ni]=b->llsq_mat+ni*b->llsq_m;
    b->llsq_b=(double*)malloc(b->llsq_m*sizeof(double));
    b->llsq_x=(double*)malloc(b->llsq_n*sizeof(double));
    b->llsq_wp=(double*)malloc(b->llsq_n*sizeof(double));
    b->llsq_zz=(double*)malloc(b->llsq_m*sizeof(double));
    b->llsq_i=(int*)malloc(b->llsq_n*sizeof(int));
    if(b->llsq_b==NULL || b->llsq_x==NULL || b->llsq_wp==NULL || b->llsq_zz==NULL || b->llsq_i==NULL) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
      freeDelayFitData(&ltdata); return(TPCERROR_OUT_OF_MEMORY);}
  }
 
  /* Integrate tissue curve */
  {
    if(verbose>5) {printf("  preparing tissue data\n"); fflush(stdout);}
    int ti=b->dtac.tacNr-1; // index of tissue curve after the input curves
    for(int i=0; i<b->dtac.sampleNr; i++) b->dtac.c[ti].y[i]=ttac->c[ci].y[i];
 
    int ret=liDerivate(b->dtac.x, b->dtac.c[ti].y, b->dtac.sampleNr, b->ddtac.c[ti].y, b->dddtac.c[ti].y, 0);
    if(b->dtac.isframe==0) {
      if(!ret) ret=liIntegrate(b->dtac.x, b->dtac.c[ti].y, b->dtac.sampleNr, b->ditac.c[ti].y, 3, 0);
      if(!ret) ret=liIntegrate(b->dtac.x, b->ditac.c[ti].y, b->dtac.sampleNr, b->diitac.c[ti].y, 3, 0);
    } else {
      if(!ret) ret=liIntegrateFE(b->dtac.x1, b->dtac.x2, b->dtac.c[ti].y, b->dtac.sampleNr, 
                                 b->ditac.c[ti].y, b->diitac.c[ti].y, 0);
    }
    if(ret) {
      freeDelayFitData(&ltdata);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return(TPCERROR_INVALID_XRANGE);
    }
 
    if(verbose>10) {
      printf("integrated tissue curves:\n");
      for(int i=0; i<b->dtac.sampleNr; i++) 
        printf("%g  %g %g %g\n", b->dtac.x[i], b->dtac.c[ti].y[i], b->ditac.c[ti].y[i], b->diitac.c[ti].y[i]);
      fflush(stdout);
    }
  }
 
 
  /* Multilinear fitting with each delayed input */
  if(verbose>5) {printf("  LLSQ\n"); fflush(stdout);}
  double ss, ssmin=nan("");
  int mini=-1;
  for(int di=0; di<b->dtac.tacNr-1; di++) {
    /* Setup data matrix A and vector B */
    for(int mi=0; mi<b->llsq_m; mi++) {
      b->llsq_b[mi]=b->dtac.c[b->dtac.tacNr-1].y[mi]; // TTAC
      b->llsq_mat[mi]=b->dtac.c[di].y[mi]; // BTAC
      if(mode==1) {
        b->llsq_mat[mi+1*b->llsq_m]=b->ddtac.c[di].y[mi]; // BTAC derivate
        b->llsq_mat[mi+2*b->llsq_m]=-b->ddtac.c[b->dtac.tacNr-1].y[mi]; // TTAC derivate
        if(b->llsq_n>3) b->llsq_mat[mi+3*b->llsq_m]=b->dddtac.c[di].y[mi]; // BTAC derivate
        if(b->llsq_n>4) b->llsq_mat[mi+4*b->llsq_m]=-b->dddtac.c[b->dtac.tacNr-1].y[mi]; // TTAC derivate
      } else if(mode==2) {
        b->llsq_mat[mi+1*b->llsq_m]=b->ditac.c[di].y[mi]; // BTAC integral
        b->llsq_mat[mi+2*b->llsq_m]=-b->ditac.c[b->dtac.tacNr-1].y[mi]; // TTAC integral
        if(b->llsq_n>3) b->llsq_mat[mi+3*b->llsq_m]=b->ddtac.c[di].y[mi]; // BTAC derivate
        if(b->llsq_n>4) b->llsq_mat[mi+4*b->llsq_m]=-b->ddtac.c[b->dtac.tacNr-1].y[mi]; // TTAC derivate
      } else { // mode 0
        b->llsq_mat[mi+1*b->llsq_m]=b->ditac.c[di].y[mi]; // BTAC integral
        b->llsq_mat[mi+2*b->llsq_m]=-b->ditac.c[b->dtac.tacNr-1].y[mi]; // TTAC integral
        if(b->llsq_n>3) b->llsq_mat[mi+3*b->llsq_m]=b->diitac.c[di].y[mi]; // BTAC 2nd integral
        if(b->llsq_n>4) b->llsq_mat[mi+4*b->llsq_m]=-b->diitac.c[b->dtac.tacNr-1].y[mi]; // TTAC 2nd integral
      }
      b->llsq_mat[mi]=b->dtac.c[di].y[mi]; // BTAC
      b->llsq_mat[mi+1*b->llsq_m]=b->ditac.c[di].y[mi]; // BTAC integral
      b->llsq_mat[mi+2*b->llsq_m]=-b->ditac.c[b->dtac.tacNr-1].y[mi]; // TTAC integral
      if(b->llsq_n>3) b->llsq_mat[mi+3*b->llsq_m]=b->diitac.c[di].y[mi]; // BTAC 2nd integral
      if(b->llsq_n>4) b->llsq_mat[mi+4*b->llsq_m]=-b->diitac.c[b->dtac.tacNr-1].y[mi]; // TTAC 2nd integral
    }
    if(verbose>3) {
      int failnr=0;
      for(int mi=0; mi<b->llsq_m; mi++) {
        if(!isfinite(b->llsq_b[mi])) failnr++;
        if(!isfinite(b->llsq_mat[mi])) failnr++;
        if(!isfinite(b->llsq_mat[mi+1*b->llsq_m])) failnr++;
        if(!isfinite(b->llsq_mat[mi+2*b->llsq_m])) failnr++;
        if(!isfinite(b->llsq_mat[mi+3*b->llsq_m])) failnr++;
        if(!isfinite(b->llsq_mat[mi+4*b->llsq_m])) failnr++;
      }
      if(failnr>0) printf("%d bad values\n", failnr);
    }
    if(verbose>9 && b->llsq_m<100) {
      printf("---- B x A ----\n");
      for(int mi=0; mi<b->llsq_m; mi++) {
        printf("  %g ", b->llsq_b[mi]);
        for(int ni=0; ni<b->llsq_n; ni++)
          printf(" %g", b->llsq_mat[mi+ni*b->llsq_m]);
         printf("\n");
      }
    }
 
    int ret=nnls(b->llsq_a, b->llsq_m, b->llsq_n, b->llsq_b, b->llsq_x, &ss, b->llsq_wp, b->llsq_zz, b->llsq_i);
    if(ret<2 && isfinite(ss)) {
      if(!isfinite(ssmin) || ssmin>ss) {ssmin=ss; mini=di;}
      if(verbose>2) printf("   %d %d %g %g\n", di, mini, ss, ssmin);
    } else {
      if(verbose>1) printf("  failed NNLS\n");
    }
  }
  if(mini<0) {
    if(verbose>1) fprintf(stderr, "Error: all fits failed.\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_BAD_FIT);
    freeDelayFitData(&ltdata); return(TPCERROR_BAD_FIT);
  }
  if(dt!=NULL && mini>=0) *dt=b->dtmin+mini*b->dtstep;
 
  freeDelayFitData(&ltdata); // the local one, not the one possibly given by caller 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/