CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TMarkov.cc
Go to the documentation of this file.
1 /*
2  * \class TMarkov
3  *
4  * $Date: 2012/02/09 10:08:10 $
5  * \author: Patrice Verrecchia - CEA/Saclay
6  */
7 
9 
10 #include <iostream>
11 #include "math.h"
12 
13 //ClassImp(TMarkov)
14 
15 
16 // Constructor...
18 {
19  // TMarkov
20  // ------ calcule les distributions invariantes de la chaine de TMarkov
21  // correspondantes au spectre original et retourne la dimension de u.
22  //
23 
24  fNPeakValues=3;
25  fNbinu=101;
26  init();
27 }
28 
29 // Destructor
31 {
32 }
33 
35 {
36  int i;
37  for(i=0;i<fNPeakValues;i++) peak[i]=0.;
38  for(i=0;i<fNbinu;i++) u[i]=0.;
39  for(i=0;i<=fNbinu;i++) binu[i]=0.;
40  return ;
41 }
42 
43 int TMarkov::computeChain(int *bing)
44 {
45  int i;int k;int nuprime;int offset=0;int m;int pass;
46  double sumUprime,sumU;
47  double jumpToNext,jumpToPrevious;
48  double chainToNext,chainToPrevious;
49  double aConst[101],uprime[101];
50 
51  pass=0;
52  for(m=3,i=1,nuprime=1;i<101;i++)
53  {
54  uprime[i]=0.;
55  for(k=1,jumpToNext=0.,jumpToPrevious=0.;k<=m;k++)
56  {
57  if(i+k < 101)
58  if(bing[i] > 0 || bing[i+k] > 0)
59  jumpToNext += exp( (double)(bing[i+k]-bing[i])
60  /sqrt((double)(bing[i+k]+bing[i])));
61  if(i-k > 0)
62  if(bing[i] > 0 || bing[i-k] > 0)
63  jumpToPrevious += exp( (double)(bing[i-k]-bing[i])
64  /sqrt((double)(bing[i-k]+bing[i])));
65  }
66  //printf(" jump %d to %d = %f\n",i,i+1,jumpToNext);
67  //printf(" jump %d to %d = %f\n",i,i-1,jumpToPrevious);
68  if(jumpToNext > 0. && jumpToPrevious > 0.)
69  {
70  aConst[i] = -log(jumpToNext+jumpToPrevious);
71  chainToNext = aConst[i]+log(jumpToNext);
72  chainToPrevious = aConst[i]+log(jumpToPrevious);
73  uprime[i]=chainToNext - chainToPrevious;
74  nuprime++; u[nuprime] = uprime[i];
75  if(pass == 0)
76  { offset=i-1; pass=1;}
77  }
78  }
79 
80  //for(i=1;i<101;i++)
81  //printf(" bin numero %d uprime = %f\n",i,uprime[i]);
82 
83  for(k=3,sumUprime=u[2],sumU=u[2];k<nuprime+1;k++)
84  {
85  sumU += u[k]; u[k] = sumU;
86  sumUprime += log(1.+exp(u[k]-u[k-1]));
87  }
88 
89  u[1] = -sumUprime;
90 
91  for(k=2;k<nuprime+1;k++)
92  u[k] += u[1];
93 
94  for(i=1;i<offset+1;i++)
95  binu[i]=0.;
96 
97  for(i=1;i<nuprime+1;i++)
98  {
99  binu[i+offset] = exp(u[i]);
100  //printf(" bin numero %d log(u) = %f\n",i+offset,u[i]);
101  //printf(" bin numero %d u = %f\n",i+offset,exp(u[i]));
102 
103  }
104 
105  return nuprime+offset;
106 }
107 
108 void TMarkov::peakFinder(int *bing)
109 {
110  int firstBin=0;int lastBin=0;
111  double barycentre=0.;
112  double sum=0.;
113  double maximum=0.;
114 
115  int nu= computeChain(&bing[0]);
116 
117  for(int i=1;i<nu+1;i++)
118  {
119  sum += binu[i];
120  barycentre += (double)i * binu[i];
121  if(binu[i] > maximum)
122  { maximum=binu[i]; imax=i; }
123  }
124 
125  maximum *= 0.75;
126  for(int i=1,pass=0;i<nu+1;i++) {
127  if(binu[i] > maximum) {
128  if(pass == 0) {
129  firstBin=i;
130  lastBin=i;
131  pass=1;
132  } else {
133  lastBin=i;
134  }
135  }
136  }
137 
138  peak[0] = (barycentre/sum);
139  peak[1]= (double)(lastBin-firstBin+1);
140  peak[2]= sum;
141 }
int i
Definition: DBlmapReader.cc:9
int imax
Definition: TMarkov.h:12
void init()
Definition: TMarkov.cc:34
virtual ~TMarkov()
Definition: TMarkov.cc:30
int computeChain(int *)
Definition: TMarkov.cc:43
TMarkov()
Definition: TMarkov.cc:17
int fNbinu
Definition: TMarkov.h:11
T sqrt(T t)
Definition: SSEVec.h:46
double binu[102]
Definition: TMarkov.h:14
unsigned int offset(bool)
void peakFinder(int *)
Definition: TMarkov.cc:108
int k[5][pyjets_maxn]
double u[101]
Definition: TMarkov.h:14
double peak[3]
Definition: TMarkov.h:13
int fNPeakValues
Definition: TMarkov.h:11