BankEfficiency.c

Go to the documentation of this file.

/*
*  Copyright (C) 2007 Thomas Cokelaer
*
*  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., 59 Temple Place, Suite 330, Boston,
*  MA  02111-1307  USA
*/


#include "BankEfficiency.h"
/* --- version information ------------------------------------------------ */
NRCSID( BANKEFFICIENCYC, "$Id: BankEfficiency.c,v 1.142 2008/07/07 17:36:46 devanka Exp $");
RCSID(  "$Id: BankEfficiency.c,v 1.142 2008/07/07 17:36:46 devanka Exp $");

#define CVS_ID_STRING_C      "$Id: BankEfficiency.c,v 1.142 2008/07/07 17:36:46 devanka Exp $"
#define CVS_REVISION_C       "$Revision: 1.142 $"
#define CVS_NAME_STRING_C    "$Name:  $"
#define CVS_SOURCE_C         "$Source: /usr/local/cvs/lscsoft/lalapps/src/findchirp/BankEfficiency.c,v $"
#define CVS_DATE_C           "$Date: 2008/07/07 17:36:46 $"
#define PROGRAM_NAME         "BankEfficiency"

/* --- global variables --------------------------------------------------- */
RandomParams *randParams = NULL;
INT4          randnStartPad = 0;  /* injections always at the same time if zero*/
INT4          ascii2xml = 0;
  
