TMarkov.cc

Go to the documentation of this file.

/* 
 *  \class TMarkov
 *
 *  \author: Patrice Verrecchia - CEA/Saclay
 */
 
#include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMarkov.h>
 
#include <iostream>
#include <cmath>
 
//ClassImp(TMarkov)
 
// Constructor...
TMarkov::TMarkov() {
  // 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();
}
 
// Destructor
TMarkov::~TMarkov() {}
 
void TMarkov::init() {
  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;
}
 
int TMarkov::computeChain(int *bing) {
  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;
}
 
void TMarkov::peakFinder(int *bing) {
  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;
  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;
}