Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
ROCCalc.cxx
Go to the documentation of this file.
1 /*! \class TMVA::ROCCalc
2 \ingroup TMVA
3 
4 */
5 
6 #include <algorithm>
7 #include <cstdlib>
8 #include <errno.h>
9 
10 #include "TObjString.h"
11 #include "TMath.h"
12 #include "TString.h"
13 #include "TTree.h"
14 #include "TLeaf.h"
15 #include "TH1.h"
16 #include "TH2.h"
17 #include "TList.h"
18 #include "TSpline.h"
19 #include "TVector.h"
20 #include "TMatrixD.h"
21 #include "TMatrixDSymEigen.h"
22 #include "TVectorD.h"
23 #include "TTreeFormula.h"
24 #include "TXMLEngine.h"
25 #include "TROOT.h"
26 #include "TMatrixDSymEigen.h"
27 #include "TColor.h"
28 #include "TMVA/Config.h"
29 
30 
31 #include "TMVA/Tools.h"
32 #include "TMVA/ROCCalc.h"
33 #include "TMVA/Config.h"
34 #include "TMVA/Event.h"
35 #include "TMVA/Version.h"
36 #include "TMVA/PDF.h"
37 #include "TMVA/MsgLogger.h"
38 
39 #include "TMVA/PDF.h"
40 #include "TMVA/TSpline1.h"
41 #include "TMVA/TSpline2.h"
42 #include "TMVA/Types.h"
43 
44 using namespace std;
45 
46 ////////////////////////////////////////////////////////////////////////////////
47 
48 TMVA::ROCCalc::ROCCalc(TH1* mvaS, TH1* mvaB) :
49  fMaxIter(100),
50  fAbsTol(0.0),
51  fStatus(kTRUE),
52  fmvaS(0),
53  fmvaB(0),
54  fmvaSpdf(0),
55  fmvaBpdf(0),
56  fSplS(0),
57  fSplB(0),
58  fSplmvaCumS(0),
59  fSplmvaCumB(0),
60  fSpleffBvsS(0),
61  fnStot(0),
62  fnBtot(0),
63  fSignificance(0),
64  fPurity(0),
65  effBvsS(0),
66  rejBvsS(0),
67  inveffBvsS(0),
68  fLogger ( new TMVA::MsgLogger("ROCCalc") )
69 {
70  fUseSplines = kTRUE;
71  fNbins = 100;
72  // fmvaS = (TH1*) mvaS->Clone("MVA Signal"); fmvaS->SetTitle("MVA Signal");
73  // fmvaB = (TH1*) mvaB->Clone("MVA Backgr"); fmvaB->SetTitle("MVA Backgr");
74  fmvaS = mvaS; fmvaS->SetTitle("MVA Signal");
75  fmvaB = mvaB; fmvaB->SetTitle("MVA Backgr");
76  fXmax = fmvaS->GetXaxis()->GetXmax();
77  fXmin = fmvaS->GetXaxis()->GetXmin();
78 
79  if (TMath::Abs(fXmax-fmvaB->GetXaxis()->GetXmax()) > 0.000001 ||
80  TMath::Abs(fXmin-fmvaB->GetXaxis()->GetXmin()) > 0.000001 ||
81  fmvaB->GetNbinsX() != fmvaS->GetNbinsX()) {
82  Log() << kERROR << "Cannot cal ROC curve etc, as in put mvaS and mvaB have differen #nbins or range "<<Endl;
83  fStatus=kFALSE;
84  }
85  if (!strcmp(fmvaS->GetXaxis()->GetTitle(),"")) fmvaS->SetXTitle("MVA-value");
86  if (!strcmp(fmvaB->GetXaxis()->GetTitle(),"")) fmvaB->SetXTitle("MVA-value");
87  if (!strcmp(fmvaS->GetYaxis()->GetTitle(),"")) fmvaS->SetYTitle("#entries");
88  if (!strcmp(fmvaB->GetYaxis()->GetTitle(),"")) fmvaB->SetYTitle("#entries");
89  ApplySignalAndBackgroundStyle(fmvaS, fmvaB);
90  // std::cout<<"mvaS->GetNbinsX()"<<mvaS->GetNbinsX()<<std::endl;
91  // std::cout<<"mvaB->GetNbinsX()"<<mvaB->GetNbinsX()<<std::endl;
92  //the output of mvaS->GetNbinsX() is about 40 and if we divide it by 100 the results is 0
93  //the I will divide it by 10 anyway doing some tests ROC integral is the same
94  fmvaSpdf = mvaS->RebinX(mvaS->GetNbinsX()/10,"MVA Signal PDF");
95  fmvaBpdf = mvaB->RebinX(mvaB->GetNbinsX()/10,"MVA Backgr PDF");
96  if(fmvaSpdf==0||fmvaBpdf==0)
97  {
98  Log() << kERROR << "Cannot Rebin Histograms mvaS and mvaB, ROC values will be calculated without Rebin histograms."<<Endl;
99  fStatus=kFALSE;
100  fmvaSpdf = (TH1*)mvaS->Clone("MVA Signal PDF");
101  fmvaBpdf = (TH1*)mvaB->Clone("MVA Backgr PDF");
102  }
103  fmvaSpdf->SetTitle("MVA Signal PDF");
104  fmvaBpdf->SetTitle("MVA Backgr PDF");
105  fmvaSpdf->Scale(1./fmvaSpdf->GetSumOfWeights());
106  fmvaBpdf->Scale(1./fmvaBpdf->GetSumOfWeights());
107  fmvaSpdf->SetMaximum(TMath::Max(fmvaSpdf->GetMaximum(), fmvaBpdf->GetMaximum()));
108  fmvaBpdf->SetMaximum(TMath::Max(fmvaSpdf->GetMaximum(), fmvaBpdf->GetMaximum()));
109  ApplySignalAndBackgroundStyle(fmvaSpdf, fmvaBpdf);
110 
111  fCutOrientation = (fmvaS->GetMean() > fmvaB->GetMean()) ? +1 : -1;
112 
113  fNevtS = 0;
114 
115 }
116 
117 ////////////////////////////////////////////////////////////////////////////////
118 
119 void TMVA::ROCCalc::ApplySignalAndBackgroundStyle( TH1* sig, TH1* bkg, TH1* any ) {
120  Int_t c_SignalLine = TColor::GetColor( "#0000ee" );
121  Int_t c_SignalFill = TColor::GetColor( "#7d99d1" );
122  Int_t c_BackgroundLine = TColor::GetColor( "#ff0000" );
123  Int_t c_BackgroundFill = TColor::GetColor( "#ff0000" );
124  // Int_t c_NovelBlue = TColor::GetColor( "#2244a5" );
125 
126  //signal
127  // const Int_t FillColor__S = 38 + 150; // change of Color Scheme in ROOT-5.16.
128  // convince yourself with gROOT->GetListOfColors()->Print()
129  Int_t FillColor__S = c_SignalFill;
130  Int_t FillStyle__S = 1001;
131  Int_t LineColor__S = c_SignalLine;
132  Int_t LineWidth__S = 2;
133 
134  // background
135  //Int_t icolor = gConfig().fVariablePlotting.fUsePaperStyle ? 2 + 100 : 2;
136  Int_t FillColor__B = c_BackgroundFill;
137  Int_t FillStyle__B = 3554;
138  Int_t LineColor__B = c_BackgroundLine;
139  Int_t LineWidth__B = 2;
140 
141  if (sig != NULL) {
142  sig->SetLineColor( LineColor__S );
143  sig->SetLineWidth( LineWidth__S );
144  sig->SetFillStyle( FillStyle__S );
145  sig->SetFillColor( FillColor__S );
146  }
147 
148  if (bkg != NULL) {
149  bkg->SetLineColor( LineColor__B );
150  bkg->SetLineWidth( LineWidth__B );
151  bkg->SetFillStyle( FillStyle__B );
152  bkg->SetFillColor( FillColor__B );
153  }
154 
155  if (any != NULL) {
156  any->SetLineColor( LineColor__S );
157  any->SetLineWidth( LineWidth__S );
158  any->SetFillStyle( FillStyle__S );
159  any->SetFillColor( FillColor__S );
160  }
161 }
162 
163 ////////////////////////////////////////////////////////////////////////////////
164 /// destructor
165 
166 TMVA::ROCCalc::~ROCCalc() {
167  // delete Splines and all histograms that were created only for internal use
168  if (fSplS) { delete fSplS; fSplS = 0; }
169  if (fSplB) { delete fSplB; fSplB = 0; }
170  if (fSpleffBvsS) { delete fSpleffBvsS; fSpleffBvsS = 0; }
171  if (fSplmvaCumS) { delete fSplmvaCumS; fSplmvaCumS = 0; }
172  if (fSplmvaCumB) { delete fSplmvaCumB; fSplmvaCumB = 0; }
173  if (fmvaScumul) { delete fmvaScumul; }
174  if (fmvaBcumul) { delete fmvaBcumul; }
175  if (effBvsS) { delete effBvsS; }
176  if (rejBvsS) { delete rejBvsS; }
177  if (inveffBvsS) { delete inveffBvsS; }
178  delete fLogger;
179 }
180 
181 ////////////////////////////////////////////////////////////////////////////////
182 /// get the ROC curve
183 
184 TH1D* TMVA::ROCCalc::GetROC(){
185  // first get the cumulative distributions of the mva distribution
186  // --> efficiencies vs cut value
187  fNevtS = fmvaS->GetSumOfWeights(); // needed to get the error on the eff.. will only be correct if the histogram is not scaled to "integral == 1" Yet;
188  if (fNevtS < 2) {
189  Log() << kERROR << "I guess the mva distributions fed into ROCCalc were already normalized, therefore the calculated error on the efficiency will be incorrect !! " << Endl;
190  fNevtS = 0; // reset to zero --> no error will be calculated on the efficiencies
191  fStatus=kFALSE;
192  }
193  fmvaScumul = gTools().GetCumulativeDist(fmvaS);
194  fmvaBcumul = gTools().GetCumulativeDist(fmvaB);
195  fmvaScumul->Scale( 1.0/TMath::Max(std::numeric_limits<double>::epsilon(),fmvaScumul->GetMaximum()) );
196  fmvaBcumul->Scale( 1.0/TMath::Max(std::numeric_limits<double>::epsilon(),fmvaBcumul->GetMaximum()) );
197  fmvaScumul->SetMinimum(0);
198  fmvaBcumul->SetMinimum(0);
199  // fmvaScumul->Draw("hist");
200  // fmvaBcumul->Draw("histsame");
201 
202  // background efficiency versus signal efficiency
203  if(effBvsS==0) effBvsS = new TH1D("effBvsS", "ROC-Curve", fNbins, 0, 1 );
204  effBvsS->SetXTitle( "Signal eff" );
205  effBvsS->SetYTitle( "Backgr eff" );
206 
207  // background rejection (=1-eff.) versus signal efficiency
208  if(rejBvsS==0) rejBvsS = new TH1D( "rejBvsS", "ROC-Curve", fNbins, 0, 1 );
209  rejBvsS->SetXTitle( "Signal eff" );
210  rejBvsS->SetYTitle( "Backgr rejection (1-eff)" );
211 
212  // inverse background eff (1/eff.) versus signal efficiency
213  if(inveffBvsS ==0) inveffBvsS = new TH1D("invBeffvsSeff", "ROC-Curve" , fNbins, 0, 1 );
214  inveffBvsS->SetXTitle( "Signal eff" );
215  inveffBvsS->SetYTitle( "Inverse backgr. eff (1/eff)" );
216 
217  // use root finder
218  // spline background efficiency plot
219  // note that there is a bin shift when going from a TH1D object to a TGraph :-(
220  if (fUseSplines) {
221  fSplmvaCumS = new TSpline1( "spline2_signal", new TGraph( fmvaScumul ) );
222  fSplmvaCumB = new TSpline1( "spline2_background", new TGraph( fmvaBcumul ) );
223  // verify spline sanity
224  gTools().CheckSplines( fmvaScumul, fSplmvaCumS );
225  gTools().CheckSplines( fmvaBcumul, fSplmvaCumB );
226  }
227 
228  Double_t effB = 0;
229  for (UInt_t bini=1; bini<=fNbins; bini++) {
230 
231  // find cut value corresponding to a given signal efficiency
232  Double_t effS = effBvsS->GetBinCenter( bini );
233  Double_t cut = Root( effS );
234 
235  // retrieve background efficiency for given cut
236  if (fUseSplines) effB = fSplmvaCumB->Eval( cut );
237  else effB = fmvaBcumul->GetBinContent( fmvaBcumul->FindBin( cut ) );
238 
239  // and fill histograms
240  effBvsS->SetBinContent( bini, effB );
241  rejBvsS->SetBinContent( bini, 1.0-effB );
242  if (effB>std::numeric_limits<double>::epsilon())
243  inveffBvsS->SetBinContent( bini, 1.0/effB );
244  }
245 
246  // create splines for histogram
247  fSpleffBvsS = new TSpline1( "effBvsS", new TGraph( effBvsS ) );
248 
249  // search for overlap point where, when cutting on it,
250  // one would obtain: eff_S = rej_B = 1 - eff_B
251 
252  Double_t effS = 0., rejB = 0., effS_ = 0., rejB_ = 0.;
253  Int_t nbins = 5000;
254  for (Int_t bini=1; bini<=nbins; bini++) {
255 
256  // get corresponding signal and background efficiencies
257  effS = (bini - 0.5)/Float_t(nbins);
258  rejB = 1.0 - fSpleffBvsS->Eval( effS );
259 
260  // find signal efficiency that corresponds to required background efficiency
261  if ((effS - rejB)*(effS_ - rejB_) < 0) break;
262  effS_ = effS;
263  rejB_ = rejB;
264  }
265  // find cut that corresponds to signal efficiency and update signal-like criterion
266  fSignalCut = Root( 0.5*(effS + effS_) );
267 
268  return rejBvsS;
269 }
270 
271 ////////////////////////////////////////////////////////////////////////////////
272 /// code to compute the area under the ROC ( rej-vs-eff ) curve
273 
274 Double_t TMVA::ROCCalc::GetROCIntegral(){
275  Double_t effS = 0, effB = 0;
276  Int_t nbins = 1000;
277  if (fSpleffBvsS == 0) this->GetROC(); // that will make the ROC calculation if not done yet
278 
279  // compute area of rej-vs-eff plot
280  Double_t integral = 0;
281  for (Int_t bini=1; bini<=nbins; bini++) {
282 
283  // get corresponding signal and background efficiencies
284  effS = (bini - 0.5)/Float_t(nbins);
285  effB = fSpleffBvsS->Eval( effS );
286  integral += (1.0 - effB);
287  }
288  integral /= nbins;
289 
290  return integral;
291 }
292 
293 ////////////////////////////////////////////////////////////////////////////////
294 /// get the signal efficiency for a particular background efficiency
295 /// that will be the value of the efficiency retured (does not affect
296 /// the efficiency-vs-bkg plot which is done anyway.
297 
298 Double_t TMVA::ROCCalc::GetEffSForEffBof(Double_t effBref, Double_t &effSerr){
299  // find precise efficiency value
300  Double_t effS=0., effB, effSOld=1., effBOld=0.;
301  Int_t nbins = 1000;
302  if (fSpleffBvsS == 0) this->GetROC(); // that will make the ROC calculation if not done yet
303 
304  Float_t step=1./nbins; // stepsize in efficiency binning
305  for (Int_t bini=1; bini<=nbins; bini++) {
306  // get corresponding signal and background efficiencies
307  effS = (bini - 0.5)*step; // efficiency goes from 0-to-1 in nbins steps of 1/nbins (take middle of the bin)
308  effB = fSpleffBvsS->Eval( effS );
309 
310  // find signal efficiency that corresponds to required background efficiency
311  if ((effB - effBref)*(effBOld - effBref) <= 0) break;
312  effSOld = effS;
313  effBOld = effB;
314  }
315 
316  // take mean between bin above and bin below
317  effS = 0.5*(effS + effSOld);
318 
319 
320  if (fNevtS > 0) effSerr = TMath::Sqrt( effS*(1.0 - effS)/fNevtS );
321  else effSerr = 0;
322 
323  return effS;
324 }
325 
326 ////////////////////////////////////////////////////////////////////////////////
327 /// returns efficiency as function of cut
328 
329 Double_t TMVA::ROCCalc::GetEffForRoot( Double_t theCut )
330 {
331  Double_t retVal=0;
332 
333  // retrieve the class object
334  if (fUseSplines) retVal = fSplmvaCumS->Eval( theCut );
335  else retVal = fmvaScumul->GetBinContent( fmvaScumul->FindBin( theCut ) );
336 
337  // caution: here we take some "forbidden" action to hide a problem:
338  // in some cases, in particular for likelihood, the binned efficiency distributions
339  // do not equal 1, at xmin, and 0 at xmax; of course, in principle we have the
340  // unbinned information available in the trees, but the unbinned minimization is
341  // too slow, and we don't need to do a precision measurement here. Hence, we force
342  // this property.
343  Double_t eps = 1.0e-5;
344  if (theCut-fXmin < eps) retVal = (fCutOrientation > 0) ? 1.0 : 0.0;
345  else if (fXmax-theCut < eps) retVal = (fCutOrientation > 0) ? 0.0 : 1.0;
346 
347 
348  return retVal;
349 }
350 
351 ////////////////////////////////////////////////////////////////////////////////
352 /// Root finding using Brents algorithm; taken from CERNLIB function RZERO
353 
354 Double_t TMVA::ROCCalc::Root( Double_t refValue )
355 {
356  Double_t a = fXmin, b = fXmax;
357  Double_t fa = GetEffForRoot( a ) - refValue;
358  Double_t fb = GetEffForRoot( b ) - refValue;
359  if (fb*fa > 0) {
360  Log() << kWARNING << "<ROCCalc::Root> initial interval w/o root: "
361  << "(a=" << a << ", b=" << b << "),"
362  << " (Eff_a=" << GetEffForRoot( a )
363  << ", Eff_b=" << GetEffForRoot( b ) << "), "
364  << "(fa=" << fa << ", fb=" << fb << "), "
365  << "refValue = " << refValue << Endl;
366  return 1;
367  }
368 
369  Bool_t ac_equal(kFALSE);
370  Double_t fc = fb;
371  Double_t c = 0, d = 0, e = 0;
372  for (Int_t iter= 0; iter <= fMaxIter; iter++) {
373  if ((fb < 0 && fc < 0) || (fb > 0 && fc > 0)) {
374 
375  // Rename a,b,c and adjust bounding interval d
376  ac_equal = kTRUE;
377  c = a; fc = fa;
378  d = b - a; e = b - a;
379  }
380 
381  if (TMath::Abs(fc) < TMath::Abs(fb)) {
382  ac_equal = kTRUE;
383  a = b; b = c; c = a;
384  fa = fb; fb = fc; fc = fa;
385  }
386 
387  Double_t tol = 0.5 * 2.2204460492503131e-16 * TMath::Abs(b);
388  Double_t m = 0.5 * (c - b);
389  if (fb == 0 || TMath::Abs(m) <= tol || TMath::Abs(fb) < fAbsTol) return b;
390 
391  // Bounds decreasing too slowly: use bisection
392  if (TMath::Abs (e) < tol || TMath::Abs (fa) <= TMath::Abs (fb)) { d = m; e = m; }
393  else {
394  // Attempt inverse cubic interpolation
395  Double_t p, q, r;
396  Double_t s = fb / fa;
397 
398  if (ac_equal) { p = 2 * m * s; q = 1 - s; }
399  else {
400  q = fa / fc; r = fb / fc;
401  p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
402  q = (q - 1) * (r - 1) * (s - 1);
403  }
404  // Check whether we are in bounds
405  if (p > 0) q = -q;
406  else p = -p;
407 
408  Double_t min1 = 3 * m * q - TMath::Abs (tol * q);
409  Double_t min2 = TMath::Abs (e * q);
410  if (2 * p < (min1 < min2 ? min1 : min2)) {
411  // Accept the interpolation
412  e = d; d = p / q;
413  }
414  else { d = m; e = m; } // Interpolation failed: use bisection.
415  }
416  // Move last best guess to a
417  a = b; fa = fb;
418  // Evaluate new trial root
419  if (TMath::Abs(d) > tol) b += d;
420  else b += (m > 0 ? +tol : -tol);
421 
422  fb = GetEffForRoot( b ) - refValue;
423 
424  }
425 
426  // Return our best guess if we run out of iterations
427  Log() << kWARNING << "<ROCCalc::Root> maximum iterations (" << fMaxIter
428  << ") reached before convergence" << Endl;
429 
430  return b;
431 }
432 
433 ////////////////////////////////////////////////////////////////////////////////
434 
435 TH1* TMVA::ROCCalc::GetPurity( Int_t nStot, Int_t nBtot)
436 {
437  if (fnStot!=nStot || fnBtot!=nBtot || !fSignificance) {
438  GetSignificance(nStot, nBtot);
439  fnStot=nStot;
440  fnBtot=nBtot;
441  }
442  return fPurity;
443 }
444 ////////////////////////////////////////////////////////////////////////////////
445 
446 TH1* TMVA::ROCCalc::GetSignificance( Int_t nStot, Int_t nBtot)
447 {
448  if (fnStot==nStot && fnBtot==nBtot && !fSignificance) return fSignificance;
449  fnStot=nStot; fnBtot=nBtot;
450 
451  fSignificance = (TH1*) fmvaScumul->Clone("Significance"); fSignificance->SetTitle("Significance");
452  fSignificance->Reset(); fSignificance->SetFillStyle(0);
453  fSignificance->SetXTitle("mva cut value");
454  fSignificance->SetYTitle("Stat. significance S/Sqrt(S+B)");
455  fSignificance->SetLineColor(2);
456  fSignificance->SetLineWidth(5);
457 
458  fPurity = (TH1*) fmvaScumul->Clone("Purity"); fPurity->SetTitle("Purity");
459  fPurity->Reset(); fPurity->SetFillStyle(0);
460  fPurity->SetXTitle("mva cut value");
461  fPurity->SetYTitle("Purity: S/(S+B)");
462  fPurity->SetLineColor(3);
463  fPurity->SetLineWidth(5);
464 
465  Double_t maxSig=0;
466  for (Int_t i=1; i<=fSignificance->GetNbinsX(); i++) {
467  Double_t S = fmvaScumul->GetBinContent( i ) * nStot;
468  Double_t B = fmvaBcumul->GetBinContent( i ) * nBtot;
469  Double_t purity;
470  Double_t sig;
471  if (S+B > 0){
472  purity = S/(S+B);
473  sig = S/TMath::Sqrt(S+B);
474  if (sig > maxSig) {
475  maxSig = sig;
476  }
477  } else {
478  purity=0;
479  sig=0;
480  }
481  cout << "S="<<S<<" B="<<B<< " purity="<<purity<< endl;
482  fPurity->SetBinContent( i, purity );
483  fSignificance->SetBinContent( i, sig );
484  }
485 
486  /*
487  TLatex* line1;
488  TLatex* line2;
489  TLatex tl;
490  tl.SetNDC();
491  tl.SetTextSize( 0.033 );
492  Int_t maxbin = fSignificance->GetMaximumBin();
493  line1 = tl.DrawLatex( 0.15, 0.23, Form("For %1.0f signal and %1.0f background", nStot, nBtot));
494  tl.DrawLatex( 0.15, 0.19, "events the maximum S/Sqrt(S+B) is");
495 
496  line2 = tl.DrawLatex( 0.15, 0.15, Form("%4.2f when cutting at %5.2f",
497  maxSig,
498  fSignificance->GetXaxis()->GetBinCenter(maxbin)) );
499  */
500  return fSignificance;
501 
502 }