Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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
35 PDF 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
59 const Int_t TMVA::PDF::fgNbin_PdfHist = 10000;
60 const Bool_t TMVA::PDF::fgManualIntegration = kTRUE;
61 const Double_t TMVA::PDF::fgEpsilon = 1.0e-12;
62 
63 ClassImp(TMVA::PDF);
64 
65 ////////////////////////////////////////////////////////////////////////////////
66 /// default constructor needed for ROOT I/O
67 
68 TMVA::PDF::PDF( const TString& name, Bool_t norm )
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 
106 TMVA::PDF::PDF( const TString& name,
107  const TH1 *hist,
108  PDF::EInterpolateMethod method,
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 
149 TMVA::PDF::PDF( const TString& name,
150  const TH1* hist,
151  KDEKernel::EKernelType ktype,
152  KDEKernel::EKernelIter kiter,
153  KDEKernel::EKernelBorder kborder,
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 
192 TMVA::PDF::PDF( const TString& name,
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;
232  fHistAvgEvtPerBin = defaultPDF->fHistAvgEvtPerBin;
233  fHistAvgEvtPerBin = defaultPDF->fHistAvgEvtPerBin;
234  fInterpolateString = defaultPDF->fInterpolateString;
235  fKDEtypeString = defaultPDF->fKDEtypeString;
236  fKDEiterString = defaultPDF->fKDEiterString;
237  fFineFactor = defaultPDF->fFineFactor;
238  fBorderMethodString = defaultPDF->fBorderMethodString;
239  fCheckHist = defaultPDF->fCheckHist;
240  fHistDefinedNBins = defaultPDF->fHistDefinedNBins;
241  }
242 }
243 
244 //___________________fNSmoothHist____________________________________________________
245 TMVA::PDF::~PDF()
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 
259 void 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 
303 Int_t TMVA::PDF::GetHistNBins ( Int_t evtNum )
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 
320 void TMVA::PDF::BuildSplinePDF()
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 
386 void TMVA::PDF::BuildKDEPDF()
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 
454 void TMVA::PDF::SmoothHistogram()
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 
513 void TMVA::PDF::FillHistToGraph()
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 
522 void TMVA::PDF::FillSplineToHist()
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 
552 void TMVA::PDF::CheckHist() const
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 
577 void 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 
634 Double_t TMVA::PDF::GetIntegral() const
635 {
636  Double_t integral = fPDFHist->GetSumOfWeights();
637  integral *= GetPdfHistBinWidth();
638 
639  return integral;
640 }
641 
642 ////////////////////////////////////////////////////////////////////////////////
643 /// static external auxiliary function (integrand)
644 
645 Double_t TMVA::PDF::IGetVal( Double_t* x, Double_t* )
646 {
647  return ThisPDF()->GetVal( x[0] );
648 }
649 
650 ////////////////////////////////////////////////////////////////////////////////
651 /// computes PDF integral within given ranges
652 
653 Double_t TMVA::PDF::GetIntegral( Double_t xmin, Double_t xmax )
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 
700 Double_t TMVA::PDF::GetVal( Double_t x ) const
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 
733 Double_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 
760 void 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 
819 void TMVA::PDF::DeclareOptions()
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 
869 void TMVA::PDF::ProcessOptions()
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 
918 void 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 
960 void 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;
974  Double_t xmin, xmax;
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 
1035 std::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 
1073 std::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 
1135 TMVA::PDF* TMVA::PDF::ThisPDF( void )
1136 {
1137  // return global "this" pointer of PDF
1138  return GetThisPdfThreadLocal();
1139 }