/* --- Main program ------------------------------------------------------- */ 
int
main (INT4 argc, CHAR **argv )
{  
  /* --- Local variables -------------------------------------------------- */  
  INT4                   i;
  INT4                   j;
  INT4                   thisTemplateIndex;
  Order                  tempOrder;  /* temporary phase order */

  /* --- the general input structure --- */
  UserParametersIn       userParam;

  /* --- signal related --- */
  REAL4Vector            signal;
  RandomInspiralSignalIn randIn;
  
  /* --- template bank related --- */
  Mybank                 mybank;  
  InspiralTemplate       insptmplt;
  InspiralCoarseBankIn   coarseBankIn;
  INT4                   sizeBank = 0;
  MetadataTable          templateBank;
  SnglInspiralTable      *tmpltHead = NULL;

  /* --- filtering related --- */
  REAL4Vector               correlation;
  BankEfficiencyBCV         bankefficiencyBCV;
  InspiralWaveOverlapIn     overlapin;      

  /* --- results and data mining --- */
  OverlapOutputIn        overlapOutputThisTemplate;
  OverlapOutputIn        overlapOutputBestTemplate;

  /* --- fft related --- */
  RealFFTPlan           *fwdp = NULL;
  RealFFTPlan           *revp = NULL;

  /* --- others ---*/
  LALStatus              status = blank_status;
  FILE                  *Foutput;

  /* --- fast option related --- */
  UINT4                   fast_index;
      
  /* --- eccentricity related to the bank */
  REAL4                  tempEccentricity = 0;

  /* --- ambiguity function and statistics --- */
  gsl_histogram         *histogramNoise = gsl_histogram_alloc (200);  
  gsl_matrix            *amb1;

  /* --- parameter for the simulation --- */
  BankEfficiencySimulation   simulation;  
  
  /* --- Some initialization --- */ 
  gsl_histogram_set_ranges_uniform (histogramNoise, 0.,20.);
  lal_errhandler = LAL_ERR_EXIT;
  lalDebugLevel = 0;
  templateBank.snglInspiralTable = NULL;
  
  /* --- Initialization of structures related to all the user parameters --- */
  BankEfficiencyParametersInitialization(&coarseBankIn, &randIn, &userParam);
  
  /* --- Read user parameters --- */
  BankEfficiencyParseParameters(&argc, argv,&coarseBankIn, &randIn, &userParam);

  /* --- Check input parameters --- */
  BankEfficiencyUpdateParams(&coarseBankIn, &randIn, &userParam);
  
  /* eccentric bank initialisation */
  BankEfficiencyEccentricBankInit(&userParam);
  
  /* --- Init a random structure using the user seed --- */  
  randParams = XLALCreateRandomParams(randIn.useed );
  
  /* --- Ascii2XML option --- */
  /* This is a call to the function that converts the ASCII output of
   * BankEfficiency into a standard XML output. Useful when several ASCII
   * output are available and needs to be converted all together into a single
   * XML file.
   * */
  if (ascii2xml == 1){
    BankEfficiencyAscii2Xml();    
  }
  
  /* --- Condor related --- */
  /* If one want to use condor script to keep track of the input parameters, 
   * we first need to create a prototype containing all parameters, which can 
   * be read later using --ascii2xml option. Use with bep.dag DAGMan script.  
  */
  if (userParam.printPrototype){
    BankEfficiencyPrintProtoXml(coarseBankIn, randIn, userParam);   
  }
  
  /* --- Allocate memory -------------------------------------------------- */
  /* --- First, estimate the size of the signal --- */
  LAL_CALL(BankEfficiencyGetMaximumSize(&status, 
      randIn, coarseBankIn, userParam, &(signal.length)), 
      &status);  
       
  /* --- Set size of the other vectors --- */
  randIn.psd.length     = signal.length/2 + 1; 
  correlation.length    = signal.length;
  /* --- Allocate memory for some vectors --- */
  signal.data      = (REAL4*)LALCalloc(1, sizeof(REAL4) * signal.length);
  correlation.data = (REAL4*)LALCalloc(1, sizeof(REAL4) * correlation.length);
  randIn.psd.data  = (REAL8*)LALCalloc(1, sizeof(REAL8) * randIn.psd.length);
  /* --- Allocate memory for BCV filtering only --- */ 
  if (userParam.template == BCV)
  {
        INT4 n;
        n = signal.length;
    bankefficiencyBCV.FilterBCV1.length = n;
    bankefficiencyBCV.FilterBCV2.length = n;
    bankefficiencyBCV.FilterBCV1.data = (REAL4*)LALCalloc(1, sizeof(REAL4) * n);
    bankefficiencyBCV.FilterBCV2.data = (REAL4*)LALCalloc(1, sizeof(REAL4) * n);
  }
  
  /* --- Create the PSD noise ----------------------------------------------- */
  LAL_CALL(BankEfficiencyCreatePsd(&status, &coarseBankIn, &randIn, userParam), 
       &status);
 
  /* --- Create the template bank ------------------------------------------- */
  LAL_CALL(BankEfficiencyCreateTemplateBank(&status, &coarseBankIn,
      &templateBank, &tmpltHead, userParam, &randIn, &sizeBank), &status);

  /* --- populate MyBank strucuture --- */
  BankEfficiencyInitMyBank(&mybank, &sizeBank, tmpltHead, userParam);    
       
  /* --- Allocate memory for the ambiguity function ------------------------- */
  /* For each bank, we have t0,t3,and the maximum */
  amb1 = gsl_matrix_alloc(4,sizeBank); /*allocated t0,t3,max*/
  for (i=0;i<4; i++)
  {
    for (j=0;j<sizeBank; j++)
    {
      gsl_matrix_set(amb1, i, j, 0.); /* set it to zero */
    }
  }      
      
  /* --- Estimate the fft's plans ----------------------------------------- */
  LALCreateForwardRealFFTPlan(&status, &fwdp, signal.length, 0);
  LALCreateReverseRealFFTPlan(&status, &revp, signal.length, 0);
  
  /* --- The overlap structure -------------------------------------------- */
  overlapin.nBegin      = 0;
  overlapin.nEnd        = 0;
  overlapin.psd         = randIn.psd;
  overlapin.fwdp        = randIn.fwdp = fwdp;
  overlapin.revp        = revp;
  overlapin.ifExtOutput = 0;  /* new option from Anand's WaveOverlap */
  overlapin.ifCorrelationOutput = 1; /* output of WaveOverlap is sqrt(x^2+y^2)*/

  /* --- Create Matrix for BCV Maximization once for all ------------------ */
  if ((userParam.template == BCV) )
  {
    LAL_CALL( BankEfficiencyCreatePowerVector(&status, 
       &(bankefficiencyBCV.powerVector), randIn, signal.length),  &status);

    BankEfficiencyCreateBCVMomentVector(
        &(bankefficiencyBCV.moments), &coarseBankIn.shf,randIn.param.tSampling, 
        randIn.param.fLower, signal.length);
  }
  
  /* ------------------------------------------------------------------------ */
  /* --- Main loop ---------------------------------------------------------- */
  /* ------------------------------------------------------------------------ */  
  /* --- ntrials is used to keep track of the simulation trials           --- */
  simulation.ntrials = 0;
  simulation.N = userParam.ntrials;
  /* --- loop over the number of simulations (ntrials) --- */  
  while (++simulation.ntrials <= userParam.ntrials) 
  {     
    if (vrbflg){
      fprintf(stdout,"Simulation number %d/%d\n", simulation.ntrials, userParam.ntrials);
    }
      
    /* --- initialisation for each simulation --- */
    BankEfficiencyInitOverlapOutputIn(&overlapOutputBestTemplate);
      
    /* --- set the fcutoff (signal,overlap) --- */
    randIn.param.fCutoff = userParam.signalfFinal; 
      
    /* --- generate the signal waveform --- */
    LAL_CALL(BankEfficiencyGenerateInputData(&status, 
        &signal, &randIn, userParam), &status);
    
    /* --- populate the main structure of the overlap ---*/
    overlapin.signal = signal;
       
    /*  --- populate the insptmplt with the signal parameter --- */ 
    insptmplt = randIn.param; /* set the sampling and other common parameters */     
    insptmplt.approximant  = userParam.template;        /* set the template   */
    insptmplt.order        = coarseBankIn.order;        /* set the phaseorder */ 
    insptmplt.ampOrder     = coarseBankIn.ampOrder;     /* set the amp order  */

    /* --- reset the bank (none of the template have been used) --- */
    for (i=0; i<sizeBank; i++){
      mybank.used[i] = 0;
    }
      
    /* --- set the general parameter for each simulation --- */
    simulation.filter_processed = 0; /* count number of templates being used  */ 
    simulation.lastEMatch       = 1; /* ematch of the last template being used*/
    simulation.filteringIndex   = 0; /* current filter number                 */
    simulation.stop             = 0; /* flag to stop this simulation          */
    simulation.bestSNR          = 0; /* value of the best SNR                 */
    simulation.bestSNRIndex     = 0; /* index of the template that has the best SNR*/
    simulation.SNRMax           = 0; /* For how long have we the best SNR     */
    simulation.fastParam1 = userParam.fastParam1;
    simulation.fastParam1 *= 1.+log10(userParam.eccentricBank.bins);
    simulation.eMatch           = userParam.eMatch; /* reset the ematch       */
    simulation.bestEMatch       = -1e16; /* reset ematch to small value*/
    simulation.randIn           = randIn;        

    /* --- loop over the bank until counter reaches size bank --- */ 
    while ( simulation.filteringIndex<sizeBank && simulation.stop == 0 )    
    {
      /* --- iterate the counter --- */
      simulation.filteringIndex++;
      
      /* -- reset the output structure --- */
      BankEfficiencyInitOverlapOutputIn(&overlapOutputThisTemplate);
    
      GetClosestValidTemplate(mybank, simulation.randIn, &fast_index);
          
      LAL_CALL(BankEfficiencyCreateListFromTmplt(&status, 
             &insptmplt, mybank, fast_index), &status); 
      thisTemplateIndex = fast_index;
                
      switch(userParam.template)
      {                 
        case BCV:
            
          insptmplt.massChoice = psi0Andpsi3;  
          LAL_CALL(LALInspiralParameterCalc(&status, &(insptmplt)),&status);
          overlapin.param = insptmplt;

          LAL_CALL(LALInspiralParameterCalc( &status,  &(overlapin.param) ), 
              &status);

          /* if faithfulness is required, the template parameters are 
           * identical to the signal except for the name of the approximant.
           *  */
          if (userParam.faithfulness)
          {
            insptmplt                    = randIn.param;
            overlapin.param              = randIn.param;
            overlapin.param.approximant  = userParam.template; 
            insptmplt.fFinal             = randIn.param.fFinal;
            /* --- the faithfulness requires only one filtering, so we stop - */
            simulation.stop              = 1;                
           }
     
          LAL_CALL(BankEfficiencyInspiralOverlapBCV(&status, 
                    &insptmplt, userParam, &randIn, 
                    &overlapin, &overlapOutputThisTemplate,
                    &correlation, &bankefficiencyBCV),
                    &status);

          /* keep track of what has been done in the overlap function */
            
          overlapOutputThisTemplate.freq =  overlapin.param.fFinal;
          overlapOutputThisTemplate.templateNumber = thisTemplateIndex;
          simulation.filter_processed++;
          break;
              
        case AmpCorPPN:
        case TaylorT1: 
        case Eccentricity:
        case TaylorT2:
        case TaylorT3:
        case TaylorF1:
        case TaylorF2:
        case EOB:
        case EOBNR:
        case PadeT1:
        case PadeF1:
        case SpinTaylor:
          if (vrbflg){
            fprintf(stderr,"closeest template is tau0= %f tau3= %f  index=%d filterintIndex= %d,",   
              mybank.tau0[fast_index], mybank.tau3[fast_index],fast_index,simulation.filteringIndex);
          }
                  
                    
          /* --- first we create the template --- */
          insptmplt.massChoice = t03;
          LAL_CALL(LALInspiralParameterCalc( &status,  &(insptmplt) ), 
              &status);
          /* --- set up the overlapin strucutre --- */
          overlapin.param = insptmplt;
          LAL_CALL(LALInspiralParameterCalc( &status,  &(overlapin.param) ), 
              &status);
          overlapin.param.fCutoff = coarseBankIn.fUpper;
          overlapin.param.fLower  = coarseBankIn.fLower;
          overlapin.param.fFinal  = randIn.param.tSampling/2. - 1;

          /* --- the fast simulation option. We compute the ematch --- */
          simulation.lastEMatch = BankEfficiencyComputeEMatch(&randIn, mybank,
              thisTemplateIndex);
          if (simulation.lastEMatch > simulation.bestEMatch)
            simulation.bestEMatch = simulation.lastEMatch;
          
          /* --- if check option is on --- */    
          if (userParam.faithfulness)
          {
            simulation.lastEMatch = 1; /* if faithfulness, masses are equal so match is 1*/                
            tempOrder = insptmplt.order;
            tempEccentricity = insptmplt.eccentricity;
            /* --- now, we overwrite insptmplt with the signal params---*/
            insptmplt = randIn.param;
            overlapin.param = randIn.param;
            /* --- but get back some parameters --- */
            LAL_CALL(LALInspiralParameterCalc( &status,  &(overlapin.param) ), &status);
            overlapin.param.fCutoff = coarseBankIn.fUpper;
            overlapin.param.fFinal = randIn.param.tSampling/2. - 1;
            overlapin.param.approximant        = userParam.template;
            overlapin.param.fLower  = coarseBankIn.fLower;
            /* we want to keep the original order except for eccentricity 
             * where order must be zero anyway*/
            overlapin.param.order              = tempOrder;
            overlapin.param.eccentricity       = tempEccentricity;
            insptmplt.eccentricity             = tempEccentricity;
            /* --- the faithfulness requires only one filtering, so we stop - */
            simulation.stop = 1;                
          }          
          /* if we want to cut integration before the fFinal */
          /*              overlapin.param.fCutoff = 1023;    */          
      
      
          if (userParam.fastSimulation == Heuristic1)
          {                        
            /* --- decrease the ematch parameter if needed --- */            
            if (simulation.lastEMatch < simulation.eMatch){
              simulation.eMatch = simulation.lastEMatch-1;
            }
            /* --- if SNRMax has reached the max iteration, we stop --- */
            if (simulation.SNRMax >= simulation.fastParam1 
               && userParam.fastSimulation){
              simulation.stop = 1;
            }
          }                        
            
          /* --- filteing!=1 ensure that at least 1 filter will be used    ---*/
          /* --- while lastEMatch is greater than the ematch --- */          
          if ((simulation.filteringIndex !=1) &&
              (userParam.fastSimulation != None) 
              && (simulation.lastEMatch < simulation.eMatch ))
          {                     
            gsl_matrix_set(amb1,2,thisTemplateIndex, 0); 
            gsl_matrix_set(amb1,3,thisTemplateIndex, 0); 
            gsl_matrix_set(amb1,0,thisTemplateIndex, insptmplt.t0); 
            gsl_matrix_set(amb1,1,thisTemplateIndex, insptmplt.t3);
            /* --- if we enter here, then we should stop the simulation --- */
            simulation.stop = 1;
          }       
          else
          {                          
            LAL_CALL(BankEfficiencyWaveOverlap(&status, 
                  &correlation,
                  &overlapin,
                  &overlapOutputThisTemplate,
                   randIn.param.nStartPad), &status);
                                
                
            overlapOutputThisTemplate.templateNumber = thisTemplateIndex;
            /* we compute the averaged ambiguity function a t=ta and 
             * the averaged maximizaed ambiguity function over time*/
            BankEfficiencyPopulateAmbiguityFunction(
                amb1,correlation,thisTemplateIndex, 
                overlapOutputThisTemplate, randIn.param.nStartPad,
                insptmplt);
                
            /* is it correct to do so ? */    
            insptmplt.fFinal = overlapin.param.fFinal;
            simulation.filter_processed++;
            
            if (vrbflg){
            fprintf(stderr, "snr = %f ematch=%f (threshold is %f) snrMax for %d\n",
                overlapOutputThisTemplate.rhoMax, simulation.lastEMatch,
                simulation.eMatch,simulation.SNRMax);
                
            fflush(stderr);
            }
            
            if (overlapOutputThisTemplate.rhoMax > simulation.bestSNR)
            {              
              simulation.bestSNR  = overlapOutputThisTemplate.rhoMax;
              simulation.bestSNRIndex = simulation.filteringIndex; 
              simulation.SNRMax = 0;
              if (userParam.fastSimulation == Heuristic1)
              {
                simulation.randIn.param.t0 = insptmplt.t0;
                simulation.randIn.param.t3 = insptmplt.t3;
              }
              simulation.bestEMatch      = simulation.lastEMatch;                            
            }
            else
            {
              simulation.SNRMax++;
            }
                                        
          }
        break;
      } /*end of the switch over template*/
      
      
      
      /* accumulates histogram of the correlations over all templates and all 
       * simulations. */
      if (userParam.printSNRHisto){      
        BankEfficiencyUpdateSNRHistogram(&correlation, histogramNoise);
      }
        
      /* for the final results, we keep only the loudest SNR over the entire
       * template bank, for each simulations.*/ 
      BankEfficiencyKeepHighestValues( overlapOutputThisTemplate, 
          &overlapOutputBestTemplate, insptmplt);
                      
    }  /* --- end of  bank process --- */
      
    /* --- print the results on stdout and xml file if requested --- */
    LAL_CALL(BankEfficiencyFinalise(&status,mybank,
        overlapOutputBestTemplate,randIn,userParam,simulation, 
        coarseBankIn), &status);
                  
    if (userParam.template == BCV) 
    {
      if (userParam.printBestOverlap || userParam.printBestTemplate) 
      {
        userParam.extraFinalPrinting = 1;
        LAL_CALL(BankEfficiencyInspiralOverlapBCV(&status, 
                    &insptmplt, userParam, &randIn,
                     &overlapin, &overlapOutputThisTemplate, &correlation, 
                     &bankefficiencyBCV),
                     &status);    
      }
    }

    /* --- print correlation of the last computed correlation --- */
    if (userParam.printBestOverlap ){
        BankEfficiencySaveVector("BankEfficiency-correlation.dat", correlation,
            randIn.param.tSampling);
    }
        
  }  /* --- end while(trial) --- */
  

  /* --- save histrogram of all the SNR output correlations --- */
  if (userParam.printSNRHisto)
  {
    Foutput = fopen(BANKEFFICIENCY_SNRHISTOGRAM,"w");
    gsl_histogram_fprintf(Foutput, histogramNoise, "%f", "%g");
    fclose(Foutput);
  }


  /* --- we save the ambiguity in a file --- */
  BankEfficiencyPrintAmbiguity(userParam,sizeBank,amb1);
  
  /* free memory */
  while ( templateBank.snglInspiralTable )
  {
    tmpltHead = templateBank.snglInspiralTable;
    templateBank.snglInspiralTable = templateBank.snglInspiralTable->next;
    LALFree( tmpltHead );
  }
  
  /* --- destroy the plans, correlation and signal --- */
  
  /*XLALDestroyRandomParams(randParams );*/
  
  if (userParam.template == BCV)
  {
    LALFree(bankefficiencyBCV.powerVector.fm5_3.data);
    LALFree(bankefficiencyBCV.powerVector.fm2_3.data);
    LALFree(bankefficiencyBCV.powerVector.fm1_2.data);
    LALFree(bankefficiencyBCV.powerVector.fm7_6.data);
    LALFree(bankefficiencyBCV.moments.a11.data);
    LALFree(bankefficiencyBCV.moments.a21.data);
    LALFree(bankefficiencyBCV.moments.a22.data);   
  
    if (userParam.template == BCV)
    {  
      LALFree(bankefficiencyBCV.FilterBCV2.data);
      LALFree(bankefficiencyBCV.FilterBCV1.data);
    }
  }
  
  LALFree(randIn.psd.data);
  LALDDestroyVector( &status, &(coarseBankIn.shf.data) );
 
  LALFree(signal.data);
  LALFree(correlation.data);

  LALDestroyRealFFTPlan(&status,&fwdp);
  LALDestroyRealFFTPlan(&status,&revp);

  gsl_matrix_free(amb1);  
  gsl_histogram_free(histogramNoise);

  free(mybank.mass1);
  free(mybank.mass2);
  /*todo  free other varaible in mybank*/
  
  LALCheckMemoryLeaks();    
    
  return 0;
}



              
/* *************************************************************************
 * Some output Results
 * ********************************************************************** */
