Logo ROOT  
Reference Guide
PDF.cxx
Go to the documentation of this file.
1// @(#)root/tmva $Id$
2// Author: Asen Christov, Andreas Hoecker, Joerg Stelzer, Helge Voss, Kai Voss, Jan Therhaag, Eckhard von Toerne
3
4/**********************************************************************************
5 * Project: TMVA - a Root-integrated toolkit for multivariate data analysis *
6 * Package: TMVA *
7 * Class : PDF *
8 * Web : http://tmva.sourceforge.net *
9 * *
10 * Description: *
11 * Implementation (see header for description) *
12 * *
13 * Authors (alphabetical): *
14 * Asen Christov <christov@physik.uni-freiburg.de> - Freiburg U., Germany *
15 * Andreas Hoecker <Andreas.Hocker@cern.ch> - CERN, Switzerland *
16 * Jan Therhaag <Jan.Therhaag@cern.ch> - U of Bonn, Germany *
17 * Eckhard von Toerne <evt@physik.uni-bonn.de> - U of Bonn, Germany *
18 * Helge Voss <Helge.Voss@cern.ch> - MPI-K Heidelberg, Germany *
19 * Kai Voss <Kai.Voss@cern.ch> - U. of Victoria, Canada *
20 * *
21 * Copyright (c) 2005-2011: *
22 * CERN, Switzerland *
23 * U. of Victoria, Canada *
24 * MPI-K Heidelberg, Germany *
25 * Freiburg U., Germany *
26 * U. of Bonn, Germany *
27 * *
28 * Redistribution and use in source and binary forms, with or without *
29 * modification, are permitted according to the terms listed in LICENSE *
30 * (http://tmva.sourceforge.net/LICENSE) *
31 **********************************************************************************/
32
33/*! \class TMVA::PDF
34\ingroup TMVA
35PDF wrapper for histograms; uses user-defined spline interpolation.
36*/
37
38#include "TMVA/PDF.h"
39
40#include "TMVA/Configurable.h"
41#include "TMVA/KDEKernel.h"
42#include "TMVA/MsgLogger.h"
43#include "TMVA/Types.h"
44#include "TMVA/Tools.h"
45#include "TMVA/TSpline1.h"
46#include "TMVA/TSpline2.h"
47#include "TMVA/Version.h"
48
49#include "Riostream.h"
50#include "TF1.h"
51#include "TH1F.h"
52#include "TMath.h"
53#include "TVectorD.h"
54
55#include <cstdlib>
56#include <iomanip>
57
58// static configuration settings
61const Double_t TMVA::PDF::fgEpsilon = 1.0e-12;
62
64
65////////////////////////////////////////////////////////////////////////////////
66/// default constructor needed for ROOT I/O
67
69: Configurable (""),
70 fUseHistogram ( kFALSE ),
71 fPDFName ( name ),
72 fNsmooth ( 0 ),
73 fMinNsmooth (-1 ),
74 fMaxNsmooth (-1 ),
75 fNSmoothHist ( 0 ),
76 fInterpolMethod( PDF::kSpline2 ),
77 fSpline ( 0 ),
78 fPDFHist ( 0 ),
79 fHist ( 0 ),
80 fHistOriginal ( 0 ),
81 fGraph ( 0 ),
82 fIGetVal ( 0 ),
83 fHistAvgEvtPerBin ( 0 ),
84 fHistDefinedNBins ( 0 ),
85 fKDEtypeString ( 0 ),
86 fKDEiterString ( 0 ),
87 fBorderMethodString( 0 ),
88 fInterpolateString ( 0 ),
89 fKDEtype ( KDEKernel::kNone ),
90 fKDEiter ( KDEKernel::kNonadaptiveKDE ),
91 fKDEborder ( KDEKernel::kNoTreatment ),
92 fFineFactor ( 0. ),
93 fReadingVersion( 0 ),
94 fCheckHist ( kFALSE ),
95 fNormalize ( norm ),
96 fSuffix ( "" ),
97 fLogger ( 0 )
98{
99 fLogger = new MsgLogger(this);
100 GetThisPdfThreadLocal() = this;
101}
102
103////////////////////////////////////////////////////////////////////////////////
104/// constructor of spline based PDF:
105
107 const TH1 *hist,
109 Int_t minnsmooth,
110 Int_t maxnsmooth,
111 Bool_t checkHist,
112 Bool_t norm) :
113 Configurable (""),
114 fUseHistogram ( kFALSE ),
115 fPDFName ( name ),
116 fMinNsmooth ( minnsmooth ),
117 fMaxNsmooth ( maxnsmooth ),
118 fNSmoothHist ( 0 ),
119 fInterpolMethod( method ),
120 fSpline ( 0 ),
121 fPDFHist ( 0 ),
122 fHist ( 0 ),
123 fHistOriginal ( 0 ),
124 fGraph ( 0 ),
125 fIGetVal ( 0 ),
126 fHistAvgEvtPerBin ( 0 ),
127 fHistDefinedNBins ( 0 ),
128 fKDEtypeString ( 0 ),
129 fKDEiterString ( 0 ),
130 fBorderMethodString( 0 ),
131 fInterpolateString ( 0 ),
132 fKDEtype ( KDEKernel::kNone ),
133 fKDEiter ( KDEKernel::kNonadaptiveKDE ),
134 fKDEborder ( KDEKernel::kNoTreatment ),
135 fFineFactor ( 0. ),
136 fReadingVersion( 0 ),
137 fCheckHist ( checkHist ),
138 fNormalize ( norm ),
139 fSuffix ( "" ),
140 fLogger ( 0 )
141{
142 fLogger = new MsgLogger(this);
143 BuildPDF( hist );
144}
145
146////////////////////////////////////////////////////////////////////////////////
147/// constructor of kernel based PDF:
148
150 const TH1* hist,
154 Float_t FineFactor,
155 Bool_t norm) :
156 Configurable (""),
157 fUseHistogram ( kFALSE ),
158 fPDFName ( name ),
159 fNsmooth ( 0 ),
160 fMinNsmooth (-1 ),
161 fMaxNsmooth (-1 ),
162 fNSmoothHist ( 0 ),
163 fInterpolMethod( PDF::kKDE ),
164 fSpline ( 0 ),
165 fPDFHist ( 0 ),
166 fHist ( 0 ),
167 fHistOriginal ( 0 ),
168 fGraph ( 0 ),
169 fIGetVal ( 0 ),
170 fHistAvgEvtPerBin ( 0 ),
171 fHistDefinedNBins ( 0 ),
172 fKDEtypeString ( 0 ),
173 fKDEiterString ( 0 ),
174 fBorderMethodString( 0 ),
175 fInterpolateString ( 0 ),
176 fKDEtype ( ktype ),
177 fKDEiter ( kiter ),
178 fKDEborder ( kborder ),
179 fFineFactor ( FineFactor),
180 fReadingVersion( 0 ),
181 fCheckHist ( kFALSE ),
182 fNormalize ( norm ),
183 fSuffix ( "" ),
184 fLogger ( 0 )
185{
186 fLogger = new MsgLogger(this);
187 BuildPDF( hist );
188}
189
190////////////////////////////////////////////////////////////////////////////////
191
193 const TString& options,
194 const TString& suffix,
195 PDF* defaultPDF,
196 Bool_t norm) :
197 Configurable (options),
198 fUseHistogram ( kFALSE ),
199 fPDFName ( name ),
200 fNsmooth ( 0 ),
201 fMinNsmooth ( -1 ),
202 fMaxNsmooth ( -1 ),
203 fNSmoothHist ( 0 ),
204 fInterpolMethod( PDF::kSpline0 ),
205 fSpline ( 0 ),
206 fPDFHist ( 0 ),
207 fHist ( 0 ),
208 fHistOriginal ( 0 ),
209 fGraph ( 0 ),
210 fIGetVal ( 0 ),
211 fHistAvgEvtPerBin ( 50 ),
212 fHistDefinedNBins ( 0 ),
213 fKDEtypeString ( "Gauss" ),
214 fKDEiterString ( "Nonadaptive" ),
215 fBorderMethodString( "None" ),
216 fInterpolateString ( "Spline2" ),
217 fKDEtype ( KDEKernel::kNone ),
218 fKDEiter ( KDEKernel::kNonadaptiveKDE ),
219 fKDEborder ( KDEKernel::kNoTreatment ),
220 fFineFactor ( 1. ),
221 fReadingVersion( 0 ),
222 fCheckHist ( kFALSE ),
223 fNormalize ( norm ),
224 fSuffix ( suffix ),
225 fLogger ( 0 )
226{
227 fLogger = new MsgLogger(this);
228 if (defaultPDF != 0) {
229 fNsmooth = defaultPDF->fNsmooth;
230 fMinNsmooth = defaultPDF->fMinNsmooth;
231 fMaxNsmooth = defaultPDF->fMaxNsmooth;
235 fKDEtypeString = defaultPDF->fKDEtypeString;
236 fKDEiterString = defaultPDF->fKDEiterString;
237 fFineFactor = defaultPDF->fFineFactor;
239 fCheckHist = defaultPDF->fCheckHist;
241 }
242}
243
244//___________________fNSmoothHist____________________________________________________
246{
247 // destructor
248 if (fSpline != NULL) delete fSpline;
249 if (fHist != NULL) delete fHist;
250 if (fPDFHist != NULL) delete fPDFHist;
251 if (fHistOriginal != NULL) delete fHistOriginal;
252 if (fIGetVal != NULL) delete fIGetVal;
253 if (fGraph != NULL) delete fGraph;
254 delete fLogger;
255}
256
257////////////////////////////////////////////////////////////////////////////////
258
259void TMVA::PDF::BuildPDF( const TH1* hist )
260{
261 GetThisPdfThreadLocal() = this;
262 // sanity check
263 if (hist == NULL) Log() << kFATAL << "Called without valid histogram pointer!" << Endl;
264
265 // histogram should be non empty
266 if (hist->GetEntries() <= 0)
267 Log() << kFATAL << "Number of entries <= 0 (" << hist->GetEntries() << " in histogram: " << hist->GetTitle() << ")" << Endl;
268
269 if (fInterpolMethod == PDF::kKDE) {
270 Log()<< kDEBUG << "Create "
271 << ((fKDEiter == KDEKernel::kNonadaptiveKDE) ? "nonadaptive " :
272 (fKDEiter == KDEKernel::kAdaptiveKDE) ? "adaptive " : "??? ")
273 << ((fKDEtype == KDEKernel::kGauss) ? "Gauss " : "??? ")
274 << "type KDE kernel for histogram: \"" << hist->GetName() << "\""
275 << Endl;
276 }
277 else {
278 // another sanity check (nsmooth<0 indicated build with KDE)
279 if (fMinNsmooth<0)
280 Log() << kFATAL << "PDF construction called with minnsmooth<0" << Endl;
281 else if (fMaxNsmooth<=0)
282 fMaxNsmooth = fMinNsmooth;
283 else if (fMaxNsmooth<fMinNsmooth)
284 Log() << kFATAL << "PDF construction called with maxnsmooth<minnsmooth" << Endl;
285 }
286
287 fHistOriginal = (TH1F*)hist->Clone( TString(hist->GetName()) + "_original" );
288 fHist = (TH1F*)hist->Clone( TString(hist->GetName()) + "_smoothed" );
289 fHistOriginal->SetTitle( fHistOriginal->GetName() ); // reset to new title as well
290 fHist ->SetTitle( fHist->GetName() );
291
292 // do not store in current target file
293 fHistOriginal->SetDirectory(0);
294 fHist ->SetDirectory(0);
295 fUseHistogram = kFALSE;
296
297 if (fInterpolMethod == PDF::kKDE) BuildKDEPDF();
298 else BuildSplinePDF();
299}
300
301////////////////////////////////////////////////////////////////////////////////
302
304{
305 Int_t ResolutionFactor = (fInterpolMethod == PDF::kKDE) ? 5 : 1;
306 if (evtNum == 0 && fHistDefinedNBins == 0)
307 Log() << kFATAL << "No number of bins set for PDF" << Endl;
308 else if (fHistDefinedNBins > 0)
309 return fHistDefinedNBins * ResolutionFactor;
310 else if ( evtNum > 0 && fHistAvgEvtPerBin > 0 )
311 return evtNum / fHistAvgEvtPerBin * ResolutionFactor;
312 else
313 Log() << kFATAL << "No number of bins or average event per bin set for PDF" << fHistAvgEvtPerBin << Endl;
314 return 0;
315}
316
317////////////////////////////////////////////////////////////////////////////////
318/// build the PDF from the original histograms
319
321{
322 // (not useful for discrete distributions, or if no splines are requested)
323 if (fInterpolMethod != PDF::kSpline0 && fCheckHist) CheckHist();
324 // use ROOT TH1 smooth method
325 fNSmoothHist = 0;
326 if (fMaxNsmooth > 0 && fMinNsmooth >=0 ) SmoothHistogram();
327
328 // fill histogramm to graph
329 FillHistToGraph();
330
331 // fGraph->Print();
332 switch (fInterpolMethod) {
333
334 case kSpline0:
335 // use original histogram as reference
336 // this is useful, eg, for discrete variables
337 fUseHistogram = kTRUE;
338 break;
339
340 case kSpline1:
341 fSpline = new TMVA::TSpline1( "spline1", new TGraph(*fGraph) );
342 break;
343
344 case kSpline2:
345 fSpline = new TMVA::TSpline2( "spline2", new TGraph(*fGraph) );
346 break;
347
348 case kSpline3:
349 fSpline = new TSpline3( "spline3", new TGraph(*fGraph) );
350 break;
351
352 case kSpline5:
353 fSpline = new TSpline5( "spline5", new TGraph(*fGraph) );
354 break;
355
356 default:
357 Log() << kWARNING << "No valid interpolation method given! Use Spline2" << Endl;
358 fSpline = new TMVA::TSpline2( "spline2", new TGraph(*fGraph) );
359 Log() << kFATAL << " Well.. .thinking about it, I better quit so you notice you are forced to fix the mistake " << Endl;
360 std::exit(1);
361 }
362 // fill into histogram
363 FillSplineToHist();
364
365 if (!UseHistogram()) {
366 fSpline->SetTitle( (TString)fHist->GetTitle() + fSpline->GetTitle() );
367 fSpline->SetName ( (TString)fHist->GetName() + fSpline->GetName() );
368 }
369
370
371 // sanity check
372 Double_t integral = GetIntegral();
373 if (integral < 0) Log() << kFATAL << "Integral: " << integral << " <= 0" << Endl;
374
375 // normalize
376 if (fNormalize)
377 if (integral>0) fPDFHist->Scale( 1.0/integral );
378
379 fPDFHist->SetDirectory(0);
380}
381
382////////////////////////////////////////////////////////////////////////////////
383/// creates high-binned reference histogram to be used instead of the
384/// PDF for speed reasons
385
387{
388 fPDFHist = new TH1F( "", "", fgNbin_PdfHist, GetXmin(), GetXmax() );
389 fPDFHist->SetTitle( (TString)fHist->GetTitle() + "_hist from_KDE" );
390 fPDFHist->SetName ( (TString)fHist->GetName() + "_hist_from_KDE" );
391
392 // create the kernel object
393 Float_t histoLowEdge = fHist->GetBinLowEdge(1);
394 Float_t histoUpperEdge = fPDFHist->GetBinLowEdge(fPDFHist->GetNbinsX()) + fPDFHist->GetBinWidth(fPDFHist->GetNbinsX());
395 KDEKernel *kern = new TMVA::KDEKernel( fKDEiter,
396 fHist, histoLowEdge, histoUpperEdge,
397 fKDEborder, fFineFactor );
398 kern->SetKernelType(fKDEtype);
399
400 for (Int_t i=1;i<fHist->GetNbinsX();i++) {
401 // loop over the bins of the original histo
402 for (Int_t j=1;j<fPDFHist->GetNbinsX();j++) {
403 // loop over the bins of the new histo and fill it
404 fPDFHist->AddBinContent(j,fHist->GetBinContent(i)*
405 kern->GetBinKernelIntegral(fPDFHist->GetBinLowEdge(j),
406 fPDFHist->GetBinLowEdge(j+1),
407 fHist->GetBinCenter(i),
408 i)
409 );
410 }
411 if (fKDEborder == 3) { // mirror the samples and fill them again
412 // in order to save time do the mirroring only for the first (the lower) 1/5 of the histo to the left;
413 // and the last (the higher) 1/5 of the histo to the right.
414 // the middle part of the histo, which is not mirrored, has no influence on the border effects anyway ...
415 if (i < fHist->GetNbinsX()/5 ) { // the first (the lower) 1/5 of the histo
416 for (Int_t j=1;j<fPDFHist->GetNbinsX();j++) {
417 // loop over the bins of the PDF histo and fill it
418 fPDFHist->AddBinContent(j,fHist->GetBinContent(i)*
419 kern->GetBinKernelIntegral(fPDFHist->GetBinLowEdge(j),
420 fPDFHist->GetBinLowEdge(j+1),
421 2*histoLowEdge-fHist->GetBinCenter(i), // mirroring to the left
422 i)
423 );
424 }
425 }
426 if (i > 4*fHist->GetNbinsX()/5) { // the last (the higher) 1/5 of the histo
427 for (Int_t j=1;j<fPDFHist->GetNbinsX();j++) {
428 // loop over the bins of the PDF histo and fill it
429 fPDFHist->AddBinContent( j,fHist->GetBinContent(i)*
430 kern->GetBinKernelIntegral(fPDFHist->GetBinLowEdge(j),
431 fPDFHist->GetBinLowEdge(j+1),
432 2*histoUpperEdge-fHist->GetBinCenter(i), i) );
433 }
434 }
435 }
436 }
437
438 fPDFHist->SetEntries(fHist->GetEntries());
439
440 delete kern;
441
442 // sanity check
443 Double_t integral = GetIntegral();
444 if (integral < 0) Log() << kFATAL << "Integral: " << integral << " <= 0" << Endl;
445
446 // normalize
447 if (fNormalize)
448 if (integral>0) fPDFHist->Scale( 1.0/integral );
449 fPDFHist->SetDirectory(0);
450}
451
452////////////////////////////////////////////////////////////////////////////////
453
455{
456 if(fHist->GetNbinsX()==1) return;
457 if (fMaxNsmooth == fMinNsmooth) {
458 fHist->Smooth( fMinNsmooth );
459 return;
460 }
461
462 //calculating Mean, RMS of the relative errors and using them to set
463 //the boundaries of the linear transformation
464 Float_t Err=0, ErrAvg=0, ErrRMS=0 ; Int_t num=0, smooth;
465 for (Int_t bin=0; bin<fHist->GetNbinsX(); bin++) {
466 if (fHist->GetBinContent(bin+1) <= fHist->GetBinError(bin+1)) continue;
467 Err = fHist->GetBinError(bin+1) / fHist->GetBinContent(bin+1);
468 ErrAvg += Err; ErrRMS += Err*Err; num++;
469 }
470 ErrAvg /= num;
471 ErrRMS = TMath::Sqrt(ErrRMS/num-ErrAvg*ErrAvg) ;
472
473 //linearly convent the histogram to a vector of smoothnum
474 Float_t MaxErr=ErrAvg+ErrRMS, MinErr=ErrAvg-ErrRMS;
475 fNSmoothHist = new TH1I("","",fHist->GetNbinsX(),0,fHist->GetNbinsX());
476 fNSmoothHist->SetTitle( (TString)fHist->GetTitle() + "_Nsmooth" );
477 fNSmoothHist->SetName ( (TString)fHist->GetName() + "_Nsmooth" );
478 for (Int_t bin=0; bin<fHist->GetNbinsX(); bin++) {
479 if (fHist->GetBinContent(bin+1) <= fHist->GetBinError(bin+1))
480 smooth=fMaxNsmooth;
481 else {
482 Err = fHist->GetBinError(bin+1) / fHist->GetBinContent(bin+1);
483 smooth=(Int_t)((Err-MinErr) /(MaxErr-MinErr) * (fMaxNsmooth-fMinNsmooth)) + fMinNsmooth;
484 }
485 smooth = TMath::Max(fMinNsmooth,TMath::Min(fMaxNsmooth,smooth));
486 fNSmoothHist->SetBinContent(bin+1,smooth);
487 }
488
489 //find regions of constant smoothnum, starting from the highest amount of
490 //smoothing. So the last iteration all the histogram will be smoothed as a whole
491 for (Int_t n=fMaxNsmooth; n>=0; n--) {
492 //all the histogram has to be smoothed at least fMinNsmooth
493 if (n <= fMinNsmooth) { fHist->Smooth(); continue; }
494 Int_t MinBin=-1,MaxBin =-1;
495 for (Int_t bin=0; bin < fHist->GetNbinsX(); bin++) {
496 if (fNSmoothHist->GetBinContent(bin+1) >= n) {
497 if (MinBin==-1) MinBin = bin;
498 else MaxBin=bin;
499 }
500 else if (MaxBin >= 0) {
501 fHist->Smooth(1,"R");
502 MinBin=MaxBin=-1;
503 }
504 else // can't smooth a single bin
505 MinBin = -1;
506 }
507 }
508}
509
510////////////////////////////////////////////////////////////////////////////////
511/// Simple conversion
512
514{
515 fGraph=new TGraph(fHist);
516}
517
518////////////////////////////////////////////////////////////////////////////////
519/// creates high-binned reference histogram to be used instead of the
520/// PDF for speed reasons
521
523{
524 if (UseHistogram()) {
525 // no spline given, use the original histogram
526 fPDFHist = (TH1*)fHist->Clone();
527 fPDFHist->SetTitle( (TString)fHist->GetTitle() + "_hist from_spline0" );
528 fPDFHist->SetName ( (TString)fHist->GetName() + "_hist_from_spline0" );
529 }
530 else {
531 // create new reference histogram
532 fPDFHist = new TH1F( "", "", fgNbin_PdfHist, GetXmin(), GetXmax() );
533 fPDFHist->SetTitle( (TString)fHist->GetTitle() + "_hist from_" + fSpline->GetTitle() );
534 fPDFHist->SetName ( (TString)fHist->GetName() + "_hist_from_" + fSpline->GetTitle() );
535
536 for (Int_t bin=1; bin <= fgNbin_PdfHist; bin++) {
537 Double_t x = fPDFHist->GetBinCenter( bin );
538 Double_t y = fSpline->Eval( x );
539 // sanity correction: in cases where strong slopes exist, accidentally, the
540 // splines can go to zero; in this case we set the corresponding bin content
541 // equal to the bin content of the original histogram
542 if (y <= fgEpsilon) y = fHist->GetBinContent( fHist->FindBin( x ) );
543 fPDFHist->SetBinContent( bin, TMath::Max(y, fgEpsilon) );
544 }
545 }
546 fPDFHist->SetDirectory(0);
547}
548
549////////////////////////////////////////////////////////////////////////////////
550/// sanity check: compare PDF with original histogram
551
553{
554 if (fHist == NULL) {
555 Log() << kFATAL << "<CheckHist> Called without valid histogram pointer!" << Endl;
556 }
557
558 Int_t nbins = fHist->GetNbinsX();
559
560 Int_t emptyBins=0;
561 // count number of empty bins
562 for (Int_t bin=1; bin<=nbins; bin++)
563 if (fHist->GetBinContent(bin) == 0) emptyBins++;
564
565 if (((Float_t)emptyBins/(Float_t)nbins) > 0.5) {
566 Log() << kWARNING << "More than 50% (" << (((Float_t)emptyBins/(Float_t)nbins)*100)
567 <<"%) of the bins in hist '"
568 << fHist->GetName() << "' are empty!" << Endl;
569 Log() << kWARNING << "X_min=" << GetXmin()
570 << " mean=" << fHist->GetMean() << " X_max= " << GetXmax() << Endl;
571 }
572}
573
574////////////////////////////////////////////////////////////////////////////////
575/// comparison of original histogram with reference PDF
576
577void TMVA::PDF::ValidatePDF( TH1* originalHist ) const
578{
579 // if no histogram is given, use the original one (the one the PDF was made from)
580 if (!originalHist) originalHist = fHistOriginal;
581
582 Int_t nbins = originalHist->GetNbinsX();
583
584 // treat errors properly
585 if (originalHist->GetSumw2()->GetSize() == 0) originalHist->Sumw2();
586
587 // ---- first validation: simple(st) possible chi2 test
588 // count number of empty bins
589 Double_t chi2 = 0;
590 Int_t ndof = 0;
591 Int_t nc1 = 0; // deviation counters
592 Int_t nc2 = 0; // deviation counters
593 Int_t nc3 = 0; // deviation counters
594 Int_t nc6 = 0; // deviation counters
595 for (Int_t bin=1; bin<=nbins; bin++) {
596 Double_t x = originalHist->GetBinCenter( bin );
597 Double_t y = originalHist->GetBinContent( bin );
598 Double_t ey = originalHist->GetBinError( bin );
599
600 Int_t binPdfHist = fPDFHist->FindBin( x );
601 if (binPdfHist<0) continue; // happens only if hist-dim>3
602
603 Double_t yref = GetVal( x );
604 Double_t rref = ( originalHist->GetSumOfWeights()/fPDFHist->GetSumOfWeights() *
605 originalHist->GetBinWidth( bin )/fPDFHist->GetBinWidth( binPdfHist ) );
606
607 if (y > 0) {
608 ndof++;
609 Double_t d = TMath::Abs( (y - yref*rref)/ey );
610 // std::cout << "bin: " << bin << " val: " << x << " data(err): " << y << "(" << ey << ") pdf: "
611 // << yref << " dev(chi2): " << d << "(" << chi2 << ") rref: " << rref << std::endl;
612 chi2 += d*d;
613 if (d > 1) { nc1++; if (d > 2) { nc2++; if (d > 3) { nc3++; if (d > 6) nc6++; } } }
614 }
615 }
616
617 Log() << kDEBUG << "Validation result for PDF \"" << originalHist->GetTitle() << "\"" << ": " << Endl;
618 Log() << kDEBUG << Form( " chi2/ndof(!=0) = %.1f/%i = %.2f (Prob = %.2f)",
619 chi2, ndof, chi2/ndof, TMath::Prob( chi2, ndof ) ) << Endl;
620 if ((1.0 - TMath::Prob( chi2, ndof )) > 0.9999994) {
621 Log() << kDEBUG << "Comparison of the original histogram \"" << originalHist->GetTitle() << "\"" << Endl;
622 Log() << kDEBUG << "with the corresponding PDF gave a chi2/ndof of " << chi2/ndof << "," << Endl;
623 Log() << kDEBUG << "which corresponds to a deviation of more than 5 sigma! Please check!" << Endl;
624 }
625 Log() << kDEBUG << Form( " #bins-found(#expected-bins) deviating > [1,2,3,6] sigmas: " \
626 "[%i(%i),%i(%i),%i(%i),%i(%i)]",
627 nc1, Int_t(TMath::Prob(1.0,1)*ndof), nc2, Int_t(TMath::Prob(4.0,1)*ndof),
628 nc3, Int_t(TMath::Prob(9.0,1)*ndof), nc6, Int_t(TMath::Prob(36.0,1)*ndof) ) << Endl;
629}
630
631////////////////////////////////////////////////////////////////////////////////
632/// computes normalisation
633
635{
636 Double_t integral = fPDFHist->GetSumOfWeights();
637 integral *= GetPdfHistBinWidth();
638
639 return integral;
640}
641
642////////////////////////////////////////////////////////////////////////////////
643/// static external auxiliary function (integrand)
644
646{
647 return ThisPDF()->GetVal( x[0] );
648}
649
650////////////////////////////////////////////////////////////////////////////////
651/// computes PDF integral within given ranges
652
654{
655 Double_t integral = 0;
656
657 if (fgManualIntegration) {
658
659 // compute integral by summing over bins
660 Int_t imin = fPDFHist->FindBin(xmin);
661 Int_t imax = fPDFHist->FindBin(xmax);
662 if (imin < 1) imin = 1;
663 if (imax > fPDFHist->GetNbinsX()) imax = fPDFHist->GetNbinsX();
664
665 for (Int_t bini = imin; bini <= imax; bini++) {
666 Float_t dx = fPDFHist->GetBinWidth(bini);
667 // correct for bin fractions in first and last bin
668 if (bini == imin) dx = fPDFHist->GetBinLowEdge(bini+1) - xmin;
669 else if (bini == imax) dx = xmax - fPDFHist->GetBinLowEdge(bini);
670 if (dx < 0 && dx > -1.0e-8) dx = 0; // protect against float->double numerical effects
671 if (dx<0) {
672 Log() << kWARNING
673 << "dx = " << dx << std::endl
674 << "bini = " << bini << std::endl
675 << "xmin = " << xmin << std::endl
676 << "xmax = " << xmax << std::endl
677 << "imin = " << imin << std::endl
678 << "imax = " << imax << std::endl
679 << "low edge of imin" << fPDFHist->GetBinLowEdge(imin) << std::endl
680 << "low edge of imin+1" << fPDFHist->GetBinLowEdge(imin+1) << Endl;
681 Log() << kFATAL << "<GetIntegral> dx = " << dx << " < 0" << Endl;
682 }
683 integral += fPDFHist->GetBinContent(bini)*dx;
684 }
685
686 }
687 else {
688
689 // compute via Gaussian quadrature (C++ version of CERNLIB function DGAUSS)
690 if (fIGetVal == 0) fIGetVal = new TF1( "IGetVal", PDF::IGetVal, GetXmin(), GetXmax(), 0 );
691 integral = fIGetVal->Integral( xmin, xmax );
692 }
693
694 return integral;
695}
696
697////////////////////////////////////////////////////////////////////////////////
698/// returns value PDF(x)
699
701{
702 // check which is filled
703 Int_t bin = fPDFHist->FindBin(x);
704 bin = TMath::Max(bin,1);
705 bin = TMath::Min(bin,fPDFHist->GetNbinsX());
706
707 Double_t retval = 0;
708
709 if (UseHistogram()) {
710 // use directly histogram bins (this is for discrete PDFs)
711 retval = fPDFHist->GetBinContent( bin );
712 }
713 else {
714 // linear interpolation
715 Int_t nextbin = bin;
716 if ((x > fPDFHist->GetBinCenter(bin) && bin != fPDFHist->GetNbinsX()) || bin == 1)
717 nextbin++;
718 else
719 nextbin--;
720
721 // linear interpolation between adjacent bins
722 Double_t dx = fPDFHist->GetBinCenter( bin ) - fPDFHist->GetBinCenter( nextbin );
723 Double_t dy = fPDFHist->GetBinContent( bin ) - fPDFHist->GetBinContent( nextbin );
724 retval = fPDFHist->GetBinContent( bin ) + (x - fPDFHist->GetBinCenter( bin ))*dy/dx;
725 }
726
727 return TMath::Max( retval, fgEpsilon );
728}
729
730////////////////////////////////////////////////////////////////////////////////
731/// returns value \f$ PDF^{-1}(y) \f$
732
733Double_t TMVA::PDF::GetValInverse( Double_t y, Bool_t isMonotonouslyIncreasingFunction ) const
734{
735 Int_t lowerBin=0, higherBin=0;
736 Double_t lowerBinValue=0, higherBinValue=0;
737 FindBinInverse(fPDFHist,lowerBin,higherBin,lowerBinValue,higherBinValue,y,isMonotonouslyIncreasingFunction);
738
739 Double_t xValueLowerBin =fPDFHist->GetBinCenter (lowerBin);
740 Double_t xValueHigherBin=fPDFHist->GetBinCenter (higherBin);
741
742 Double_t length =(higherBinValue-lowerBinValue);
743 Double_t fraction=lowerBinValue;
744 if (length>1.0e-10)
745 fraction=(y-lowerBinValue)/length;
746
747 Double_t lengthX =xValueHigherBin-xValueLowerBin;
748 Double_t x =xValueLowerBin+lengthX*fraction;
749
750 // comparison for test purposes
751 // std::cout << "lb " << lowerBin << " hb " << higherBin << " lbv " << lowerBinValue << " hbv " << higherBinValue << " frac " << fraction << std::endl;
752 // std::cout << "y " << y << " inv x " << x << " straight y " << GetVal(x) << std::endl;
753
754 return x;
755}
756
757////////////////////////////////////////////////////////////////////////////////
758/// find bin from value on ordinate
759
760void TMVA::PDF::FindBinInverse( const TH1* histogram, Int_t& lowerBin, Int_t& higherBin, Double_t& lowerBinValue, Double_t& higherBinValue,
761 Double_t y, Bool_t isMonotonouslyIncreasingFunction ) const
762{
763 if (isMonotonouslyIncreasingFunction) {
764 higherBin=histogram->GetNbinsX();
765 lowerBin =0;
766
767 Int_t bin=higherBin/2;
768
769 while (bin>lowerBin && bin<higherBin) {
770 Double_t binContent=histogram->GetBinContent(bin);
771
772 if (y<binContent) {
773 higherBin =bin;
774 higherBinValue=binContent;
775 }
776 else if (y>=binContent){
777 lowerBin =bin;
778 lowerBinValue =binContent;
779 }
780 bin=lowerBin+(higherBin-lowerBin)/2;
781 }
782 return;
783 }
784 // search sequentially
785 for (Int_t bin=0, binEnd=histogram->GetNbinsX(); bin<binEnd; ++bin) {
786 Double_t binContent=histogram->GetBinContent(bin);
787 if (binContent<y) {
788 lowerBin =bin;
789 higherBin=bin;
790 lowerBinValue =binContent;
791 higherBinValue=binContent;
792 }
793 else {
794 higherBin=bin;
795 higherBinValue=binContent;
796 break;
797 }
798 }
799}
800
801
802////////////////////////////////////////////////////////////////////////////////
803/// define the options (their key words) that can be set in the option string
804///
805/// know options:
806///
807/// PDFInterpol[ivar] <string> Spline0, Spline1, Spline2 <default>, Spline3, Spline5, KDE used to interpolate reference histograms
808/// if no variable index is given, it is valid for ALL the variables
809///
810/// - NSmooth <int> how often the input histos are smoothed
811/// - MinNSmooth <int> min number of smoothing iterations, for bins with most data
812/// - MaxNSmooth <int> max number of smoothing iterations, for bins with least data
813/// - NAvEvtPerBin <int> minimum average number of events per PDF bin
814/// - TransformOutput <bool> transform (often strongly peaked) likelihood output through sigmoid inversion
815/// - fKDEtype <KernelType> type of the Kernel to use (1 is Gaussian)
816/// - fKDEiter <KerneIter> number of iterations (1 --> "static KDE", 2 --> "adaptive KDE")
817/// - fBorderMethod <KernelBorder> the method to take care about "border" effects (1=no treatment , 2=kernel renormalization, 3=sample mirroring)
818
820{
821 DeclareOptionRef( fNsmooth, Form("NSmooth%s",fSuffix.Data()),
822 "Number of smoothing iterations for the input histograms" );
823 DeclareOptionRef( fMinNsmooth, Form("MinNSmooth%s",fSuffix.Data()),
824 "Min number of smoothing iterations, for bins with most data" );
825
826 DeclareOptionRef( fMaxNsmooth, Form("MaxNSmooth%s",fSuffix.Data()),
827 "Max number of smoothing iterations, for bins with least data" );
828
829 DeclareOptionRef( fHistAvgEvtPerBin, Form("NAvEvtPerBin%s",fSuffix.Data()),
830 "Average number of events per PDF bin" );
831
832 DeclareOptionRef( fHistDefinedNBins, Form("Nbins%s",fSuffix.Data()),
833 "Defined number of bins for the histogram from which the PDF is created" );
834
835 DeclareOptionRef( fCheckHist, Form("CheckHist%s",fSuffix.Data()),
836 "Whether or not to check the source histogram of the PDF" );
837
838 DeclareOptionRef( fInterpolateString, Form("PDFInterpol%s",fSuffix.Data()),
839 "Interpolation method for reference histograms (e.g. Spline2 or KDE)" );
840 AddPreDefVal(TString("Spline0")); // take histogram
841 AddPreDefVal(TString("Spline1")); // linear interpolation between bins
842 AddPreDefVal(TString("Spline2")); // quadratic interpolation
843 AddPreDefVal(TString("Spline3")); // cubic interpolation
844 AddPreDefVal(TString("Spline5")); // fifth order polynom interpolation
845 AddPreDefVal(TString("KDE")); // use kernel density estimator
846
847 DeclareOptionRef( fKDEtypeString, Form("KDEtype%s",fSuffix.Data()), "KDE kernel type (1=Gauss)" );
848 AddPreDefVal(TString("Gauss"));
849
850 DeclareOptionRef( fKDEiterString, Form("KDEiter%s",fSuffix.Data()), "Number of iterations (1=non-adaptive, 2=adaptive)" );
851 AddPreDefVal(TString("Nonadaptive"));
852 AddPreDefVal(TString("Adaptive"));
853
854 DeclareOptionRef( fFineFactor , Form("KDEFineFactor%s",fSuffix.Data()),
855 "Fine tuning factor for Adaptive KDE: Factor to multiply the width of the kernel");
856
857 DeclareOptionRef( fBorderMethodString, Form("KDEborder%s",fSuffix.Data()),
858 "Border effects treatment (1=no treatment , 2=kernel renormalization, 3=sample mirroring)" );
859 AddPreDefVal(TString("None"));
860 AddPreDefVal(TString("Renorm"));
861 AddPreDefVal(TString("Mirror"));
862
863 SetConfigName( GetName() );
864 SetConfigDescription( "Configuration options for the PDF class" );
865}
866
867////////////////////////////////////////////////////////////////////////////////
868
870{
871 if (fNsmooth < 0) fNsmooth = 0; // set back to default
872
873 if (fMaxNsmooth < 0 || fMinNsmooth < 0) { // use "Nsmooth" variable
874 fMinNsmooth = fMaxNsmooth = fNsmooth;
875 }
876
877 if (fMaxNsmooth < fMinNsmooth && fMinNsmooth >= 0) { // sanity check
878 Log() << kFATAL << "ERROR: MaxNsmooth = "
879 << fMaxNsmooth << " < MinNsmooth = " << fMinNsmooth << Endl;
880 }
881
882 if (fMaxNsmooth < 0 || fMinNsmooth < 0) {
883 Log() << kFATAL << "ERROR: MaxNsmooth = "
884 << fMaxNsmooth << " or MinNsmooth = " << fMinNsmooth << " smaller than zero" << Endl;
885 }
886
887 //processing the options
888 if (fInterpolateString == "Spline0") fInterpolMethod = TMVA::PDF::kSpline0;
889 else if (fInterpolateString == "Spline1") fInterpolMethod = TMVA::PDF::kSpline1;
890 else if (fInterpolateString == "Spline2") fInterpolMethod = TMVA::PDF::kSpline2;
891 else if (fInterpolateString == "Spline3") fInterpolMethod = TMVA::PDF::kSpline3;
892 else if (fInterpolateString == "Spline5") fInterpolMethod = TMVA::PDF::kSpline5;
893 else if (fInterpolateString == "KDE" ) fInterpolMethod = TMVA::PDF::kKDE;
894 else if (fInterpolateString != "" ) {
895 Log() << kFATAL << "unknown setting for option 'InterpolateMethod': " << fKDEtypeString << ((fSuffix=="")?"":Form(" for pdf with suffix %s",fSuffix.Data())) << Endl;
896 }
897
898 // init KDE options
899 if (fKDEtypeString == "Gauss" ) fKDEtype = KDEKernel::kGauss;
900 else if (fKDEtypeString != "" )
901 Log() << kFATAL << "unknown setting for option 'KDEtype': " << fKDEtypeString << ((fSuffix=="")?"":Form(" for pdf with suffix %s",fSuffix.Data())) << Endl;
902 if (fKDEiterString == "Nonadaptive") fKDEiter = KDEKernel::kNonadaptiveKDE;
903 else if (fKDEiterString == "Adaptive" ) fKDEiter = KDEKernel::kAdaptiveKDE;
904 else if (fKDEiterString != "" )// nothing more known
905 Log() << kFATAL << "unknown setting for option 'KDEiter': " << fKDEtypeString << ((fSuffix=="")?"":Form(" for pdf with suffix %s",fSuffix.Data())) << Endl;
906
907 if ( fBorderMethodString == "None" ) fKDEborder= KDEKernel::kNoTreatment;
908 else if ( fBorderMethodString == "Renorm" ) fKDEborder= KDEKernel::kKernelRenorm;
909 else if ( fBorderMethodString == "Mirror" ) fKDEborder= KDEKernel::kSampleMirror;
910 else if ( fKDEiterString != "" ) { // nothing more known
911 Log() << kFATAL << "unknown setting for option 'KDEBorder': " << fKDEtypeString << ((fSuffix=="")?"":Form(" for pdf with suffix %s",fSuffix.Data())) << Endl;
912 }
913}
914
915////////////////////////////////////////////////////////////////////////////////
916/// XML file writing
917
918void TMVA::PDF::AddXMLTo( void* parent )
919{
920 void* pdfxml = gTools().AddChild(parent, "PDF");
921 gTools().AddAttr(pdfxml, "Name", fPDFName );
922 gTools().AddAttr(pdfxml, "MinNSmooth", fMinNsmooth );
923 gTools().AddAttr(pdfxml, "MaxNSmooth", fMaxNsmooth );
924 gTools().AddAttr(pdfxml, "InterpolMethod", fInterpolMethod );
925 gTools().AddAttr(pdfxml, "KDE_type", fKDEtype );
926 gTools().AddAttr(pdfxml, "KDE_iter", fKDEiter );
927 gTools().AddAttr(pdfxml, "KDE_border", fKDEborder );
928 gTools().AddAttr(pdfxml, "KDE_finefactor", fFineFactor );
929 void* pdfhist = gTools().AddChild(pdfxml,"Histogram" );
930 TH1* histToWrite = GetOriginalHist();
931 Bool_t hasEquidistantBinning = gTools().HistoHasEquidistantBins(*histToWrite);
932 gTools().AddAttr(pdfhist, "Name", histToWrite->GetName() );
933 gTools().AddAttr(pdfhist, "NBins", histToWrite->GetNbinsX() );
934 gTools().AddAttr(pdfhist, "XMin", histToWrite->GetXaxis()->GetXmin() );
935 gTools().AddAttr(pdfhist, "XMax", histToWrite->GetXaxis()->GetXmax() );
936 gTools().AddAttr(pdfhist, "HasEquidistantBins", hasEquidistantBinning );
937
938 TString bincontent("");
939 for (Int_t i=0; i<histToWrite->GetNbinsX(); i++) {
940 bincontent += gTools().StringFromDouble(histToWrite->GetBinContent(i+1));
941 bincontent += " ";
942 }
943 gTools().AddRawLine(pdfhist, bincontent );
944
945 if (!hasEquidistantBinning) {
946 void* pdfhistbins = gTools().AddChild(pdfxml,"HistogramBinning" );
947 gTools().AddAttr(pdfhistbins, "NBins", histToWrite->GetNbinsX() );
948 TString binns("");
949 for (Int_t i=1; i<=histToWrite->GetNbinsX()+1; i++) {
950 binns += gTools().StringFromDouble(histToWrite->GetXaxis()->GetBinLowEdge(i));
951 binns += " ";
952 }
953 gTools().AddRawLine(pdfhistbins, binns );
954 }
955}
956
957////////////////////////////////////////////////////////////////////////////////
958/// XML file reading
959
960void TMVA::PDF::ReadXML( void* pdfnode )
961{
962 UInt_t enumint;
963
964 gTools().ReadAttr(pdfnode, "MinNSmooth", fMinNsmooth );
965 gTools().ReadAttr(pdfnode, "MaxNSmooth", fMaxNsmooth );
966 gTools().ReadAttr(pdfnode, "InterpolMethod", enumint ); fInterpolMethod = EInterpolateMethod(enumint);
967 gTools().ReadAttr(pdfnode, "KDE_type", enumint ); fKDEtype = KDEKernel::EKernelType(enumint);
968 gTools().ReadAttr(pdfnode, "KDE_iter", enumint ); fKDEiter = KDEKernel::EKernelIter(enumint);
969 gTools().ReadAttr(pdfnode, "KDE_border", enumint ); fKDEborder = KDEKernel::EKernelBorder(enumint);
970 gTools().ReadAttr(pdfnode, "KDE_finefactor", fFineFactor );
971
972 TString hname;
973 UInt_t nbins;
975 Bool_t hasEquidistantBinning;
976
977 void* histch = gTools().GetChild(pdfnode);
978 gTools().ReadAttr( histch, "Name", hname );
979 gTools().ReadAttr( histch, "NBins", nbins );
980 gTools().ReadAttr( histch, "XMin", xmin );
981 gTools().ReadAttr( histch, "XMax", xmax );
982 gTools().ReadAttr( histch, "HasEquidistantBins", hasEquidistantBinning );
983
984 // recreate the original hist
985 TH1* newhist = 0;
986 if (hasEquidistantBinning) {
987 newhist = new TH1F( hname, hname, nbins, xmin, xmax );
988 newhist->SetDirectory(0);
989 const char* content = gTools().GetContent(histch);
990 std::stringstream s(content);
991 Double_t val;
992 for (UInt_t i=0; i<nbins; i++) {
993 s >> val;
994 newhist->SetBinContent(i+1,val);
995 }
996 }
997 else {
998 const char* content = gTools().GetContent(histch);
999 std::stringstream s(content);
1000 Double_t val;
1001 void* binch = gTools().GetNextChild(histch);
1002 UInt_t nbinning;
1003 gTools().ReadAttr( binch, "NBins", nbinning );
1004 TVectorD binns(nbinning+1);
1005 if (nbinning != nbins) {
1006 Log() << kFATAL << "Number of bins in content and binning array differs"<<Endl;
1007 }
1008 const char* binString = gTools().GetContent(binch);
1009 std::stringstream sb(binString);
1010 for (UInt_t i=0; i<=nbins; i++) sb >> binns[i];
1011 newhist = new TH1F( hname, hname, nbins, binns.GetMatrixArray() );
1012 newhist->SetDirectory(0);
1013 for (UInt_t i=0; i<nbins; i++) {
1014 s >> val;
1015 newhist->SetBinContent(i+1,val);
1016 }
1017 }
1018
1019 TString hnameSmooth = hname;
1020 hnameSmooth.ReplaceAll( "_original", "_smoothed" );
1021
1022 if (fHistOriginal != 0) delete fHistOriginal;
1023 fHistOriginal = newhist;
1024 fHist = (TH1F*)fHistOriginal->Clone( hnameSmooth );
1025 fHist->SetTitle( hnameSmooth );
1026 fHist->SetDirectory(0);
1027
1028 if (fInterpolMethod == PDF::kKDE) BuildKDEPDF();
1029 else BuildSplinePDF();
1030}
1031
1032////////////////////////////////////////////////////////////////////////////////
1033/// write the pdf
1034
1035std::ostream& TMVA::operator<< ( std::ostream& os, const PDF& pdf )
1036{
1037 Int_t dp = os.precision();
1038 os << "MinNSmooth " << pdf.fMinNsmooth << std::endl;
1039 os << "MaxNSmooth " << pdf.fMaxNsmooth << std::endl;
1040 os << "InterpolMethod " << pdf.fInterpolMethod << std::endl;
1041 os << "KDE_type " << pdf.fKDEtype << std::endl;
1042 os << "KDE_iter " << pdf.fKDEiter << std::endl;
1043 os << "KDE_border " << pdf.fKDEborder << std::endl;
1044 os << "KDE_finefactor " << pdf.fFineFactor << std::endl;
1045
1046 TH1* histToWrite = pdf.GetOriginalHist();
1047
1048 const Int_t nBins = histToWrite->GetNbinsX();
1049
1050 // note: this is a schema change introduced for v3.7.3
1051 os << "Histogram "
1052 << histToWrite->GetName()
1053 << " " << nBins // nbins
1054 << " " << std::setprecision(12) << histToWrite->GetXaxis()->GetXmin() // x_min
1055 << " " << std::setprecision(12) << histToWrite->GetXaxis()->GetXmax() // x_max
1056 << std::endl;
1057
1058 // write the smoothed hist
1059 os << "Weights " << std::endl;
1060 os << std::setprecision(8);
1061 for (Int_t i=0; i<nBins; i++) {
1062 os << std::setw(15) << std::left << histToWrite->GetBinContent(i+1) << std::right << " ";
1063 if ((i+1)%5==0) os << std::endl;
1064 }
1065
1066 os << std::setprecision(dp);
1067 return os; // Return the output stream.
1068}
1069
1070////////////////////////////////////////////////////////////////////////////////
1071/// read the tree from an std::istream
1072
1073std::istream& TMVA::operator>> ( std::istream& istr, PDF& pdf )
1074{
1075 TString devnullS;
1076 Int_t valI;
1077 Int_t nbins=-1; // default binning will cause an exit
1078 Float_t xmin=-1., xmax=-1.;
1079 TString hname="_original";
1080 Bool_t doneReading = kFALSE;
1081 while (!doneReading) {
1082 istr >> devnullS;
1083 if (devnullS=="NSmooth")
1084 {istr >> pdf.fMinNsmooth; pdf.fMaxNsmooth=pdf.fMinNsmooth;}
1085 else if (devnullS=="MinNSmooth") istr >> pdf.fMinNsmooth;
1086 else if (devnullS=="MaxNSmooth") istr >> pdf.fMaxNsmooth;
1087 // have to do this with strings to be more stable with developing code
1088 else if (devnullS == "InterpolMethod") { istr >> valI; pdf.fInterpolMethod = PDF::EInterpolateMethod(valI);}
1089 else if (devnullS == "KDE_type") { istr >> valI; pdf.fKDEtype = KDEKernel::EKernelType(valI); }
1090 else if (devnullS == "KDE_iter") { istr >> valI; pdf.fKDEiter = KDEKernel::EKernelIter(valI);}
1091 else if (devnullS == "KDE_border") { istr >> valI; pdf.fKDEborder = KDEKernel::EKernelBorder(valI);}
1092 else if (devnullS == "KDE_finefactor") {
1093 istr >> pdf.fFineFactor;
1094 if (pdf.GetReadingVersion() != 0 && pdf.GetReadingVersion() < TMVA_VERSION(3,7,3)) { // here we expect the histogram limits if the version is below 3.7.3. When version == 0, the newest TMVA version is assumed.
1095 // coverity[tainted_data_argument]
1096 istr >> nbins >> xmin >> xmax;
1097 doneReading = kTRUE;
1098 }
1099 }
1100 else if (devnullS == "Histogram") { istr >> hname >> nbins >> xmin >> xmax; }
1101 else if (devnullS == "Weights") { doneReading = kTRUE; }
1102 }
1103
1104 TString hnameSmooth = hname;
1105 hnameSmooth.ReplaceAll( "_original", "_smoothed" );
1106
1107 // recreate the original hist
1108 if (nbins==-1) {
1109 std::cout << "PDF, trying to create a histogram without defined binning"<< std::endl;
1110 std::exit(1);
1111 }
1112 TH1* newhist = new TH1F( hname,hname, nbins, xmin, xmax );
1113 newhist->SetDirectory(0);
1114 Float_t val;
1115 for (Int_t i=0; i<nbins; i++) {
1116 istr >> val;
1117 newhist->SetBinContent(i+1,val);
1118 }
1119
1120 if (pdf.fHistOriginal != 0) delete pdf.fHistOriginal;
1121 pdf.fHistOriginal = newhist;
1122 pdf.fHist = (TH1F*)pdf.fHistOriginal->Clone( hnameSmooth );
1123 pdf.fHist->SetTitle( hnameSmooth );
1124 pdf.fHist->SetDirectory(0);
1125
1126 if (pdf.fMinNsmooth>=0) pdf.BuildSplinePDF();
1127 else {
1128 pdf.fInterpolMethod = PDF::kKDE;
1129 pdf.BuildKDEPDF();
1130 }
1131
1132 return istr;
1133}
1134
1136{
1137 // return global "this" pointer of PDF
1138 return GetThisPdfThreadLocal();
1139}
const Handle_t kNone
Definition: GuiTypes.h:87
#define d(i)
Definition: RSha256.hxx:102
int Int_t
Definition: RtypesCore.h:43
unsigned int UInt_t
Definition: RtypesCore.h:44
const Bool_t kFALSE
Definition: RtypesCore.h:90
bool Bool_t
Definition: RtypesCore.h:61
double Double_t
Definition: RtypesCore.h:57
float Float_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
char name[80]
Definition: TGX11.cxx:109
float xmin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
char * Form(const char *fmt,...)
#define TMVA_VERSION(a, b, c)
Definition: Version.h:48
Int_t GetSize() const
Definition: TArray.h:47
Double_t GetXmax() const
Definition: TAxis.h:134
virtual Double_t GetBinLowEdge(Int_t bin) const
Return low edge of bin.
Definition: TAxis.cxx:515
Double_t GetXmin() const
Definition: TAxis.h:133
1-Dim function class
Definition: TF1.h:210
A TGraph is an object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:571
1-D histogram with an int per channel (see TH1 documentation)}
Definition: TH1.h:530
The TH1 histogram class.
Definition: TH1.h:56
virtual void SetDirectory(TDirectory *dir)
By default when an histogram is created, it is added to the list of histogram objects in the current ...
Definition: TH1.cxx:8393
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition: TH1.cxx:8597
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8519
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2665
virtual Int_t GetNbinsX() const
Definition: TH1.h:292
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:8678
virtual Double_t GetEntries() const
Return the current number of entries.
Definition: TH1.cxx:4302
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4907
virtual TArrayD * GetSumw2()
Definition: TH1.h:308
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width for 1D histogram.
Definition: TH1.cxx:8619
virtual Double_t GetSumOfWeights() const
Return the sum of weights excluding under/overflows.
Definition: TH1.cxx:7426
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition: TH1.cxx:8476
KDE Kernel for "smoothing" the PDFs.
Definition: KDEKernel.h:50
void SetKernelType(EKernelType ktype=kGauss)
fIter == 1 —> nonadaptive KDE fIter == 2 —> adaptive KDE
Definition: KDEKernel.cxx:112
Float_t GetBinKernelIntegral(Float_t lowr, Float_t highr, Float_t mean, Int_t binnum)
calculates the integral of the Kernel
Definition: KDEKernel.cxx:217
ostringstream derivative to redirect and format output
Definition: MsgLogger.h:59
PDF wrapper for histograms; uses user-defined spline interpolation.
Definition: PDF.h:63
static const Bool_t fgManualIntegration
Definition: PDF.h:162
KDEKernel::EKernelType fKDEtype
Definition: PDF.h:188
void FillHistToGraph()
Simple conversion.
Definition: PDF.cxx:513
void BuildPDF(const TH1 *theHist)
Definition: PDF.cxx:259
KDEKernel::EKernelBorder fKDEborder
Definition: PDF.h:190
static Double_t IGetVal(Double_t *, Double_t *)
static external auxiliary function (integrand)
Definition: PDF.cxx:645
TString fBorderMethodString
Definition: PDF.h:185
Int_t fMaxNsmooth
Definition: PDF.h:169
TMVA::PDF::EInterpolateMethod fInterpolMethod
Definition: PDF.h:172
void SmoothHistogram()
Definition: PDF.cxx:454
Double_t GetValInverse(Double_t y, Bool_t isMonotonouslyIncreasingFunction=kFALSE) const
returns value
Definition: PDF.cxx:733
void CheckHist() const
sanity check: compare PDF with original histogram
Definition: PDF.cxx:552
void FindBinInverse(const TH1 *histogram, Int_t &lowerBin, Int_t &higherBin, Double_t &lowerBinValue, Double_t &higherBinValue, Double_t y, Bool_t isMonotonouslyIncreasingFunction=kFALSE) const
find bin from value on ordinate
Definition: PDF.cxx:760
Int_t fHistDefinedNBins
Definition: PDF.h:181
TString fKDEiterString
Definition: PDF.h:184
void ReadXML(void *pdfnode)
XML file reading.
Definition: PDF.cxx:960
TString fInterpolateString
Definition: PDF.h:186
Int_t fNsmooth
Definition: PDF.h:167
void DeclareOptions()
define the options (their key words) that can be set in the option string
Definition: PDF.cxx:819
void AddXMLTo(void *parent)
XML file writing.
Definition: PDF.cxx:918
void BuildSplinePDF()
build the PDF from the original histograms
Definition: PDF.cxx:320
void FillSplineToHist()
creates high-binned reference histogram to be used instead of the PDF for speed reasons
Definition: PDF.cxx:522
void BuildKDEPDF()
creates high-binned reference histogram to be used instead of the PDF for speed reasons
Definition: PDF.cxx:386
Int_t fMinNsmooth
Definition: PDF.h:168
Int_t GetHistNBins(Int_t evtNum=0)
Definition: PDF.cxx:303
static const Double_t fgEpsilon
Definition: PDF.h:163
TH1 * GetOriginalHist() const
Definition: PDF.h:94
Float_t fFineFactor
Definition: PDF.h:191
static PDF *& GetThisPdfThreadLocal()
Definition: PDF.h:205
virtual ~PDF()
Definition: PDF.cxx:245
void ValidatePDF(TH1 *original=0) const
comparison of original histogram with reference PDF
Definition: PDF.cxx:577
Double_t GetVal(Double_t x) const
returns value PDF(x)
Definition: PDF.cxx:700
Int_t fHistAvgEvtPerBin
Definition: PDF.h:180
Double_t GetIntegral() const
computes normalisation
Definition: PDF.cxx:634
Bool_t fCheckHist
Definition: PDF.h:195
MsgLogger * fLogger
the suffix for options
Definition: PDF.h:199
TString fKDEtypeString
Definition: PDF.h:183
void ProcessOptions()
Definition: PDF.cxx:869
EInterpolateMethod
Definition: PDF.h:70
@ kKDE
Definition: PDF.h:70
@ kSpline1
Definition: PDF.h:70
@ kSpline0
Definition: PDF.h:70
@ kSpline3
Definition: PDF.h:70
@ kSpline2
Definition: PDF.h:70
@ kSpline5
Definition: PDF.h:70
static const Int_t fgNbin_PdfHist
Definition: PDF.h:161
PDF(const TString &name, Bool_t norm=kTRUE)
default constructor needed for ROOT I/O
Definition: PDF.cxx:68
KDEKernel::EKernelIter fKDEiter
Definition: PDF.h:189
static PDF * ThisPDF(void)
Definition: PDF.cxx:1135
Linear interpolation of TGraph.
Definition: TSpline1.h:43
Quadratic interpolation of TGraph.
Definition: TSpline2.h:43
Bool_t HistoHasEquidistantBins(const TH1 &h)
Definition: Tools.cxx:1498
void * GetNextChild(void *prevchild, const char *childname=0)
XML helpers.
Definition: Tools.cxx:1173
void * AddChild(void *parent, const char *childname, const char *content=0, bool isRootNode=false)
add child node
Definition: Tools.cxx:1135
Bool_t AddRawLine(void *node, const char *raw)
XML helpers.
Definition: Tools.cxx:1201
TString StringFromDouble(Double_t d)
string tools
Definition: Tools.cxx:1244
const char * GetContent(void *node)
XML helpers.
Definition: Tools.cxx:1185
void * GetChild(void *parent, const char *childname=0)
get child node
Definition: Tools.cxx:1161
void ReadAttr(void *node, const char *, T &value)
read attribute from xml
Definition: Tools.h:335
void AddAttr(void *node, const char *, const T &value, Int_t precision=16)
add attribute to xml
Definition: Tools.h:353
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Class to create third splines to interpolate knots Arbitrary conditions can be introduced for first a...
Definition: TSpline.h:192
Class to create quintic natural splines to interpolate knots Arbitrary conditions can be introduced f...
Definition: TSpline.h:248
Basic string class.
Definition: TString.h:131
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:687
Element * GetMatrixArray()
Definition: TVectorT.h:78
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Double_t ey[n]
Definition: legend1.C:17
static constexpr double s
std::istream & operator>>(std::istream &istr, BinaryTree &tree)
Tools & gTools()
std::ostream & operator<<(std::ostream &os, const BinaryTree &tree)
MsgLogger & Endl(MsgLogger &ml)
Definition: MsgLogger.h:158
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
Double_t Prob(Double_t chi2, Int_t ndf)
Computation of the probability for a certain Chi-squared (chi2) and number of degrees of freedom (ndf...
Definition: TMath.cxx:614
Double_t Log(Double_t x)
Definition: TMath.h:750
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
Short_t Abs(Short_t d)
Definition: TMathBase.h:120