LALInferenceInitCBC.c

/*
 *  LALInferenceCBCInit.c:  Bayesian Followup initialisation routines.
 *
 *  Copyright (C) 2012 Vivien Raymond, John Veitch, Salvatore Vitale
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with with program; see the file COPYING. If not, write to the
 *  Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 *  MA  02110-1301  USA
 */


#include <stdio.h>
#include <assert.h>
#include <errno.h>
#include <lal/Date.h>
#include <lal/GenerateInspiral.h>
#include <lal/LALInference.h>
#include <lal/FrequencySeries.h>
#include <lal/Units.h>
#include <lal/StringInput.h>
#include <lal/LIGOLwXMLInspiralRead.h>
#include <lal/TimeSeries.h>
#include <lal/LALInferencePrior.h>
#include <lal/LALInferenceTemplate.h>
#include <lal/LALInferenceProposal.h>
#include <lal/LALInferenceLikelihood.h>
#include <lal/LALInferenceReadData.h>
#include <lal/LALInferenceInit.h>
#include <lal/LALInferenceCalibrationErrors.h>
#include <lal/LALSimNeutronStar.h>

static int checkParamInList(const char *list, const char *param);
static int checkParamInList(const char *list, const char *param)
{
  /* Check for param in comma-seperated list */
  char *post=NULL,*pos=NULL;
  if (list==NULL) return 0;
  if (param==NULL) return 0;

  if(!(pos=strstr(list,param))) return 0;

  /* The string is a substring. Check that it is a token */
  /* Check the character before and after */
  if(pos!=list)
  if(*(pos-1)!=',')
  return 0;

  post=&(pos[strlen(param)]);
  if(*post!='\0')
  if(*post!=',')
  return 0;
  return 1;
}

static void print_flags_orders_warning(SimInspiralTable *injt, ProcessParamsTable *commline);
static void LALInferenceInitSpinVariables(LALInferenceRunState *state, LALInferenceModel *model);
static void LALInferenceInitMassVariables(LALInferenceRunState *state);
static void LALInferenceInitNonGRParams(LALInferenceRunState *state, LALInferenceModel *model);
static void LALInferenceCheckApproximantNeeds(LALInferenceRunState *state,Approximant approx);


/* Initialize a bare-bones run-state,
   calling the "ReadData()" function to gather data & PSD from files,
   sets up the random seed and rng, and initializes other variables accordingly.
*/
LALInferenceRunState *LALInferenceInitRunState(ProcessParamsTable *command_line) {
    ProcessParamsTable *ppt=NULL;
    INT4 randomseed;
    FILE *devrandom;
    struct timeval tv;
    LALInferenceRunState *run_state = XLALCalloc(1, sizeof(LALInferenceRunState));

    /* Check that command line is consistent first */
    LALInferenceCheckOptionsConsistency(command_line);
    run_state->commandLine = command_line;

    /* Initialize parameters structure */
    run_state->algorithmParams = XLALCalloc(1, sizeof(LALInferenceVariables));
    run_state->priorArgs = XLALCalloc(1, sizeof(LALInferenceVariables));
    run_state->proposalArgs = XLALCalloc(1, sizeof(LALInferenceVariables));

    /* Read data from files or generate fake data */
    run_state->data = LALInferenceReadData(command_line);
    if (run_state->data == NULL)
        return(NULL);

    /* Setup the random number generator */
    gsl_rng_env_setup();
    run_state->GSLrandom = gsl_rng_alloc(gsl_rng_mt19937);

    /* Use clocktime if seed isn't provided */
    ppt = LALInferenceGetProcParamVal(command_line, "--randomseed");
    if (ppt)
        randomseed = atoi(ppt->value);
    else {
        if ((devrandom = fopen("/dev/urandom","r")) == NULL) {
            gettimeofday(&tv, 0);
            randomseed = tv.tv_sec + tv.tv_usec;
        } else {
            fread(&randomseed, sizeof(randomseed), 1, devrandom);
            fclose(devrandom);
        }
    }
    gsl_rng_set(run_state->GSLrandom, randomseed);

    /* Save the random seed */
    LALInferenceAddVariable(run_state->algorithmParams, "random_seed", &randomseed,
                            LALINFERENCE_INT4_t, LALINFERENCE_PARAM_OUTPUT);

    return(run_state);
}

/* Draw initial parameters for each of the threads in run state */
void LALInferenceDrawThreads(LALInferenceRunState *run_state) {
    if (run_state == NULL)
        return;

    if (run_state->threads == NULL) {
        XLALPrintError("Error: LALInferenceDrawThreads expects initialized run_state->threads\n");
        XLAL_ERROR_VOID(XLAL_EINVAL);
    }

    LALInferenceThreadState *thread = &(run_state->threads[0]);
    INT4 t;

    /* If using a malmquist prior, force a strict prior window on distance for starting point, otherwise
     * the approximate prior draws are very unlikely to be within the malmquist prior */
    REAL8 dist_low, dist_high;
    REAL8 restricted_dist_low = log(10.0);
    REAL8 restricted_dist_high = log(100.0);
    INT4 changed_dist = 0;
    if (LALInferenceCheckVariable(run_state->priorArgs, "malmquist") &&
        LALInferenceCheckVariableNonFixed(thread->currentParams, "logdistance")) {
        changed_dist = 1;
        LALInferenceGetMinMaxPrior(run_state->priorArgs, "logdistance", &dist_low, &dist_high);
        LALInferenceRemoveMinMaxPrior(run_state->priorArgs, "logdistance");
        LALInferenceAddMinMaxPrior(run_state->priorArgs, "logdistance",
                                   &restricted_dist_low, &restricted_dist_high, LALINFERENCE_REAL8_t);
    }

    /* If the currentParams are not in the prior, overwrite and pick paramaters
     *   from the priors. OVERWRITE EVEN USER CHOICES.
     *   (necessary for complicated prior shapes where
     *   LALInferenceCyclicReflectiveBound() is not enough) */
    #pragma omp parallel for private(thread)
    for (t = 0; t < run_state->nthreads; t++) {
        LALInferenceVariables *priorDraw = XLALCalloc(1, sizeof(LALInferenceVariables));

        thread = &(run_state->threads[t]);

        /* Try not to clobber values given on the command line */
        LALInferenceCopyVariables(thread->currentParams, priorDraw);
        LALInferenceDrawApproxPrior(thread, priorDraw, priorDraw);
        LALInferenceCopyUnsetREAL8Variables(priorDraw, thread->currentParams,
                                            run_state->commandLine);

        while(isinf(run_state->prior(run_state,
                                     thread->currentParams,
                                     thread->model))) {
            LALInferenceDrawApproxPrior(thread,
                                        thread->currentParams,
                                        thread->currentParams);
        }

        /* Make sure that our initial value is within the
        *     prior-supported volume. */
        LALInferenceCyclicReflectiveBound(thread->currentParams, run_state->priorArgs);

        /* Initialize starting likelihood and prior */
        thread->currentPrior = run_state->prior(run_state,
                                                thread->currentParams,
                                                thread->model);

        thread->currentLikelihood = run_state->likelihood(thread->currentParams,
                                                          run_state->data,
                                                          thread->model);

        LALInferenceClearVariables(priorDraw);
        XLALFree(priorDraw);
    }

    /* Replace distance prior if changed for initial sample draw */
    if (changed_dist) {
        LALInferenceRemoveMinMaxPrior(run_state->priorArgs, "logdistance");
        LALInferenceAddMinMaxPrior(run_state->priorArgs, "logdistance",
                                   &dist_low, &dist_high, LALINFERENCE_REAL8_t);
    }
}