void BankEfficiencyKeepHighestValues(
  OverlapOutputIn  resultThisTemplate,
  OverlapOutputIn *resultBestTemplate,
  InspiralTemplate insptmplt)
{
  
  if (resultThisTemplate.rhoMax > resultBestTemplate->rhoMax)
  {        
    resultBestTemplate->rhoMax         = resultThisTemplate.rhoMax;
    resultBestTemplate->phase          = resultThisTemplate.phase;
    resultBestTemplate->alpha          = resultThisTemplate.alpha;
    resultBestTemplate->rhoBin         = resultThisTemplate.rhoBin;
    resultBestTemplate->freq           = resultThisTemplate.freq;
    resultBestTemplate->templateNumber = resultThisTemplate.templateNumber;    
    resultBestTemplate->snrAtCoaTime   = resultThisTemplate.snrAtCoaTime;
    resultBestTemplate->eccentricity   = insptmplt.eccentricity;
  }
}


/* ************************************************************************
 * Some output Results
 * ********************************************************************** */
void BankEfficiencyGetResult(
  LALStatus        *status,
  InspiralTemplate *list,
  InspiralTemplate  injected,
  OverlapOutputIn   bestOverlap, 
  ResultIn         *result,
  UserParametersIn  userParam
)
{
  INT4 templateNumber;
  InspiralTemplate trigger;

  INITSTATUS (status, "GetResult", BANKEFFICIENCYC);
  ATTATCHSTATUSPTR(status);
     
  templateNumber = bestOverlap.templateNumber;

  trigger = *list;

  if (userParam.template == BCV ){

    /*    LALInspiralParameterCalc( status->statusPtr,  &trigger );
    CHECKSTATUSPTR(status);                         
    */
    result->psi0_inject  = injected.psi0;
    result->psi3_inject  = injected.psi3;     
    result->psi0_trigger = trigger.psi0;
    result->psi3_trigger = trigger.psi3;     
    result->tau0_trigger = 0.;
    result->tau3_trigger = 0.;
    result->tau0_inject  = 0.;
    result->tau3_inject  = 0.; 
  }
  else
  {
    LALInspiralParameterCalc( status->statusPtr,  &trigger );
    CHECKSTATUSPTR(status);                         
    result->psi0_inject  = 0.;
    result->psi3_inject  = 0.;     
    result->psi0_trigger = 0.;
    result->psi3_trigger = 0.;     
    result->tau0_trigger = trigger.t0;
    result->tau3_trigger = trigger.t3;
    result->tau0_inject  = injected.t0;
    result->tau3_inject  = injected.t3; 
    result->polarisationAngle  = injected.polarisationAngle;
    result->inclination  = injected.inclination;     
  }

  result->mass1_inject = injected.mass1;
  result->mass2_inject = injected.mass2;
  result->fend_inject  = injected.fFinal;
  result->fend_trigger = bestOverlap.freq;
  result->rho_final    = bestOverlap.rhoMax;
  result->alphaF       = bestOverlap.alpha*pow(bestOverlap.freq,2./3.);
  result->bin          = bestOverlap.rhoBin;
  result->phase        = bestOverlap.phase;
  result->layer        = bestOverlap.layer;
  result->phase        = bestOverlap.phase;
  result->snrAtCoaTime = bestOverlap.snrAtCoaTime;
  result->eccentricity = bestOverlap.eccentricity;
  
  DETATCHSTATUSPTR(status);
  RETURN (status);
}


