Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 * *
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 * (see tmva/doc/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
59
60////////////////////////////////////////////////////////////////////////////////
61/// Default constructor for streamer, user should not use it.
62
68
69////////////////////////////////////////////////////////////////////////////////
70
75
76////////////////////////////////////////////////////////////////////////////////
77/// Copy Constructor NOT IMPLEMENTED (NEVER USED)
78
80 : PDEFoam(from)
81 , fClass(from.fClass)
82{
83 Log() << kFATAL << "COPY CONSTRUCTOR NOT IMPLEMENTED" << Endl;
84}
85
86////////////////////////////////////////////////////////////////////////////////
87/// This function fills an event into the discriminant PDEFoam. The
88/// event weight 'wt' is filled into cell element 0 if the event is
89/// of class fClass, and filled into cell element 1 otherwise.
90
92{
93 // find corresponding foam cell
94 std::vector<Float_t> values = ev->GetValues();
95 std::vector<Float_t> tvalues = VarTransform(values);
96 PDEFoamCell *cell = FindCell(tvalues);
97
98 // 0. Element: Number of signal events (even class == fClass)
99 // 1. Element: Number of background events times normalization
100 if (ev->GetClass() == fClass)
101 SetCellElement(cell, 0, GetCellElement(cell, 0) + wt);
102 else
103 SetCellElement(cell, 1, GetCellElement(cell, 1) + wt);
104}
105
106////////////////////////////////////////////////////////////////////////////////
107/// Calc discriminator and its error for every cell and save it to
108/// the cell.
109
111{
112 // loop over cells
113 for (Long_t iCell = 0; iCell <= fLastCe; iCell++) {
114 if (!(fCells[iCell]->GetStat()))
115 continue;
116
117 Double_t n_sig = GetCellElement(fCells[iCell], 0); // get number of signal events
118 Double_t n_bg = GetCellElement(fCells[iCell], 1); // get number of bg events
119
120 if (n_sig < 0.) {
121 Log() << kWARNING << "Negative number of signal events in cell " << iCell
122 << ": " << n_sig << ". Set to 0." << Endl;
123 n_sig = 0.;
124 }
125 if (n_bg < 0.) {
126 Log() << kWARNING << "Negative number of background events in cell " << iCell
127 << ": " << n_bg << ". Set to 0." << Endl;
128 n_bg = 0.;
129 }
130
131 // calculate discriminant
132 if (n_sig + n_bg > 0) {
133 // discriminant
134 SetCellElement(fCells[iCell], 0, n_sig / (n_sig + n_bg));
135 // discriminant error
136 SetCellElement(fCells[iCell], 1, TMath::Sqrt(Sqr(n_sig / Sqr(n_sig + n_bg))*n_sig +
137 Sqr(n_bg / Sqr(n_sig + n_bg))*n_bg));
138
139 } else {
140 SetCellElement(fCells[iCell], 0, 0.5); // set discriminator
141 SetCellElement(fCells[iCell], 1, 1.); // set discriminator error
142 }
143 }
144}
145
146////////////////////////////////////////////////////////////////////////////////
147/// Project foam variable idim1 and variable idim2 to histogram.
148/// The projection algorithm is modified such that the z axis range
149/// of the returned histogram is [0, 1], as necessary for the
150/// interpretation as a discriminator. This is done by weighting
151/// the cell values (in case of cell_value = kValue) by the cell
152/// volume in all dimensions, excluding 'idim1' and 'idim2'.
153///
154/// Parameters:
155///
156/// - idim1, idim2 - dimensions to project to
157///
158/// - cell_value - the cell value to draw
159///
160/// - kernel - a PDEFoam kernel (optional). If NULL is given, the
161/// kernel is ignored and the pure cell values are
162/// plotted.
163///
164/// - nbin - number of bins in x and y direction of result histogram
165/// (optional, default is 50).
166///
167/// Returns:
168/// a 2-dimensional histogram
169
171{
172 // avoid plotting of wrong dimensions
173 if ((idim1 >= GetTotDim()) || (idim1 < 0) ||
174 (idim2 >= GetTotDim()) || (idim2 < 0) ||
175 (idim1 == idim2))
176 Log() << kFATAL << "<Project2>: wrong dimensions given: "
177 << idim1 << ", " << idim2 << Endl;
178
179 // root can not handle too many bins in one histogram --> catch this
180 // Furthermore, to have more than 1000 bins in the histogram doesn't make
181 // sense.
182 if (nbin > 1000) {
183 Log() << kWARNING << "Warning: number of bins too big: " << nbin
184 << " Using 1000 bins for each dimension instead." << Endl;
185 nbin = 1000;
186 } else if (nbin < 1) {
187 Log() << kWARNING << "Wrong bin number: " << nbin
188 << "; set nbin=50" << Endl;
189 nbin = 50;
190 }
191
192 // create result histogram
193 TString hname = TString::Format("h_%d_vs_%d", idim1, idim2);
194
195 // if histogram with this name already exists, delete it
196 TH2D* h1 = (TH2D*)gDirectory->Get(hname.Data());
197 if (h1) delete h1;
198 h1 = new TH2D(hname.Data(), TString::Format("var%d vs var%d", idim1, idim2).Data(), nbin, fXmin[idim1], fXmax[idim1], nbin, fXmin[idim2], fXmax[idim2]);
199
200 if (!h1) Log() << kFATAL << "ERROR: Can not create histo" << hname << Endl;
201 if (cell_value == kValue)
202 h1->GetZaxis()->SetRangeUser(-std::numeric_limits<float>::epsilon(),
203 1. + std::numeric_limits<float>::epsilon());
204
205 // ============== start projection algorithm ================
206 // loop over all histogram bins (2-dim)
207 for (Int_t xbin = 1; xbin <= h1->GetNbinsX(); ++xbin) {
208 for (Int_t ybin = 1; ybin <= h1->GetNbinsY(); ++ybin) {
209 // calculate the phase space point, which corresponds to this
210 // bin combination
211 std::map<Int_t, Float_t> txvec;
212 txvec[idim1] = VarTransform(idim1, h1->GetXaxis()->GetBinCenter(xbin));
213 txvec[idim2] = VarTransform(idim2, h1->GetYaxis()->GetBinCenter(ybin));
214
215 // find the cells, which corresponds to this phase space
216 // point
217 std::vector<TMVA::PDEFoamCell*> cells = FindCells(txvec);
218
219 // loop over cells and fill the histogram with the cell
220 // values
221 Float_t sum_cv = 0; // sum of the cell values
222 for (std::vector<TMVA::PDEFoamCell*>::const_iterator it = cells.begin();
223 it != cells.end(); ++it) {
224 // get cell position and size
225 PDEFoamVect cellPosi(GetTotDim()), cellSize(GetTotDim());
226 (*it)->GetHcub(cellPosi, cellSize);
227 // Create complete event vector from txvec. The missing
228 // coordinates of txvec are set to the cell center.
229 std::vector<Float_t> tvec;
230 for (Int_t i = 0; i < GetTotDim(); ++i) {
231 if (i != idim1 && i != idim2)
232 tvec.push_back(cellPosi[i] + 0.5 * cellSize[i]);
233 else
234 tvec.push_back(txvec[i]);
235 }
236 // get the cell value using the kernel
237 Float_t cv = 0;
238 if (kernel != NULL) {
239 cv = kernel->Estimate(this, tvec, cell_value);
240 } else {
241 cv = GetCellValue(FindCell(tvec), cell_value);
242 }
243 if (cell_value == kValue) {
244 // calculate cell volume in other dimensions (not
245 // including idim1 and idim2)
246 Float_t area_cell = 1.;
247 for (Int_t d1 = 0; d1 < GetTotDim(); ++d1) {
248 if ((d1 != idim1) && (d1 != idim2))
249 area_cell *= cellSize[d1];
250 }
251 // calc discriminator * (cell area times foam area)
252 // foam is normalized -> length of foam = 1.0
253 cv *= area_cell;
254 }
255 sum_cv += cv;
256 }
257
258 // fill the bin content
260 }
261 }
262
263 return h1;
264}
Cppyy::TCppType_t fClass
long Long_t
Signed long integer 4 bytes (long). Size depends on architecture.
Definition RtypesCore.h:68
float Float_t
Float 4 bytes (float)
Definition RtypesCore.h:71
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
#define gDirectory
Definition TDirectory.h:385
char name[80]
Definition TGX11.cxx:110
const_iterator begin() const
const_iterator end() const
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition TAxis.cxx:481
virtual void SetRangeUser(Double_t ufirst, Double_t ulast)
Set the viewing range for the axis from ufirst to ulast (in user coordinates, that is,...
Definition TAxis.cxx:1074
TAxis * GetZaxis()
Definition TH1.h:574
virtual Int_t GetNbinsY() const
Definition TH1.h:543
TAxis * GetXaxis()
Definition TH1.h:572
virtual Int_t GetNbinsX() const
Definition TH1.h:542
TAxis * GetYaxis()
Definition TH1.h:573
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:9260
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:5076
2-D histogram with a double per channel (see TH1 documentation)
Definition TH2.h:356
This PDEFoam variant stores in every cell the discriminant.
void Finalize() override
Calc discriminator and its error for every cell and save it to the cell.
void FillFoamCells(const Event *ev, Float_t wt) override
This function fills an event into the discriminant PDEFoam.
TH2D * Project2(Int_t, Int_t, ECellValue, PDEFoamKernelBase *, UInt_t) override
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.
Implementation of PDEFoam.
Definition PDEFoam.h:79
MsgLogger & Log() const
Definition PDEFoam.h:240
Basic string class.
Definition TString.h:138
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition TString.cxx:2384
TH1F * h1
Definition legend1.C:5
MsgLogger & Endl(MsgLogger &ml)
Definition MsgLogger.h:148
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:673