/*
 * Initialize threads in memory, using LALInferenceInitCBCModel() to init models.
 */
void LALInferenceInitCBCThreads(LALInferenceRunState *run_state, INT4 nthreads) {
  if (run_state == NULL){
    LALInferenceInitCBCModel(run_state);
    return;
  }

  ProcessParamsTable *commandLine=run_state->commandLine;

  LALInferenceThreadState *thread;
  INT4 t, nifo;
  INT4 randomseed;
  LALInferenceIFOData *data;
  run_state->nthreads = nthreads;
  run_state->threads = LALInferenceInitThreads(nthreads);

  for (t = 0; t < nthreads; t++) {
    thread = &(run_state->threads[t]);

    /* Link back to run-state */
    thread->parent = run_state;

    /* Set up CBC model and parameter array */
    thread->model = LALInferenceInitCBCModel(run_state);
    thread->model->roq_flag = 0;

    /* Allocate IFO likelihood holders */
    nifo = 0;
    data = run_state->data;
    while (data != NULL) {
        data = data->next;
        nifo++;
    }
    thread->currentIFOLikelihoods = XLALCalloc(nifo, sizeof(REAL8));

    /* Setup ROQ */
    if (LALInferenceGetProcParamVal(commandLine, "--roqtime_steps")){

        LALInferenceSetupROQmodel(thread->model, commandLine);
        fprintf(stderr, "done LALInferenceSetupROQmodel\n");

    }else{
      thread->model->roq_flag=0;
    }

    LALInferenceCopyVariables(thread->model->params, thread->currentParams);
    LALInferenceCopyVariables(run_state->proposalArgs, thread->proposalArgs);

    /* Link thread-state prior-args to the parent run-state's */
    thread->priorArgs = run_state->priorArgs;

    /* Use clocktime if seed isn't provided */
    thread->GSLrandom = gsl_rng_alloc(gsl_rng_mt19937);
    randomseed = gsl_rng_get(run_state->GSLrandom);
    gsl_rng_set(thread->GSLrandom, randomseed);
  }

  return;
}


/* Setup the template generation */
/* Defaults to using LALSimulation */
LALInferenceTemplateFunction LALInferenceInitCBCTemplate(LALInferenceRunState *runState)
{
  char help[]="(--template [LAL,PhenSpin,LALGenerateInspiral,LALSim,multiband]\tSpecify template (default LAL)\n";
  ProcessParamsTable *ppt=NULL;
  ProcessParamsTable *commandLine=runState->commandLine;
  /* Print command line arguments if help requested */
  //Help is taken care of in LALInferenceInitCBCVariables
  //ppt=LALInferenceGetProcParamVal(commandLine,"--help");
  //if(ppt)
  //{
  //	fprintf(stdout,"%s",help);
  //	return;
  //}
  /* This is the LAL template generator for inspiral signals */
  LALInferenceTemplateFunction templt = &LALInferenceTemplateXLALSimInspiralChooseWaveform;
  ppt=LALInferenceGetProcParamVal(commandLine,"--template");
  if(ppt) {
    if(!strcmp("LALSim",ppt->value))
      templt=&LALInferenceTemplateXLALSimInspiralChooseWaveform;
    else if(!strcmp("null",ppt->value))
        templt=&LALInferenceTemplateNullFreqdomain;
	else if(!strcmp("multiband",ppt->value)){
        templt=&LALInferenceTemplateXLALSimInspiralChooseWaveformPhaseInterpolated;
        fprintf(stdout,"Template function called is \"LALInferenceTemplateXLALSimInspiralChooseWaveformPhaseInterpolated\"\n");
    }
    else {
        XLALPrintError("Error: unknown template %s\n",ppt->value);
        XLALPrintError("%s", help);
        XLAL_ERROR_NULL(XLAL_EINVAL);
    }
  }
  else if(LALInferenceGetProcParamVal(commandLine,"--LALSimulation")){
    fprintf(stderr,"Warning: --LALSimulation is deprecated, the LALSimulation package is now the default. To use LALInspiral specify:\n\
                    --template LALGenerateInspiral (for time-domain templates)\n\
                    --template LAL (for frequency-domain templates)\n");
  }
  else if(LALInferenceGetProcParamVal(commandLine,"--roqtime_steps")){
  templt=&LALInferenceROQWrapperForXLALSimInspiralChooseFDWaveformSequence;
        fprintf(stderr, "template is \"LALInferenceROQWrapperForXLALSimInspiralChooseFDWaveformSequence\"\n");
  }
  else {
    fprintf(stdout,"Template function called is \"LALInferenceTemplateXLALSimInspiralChooseWaveform\"\n");
  }
  return templt;
}