void BankEfficiencyPrintResults(
  ResultIn                 result, 
  RandomInspiralSignalIn   randIn,
  BankEfficiencySimulation simulation)
{
  fprintf(stdout, 
  "%e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e  %e %e %e %e %d %d %d %d %f %d/%d\n",  
      result.psi0_trigger, result.psi3_trigger,
      randIn.param.psi0, randIn.param.psi3, 
      result.tau0_trigger, result.tau3_trigger,
      randIn.param.t0, randIn.param.t3, 
      result.eccentricity, randIn.param.eccentricity,randIn.param.alpha1,
      result.fend_trigger, randIn.param.fFinal,
      randIn.param.mass1,randIn.param.mass2,
      randIn.param.inclination,randIn.param.polarisationAngle,
      randIn.param.startPhase, result.rho_final, 
      result.snrAtCoaTime, result.phase, 
      result.alphaF,result.bin, randIn.param.nStartPad, result.nfast, 
      simulation.bestSNRIndex,simulation.bestEMatch,simulation.ntrials, simulation.N);

  fflush(stdout);
}


/* ****************************************************************************
 * Function to generate and stored the moments in  REAL4vectors. 
 * ************************************************************************* */
void BankEfficiencyCreateBCVMomentVector(
  BankEfficiencyMoments *moments,
  REAL8FrequencySeries  *psd,
  REAL8                 sampling,
  REAL8                 fLower,
  INT4                  length)
{
  REAL8 m7 = 0.;                            /* the moments */
  REAL8 m5 = 0.;    
  REAL8 m3 = 0.;
  REAL8 f;
  
  UINT4 kMin;
  UINT4 kMax;
  UINT4 k;
  
  InspiralMomentsIn     in;

  moments->a11.length =  length / 2 ;
  moments->a22.length =  length / 2 ;
  moments->a21.length =  length / 2 ;
  moments->a11.data   = (REAL4*) LALCalloc(1, sizeof(REAL4) * moments->a11.length);
  moments->a22.data   = (REAL4*) LALCalloc(1, sizeof(REAL4) * moments->a22.length);
  moments->a21.data   = (REAL4*) LALCalloc(1, sizeof(REAL4) * moments->a21.length);

  /* inspiral structure */
  in.shf    = psd;          /* The spectrum             */
  in.xmin   = fLower;       /* Lower frequency of integral? is it correct or 
                               should we put zero? */
  in.xmax   = sampling/2.;  /* We dont' need to cut here */
  in.norm   = 1./4.;        /* Some normalization */

  in.norm *= sampling * sampling;

  /* The minimum and maximum frequency */
  kMin = (UINT8) floor( (in.xmin - psd->f0) / psd->deltaF );
  kMax = (UINT8) floor( (in.xmax - psd->f0) / psd->deltaF );
     
  /* Skip frequency less than fLower.*/
  for (k = 0; k < kMin ; k++)
  {
    moments->a11.data[k] = 0.;
    moments->a21.data[k] = 0.;
    moments->a22.data[k] = 0.;
  }
  for ( k = kMin; k < kMax; ++k )
  {
    f = psd->f0 + (REAL8) k * psd->deltaF;
    if (psd->data->data[k])
    {
      m7 +=  pow( f, -(7./3.) ) / psd->data->data[k] * psd->deltaF / in.norm;
      m5 +=  pow( f, -(5./3) )  / psd->data->data[k] * psd->deltaF / in.norm;
      m3 +=  pow( f, -(1.) )    / psd->data->data[k] * psd->deltaF / in.norm;
      moments->a11.data[k] = 1./sqrt(m7);
      moments->a22.data[k] = 1. / sqrt(m3 - m5*m5/m7);
      moments->a21.data[k] = -m5 / m7 * moments->a22.data[k];
    }
  }
  BankEfficiencySaveVector("moment1.dat",moments->a11,sampling);
}


