Logo ROOT   6.18/05
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
32This PDEFoam variant stores in every cell the discriminant
33
34 D = #events with given class / total number of events
35
36as well as the statistical error on the discriminant. It therefore
37acts as a discriminant estimator. It should be booked together
38with the PDEFoamDiscriminantDensity density estimator, which
39returns the discriminant density at a given phase space point
40during 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"
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
171TH2D* 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}
int Int_t
Definition: RtypesCore.h:41
unsigned int UInt_t
Definition: RtypesCore.h:42
long Long_t
Definition: RtypesCore.h:50
double Double_t
Definition: RtypesCore.h:55
float Float_t
Definition: RtypesCore.h:53
#define ClassImp(name)
Definition: Rtypes.h:365
#define gDirectory
Definition: TDirectory.h:218
char name[80]
Definition: TGX11.cxx:109
char * Form(const char *fmt,...)
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:464
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:928
TAxis * GetZaxis()
Definition: TH1.h:318
virtual Int_t GetNbinsY() const
Definition: TH1.h:293
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
virtual Int_t GetNbinsX() const
Definition: TH1.h:292
TAxis * GetYaxis()
Definition: TH1.h:317
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:8635
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4882
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:289
UInt_t GetClass() const
Definition: Event.h:87
std::vector< Float_t > & GetValues()
Definition: Event.h:95
This PDEFoam variant stores in every cell the discriminant.
virtual void FillFoamCells(const Event *ev, Float_t wt)
This function fills an event into the discriminant PDEFoam.
virtual void Finalize()
Calc discriminator and its error for every cell and save it to the cell.
virtual TH2D * Project2(Int_t, Int_t, ECellValue, PDEFoamKernelBase *, UInt_t)
Project foam variable idim1 and variable idim2 to histogram.
PDEFoamDiscriminant()
Default constructor for streamer, user should not use it.
This class is the abstract kernel interface for PDEFoam.
virtual Float_t Estimate(PDEFoam *, std::vector< Float_t > &, ECellValue)=0
Implementation of PDEFoam.
Definition: PDEFoam.h:77
MsgLogger & Log() const
Definition: PDEFoam.h:238
Basic string class.
Definition: TString.h:131
const char * Data() const
Definition: TString.h:364
TH1F * h1
Definition: legend1.C:5
MsgLogger & Endl(MsgLogger &ml)
Definition: MsgLogger.h:158
Double_t Log(Double_t x)
Definition: TMath.h:748
Double_t Sqrt(Double_t x)
Definition: TMath.h:679
REAL epsilon
Definition: triangle.c:617