#include <TMarkov.h>

Public Member Functions
int	getBinMax () const
double	getPeakValue (int i) const
void	peakFinder (int *)
	TMarkov ()
virtual	~TMarkov ()
Private Member Functions
int	computeChain (int *)
void	init ()
Private Attributes
double	binu [102]
int	fNbinu
int	fNPeakValues
int	imax
double	peak [3]
double	u [101]

Detailed Description

Definition at line 6 of file TMarkov.h.

Constructor & Destructor Documentation

TMarkov::TMarkov ( )

Definition at line 17 of file TMarkov.cc.

References fNbinu, fNPeakValues, and init().

{ 
  // TMarkov
  // ------ calcule les distributions invariantes de la chaine de TMarkov
  //   correspondantes au spectre original et retourne la dimension de u.
  //
 
  fNPeakValues=3;
  fNbinu=101;
  init();
}

TMarkov::~TMarkov ( ) [virtual]

Definition at line 30 of file TMarkov.cc.

{ 
}

Member Function Documentation

int TMarkov::computeChain ( int * bing ) [private]

Definition at line 43 of file TMarkov.cc.

References binu, funct::exp(), i, gen::k, funct::log(), m, evf::evtn::offset(), mathSSE::sqrt(), and u.

Referenced by peakFinder().

{
  int i;int k;int nuprime;int offset=0;int m;int pass;
  double sumUprime,sumU;
  double jumpToNext,jumpToPrevious;
  double chainToNext,chainToPrevious;
  double aConst[101],uprime[101];  

  pass=0;
  for(m=3,i=1,nuprime=1;i<101;i++)
  {
       uprime[i]=0.;
       for(k=1,jumpToNext=0.,jumpToPrevious=0.;k<=m;k++)
       {
          if(i+k < 101)
            if(bing[i] > 0 || bing[i+k] > 0)        
              jumpToNext += exp( (double)(bing[i+k]-bing[i])
                                /sqrt((double)(bing[i+k]+bing[i])));
          if(i-k > 0)
            if(bing[i] > 0 || bing[i-k] > 0)
              jumpToPrevious += exp( (double)(bing[i-k]-bing[i])
                                    /sqrt((double)(bing[i-k]+bing[i])));
        }
       //printf(" jump %d to %d = %f\n",i,i+1,jumpToNext);
       //printf(" jump %d to %d = %f\n",i,i-1,jumpToPrevious);
        if(jumpToNext > 0. && jumpToPrevious > 0.)
        {
          aConst[i] = -log(jumpToNext+jumpToPrevious);
          chainToNext = aConst[i]+log(jumpToNext);
          chainToPrevious = aConst[i]+log(jumpToPrevious);
          uprime[i]=chainToNext - chainToPrevious;
          nuprime++; u[nuprime] = uprime[i];
          if(pass == 0)
          {  offset=i-1; pass=1;}
        }                                                    
    }
    
  //for(i=1;i<101;i++)
  //printf(" bin numero %d   uprime = %f\n",i,uprime[i]);
    
    for(k=3,sumUprime=u[2],sumU=u[2];k<nuprime+1;k++)
    {
       sumU += u[k];  u[k] = sumU;
       sumUprime += log(1.+exp(u[k]-u[k-1]));
    }      
        
    u[1] = -sumUprime;

    for(k=2;k<nuprime+1;k++)
        u[k] += u[1];

    for(i=1;i<offset+1;i++)
       binu[i]=0.;

    for(i=1;i<nuprime+1;i++)
    {
       binu[i+offset] = exp(u[i]);
       //printf(" bin numero %d   log(u) = %f\n",i+offset,u[i]); 
       //printf(" bin numero %d   u = %f\n",i+offset,exp(u[i])); 
       
    }
  
    return nuprime+offset;
}

int TMarkov::getBinMax ( ) const [inline]

Definition at line 28 of file TMarkov.h.

References imax.

{ return imax; }

double TMarkov::getPeakValue ( int i ) const [inline]

Definition at line 27 of file TMarkov.h.

References i, and peak.

{ return peak[i]; }

void TMarkov::init ( ) [private]

Definition at line 34 of file TMarkov.cc.

References binu, fNbinu, fNPeakValues, i, peak, and u.

Referenced by TMarkov().

{
  int i;
  for(i=0;i<fNPeakValues;i++) peak[i]=0.;
  for(i=0;i<fNbinu;i++) u[i]=0.;
  for(i=0;i<=fNbinu;i++) binu[i]=0.;
  return ;
}

void TMarkov::peakFinder ( int * bing )

Definition at line 108 of file TMarkov.cc.

References binu, computeChain(), i, imax, and peak.

{
    int pass;int firstBin=0;int lastBin=0;
    double barycentre=0.;
    double sum=0.;
    double maximum=0.;

    int nu= computeChain(&bing[0]);

    for(int i=1;i<nu+1;i++)
    {
       sum += binu[i];
       barycentre += (double)i * binu[i];
       if(binu[i] > maximum)
       { maximum=binu[i]; imax=i; }
    }
    
    maximum *= 0.75;
    pass=0;
    for(int i=1,pass=0;i<nu+1;i++) {
       if(binu[i] > maximum) {
         if(pass == 0) {
                firstBin=i;
                lastBin=i;
                pass=1;
         } else {
                 lastBin=i;
         }
       }
    }

    peak[0] = (barycentre/sum);
    peak[1]= (double)(lastBin-firstBin+1);
    peak[2]= sum;
}