Logo ROOT   6.10/09
Reference Guide
PDEFoamDiscriminant.cxx
Go to the documentation of this file.
1 // @(#)root/tmva $Id$
2 // Author: Tancredi Carli, Dominik Dannheim, Alexander Voigt
3 
4 /**********************************************************************************
5  * Project: TMVA - a Root-integrated toolkit for multivariate data analysis *
6  * Package: TMVA *
7  * Classes: PDEFoamDiscriminant *
8  * Web : http://tmva.sourceforge.net *
9  * *
10  * Description: *
11  * Implementation. *
12  * *
13  * Authors (alphabetical): *
14  * Tancredi Carli - CERN, Switzerland *
15  * Dominik Dannheim - CERN, Switzerland *
16  * S. Jadach - Institute of Nuclear Physics, Cracow, Poland *
17  * Alexander Voigt - TU Dresden, Germany *
18  * Peter Speckmayer - CERN, Switzerland *
19  * *
20  * Copyright (c) 2008, 2010: *
21  * CERN, Switzerland *
22  * MPI-K Heidelberg, Germany *
23  * *
24  * Redistribution and use in source and binary forms, with or without *
25  * modification, are permitted according to the terms listed in LICENSE *
26  * (http://tmva.sourceforge.net/LICENSE) *
27  **********************************************************************************/
28 
29 /*! \class TMVA::PDEFoamDiscriminant
30 \ingroup TMVA
31 
32 This PDEFoam variant stores in every cell the discriminant
33 
34  D = #events with given class / total number of events
35 
36 as well as the statistical error on the discriminant. It therefore
37 acts as a discriminant estimator. It should be booked together
38 with the PDEFoamDiscriminantDensity density estimator, which
39 returns the discriminant density at a given phase space point
40 during the foam build-up.
41 */
42 
44 
45 
46 #include "TMVA/Event.h"
47 #include "TMVA/MsgLogger.h"
48 #include "TMVA/PDEFoam.h"
49 #include "TMVA/PDEFoamCell.h"
50 #include "TMVA/PDEFoamKernelBase.h"
51 #include "TMVA/Types.h"
52 
53 #include "TDirectory.h"
54 #include "TMath.h"
55 #include "TH2D.h"
56 
57 #include <climits>
58 
60 
61 ////////////////////////////////////////////////////////////////////////////////
62 /// Default constructor for streamer, user should not use it.
63 
65 : PDEFoam()
66  , fClass(0)
67 {
68 }
69 
70 ////////////////////////////////////////////////////////////////////////////////
71 
73  : PDEFoam(name)
74  , fClass(cls)
75 {}
76 
77 ////////////////////////////////////////////////////////////////////////////////
78 /// Copy Constructor NOT IMPLEMENTED (NEVER USED)
79 
81  : PDEFoam(from)
82  , fClass(from.fClass)
83 {
84  Log() << kFATAL << "COPY CONSTRUCTOR NOT IMPLEMENTED" << Endl;
85 }
86 
87 ////////////////////////////////////////////////////////////////////////////////
88 /// This function fills an event into the discriminant PDEFoam. The
89 /// event weight 'wt' is filled into cell element 0 if the event is
90 /// of class fClass, and filled into cell element 1 otherwise.
91 
93 {
94  // find corresponding foam cell
95  std::vector<Float_t> values = ev->GetValues();
96  std::vector<Float_t> tvalues = VarTransform(values);
97  PDEFoamCell *cell = FindCell(tvalues);
98 
99  // 0. Element: Number of signal events (even class == fClass)
100  // 1. Element: Number of background events times normalization
101  if (ev->GetClass() == fClass)
102  SetCellElement(cell, 0, GetCellElement(cell, 0) + wt);
103  else
104  SetCellElement(cell, 1, GetCellElement(cell, 1) + wt);
105 }
106 
107 ////////////////////////////////////////////////////////////////////////////////
108 /// Calc discriminator and its error for every cell and save it to
109 /// the cell.
110 
112 {
113  // loop over cells
114  for (Long_t iCell = 0; iCell <= fLastCe; iCell++) {
115  if (!(fCells[iCell]->GetStat()))
116  continue;
117 
118  Double_t n_sig = GetCellElement(fCells[iCell], 0); // get number of signal events
119  Double_t n_bg = GetCellElement(fCells[iCell], 1); // get number of bg events
120 
121  if (n_sig < 0.) {
122  Log() << kWARNING << "Negative number of signal events in cell " << iCell
123  << ": " << n_sig << ". Set to 0." << Endl;
124  n_sig = 0.;
125  }
126  if (n_bg < 0.) {
127  Log() << kWARNING << "Negative number of background events in cell " << iCell
128  << ": " << n_bg << ". Set to 0." << Endl;
129  n_bg = 0.;
130  }
131 
132  // calculate discriminant
133  if (n_sig + n_bg > 0) {
134  // discriminant
135  SetCellElement(fCells[iCell], 0, n_sig / (n_sig + n_bg));
136  // discriminant error
137  SetCellElement(fCells[iCell], 1, TMath::Sqrt(Sqr(n_sig / Sqr(n_sig + n_bg))*n_sig +
138  Sqr(n_bg / Sqr(n_sig + n_bg))*n_bg));
139 
140  } else {
141  SetCellElement(fCells[iCell], 0, 0.5); // set discriminator
142  SetCellElement(fCells[iCell], 1, 1.); // set discriminator error
143  }
144  }
145 }
146 
147 ////////////////////////////////////////////////////////////////////////////////
148 /// Project foam variable idim1 and variable idim2 to histogram.
149 /// The projection algorithm is modified such that the z axis range
150 /// of the returned histogram is [0, 1], as necessary for the
151 /// interpretation as a discriminator. This is done by weighting
152 /// the cell values (in case of cell_value = kValue) by the cell
153 /// volume in all dimensions, excluding 'idim1' and 'idim2'.
154 ///
155 /// Parameters:
156 ///
157 /// - idim1, idim2 - dimensions to project to
158 ///
159 /// - cell_value - the cell value to draw
160 ///
161 /// - kernel - a PDEFoam kernel (optional). If NULL is given, the
162 /// kernel is ignored and the pure cell values are
163 /// plotted.
164 ///
165 /// - nbin - number of bins in x and y direction of result histogram
166 /// (optional, default is 50).
167 ///
168 /// Returns:
169 /// a 2-dimensional histogram
170 
171 TH2D* TMVA::PDEFoamDiscriminant::Project2(Int_t idim1, Int_t idim2, ECellValue cell_value, PDEFoamKernelBase *kernel, UInt_t nbin)
172 {
173  // avoid plotting of wrong dimensions
174  if ((idim1 >= GetTotDim()) || (idim1 < 0) ||
175  (idim2 >= GetTotDim()) || (idim2 < 0) ||
176  (idim1 == idim2))
177  Log() << kFATAL << "<Project2>: wrong dimensions given: "
178  << idim1 << ", " << idim2 << Endl;
179 
180  // root can not handle too many bins in one histogram --> catch this
181  // Furthermore, to have more than 1000 bins in the histogram doesn't make
182  // sense.
183  if (nbin > 1000) {
184  Log() << kWARNING << "Warning: number of bins too big: " << nbin
185  << " Using 1000 bins for each dimension instead." << Endl;
186  nbin = 1000;
187  } else if (nbin < 1) {
188  Log() << kWARNING << "Wrong bin number: " << nbin
189  << "; set nbin=50" << Endl;
190  nbin = 50;
191  }
192 
193  // create result histogram
194  TString hname(Form("h_%d_vs_%d", idim1, idim2));
195 
196  // if histogram with this name already exists, delete it
197  TH2D* h1 = (TH2D*)gDirectory->Get(hname.Data());
198  if (h1) delete h1;
199  h1 = new TH2D(hname.Data(), Form("var%d vs var%d", idim1, idim2), nbin, fXmin[idim1], fXmax[idim1], nbin, fXmin[idim2], fXmax[idim2]);
200 
201  if (!h1) Log() << kFATAL << "ERROR: Can not create histo" << hname << Endl;
202  if (cell_value == kValue)
205 
206  // ============== start projection algorithm ================
207  // loop over all histogram bins (2-dim)
208  for (Int_t xbin = 1; xbin <= h1->GetNbinsX(); ++xbin) {
209  for (Int_t ybin = 1; ybin <= h1->GetNbinsY(); ++ybin) {
210  // calculate the phase space point, which corresponds to this
211  // bin combination
212  std::map<Int_t, Float_t> txvec;
213  txvec[idim1] = VarTransform(idim1, h1->GetXaxis()->GetBinCenter(xbin));
214  txvec[idim2] = VarTransform(idim2, h1->GetYaxis()->GetBinCenter(ybin));
215 
216  // find the cells, which corresponds to this phase space
217  // point
218  std::vector<TMVA::PDEFoamCell*> cells = FindCells(txvec);
219 
220  // loop over cells and fill the histogram with the cell
221  // values
222  Float_t sum_cv = 0; // sum of the cell values
223  for (std::vector<TMVA::PDEFoamCell*>::const_iterator it = cells.begin();
224  it != cells.end(); ++it) {
225  // get cell position and size
226  PDEFoamVect cellPosi(GetTotDim()), cellSize(GetTotDim());
227  (*it)->GetHcub(cellPosi, cellSize);
228  // Create complete event vector from txvec. The missing
229  // coordinates of txvec are set to the cell center.
230  std::vector<Float_t> tvec;
231  for (Int_t i = 0; i < GetTotDim(); ++i) {
232  if (i != idim1 && i != idim2)
233  tvec.push_back(cellPosi[i] + 0.5 * cellSize[i]);
234  else
235  tvec.push_back(txvec[i]);
236  }
237  // get the cell value using the kernel
238  Float_t cv = 0;
239  if (kernel != NULL) {
240  cv = kernel->Estimate(this, tvec, cell_value);
241  } else {
242  cv = GetCellValue(FindCell(tvec), cell_value);
243  }
244  if (cell_value == kValue) {
245  // calculate cell volume in other dimensions (not
246  // including idim1 and idim2)
247  Float_t area_cell = 1.;
248  for (Int_t d1 = 0; d1 < GetTotDim(); ++d1) {
249  if ((d1 != idim1) && (d1 != idim2))
250  area_cell *= cellSize[d1];
251  }
252  // calc discriminator * (cell area times foam area)
253  // foam is normalized -> length of foam = 1.0
254  cv *= area_cell;
255  }
256  sum_cv += cv;
257  }
258 
259  // fill the bin content
260  h1->SetBinContent(xbin, ybin, sum_cv + h1->GetBinContent(xbin, ybin));
261  }
262  }
263 
264  return h1;
265 }
MsgLogger & Endl(MsgLogger &ml)
Definition: MsgLogger.h:158
This class is the abstract kernel interface for PDEFoam.
PDEFoamDiscriminant()
Default constructor for streamer, user should not use it.
float Float_t
Definition: RtypesCore.h:53
Float_t VarTransform(Int_t idim, Float_t x) const
Definition: PDEFoam.h:279
std::vector< TMVA::PDEFoamCell * > FindCells(const std::vector< Float_t > &) const
Find all cells, that contain txvec.
Definition: PDEFoam.cxx:1166
Double_t * fXmax
Definition: PDEFoam.h:105
Basic string class.
Definition: TString.h:129
int Int_t
Definition: RtypesCore.h:41
virtual TH2D * Project2(Int_t, Int_t, ECellValue, PDEFoamKernelBase *, UInt_t)
Project foam variable idim1 and variable idim2 to histogram.
Double_t GetCellElement(const PDEFoamCell *cell, UInt_t i) const
Returns cell element i of cell &#39;cell&#39;.
Definition: PDEFoam.cxx:1417
void SetCellElement(PDEFoamCell *cell, UInt_t i, Double_t value)
Set cell element i of cell to value.
Definition: PDEFoam.cxx:1433
#define NULL
Definition: RtypesCore.h:88
Int_t GetTotDim() const
Definition: PDEFoam.h:195
virtual void SetRangeUser(Double_t ufirst, Double_t ulast)
Set the viewing range for the axis from ufirst to ulast (in user coordinates).
Definition: TAxis.cxx:927
virtual Float_t Estimate(PDEFoam *, std::vector< Float_t > &, ECellValue)=0
PDEFoamCell * FindCell(const std::vector< Float_t > &) const
Find cell that contains &#39;xvec&#39; (in foam coordinates [0,1]).
Definition: PDEFoam.cxx:1081
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:464
UInt_t GetClass() const
Definition: Event.h:81
TH1F * h1
Definition: legend1.C:5
Implementation of PDEFoam.
Definition: PDEFoam.h:77
unsigned int UInt_t
Definition: RtypesCore.h:42
char * Form(const char *fmt,...)
TAxis * GetYaxis()
Definition: TH1.h:301
MsgLogger & Log() const
Definition: PDEFoam.h:238
Double_t * fXmin
Definition: PDEFoam.h:104
REAL epsilon
Definition: triangle.c:617
long Long_t
Definition: RtypesCore.h:50
#define ClassImp(name)
Definition: Rtypes.h:336
double Double_t
Definition: RtypesCore.h:55
PDEFoamCell ** fCells
Definition: PDEFoam.h:94
THist< 2, double, THistStatContent, THistStatUncertainty > TH2D
Definition: THist.hxx:316
virtual Float_t GetCellValue(const std::vector< Float_t > &xvec, ECellValue cv, PDEFoamKernelBase *)
This function finds the cell, which corresponds to the given untransformed event vector &#39;xvec&#39; and re...
Definition: PDEFoam.cxx:1017
virtual void Finalize()
Calc discriminator and its error for every cell and save it to the cell.
TAxis * GetZaxis()
Definition: TH1.h:302
virtual void FillFoamCells(const Event *ev, Float_t wt)
This function fills an event into the discriminant PDEFoam.
Abstract ClassifierFactory template that handles arbitrary types.
std::vector< Float_t > & GetValues()
Definition: Event.h:89
T Sqr(T x) const
Definition: PDEFoam.h:159
Int_t fLastCe
Definition: PDEFoam.h:93
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH2.h:84
#define gDirectory
Definition: TDirectory.h:211
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content.
Definition: TH2.cxx:2471
virtual Int_t GetNbinsX() const
Definition: TH1.h:277
Double_t Sqrt(Double_t x)
Definition: TMath.h:591
This PDEFoam variant stores in every cell the discriminant.
TAxis * GetXaxis()
Definition: TH1.h:300
virtual Int_t GetNbinsY() const
Definition: TH1.h:278
tomato 2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:290
const char * Data() const
Definition: TString.h:347