meerkat is hosted by Hepforge, IPPP Durham
Meerkat  v1r3
Multidimensional kernel density estimation package
PolynomialDensity.cpp
Go to the documentation of this file.
1 #include "PolynomialDensity.hh"
2 #include "AbsPhaseSpace.hh"
3 #include "OneDimPhaseSpace.hh"
4 #include "UniformDensity.hh"
5 
6 #include "TMath.h"
7 #include "TVirtualFitter.h"
8 #include "TTree.h"
9 
10 #include "Logger.hh"
11 
12 static std::vector<std::vector<Double_t> > gData;
13 static Int_t gPower;
15 static std::vector<Double_t> gMiddle;
16 static Double_t gMaximum;
17 static std::vector<std::pair<Double_t, Double_t> > gIntegVector2D;
18 static Int_t gIter;
19 static const char* gName;
20 
21 Double_t gPoly1D(Double_t* pars, Double_t x, Int_t pow) {
22  // Number of parameters = pow
23  Int_t psum, pnum;
24  Double_t sum = 1.;
25  pnum = 0;
26  for (psum=1; psum<=pow; psum++) {
27  sum += pars[pnum]*TMath::Power(x - gMiddle[0], psum);
28  pnum ++;
29  }
30 
31  return sum;
32 
33 }
34 
35 Double_t gPoly1DInt(Double_t* pars, Int_t pow, AbsPhaseSpace* phsp) {
36  Double_t low = phsp->lowerLimit(0);
37  Double_t up = phsp->upperLimit(0);
38 
39  Int_t psum, pnum;
40  Double_t sum = up-low;
41 
42  pnum = 0;
43  for (psum=1; psum<=pow; psum++) {
44  sum += pars[pnum]*(TMath::Power(up - gMiddle[0], psum+1) - TMath::Power(low - gMiddle[0], psum+1))/(psum+1.);
45  pnum ++;
46  }
47 
48  return sum;
49 
50 }
51 
52 void gPoly1DLH(Int_t &, Double_t *, Double_t &f, Double_t *pars, Int_t) {
53  std::vector<std::vector<Double_t> >::iterator i;
54  Double_t lhsum = 0.;
55  Int_t error = 0;
56  gIter++;
57 
58  Double_t integ = gPoly1DInt(pars, gPower, gPhsp);
59  for (i=gData.begin(); i != gData.end(); i++) {
60  Double_t x = (*i)[0];
61  Double_t pdf = gPoly1D(pars, x, gPower);
62  if (pdf < 0) {
63  error = 1;
64  } else {
65  lhsum += TMath::Log( pdf/integ );
66  }
67  }
68  if (error) {
69  f = gMaximum;
70  } else {
71  f = -2.*lhsum;
72  if (f > gMaximum) gMaximum = f;
73  }
74 
75  if (Logger::timer(4)) {
76  Logger::print(0,"%20.20s INFO: Iteration %d\n", gName, gIter);
77  for (Int_t v=0; v<gPower; v++)
78  Logger::print(0,"%20.20s INFO: P%d = %f\n", gName, v+1, pars[v]);
79  }
80 }
81 
82 Double_t gPoly2D(Double_t* pars, Double_t x, Double_t y, Int_t pow) {
83  // Number of parameters: 2 + 3 + ... + (pow+1) = (pow+1)*pow/2-1
84 
85  Int_t psum, pnum;
86  Double_t sum = 1.;
87 
88  pnum = 0;
89  for (psum=1; psum<=pow; psum++) {
90  Int_t px;
91  for (px = 0; px <= psum; px++) {
92  Int_t py = psum-px;
93  sum += pars[pnum]*TMath::Power(x - gMiddle[0], px)*TMath::Power(y - gMiddle[1], py);
94  pnum ++;
95  }
96  }
97 
98  return sum;
99 }
100 
101 Double_t gPoly2DInt(Double_t* pars, Int_t pow, __attribute__((unused)) AbsPhaseSpace* phsp) {
102  Double_t sum = 0;
103 
104  std::vector<std::pair<Double_t, Double_t> >::iterator i;
105  for (i = gIntegVector2D.begin(); i != gIntegVector2D.end(); i++) {
106  Double_t x = (*i).first;
107  Double_t y = (*i).second;
108 
109  Double_t pdf = gPoly2D(pars, x, y, pow);
110  sum += pdf;
111  }
112 // sum /= (Double_t) gIntegVector2D.size();
113 
114  return sum;
115 }
116 
117 void gPoly2DLH(Int_t &, Double_t *, Double_t &f, Double_t *pars, Int_t) {
118  std::vector<std::vector<Double_t> >::iterator i;
119  Double_t lhsum = 0.;
120  Int_t error = 0;
121  gIter++;
122 
123  Double_t integ = gPoly2DInt(pars, gPower, gPhsp);
124  for (i=gData.begin(); i != gData.end(); i++) {
125  Double_t x = (*i)[0];
126  Double_t y = (*i)[1];
127  Double_t pdf = gPoly2D(pars, x, y, gPower);
128  if (pdf < 0) {
129  error = 1;
130  } else {
131  lhsum += TMath::Log( pdf/integ );
132  }
133  }
134  if (error) {
135  f = gMaximum;
136  } else {
137  f = -2.*lhsum;
138  if (f > gMaximum) gMaximum = f;
139  }
140 
141  if (Logger::timer(4)) {
142  Logger::print(0,"%20.20s INFO: Iteration %d\n", gName, gIter);
143  Int_t pnum = 0;
144  Int_t psum;
145  for (psum=1; psum<=(Int_t)gPower; psum++) {
146  Int_t px;
147  for (px = 0; px <= psum; px++) {
148  Int_t py = psum-px;
149  pnum++;
150  Logger::print(0,"%20.20s INFO: P%d_%d = %f\n", gName, px, py, pars[pnum-1]);
151  }
152  }
153  }
154 }
155 
157  OneDimPhaseSpace* thePhaseSpace,
158  UInt_t maxPower,
159  TTree* tree,
160  const char* var,
161  UInt_t maxEvents) : AbsDensity(pdfName) {
162  m_power = maxPower;
163  m_phaseSpace = thePhaseSpace;
165  gPhsp = thePhaseSpace;
166  gPower = maxPower;
167  gMaximum = -1.e50;
168  gIter = 0;
169  gName = pdfName;
170 
171  if (m_dim != 1) {
172  Logger::print(0,"%20.20s INFO: Dimensionality of the phase space is %d for constructor with 1 parameter\n", m_name, m_dim);
173  abort();
174  }
175 
176 
177  gMiddle.resize(1);
178  gMiddle[0] = (gPhsp->upperLimit(0) + gPhsp->lowerLimit(0))/2.;
179 
180  std::vector<TString> vars(1);
181  vars[0] = TString(var);
182  readTuple(tree, vars, maxEvents);
183 
184  // Create fitter and parameters
185  TVirtualFitter *minuit = TVirtualFitter::Fitter(0, gPower);
186  minuit->SetFCN(gPoly1DLH);
187 
188  double arglist[100];
189 
190  Double_t halfSize = (gPhsp->upperLimit(0) - gPhsp->lowerLimit(0))/2.;
191  for (Int_t v=1; v<=gPower; v++) {
192  char parname[32];
193  sprintf(parname, "P%d", v);
194  minuit->SetParameter(v-1, parname, 0., 0.01*TMath::Power(halfSize, 1./(Double_t)v), 0, 0);
195  }
196 
197  Logger::setTimer();
198 
199  // Fit polynomial
200  arglist[0] = 1;
201  minuit->ExecuteCommand("MIGRAD", arglist, 0);
202 
203  for (Int_t v=0; v<gPower; v++) {
204  m_par[v] = minuit->GetParameter(v);
205  Logger::print(0,"%20.20s INFO: Parameter %d = %f\n", m_name, v, m_par[v]);
206  }
207 
208 }
209 
211  AbsPhaseSpace* thePhaseSpace,
212  UInt_t maxPower,
213  TTree* tree,
214  const char* var1,
215  const char* var2,
216  UInt_t integEvents,
217  UInt_t maxEvents) : AbsDensity(pdfName) {
218  gPhsp = thePhaseSpace;
219  gPower = maxPower;
220  m_power = maxPower;
221  m_phaseSpace = thePhaseSpace;
223  gMaximum = -1.e50;
224  gIter = 0;
225  gName = pdfName;
226 
227  if (m_dim != 2) {
228  Logger::print(0,"%20.20s INFO: Dimensionality of the phase space is %d for constructor with 2 parameters\n", m_name, m_dim);
229  abort();
230  }
231 
232  gMiddle.resize(2);
233  gMiddle[0] = (gPhsp->upperLimit(0) + gPhsp->lowerLimit(0))/2.;
234  gMiddle[1] = (gPhsp->upperLimit(1) + gPhsp->lowerLimit(1))/2.;
235 
236  std::vector<TString> vars(2);
237  vars[0] = TString(var1);
238  vars[1] = TString(var2);
239  readTuple(tree, vars, maxEvents);
240 
241  UInt_t i;
242 
243  UniformDensity uniform(pdfName, m_phaseSpace);
244  std::vector<Double_t> genVector(2);
245  gIntegVector2D.clear();
246  for (i=0; i<integEvents; i++) {
247  uniform.generate(genVector);
248  Double_t x = genVector[0];
249  Double_t y = genVector[1];
250  std::pair<Double_t, Double_t> point(x, y);
251  gIntegVector2D.push_back(point);
252  }
253 
254  // Create fitter and parameters
255 
256  Int_t npars = (gPower+2)*(gPower+1)/2 - 1;
257 
258  TVirtualFitter *minuit = TVirtualFitter::Fitter(0, npars);
259  minuit->SetFCN(gPoly2DLH);
260 
261  double arglist[100];
262 
263  Int_t psum, pnum;
264  pnum = 0;
265  Double_t halfSizeX = (gPhsp->upperLimit(0) - gPhsp->lowerLimit(0))/2.;
266  Double_t halfSizeY = (gPhsp->upperLimit(1) - gPhsp->lowerLimit(1))/2.;
267 
268  for (psum=1; psum<=(Int_t)m_power; psum++) {
269  Int_t px;
270  for (px = 0; px <= psum; px++) {
271  Int_t py = psum-px;
272  pnum ++;
273  char parname[32];
274  sprintf(parname, "P%d_%d", px, py);
275  Double_t err = 0.01;
276  if (px > 0) err /= TMath::Power(halfSizeX, 1./(Double_t)px);
277  if (py > 0) err /= TMath::Power(halfSizeY, 1./(Double_t)py);
278  minuit->SetParameter(pnum-1, parname, 0., err, 0, 0);
279  }
280  }
281 
282  Logger::setTimer();
283 
284  // Fit polynomial
285  arglist[0] = 1;
286  minuit->ExecuteCommand("MIGRAD", arglist, 0);
287 
288  for (Int_t v=0; v<npars; v++) {
289  m_par[v] = minuit->GetParameter(v);
290  Logger::print(0,"%20.20s INFO: Parameter %d = %f\n", m_name, v, m_par[v]);
291  }
292 
293 }
294 
296 
297 }
298 
299 void PolynomialDensity::readTuple(TTree* tree, std::vector<TString> &vars, UInt_t maxEvents) {
300 
301  if (vars.size() != m_phaseSpace->dimensionality() ) {
302  Logger::print(2,"%20.20s ERROR: Number of TTree variables (%d) in tree \"%s\" does not correspond to phase space dimensionality (%d)\n",
303  m_name, (UInt_t)vars.size(), tree->GetName(), m_phaseSpace->dimensionality() );
304  abort();
305  }
306 
307  UInt_t nvars = vars.size();
308 
309  tree->ResetBranchAddresses();
310 
311  Long64_t nentries = tree->GetEntries();
312 
313  if (maxEvents > 0 && maxEvents < nentries) nentries = maxEvents;
314 
315  Long64_t i;
316 
317  std::vector<Float_t> varArray(nvars);
318 
319  UInt_t n;
320  for (n=0; n < nvars; n++) {
321  Logger::print(0,"%20.20s INFO: Will read branch \"%s\" from tree \"%s\"\n", m_name, vars[n].Data(), tree->GetName());
322  Int_t status = tree->SetBranchAddress(vars[n], &( varArray[n] ));
323  if (status < 0) {
324  Logger::print(2,"%20.20s ERROR: Error setting branch, status=%d\n", m_name, status);
325  abort();
326  }
327  }
328 
329  std::vector<Double_t> point(nvars);
330 
331  gData.clear();
332  gData.reserve(nentries);
333 
334  Logger::setTimer();
335 
336  UInt_t nout = 0;
337  for(i=0; i<nentries; i++) {
338  tree->GetEntry(i);
339 
340  for (n=0; n<nvars; n++) {
341  point[n] = varArray[n];
342  }
343 
344  if (!m_phaseSpace->withinLimits( point )) {
345  nout ++;
346  Logger::print(1,"%20.20s WARNING: Ntuple point (", m_name);
347  for (n=0; n<nvars; n++) {
348  Logger::print(1,"%f ", point[n]);
349  }
350  Logger::print(1,", %f%%) outside phase space\n", 100.*float(nout)/float(i));
351  } else {
352  for (UInt_t v = 0; v < nvars; v++) {
353  gData.push_back(point);
354  }
355  }
356 
357  if (i % 100 == 0 && Logger::timer(2)) {
358  Logger::print(0,"%20.20s INFO: Read %lld/%lld events (%f%%)\n", m_name, i, nentries, 100.*float(i)/float(nentries));
359  }
360  }
361 
362  Logger::print(0,"%20.20s INFO: %lld events read in from \"%s\"\n", m_name, nentries-nout, tree->GetName() );
363 
364 }
365 
366 Double_t PolynomialDensity::density(std::vector<Double_t> &x) {
367  gPhsp = m_phaseSpace;
368  gPower = m_power;
369  gMaximum = -1.e50;
370  if (m_dim == 1) {
371  gMiddle.resize(1);
372  gMiddle[0] = (gPhsp->upperLimit(0) + gPhsp->lowerLimit(0))/2.;
373  return gPoly1D(m_par, x[0], m_power);
374  } else if (m_dim == 2) {
375  gMiddle.resize(2);
376  gMiddle[0] = (gPhsp->upperLimit(0) + gPhsp->lowerLimit(0))/2.;
377  gMiddle[1] = (gPhsp->upperLimit(1) + gPhsp->lowerLimit(1))/2.;
378  return gPoly2D(m_par, x[0], x[1], m_power);
379  }
380  return 0.;
381 }
static Double_t gMaximum
Double_t gPoly2DInt(Double_t *pars, Int_t pow, __attribute__((unused)) AbsPhaseSpace *phsp)
void readTuple(TTree *tree, std::vector< TString > &vars, UInt_t maxEvents=0)
Read Ntuple contaning data points to be used for the fit.
static Int_t gIter
virtual Double_t upperLimit(UInt_t var)=0
Return upper allowed limit of the variable.
Abstract class which defines probability density interface.
Definition: AbsDensity.hh:16
virtual ~PolynomialDensity()
Destructor.
void gPoly2DLH(Int_t &, Double_t *, Double_t &f, Double_t *pars, Int_t)
void gPoly1DLH(Int_t &, Double_t *, Double_t &f, Double_t *pars, Int_t)
void print(int level, const char *format,...)
Definition: Logger.cpp:27
char m_name[256]
PDF name.
Definition: AbsDensity.hh:118
UInt_t m_dim
Cached dimensionality of the phase space.
Class that describes the uniform density over any phase space.
static std::vector< Double_t > gMiddle
Abstract class which defines phase space interface.
virtual UInt_t dimensionality()=0
Get dimensionality of phase space.
void setTimer(void)
Reset time counter.
Definition: Logger.cpp:9
virtual Double_t lowerLimit(UInt_t var)=0
Return lower allowed limit of the variable.
static std::vector< std::vector< Double_t > > gData
static const char * gName
int timer(int secs)
Return 1 if more than a certain number of seconds have passed since the last call of this function or...
Definition: Logger.cpp:13
Double_t generate(std::vector< Double_t > &x)
Generate a single point within the phase space according to the PDF using accept-reject method...
Definition: AbsDensity.cpp:99
AbsPhaseSpace * m_phaseSpace
Reference to phase space.
Double_t gPoly2D(Double_t *pars, Double_t x, Double_t y, Int_t pow)
static AbsPhaseSpace * gPhsp
static Int_t gPower
virtual Bool_t withinLimits(std::vector< Double_t > &x)=0
Check if the point is within phase space limits.
Double_t gPoly1DInt(Double_t *pars, Int_t pow, AbsPhaseSpace *phsp)
PolynomialDensity(const char *pdfName, OneDimPhaseSpace *thePhaseSpace, UInt_t maxPower, TTree *tree, const char *var, UInt_t maxEvents=0)
static std::vector< std::pair< Double_t, Double_t > > gIntegVector2D
Double_t m_par[100]
Vector of fitted parameters.
UInt_t m_power
Power of the polynomial.
Double_t density(std::vector< Double_t > &x)
Calculate PDF value at the given point.
Double_t gPoly1D(Double_t *pars, Double_t x, Int_t pow)