EcalLaserAnalyzer/src/TMarkov.cc

0001 /*
0002  *  \class TMarkov
0003  *
0004  *  \author: Patrice Verrecchia - CEA/Saclay
0005  */
0006
0007 #include <CalibCalorimetry/EcalLaserAnalyzer/interface/TMarkov.h>
0008
0009 #include <iostream>
0010 #include <cmath>
0011
0012 //ClassImp(TMarkov)
0013
0014 // Constructor...
0015 TMarkov::TMarkov() {
0016   // TMarkov
0017   // ------ calcule les distributions invariantes de la chaine de TMarkov
0018   //   correspondantes au spectre original et retourne la dimension de u.
0019   //
0020
0021   fNPeakValues = 3;
0022   fNbinu = 101;
0023   init();
0024 }
0025
0026 // Destructor
0027 TMarkov::~TMarkov() {}
0028
0029 void TMarkov::init() {
0030   int i;
0031   for (i = 0; i < fNPeakValues; i++)
0032     peak[i] = 0.;
0033   for (i = 0; i < fNbinu; i++)
0034     u[i] = 0.;
0035   for (i = 0; i <= fNbinu; i++)
0036     binu[i] = 0.;
0037   return;
0038 }
0039
0040 int TMarkov::computeChain(int *bing) {
0041   int i;
0042   int k;
0043   int nuprime;
0044   int offset = 0;
0045   int m;
0046   int pass;
0047   double sumUprime, sumU;
0048   double jumpToNext, jumpToPrevious;
0049   double chainToNext, chainToPrevious;
0050   double aConst[101], uprime[101];
0051
0052   pass = 0;
0053   for (m = 3, i = 1, nuprime = 1; i < 101; i++) {
0054     uprime[i] = 0.;
0055     for (k = 1, jumpToNext = 0., jumpToPrevious = 0.; k <= m; k++) {
0056       if (i + k < 101)
0057         if (bing[i] > 0 || bing[i + k] > 0)
0058           jumpToNext += exp((double)(bing[i + k] - bing[i]) / sqrt((double)(bing[i + k] + bing[i])));
0059       if (i - k > 0)
0060         if (bing[i] > 0 || bing[i - k] > 0)
0061           jumpToPrevious += exp((double)(bing[i - k] - bing[i]) / sqrt((double)(bing[i - k] + bing[i])));
0062     }
0063     //printf(" jump %d to %d = %f\n",i,i+1,jumpToNext);
0064     //printf(" jump %d to %d = %f\n",i,i-1,jumpToPrevious);
0065     if (jumpToNext > 0. && jumpToPrevious > 0.) {
0066       aConst[i] = -log(jumpToNext + jumpToPrevious);
0067       chainToNext = aConst[i] + log(jumpToNext);
0068       chainToPrevious = aConst[i] + log(jumpToPrevious);
0069       uprime[i] = chainToNext - chainToPrevious;
0070       nuprime++;
0071       u[nuprime] = uprime[i];
0072       if (pass == 0) {
0073         offset = i - 1;
0074         pass = 1;
0075       }
0076     }
0077   }
0078
0079   //for(i=1;i<101;i++)
0080   //printf(" bin numero %d   uprime = %f\n",i,uprime[i]);
0081
0082   for (k = 3, sumUprime = u[2], sumU = u[2]; k < nuprime + 1; k++) {
0083     sumU += u[k];
0084     u[k] = sumU;
0085     sumUprime += log(1. + exp(u[k] - u[k - 1]));
0086   }
0087
0088   u[1] = -sumUprime;
0089
0090   for (k = 2; k < nuprime + 1; k++)
0091     u[k] += u[1];
0092
0093   for (i = 1; i < offset + 1; i++)
0094     binu[i] = 0.;
0095
0096   for (i = 1; i < nuprime + 1; i++) {
0097     binu[i + offset] = exp(u[i]);
0098     //printf(" bin numero %d   log(u) = %f\n",i+offset,u[i]);
0099     //printf(" bin numero %d   u = %f\n",i+offset,exp(u[i]));
0100   }
0101
0102   return nuprime + offset;
0103 }
0104
0105 void TMarkov::peakFinder(int *bing) {
0106   int firstBin = 0;
0107   int lastBin = 0;
0108   double barycentre = 0.;
0109   double sum = 0.;
0110   double maximum = 0.;
0111
0112   int nu = computeChain(&bing[0]);
0113
0114   for (int i = 1; i < nu + 1; i++) {
0115     sum += binu[i];
0116     barycentre += (double)i * binu[i];
0117     if (binu[i] > maximum) {
0118       maximum = binu[i];
0119       imax = i;
0120     }
0121   }
0122
0123   maximum *= 0.75;
0124   for (int i = 1, pass = 0; i < nu + 1; i++) {
0125     if (binu[i] > maximum) {
0126       if (pass == 0) {
0127         firstBin = i;
0128         lastBin = i;
0129         pass = 1;
0130       } else {
0131         lastBin = i;
0132       }
0133     }
0134   }
0135
0136   peak[0] = (barycentre / sum);
0137   peak[1] = (double)(lastBin - firstBin + 1);
0138   peak[2] = sum;
0139 }