/* Setup the glitch model */
void LALInferenceInitGlitchVariables(LALInferenceRunState *runState, LALInferenceVariables *currentParams)
{
  ProcessParamsTable    *commandLine   = runState->commandLine;
  LALInferenceIFOData   *dataPtr       = runState->data;
  LALInferenceVariables *priorArgs     = runState->priorArgs;

  UINT4 i,nifo;
  UINT4 n = (UINT4)dataPtr->timeData->data->length;
  UINT4 gflag  = 1;
  REAL8 gmin   = 0.0;
  REAL8 gmax   = 20.0;

  //over-ride default gmax from command line
  if(LALInferenceGetProcParamVal(commandLine, "--glitchNmax"))
    gmax = (REAL8)atoi(LALInferenceGetProcParamVal(commandLine, "--glitchNmax")->value);

  //count interferometers in network before allocating memory
  //compute imin,imax for each IFO -- may be different
  nifo=0;
  dataPtr = runState->data;
  while (dataPtr != NULL)
  {
    dataPtr = dataPtr->next;
    nifo++;
  }
  dataPtr = runState->data;

  UINT4Vector *gsize  = XLALCreateUINT4Vector(nifo);
  //Meyer?? REAL8Vector *gprior = XLALCreateREAL8Vector((int)gmax+1);

  //Morlet??
  gsl_matrix *mAmp = gsl_matrix_alloc(nifo,(int)(gmax));
  gsl_matrix *mf0  = gsl_matrix_alloc(nifo,(int)(gmax));
  gsl_matrix *mQ   = gsl_matrix_alloc(nifo,(int)(gmax));
  gsl_matrix *mt0  = gsl_matrix_alloc(nifo,(int)(gmax));
  gsl_matrix *mphi = gsl_matrix_alloc(nifo,(int)(gmax));

  double Amin,Amax;
  double Qmin,Qmax;
  double f_min,f_max;
  double tmin,tmax;
  double pmin,pmax;
  double Anorm;

  REAL8 TwoDeltaToverN = 2.0 * dataPtr->timeData->deltaT / ((double) dataPtr->timeData->data->length);

  Anorm = sqrt(TwoDeltaToverN);
  Amin = 10.0/Anorm;
  Amax = 10000.0/Anorm;

  Qmin = 3.0;
  Qmax = 30.0;
  tmin = 0.0;
  tmax = dataPtr->timeData->data->length*dataPtr->timeData->deltaT;
  f_min = dataPtr->fLow;
  f_max = dataPtr->fHigh;
  pmin = 0.0;
  pmax = LAL_TWOPI;

  gsl_matrix_set_all(mAmp, Amin);
  gsl_matrix_set_all(mf0,  f_min);
  gsl_matrix_set_all(mQ,   Qmin);
  gsl_matrix_set_all(mt0,  tmin);
  gsl_matrix_set_all(mphi, pmin);

  gsl_matrix  *gFD       = gsl_matrix_alloc(nifo,(int)n); //store the Fourier-domain glitch signal

  for(i=0; i<nifo; i++) gsize->data[i]=0;

  //Morlet wavelet parameters
  LALInferenceAddVariable(currentParams, "morlet_FD",  &gFD,  LALINFERENCE_gslMatrix_t, LALINFERENCE_PARAM_LINEAR);
  LALInferenceAddVariable(currentParams, "morlet_Amp", &mAmp, LALINFERENCE_gslMatrix_t, LALINFERENCE_PARAM_LINEAR);
  LALInferenceAddVariable(currentParams, "morlet_f0" , &mf0,  LALINFERENCE_gslMatrix_t, LALINFERENCE_PARAM_LINEAR);
  LALInferenceAddVariable(currentParams, "morlet_Q"  , &mQ,   LALINFERENCE_gslMatrix_t, LALINFERENCE_PARAM_LINEAR);
  LALInferenceAddVariable(currentParams, "morlet_t0" , &mt0,  LALINFERENCE_gslMatrix_t, LALINFERENCE_PARAM_LINEAR);
  LALInferenceAddVariable(currentParams, "morlet_phi", &mphi, LALINFERENCE_gslMatrix_t, LALINFERENCE_PARAM_LINEAR);

  LALInferenceAddVariable(currentParams, "glitch_size", &gsize, LALINFERENCE_UINT4Vector_t, LALINFERENCE_PARAM_LINEAR);
  LALInferenceAddVariable(currentParams, "glitchFitFlag", &gflag, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED);

  LALInferenceAddMinMaxPrior(priorArgs, "morlet_Amp_prior", &Amin, &Amax, LALINFERENCE_REAL8_t);
  LALInferenceAddMinMaxPrior(priorArgs, "morlet_f0_prior" , &f_min, &f_max, LALINFERENCE_REAL8_t);
  LALInferenceAddMinMaxPrior(priorArgs, "morlet_Q_prior"  , &Qmin, &Qmax, LALINFERENCE_REAL8_t);
  LALInferenceAddMinMaxPrior(priorArgs, "morlet_t0_prior" , &tmin, &tmax, LALINFERENCE_REAL8_t);
  LALInferenceAddMinMaxPrior(priorArgs, "morlet_phi_prior", &pmin, &pmax, LALINFERENCE_REAL8_t);

  LALInferenceAddMinMaxPrior(priorArgs, "glitch_size", &gmin, &gmax, LALINFERENCE_REAL8_t);
  LALInferenceAddMinMaxPrior(priorArgs, "glitch_dim", &gmin, &gmax, LALINFERENCE_REAL8_t);

  LALInferenceAddVariable(priorArgs, "glitch_norm", &Anorm, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
}

struct spcal_envelope
{
    gsl_spline  *amp_median,*amp_std,
                *phase_median,*phase_std;
};

/* Format string for the calibratino envelope file */
/* Frequency    Median Mag     Phase (Rad)    -1 Sigma Mag   -1 Sigma Phase +1 Sigma Mag   +1 Sigma Phase */

#define CAL_ENV_FORMAT "%lf %lf %lf %lf %lf %lf %lf\n"

static struct spcal_envelope *initCalibrationEnvelope(char *filename);

static struct spcal_envelope *initCalibrationEnvelope(char *filename)
{
    FILE *fp=fopen(filename,"r");
    char tmpline[1024];
    if(!fp) {fprintf(stderr,"Unable to open %s: Error %i %s\n",filename,errno,strerror(errno)); exit(1);}
    int Nlines=0;
    REAL8 freq, *logfreq=NULL, *mag_med=NULL, mag_low, mag_hi, *mag_std=NULL, *phase_med=NULL, phase_low, phase_hi, *phase_std=NULL;
    for(Nlines=0;fgets(tmpline,1024,fp); )
    {
        /* Skip header */
        if(tmpline[0]=='#') continue;
        /* Grow arrays */
        logfreq=realloc(logfreq,sizeof(*logfreq)*(Nlines+1));
        mag_med=realloc(mag_med,sizeof(*mag_med)*(Nlines+1));
        mag_std=realloc(mag_std,sizeof(*mag_std)*(Nlines+1));
        phase_med=realloc(phase_med,sizeof(*phase_med)*(Nlines+1));
        phase_std=realloc(phase_std,sizeof(*phase_std)*(Nlines+1));

        if((7!=sscanf(tmpline,CAL_ENV_FORMAT, &freq, &(mag_med[Nlines]), &(phase_med[Nlines]), &mag_low, &phase_low, &mag_hi, &phase_hi)))
        {
            fprintf(stderr,"Malformed input line in file %s: %s\n",filename,tmpline);
            exit(1);
        }
		mag_med[Nlines]-=1.0; /* Subtract off 1 to get delta */
        logfreq[Nlines]=log(freq);
        mag_std[Nlines]=(mag_hi - mag_low ) /2.0;
        phase_std[Nlines]=(phase_hi - phase_low) /2.0;
		Nlines++;
    }
    fprintf(stdout,"Read %i lines from calibration envelope %s\n",Nlines,filename);
    fclose(fp);

    struct spcal_envelope *env=XLALCalloc(1,sizeof(*env));
    env->amp_median = gsl_spline_alloc ( gsl_interp_cspline, Nlines);
    env->amp_std = gsl_spline_alloc ( gsl_interp_cspline, Nlines);
    env->phase_median = gsl_spline_alloc ( gsl_interp_cspline, Nlines);
    env->phase_std = gsl_spline_alloc ( gsl_interp_cspline, Nlines);

    gsl_spline_init(env->amp_median, logfreq, mag_med, Nlines);
    gsl_spline_init(env->amp_std, logfreq, mag_std, Nlines);
    gsl_spline_init(env->phase_median, logfreq, phase_med, Nlines);
    gsl_spline_init(env->phase_std, logfreq, phase_std, Nlines);

    free(logfreq); free(mag_med); free(mag_std); free(phase_med); free(phase_std);

    return(env);
}

static int destroyCalibrationEnvelope(struct spcal_envelope *env);
static int destroyCalibrationEnvelope(struct spcal_envelope *env)
{
    if(!env) XLAL_ERROR(XLAL_EINVAL);
    if(env->amp_median) gsl_spline_free(env->amp_median);
    if(env->amp_std) gsl_spline_free(env->amp_std);
    if(env->phase_median) gsl_spline_free(env->phase_median);
    if(env->phase_std) gsl_spline_free(env->phase_std);
    XLALFree(env);
    return XLAL_SUCCESS;
}

void LALInferenceInitCalibrationVariables(LALInferenceRunState *runState, LALInferenceVariables *currentParams) {
  ProcessParamsTable *ppt = NULL;
  LALInferenceIFOData *ifo = NULL;
  LALInferenceIFOData *dataPtr=NULL;
  UINT4 calOn = 1;
  if ((ppt = LALInferenceGetProcParamVal(runState->commandLine, "--enable-spline-calibration"))){
    /* Use spline to marginalize*/
    UINT4 ncal = 5; /* Number of calibration nodes, log-distributed
		between fmin and fmax. */
    REAL8 ampUncertaintyPrior = 0.1; /* 10% amplitude */
    REAL8 phaseUncertaintyPrior = 5*M_PI/180.0; /* 5 degrees phase */
    if ((ppt = LALInferenceGetProcParamVal(runState->commandLine, "--spcal-nodes"))) {
      ncal = atoi(ppt->value);
      if (ncal < 3) { /* Cannot do spline with fewer than 3 points! */
	fprintf(stderr, "ERROR: given '--spcal-nodes %d', but cannot spline with fewer than 3\n", ncal);
	exit(1);
      }
    }

    LALInferenceAddVariable(currentParams, "spcal_active", &calOn, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED);
    LALInferenceAddVariable(currentParams, "spcal_npts", &ncal, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED);

    for(ifo=runState->data;ifo;ifo=ifo->next) {
      UINT4 i;

      char freqVarName[VARNAME_MAX];
      char ampVarName[VARNAME_MAX];
      char phaseVarName[VARNAME_MAX];

      REAL8 fMin = ifo->fLow;
      REAL8 fMax = ifo->fHigh;
      REAL8 logFMin = log(fMin);
      REAL8 logFMax = log(fMax);
      REAL8 dLogF = (logFMax - logFMin)/(ncal-1);

      char amp_uncert_op[VARNAME_MAX];
      char pha_uncert_op[VARNAME_MAX];
      char env_uncert_op[VARNAME_MAX];
      struct spcal_envelope *env=NULL;

      if((VARNAME_MAX <= snprintf(amp_uncert_op, VARNAME_MAX, "--%s-spcal-amp-uncertainty", ifo->name))
          || (VARNAME_MAX <= snprintf(pha_uncert_op, VARNAME_MAX, "--%s-spcal-phase-uncertainty", ifo->name))
          || (VARNAME_MAX <= snprintf(env_uncert_op, VARNAME_MAX, "--%s-spcal-envelope",ifo->name)) )
      {
        fprintf(stderr,"variable name too long\n"); exit(1);
      }

      if( (ppt=LALInferenceGetProcParamVal(runState->commandLine, env_uncert_op)))
          env = initCalibrationEnvelope(ppt->value);
      else
      {
        if ((ppt = LALInferenceGetProcParamVal(runState->commandLine, amp_uncert_op))) {
            ampUncertaintyPrior = atof(ppt->value);
        }
        else{
            fprintf(stderr,"Error, missing %s or %s\n",amp_uncert_op, env_uncert_op);
            exit(1);
        }

        if ((ppt = LALInferenceGetProcParamVal(runState->commandLine, pha_uncert_op))) {
            phaseUncertaintyPrior = M_PI/180.0*atof(ppt->value); /* CL arg in degrees, variable in radians */
        }
        else{
            fprintf(stderr,"Error, missing %s or %s\n",pha_uncert_op,env_uncert_op);
            exit(1);
        }
      }
      /* Now add each spline node */
      for(i=0;i<ncal;i++)
	  {
			  if((VARNAME_MAX <= snprintf(freqVarName, VARNAME_MAX, "%s_spcal_logfreq_%i",ifo->name,i))
                  || (VARNAME_MAX <= snprintf(ampVarName, VARNAME_MAX, "%s_spcal_amp_%i", ifo->name,i))
                  || (VARNAME_MAX <= snprintf(phaseVarName, VARNAME_MAX, "%s_spcal_phase_%i", ifo->name,i)) )
              {
                  fprintf(stderr,"Variable name too long\n"); exit(1);
              }
              REAL8 amp_std=ampUncertaintyPrior,amp_mean=0.0;
			  REAL8 phase_std=phaseUncertaintyPrior,phase_mean=0.0;
			  REAL8 logFreq = logFMin + i*dLogF;
			  LALInferenceAddREAL8Variable(currentParams,freqVarName,logFreq,LALINFERENCE_PARAM_FIXED);
			  if(env)
			  {
					  amp_std = gsl_spline_eval(env->amp_std, logFreq, NULL);
					  amp_mean = gsl_spline_eval(env->amp_median, logFreq, NULL);
					  phase_std = gsl_spline_eval(env->phase_std, logFreq, NULL);
					  phase_mean = gsl_spline_eval(env->phase_median, logFreq, NULL);
			  }
			  LALInferenceRegisterGaussianVariableREAL8(runState, currentParams, ampVarName, 0, amp_mean, amp_std, LALINFERENCE_PARAM_LINEAR);
			  LALInferenceRegisterGaussianVariableREAL8(runState, currentParams, phaseVarName, 0, phase_mean, phase_std, LALINFERENCE_PARAM_LINEAR);
	  } /* End loop over spline nodes */

	  if(env) destroyCalibrationEnvelope(env);
	} /* End loop over IFOs */
  } /* End case of spline calibration error */
  else if(LALInferenceGetProcParamVal(runState->commandLine, "--MarginalizeConstantCalAmp") ||LALInferenceGetProcParamVal(runState->commandLine, "--MarginalizeConstantCalPha")){
    /* Use constant (in frequency) approximation for the errors */
    if (LALInferenceGetProcParamVal(runState->commandLine, "--MarginalizeConstantCalAmp")){
      /*For the moment the prior ranges are the same for the three IFOs */
      REAL8 camp_max_A=0.25; /* plus minus 25% amplitude errors*/
      REAL8 camp_min_A=-0.25;
      REAL8 zero=0.0;
      dataPtr = runState->data;
      while (dataPtr != NULL){
        char CA_A[320];
        sprintf(CA_A,"%s_%s","calamp",dataPtr->name);
        LALInferenceRegisterUniformVariableREAL8(runState, currentParams, CA_A, zero, camp_min_A, camp_max_A, LALINFERENCE_PARAM_LINEAR);
        dataPtr = dataPtr->next;
      }

      LALInferenceAddVariable(currentParams, "constantcal_active", &calOn, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED);
      /*If user specifies a width for the error prior, a gaussian prior will be used, otherwise a flat prior will be used*/
      REAL8 ampUncertaintyPrior=-1.0;
      ppt = LALInferenceGetProcParamVal(runState->commandLine, "--constcal_ampsigma");
      if (ppt) {
        ampUncertaintyPrior = atof(ppt->value);
      }
      LALInferenceAddVariable(runState->priorArgs, "constcal_amp_uncertainty", &ampUncertaintyPrior, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);
    }
    if (LALInferenceGetProcParamVal(runState->commandLine, "--MarginalizeConstantCalPha")){
      /* Add linear calibration phase errors to the measurement. For the moment the prior ranges are the same for the three IFOs */
      REAL8 cpha_max_A=0.349;  /* plus/minus 20 degs*/
      REAL8 cpha_min_A=-0.349;
      REAL8 zero=0.0;
      dataPtr = runState->data;
      while (dataPtr != NULL)
      {
        char CP_A[320];
        sprintf(CP_A,"%s_%s","calpha",dataPtr->name);
        LALInferenceRegisterUniformVariableREAL8(runState, currentParams, CP_A, zero, cpha_min_A, cpha_max_A, LALINFERENCE_PARAM_LINEAR);
        dataPtr = dataPtr->next;
      }
      LALInferenceAddVariable(currentParams, "constantcal_active", &calOn, LALINFERENCE_UINT4_t, LALINFERENCE_PARAM_FIXED);

     /*If user specifies a width for the error prior, a gaussian prior will be used, otherwise a flat prior will be used*/
      REAL8 phaseUncertaintyPrior=-1.0;
      ppt = LALInferenceGetProcParamVal(runState->commandLine, "--constcal_phasigma");
      if (ppt) {
        phaseUncertaintyPrior = M_PI/180.0*atof(ppt->value); /* CL arg in degrees, variable in radians */
      }
      LALInferenceAddVariable(runState->priorArgs, "constcal_phase_uncertainty", &phaseUncertaintyPrior, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_FIXED);

    }
  }
  else{
    /* No calibration marginalization asked. Just exit */
    return;
  }
}

void LALInferenceRegisterGaussianVariableREAL8(LALInferenceRunState *state, LALInferenceVariables *var, const char name[VARNAME_MAX], REAL8 startval, REAL8 mean, REAL8 stdev, LALInferenceParamVaryType varytype)
{
  char meanopt[VARNAME_MAX+8];
  char sigmaopt[VARNAME_MAX+9];
  char valopt[VARNAME_MAX+3];
  char fixopt[VARNAME_MAX+7];
  ProcessParamsTable *ppt=NULL;

  sprintf(meanopt,"--%s-mean",name);
  sprintf(sigmaopt,"--%s-sigma",name);
  sprintf(valopt,"--%s",name);
  sprintf(fixopt,"--fix-%s",name);

  if((ppt=LALInferenceGetProcParamVal(state->commandLine,meanopt))) mean=atof(ppt->value);
  if((ppt=LALInferenceGetProcParamVal(state->commandLine,sigmaopt))) stdev=atof(ppt->value);
  if((ppt=LALInferenceGetProcParamVal(state->commandLine,fixopt)))
  {
    varytype = LALINFERENCE_PARAM_FIXED;
    startval = atof(ppt->value);
  }
  if((ppt=LALInferenceGetProcParamVal(state->commandLine,valopt))) startval=atof(ppt->value);

  assert(stdev>0);
  LALInferenceAddVariable(var,name,&startval,LALINFERENCE_REAL8_t,varytype);
  LALInferenceAddGaussianPrior(state->priorArgs, name, &mean, &stdev, LALINFERENCE_REAL8_t);

}

void LALInferenceRegisterUniformVariableREAL8(LALInferenceRunState *state, LALInferenceVariables *var, const char name[VARNAME_MAX], REAL8 startval, REAL8 min, REAL8 max, LALInferenceParamVaryType varytype)
{
  char minopt[VARNAME_MAX+7];
  char maxopt[VARNAME_MAX+7];
  char valopt[VARNAME_MAX+3];
  char fixopt[VARNAME_MAX+7];
  ProcessParamsTable *ppt=NULL;

  sprintf(minopt,"--%s-min",name);
  sprintf(maxopt,"--%s-max",name);
  sprintf(valopt,"--%s",name);
  sprintf(fixopt,"--fix-%s",name);

  if((ppt=LALInferenceGetProcParamVal(state->commandLine,minopt))) min=atof(ppt->value);
  if((ppt=LALInferenceGetProcParamVal(state->commandLine,maxopt))) max=atof(ppt->value);
  if((ppt=LALInferenceGetProcParamVal(state->commandLine,fixopt))) {
		  varytype=LALINFERENCE_PARAM_FIXED;
		  startval=atof(ppt->value);
  }
  if((ppt=LALInferenceGetProcParamVal(state->commandLine,valopt))) startval=atof(ppt->value);
  else if(varytype!=LALINFERENCE_PARAM_FIXED) startval=min+(max-min)*gsl_rng_uniform(state->GSLrandom);

  /* Error checking */
  if(min>max) {
    fprintf(stderr,"ERROR: Prior for %s has min(%lf) > max(%lf)\n",name,min,max);
    exit(1);
  }
  if(startval<min || startval>max){
    fprintf(stderr,"ERROR: Initial value %lf for %s lies outwith prior (%lf,%lf)\n",startval,name,min,max);
    exit(1);
  }
  /* Mass parameters checks*/
  if (!strcmp(name,"eta"))
    if (max>0.25){
      fprintf(stderr,"ERROR: maximum of eta cannot be larger than 0.25. Check --eta-max\n");
      exit(1);
    }
  if (!strcmp(name,"q")){
    REAL8 qMin=min;
    REAL8 qMax=max;

    if (qMin <= 0.0 || qMin > 1.0)
    {
        fprintf(stderr,"ERROR: qMin must be between 0 and 1, got value qMin=%f\n",qMin);
		exit(1);
    }
    if (qMax > 1.0 || qMax <0.0 || qMax < qMin)
    {
      fprintf(stderr,"ERROR: qMax must be between 0 and 1, and qMax > qMin. Got value qMax=%f, qMin=%f\n",qMax,qMin);
	  exit(1);
    }
  }
  /*End of mass parameters check */

  LALInferenceAddVariable(var,name,&startval,LALINFERENCE_REAL8_t,varytype);
  LALInferenceAddMinMaxPrior(state->priorArgs, name, &min, &max, LALINFERENCE_REAL8_t);

}


/* Setup the variables to control template generation for the CBC model */
/* Includes specification of prior ranges. Returns address of new LALInferenceVariables */

LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
  char help[]="\
    ----------------------------------------------\n\
    --- Injection Arguments ----------------------\n\
    ----------------------------------------------\n\
    (--inj injections.xml) Injection XML file to use\n\
    (--event N)            Event number from Injection XML file to use\n\
\n\
    ----------------------------------------------\n\
    --- Template Arguments -----------------------\n\
    ----------------------------------------------\n\
    (--use-eta)            Jump in symmetric mass ratio eta, instead of q=m1/m2 (m1>m2)\n\
    (--approx)             Specify a template approximant and phase order to use\n\
                         (default TaylorF2threePointFivePN). Available approximants:\n\
                         default modeldomain=\"time\": GeneratePPN, TaylorT1, TaylorT2,\n\
                                                       TaylorT3, TaylorT4, EOB, EOBNR,\n\
                                                       EOBNRv2, EOBNRv2HM, SEOBNRv1,\n\
                                                       SpinTaylor, SpinQuadTaylor, \n\
                                                       SpinTaylorFrameless, SpinTaylorT4,\n\
                                                       PhenSpinTaylorRD, NumRel.\n\
                         default modeldomain=\"frequency\": TaylorF1, TaylorF2, TaylorF2RedSpin,\n\
                                                       TaylorF2RedSpinTidal, IMRPhenomA,\n\
                                                       IMRPhenomB, IMRPhenomP, IMRPhenomPv2.\n\
    (--amporder PNorder)            Specify a PN order in amplitude to use (defaults: LALSimulation: max available; LALInspiral: newtownian).\n\
    (--fref f_ref)                  Specify a reference frequency at which parameters are defined (default 100).\n\
    (--tidal)                   Enables tidal corrections, only with LALSimulation.\n\
    (--tidalT)                  Enables reparmeterized tidal corrections, only with LALSimulation.\n\
    (--4PolyEOS)                Enables 4-piece polytropic EOS parmeterization, only with LALSimulation.\n\
    (--4SpectralDecomp)         Enables 4-coeff. spectral decomposition EOS parmeterization, only with LALSimulation.\n\
    (--eos   EOS)               Fix the neutron star EOS. Use \"--eos help\" for allowed names\n\
    (--BinaryLove)              Enable the Binary Neutron Star common EoS tidal model, expressed through EOS-insensitive relations\n\
    (--spinOrder PNorder)           Specify twice the PN order (e.g. 5 <==> 2.5PN) of spin effects to use, only for LALSimulation (default: -1 <==> Use all spin effects).\n\
    (--tidalOrder PNorder)          Specify twice the PN order (e.g. 10 <==> 5PN) of tidal effects to use, only for LALSimulation (default: -1 <==> Use all tidal effects).\n\
    (--numreldata FileName)         Location of NR data file for NR waveforms (with NR_hdf5 approx).\n\
    (--modeldomain)                 domain the waveform template will be computed in (\"time\" or \"frequency\"). If not given will use LALSim to decide\n\
    (--spinAligned or --aligned-spin)  template will assume spins aligned with the orbital angular momentum.\n\
    (--singleSpin)                  template will assume only the spin of the most massive binary component exists.\n\
    (--noSpin, --disable-spin)      template will assume no spins (giving this will void spinOrder!=0) \n\
    (--no-detector-frame)              model will NOT use detector-centred coordinates and instead RA,dec\n\
    (--nonGR_alpha value) this is a LIV parameter which should only be passed when log10lambda_eff/lambda_eff is passed as a grtest-parameter for LIV test\n\
    (--LIV_A_sign) this is a LIV parameter determining if +A or -A is being tested; A occurs in the modified dispersion relation. LIV_A_sign has to be either +1 or -1 \n\
    (--liv) this flag is now needed for launching a LIV run\n\
    (--grtest-parameters dchi0,..,dxi1,..,dalpha1,..) template will assume deformations in the corresponding phase coefficients.\n\
    (--ppe-parameters aPPE1,....     template will assume the presence of an arbitrary number of PPE parameters. They must be paired correctly.\n\
    (--modeList lm,l-m...,lm,l-m)           List of modes to be used by the model. The chosen modes ('lm') should be passed as a ',' seperated list.\n\
    (--phenomXHMMband float)     Threshold parameter for the Multibanding of the non-precessing hlm modes in IMRPhenomXHM and IMRPhenomXPHM. If set to 0 then do not use multibanding.\n\
                                 Options and default values can be found in https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom_x__c.html\n\
    (--phenomXPHMMband float)    Threshold parameter for the Multibanding of the Euler angles in IMRPhenomXPHM. If set to 0 then do not use multibanding.\n\
                                 Options and default values can be found in https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom_x__c.html\n\
    (--phenomXPFinalSpinMod int) Change version for the final spin model used in IMRPhenomXP/IMRPhenomXPHM.\n\
                                 Options and default values can be found in https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom_x__c.html\n\
    (--phenomXPrecVersion int)   Change version of the Euler angles for the twisting-up of IMRPhenomXP/IMRPhenomXPHM.\n\
                                 Options and default values can be found in https://lscsoft.docs.ligo.org/lalsuite/lalsimulation/group___l_a_l_sim_i_m_r_phenom_x__c.html\n\
\n\
    ----------------------------------------------\n\
    --- Starting Parameters ----------------------\n\
    ----------------------------------------------\n\
    You can generally have MCMC chains to start from a given parameter value by using --parname VALUE. Names currently known to the code are:\n\
     time                         Waveform time (overrides random about trigtime).\n\
     chirpmass                    Chirpmass\n\
     eta                          Symmetric massratio (needs --use-eta)\n\
     q                            Asymmetric massratio (a.k.a. q=m2/m1 with m1>m2)\n\
     phase                        Coalescence phase.\n\
     costheta_jn                  Cosine of angle between J and line of sight [rads]\n\
     logdistance                  Log Distance (requires --use-logdistance)\n\
     rightascension               Rightascensions\n\
     declination                  Declination.\n\
     polarisation                 Polarisation angle.\n\
    * Spin Parameters:\n\
     a_spin1                      Spin1 magnitude\n\
     a_spin2                      Spin2 magnitude\n\
     tilt_spin1                   Angle between spin1 and orbital angular momentum\n\
     tilt_spin2                   Angle between spin2 and orbital angular momentum \n\
     phi_12                       Difference between spins' azimuthal angles \n\
     phi_jl                       Difference between total and orbital angular momentum azimuthal angles\n\
    * Equation of State parameters:\n\
     (requires --tidal)\n\
     lambda1                      lambda1.\n\
     lambda2                      lambda2.\n\
     (requires --tidalT)\n\
     lambdaT                      lambdaT.\n\
     dLambdaT                     dLambdaT.\n\
     (requires --4PolyEOS)\n\
     logp1                        logp1.\n\
     gamma1                       gamma1.\n\
     gamma2                       gamma2.\n\
     gamma3                       gamma3.\n\
     (requires --BinaryLove)\n\
     lambdaS                      Symmetric tidal deformability.\n\
     (requires --4SpectralDecomp):\n\
     SDgamma0                     SDgamma0.\n\
     SDgamma1                     SDgamma1.\n\
     SDgamma2                     SDgamma2.\n\
     SDgamma3                     SDgamma3.\n\
    * \n\
    * Spin-induced quadrupole moment test parameters:\n\
     (requires --dQuadMon12)\n\
     dQuadMon1                      dQuadMon1.\n\
     dQuadMon2                      dQuadMon2.\n\
     (requires --dQuadMonSA) \n\
     dQuadMonS                      dQuadMonS.\n\
     dQuadMonA                      dQuadMonA.\n\
    (dQuadMonS and dQuadMonA are the symmetric and antisymmetric combinations of dQuadMon1 and dQuadMon2).\n\
    ----------------------------------------------\n\
    --- Prior Ranges -----------------------------\n\
    ----------------------------------------------\n\
    You can generally use --paramname-min MIN --paramname-max MAX to set the prior range for the parameter paramname\n\
    The names known to the code are listed below.\n\
    Component masses, total mass and time have dedicated options listed here:\n\n\
    (--trigtime time)                       Center of the prior for the time variable.\n\
    (--comp-min min)                        Minimum component mass (1.0).\n\
    (--comp-max max)                        Maximum component mass (100.0).\n\
    (--mass1-min min, --mass1-max max)      Min and max for mass1 (default: same as comp-min,comp-max, will over-ride these.\n\
    (--mass2-min min, --mass2-max max)      Min and max for mass2 (default: same as comp-min,comp-max, will over-ride these.\n\
    (--mtotal-min min)                      Minimum total mass (2.0).\n\
    (--mtotal-max max)                      Maximum total mass (200.0).\n\
    (--dt time)                             Width of time prior, centred around trigger (0.2s).\n\
\n\
    (--ns-max-mass)                          Maximum observed NS mass used for EOS prior (none).\n\
    (--varyFlow, --flowMin, --flowMax)       Allow the lower frequency bound of integration to vary in given range.\n\
    (--pinparams)                            List of parameters to set to injected values [mchirp,asym_massratio,etc].\n\
    ----------------------------------------------\n\
    --- Fix Parameters ---------------------------\n\
    ----------------------------------------------\n\
    You can generally fix a parameter to be fixed to a given values by using --fix-paramname VALUE\n\
    where the known names have been listed above.\n\
\n";

  /* Print command line arguments if state was not allocated */
  if(state==NULL)
    {
      fprintf(stdout,"%s",help);
      return(NULL);
    }

  /* Print command line arguments if help requested */
  if(LALInferenceGetProcParamVal(state->commandLine,"--help"))
    {
      fprintf(stdout,"%s",help);
      return(NULL);
    }

  LALStatus status;
  memset(&status,0,sizeof(status));
  int errnum;
  SimInspiralTable *injTable=NULL;
  LALInferenceVariables *priorArgs=state->priorArgs;
  LALInferenceVariables *proposalArgs=state->proposalArgs;
  ProcessParamsTable *commandLine=state->commandLine;
  ProcessParamsTable *ppt=NULL;
  ProcessParamsTable *ppt_order=NULL;
  LALPNOrder PhaseOrder=-1;
  LALPNOrder AmpOrder=-1;
  Approximant approx=NumApproximants;
  REAL8 f_ref = 100.0;
  LALInferenceIFOData *dataPtr;
  UINT4 event=0;
  UINT4 i=0;
  /* Default priors */
  REAL8 Dmin=1.0;
  REAL8 Dmax=2000.0;
  REAL8 Dinitial = (Dmax + Dmin)/2.0;
  REAL8 mcMin=1.0;
  REAL8 mcMax=15.3;
  REAL8 etaMin=0.0312;
  REAL8 etaMax=0.25;
  REAL8 qMin=1./30.; // The ratio between min and max component mass (see InitMassVariables)
  REAL8 qMax=1.0;
  REAL8 psiMin=0.0,psiMax=LAL_PI;
  REAL8 decMin=-LAL_PI/2.0,decMax=LAL_PI/2.0;
  REAL8 raMin=0.0,raMax=LAL_TWOPI;
  REAL8 phiMin=0.0,phiMax=LAL_TWOPI;
  REAL8 costhetaJNmin=-1.0 , costhetaJNmax=1.0;
  REAL8 dt=0.1;  /* Half the width of time prior */
  REAL8 lambda1Min=0.0;
  REAL8 lambda1Max=3000.0;
  REAL8 lambda2Min=0.0;
  REAL8 lambda2Max=3000.0;
  REAL8 lambdaTMin=0.0;
  REAL8 lambdaTMax=3000.0;
  REAL8 dLambdaTMin=-500.0;
  REAL8 dLambdaTMax=500.0;
  REAL8 dQuadMonMin=-200.0;
  REAL8 dQuadMonMax=200.0;
  REAL8 logp1Min=33.6;
  REAL8 logp1Max=35.4;
  REAL8 gamma1Min=2.0;
  REAL8 gamma1Max=4.5;
  REAL8 gamma2Min=1.1;
  REAL8 gamma2Max=4.5;
  REAL8 gamma3Min=1.1;
  REAL8 gamma3Max=4.5;
  REAL8 lambdaSMin=0.0;
  REAL8 lambdaSMax=5000.0;
  REAL8 SDgamma0Min=0.2;
  REAL8 SDgamma0Max=2.0;
  REAL8 SDgamma1Min=-1.6;
  REAL8 SDgamma1Max=1.7;
  REAL8 SDgamma2Min=-0.6;
  REAL8 SDgamma2Max=0.6;
  REAL8 SDgamma3Min=-0.02;
  REAL8 SDgamma3Max=0.02;
  gsl_rng *GSLrandom=state->GSLrandom;
  REAL8 endtime=0.0, timeParam=0.0;
  REAL8 timeMin=endtime-dt,timeMax=endtime+dt;
  REAL8 zero=0.0; /* just a number that will be overwritten anyway*/

  /* Over-ride prior bounds if analytic test */
  if (LALInferenceGetProcParamVal(commandLine, "--correlatedGaussianLikelihood"))
  {
    return(LALInferenceInitModelReviewEvidence(state));
  }
  else if (LALInferenceGetProcParamVal(commandLine, "--bimodalGaussianLikelihood"))
  {
    return(LALInferenceInitModelReviewEvidence_bimod(state));
  }
  else if (LALInferenceGetProcParamVal(commandLine, "--rosenbrockLikelihood"))
  {
    return(LALInferenceInitModelReviewEvidence(state)); /* CHECKME: Use the default prior for unimodal */
  }

  LALInferenceModel *model = XLALMalloc(sizeof(LALInferenceModel));
  model->params = XLALCalloc(1, sizeof(LALInferenceVariables));
  memset(model->params, 0, sizeof(LALInferenceVariables));
  model->eos_fam = NULL;

  UINT4 signal_flag=1;
  ppt = LALInferenceGetProcParamVal(commandLine, "--noiseonly");
  if(ppt)signal_flag=0;
  LALInferenceAddVariable(model->params, "signalModelFlag", &signal_flag,  LALINFERENCE_INT4_t,  LALINFERENCE_PARAM_FIXED);

  /* Read injection XML file for parameters if specified */
  ppt=LALInferenceGetProcParamVal(commandLine,"--inj");
  if(ppt){
    SimInspiralTableFromLIGOLw(&injTable,ppt->value,0,0);
    if(!injTable){
      fprintf(stderr,"Unable to open injection file %s\n",ppt->value);
      exit(1);
    }
    ppt=LALInferenceGetProcParamVal(commandLine,"--event");
    if(ppt){
      event= atoi(ppt->value);
      fprintf(stderr,"Reading event %d from file\n",event);
      i=0;
      while(i<event) {i++; injTable=injTable->next;} /* select event */
    }
  }

  /* See if there are any parameters pinned to injection values */
  if((ppt=LALInferenceGetProcParamVal(commandLine,"--pinparams"))){
    char *pinned_params=ppt->value;
    LALInferenceVariables tempParams;
    memset(&tempParams,0,sizeof(tempParams));
    char **strings=NULL;
    UINT4 N;
    LALInferenceParseCharacterOptionString(pinned_params,&strings,&N);
    LALInferenceInjectionToVariables(injTable,&tempParams);
    LALInferenceVariableItem *node=NULL;
    while(N>0){
      N--;
      char *name=strings[N];
      fprintf(stdout,"Pinning parameter %s\n",node->name);
      node=LALInferenceGetItem(&tempParams,name);
      if(node) LALInferenceAddVariable(model->params,node->name,node->value,node->type,node->vary);
      else {fprintf(stderr,"Error: Cannot pin parameter %s. No such parameter found in injection!\n",node->name);}
    }
  }

  /* Over-ride approximant if user specifies */
  ppt=LALInferenceGetProcParamVal(commandLine,"--approximant");
  if(ppt){
    approx = XLALGetApproximantFromString(ppt->value);
    ppt_order=LALInferenceGetProcParamVal(commandLine,"--order");
    if(ppt_order) PhaseOrder = XLALGetOrderFromString(ppt_order->value);
  }
  ppt=LALInferenceGetProcParamVal(commandLine,"--approx");
  if(ppt){
    approx = XLALGetApproximantFromString(ppt->value);
    XLAL_TRY(PhaseOrder = XLALGetOrderFromString(ppt->value),errnum);
    if( (int) PhaseOrder == XLAL_FAILURE || errnum) {
      PhaseOrder=-1;
    }
  }