/* ****************************************************************************
 * Create an orthogonal filter. by taking its complex conjuguate
 * ************************************************************************** */
void BankEfficiencyGetOrthogonalFilter(
  REAL4Vector *filter)
{
  UINT4   i;
  UINT4   n      = filter->length;
  UINT4   nby2   = filter->length / 2;
  REAL4   temp;
  
  for (i = 1; i < nby2; i++)
  {
    temp              =  filter->data[i];
    filter->data[i]   = -filter->data[n-i];
    filter->data[n-i] = temp;
  }
}

/* *****************************************************************************
 * Create a vector f^{a/b} 
 * ************************************************************************** */
void BankEfficiencyCreateVectorFreqPower(
  REAL4Vector      *vector,
  InspiralTemplate  params,
  INT4              a,
  INT4              b)
{
  INT4  i;
  INT4  n = vector->length; /* Length of the vector to create */
      
  REAL8 power = (REAL8)a / (REAL8)b;       /* the frequency power */
  REAL8 f;                                 /* frequency */
  REAL8 df = params.tSampling/(REAL8)n/2.; /* sampling frequency */

  /* First value equal to zero */
  vector->data[0] = 0.;
  for( i = 1; i < n; i++)
  {
    f = i * df; 
    /* Probably not efficient but this function is 
     * just called once at the beginning */
    vector->data[i] = pow(f, power);
  }
}


/* *****************************************************************************
 * Create the filters for BCV correlation                   
 * ************************************************************************** */
void BankEfficiencyCreateBCVFilters(
  BankEfficiencyBCV *bankefficiencyBCV,
  UINT4              kMin,           /* position of the integration   */
  UINT4              kMax,
  REAL4              psi0,