Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooEllipse.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitCore *
4 * @(#)root/roofitcore:$Id$
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2000-2005, Regents of the University of California *
10 * and Stanford University. All rights reserved. *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17/**
18\file RooEllipse.cxx
19\class RooEllipse
20\ingroup Roofitcore
21
22Two-dimensional ellipse that can be used to represent an error contour.
23**/
24
25#include "RooEllipse.h"
26#include "TMath.h"
27#include "RooMsgService.h"
28
29#include "Riostream.h"
30#include "TClass.h"
31#include <cmath>
32
33
34////////////////////////////////////////////////////////////////////////////////
35/// Create a 2-dimensional ellipse centered at `(x1,x2)` that represents the confidence
36/// level contour for a measurement with standard deviations`(s1,s2)` and correlation coefficient rho,
37/// and semiaxes corresponding to `k` standard deviations.
38/// It is assumed that the data are distributed according to a bivariate normal distribution:
39/// \f[ p(x,y) = {1 \over 2 \pi s_1 s_2 \sqrt{1-\rho^2}} \exp (-((x-x_1)^2/s_1^2 + (y-x_2)^2/s_2^2 - 2 \rho x y/(s_1s_2))/2(1-\rho^2)) \f]
40/// \see ROOT::Math::bigaussian_pdf
41/// Or in a split form:
42/// \f[ p(z(x,y)) = {1 \over 2 \pi s_1 s_2 \sqrt{1-\rho^2}} \exp (-z/2(1-\rho^2)) \f]
43/// with:
44/// \f[ z = ((x-x_1)^2/s_1^2 + (y-x_2)^2/s_2^2 - 2 \rho x y/(s_1s_2)) \f]
45/// As demonstrated in https://root-forum.cern.ch/t/drawing-convergence-ellipse-in-2d/61936/10, the
46/// "oriented" ellipse with semi-axis = (k * s_1, k * s_2) includes 39% (for k = 1) of
47/// all data, this confidence ellipse can thus be interpreted as a confidence contour level.
48/// This ellipse can be described by the implicit equation `z/(1-rho*rho) = k^2 = - 2 * nll ratio`, or explictly:
49///
50/// x*x 2*rho*x*y y*y
51/// ----- - --------- + ----- = k*k * (1 - rho*rho)
52/// s1*s1 s1*s2 s2*s2
53///
54/// The input parameters s1,s2,k must be > 0 and also |rho| <= 1.
55/// The degenerate case |rho|=1 corresponds to a straight line and
56/// is handled as a special case.
57/// \warning The default value of k = 1 (semiaxes at +/-sigma) corresponds to the 39% CL contour. For a 1-sigma confidence
58/// level (68%), set instead k ~ 1.51 (semiaxes at +/-1.51sigma).
59
60RooEllipse::RooEllipse(const char *name, double x1, double x2, double s1, double s2, double rho, Int_t points, double k)
61{
64
65 if(s1 <= 0 || s2 <= 0 || k <= 0) {
66 coutE(InputArguments) << "RooEllipse::RooEllipse: bad parameter s1, s2 or k <= 0" << std::endl;
67 return;
68 }
69 double tmp = k*k*(1-rho*rho);
70 if(tmp < 0) {
71 coutE(InputArguments) << "RooEllipse::RooEllipse: bad parameter |rho| > 1" << std::endl;
72 return;
73 }
74
75 if(tmp == 0) {
76 // handle the degenerate case of |rho|=1
77 SetPoint(0,x1-s1,x2-s2);
78 SetPoint(1,x1+s1,x2+s2);
80 }
81 else {
82 double r;
83 double psi;
84 double phi;
85 double u1;
86 double u2;
87 double xx1;
88 double xx2;
89 double dphi(2 * TMath::Pi() / points);
90 for(Int_t index= 0; index < points; index++) {
91 phi= index*dphi;
92 // adjust the angular spacing of the points for the aspect ratio
93 psi= atan2(s2*sin(phi),s1*cos(phi));
94 u1= cos(psi)/s1;
95 u2= sin(psi)/s2;
96 r= sqrt(tmp/(u1*u1 - 2*rho*u1*u2 + u2*u2));
97 xx1= x1 + r*u1*s1;
98 xx2= x2 + r*u2*s2;
100 if(index == 0) {
102 // add an extra segment to close the curve
104 }
105 else {
107 }
108 }
109 }
110}
111
112
113
114////////////////////////////////////////////////////////////////////////////////
115/// Print name of ellipse on ostream
116
117void RooEllipse::printName(std::ostream& os) const
118{
119 os << GetName() ;
120}
121
122
123////////////////////////////////////////////////////////////////////////////////
124/// Print title of ellipse on ostream
125
126void RooEllipse::printTitle(std::ostream& os) const
127{
128 os << GetName() ;
129}
130
131
132////////////////////////////////////////////////////////////////////////////////
133/// Print class name of ellipse on ostream
134
135void RooEllipse::printClassName(std::ostream& os) const
136{
137 os << ClassName() ;
138}
139
140
141////////////////////////////////////////////////////////////////////////////////
142/// Print detailed multi line information on ellipse on ostreamx
143
144void RooEllipse::printMultiline(std::ostream& os, Int_t contents, bool verbose, TString indent) const
145{
146 RooPlotable::printMultiline(os,contents,verbose,indent);
147 for(Int_t index=0; index < fNpoints; index++) {
148 os << indent << "Point [" << index << "] is at (" << fX[index] << "," << fY[index] << ")" << std::endl;
149 }
150}
#define s1(x)
Definition RSha256.hxx:91
#define coutE(a)
static void indent(ostringstream &buf, int indent_level)
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Option_t Option_t TPoint TPoint const char x2
Option_t Option_t TPoint TPoint const char x1
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t points
char name[80]
Definition TGX11.cxx:110
void printTitle(std::ostream &os) const override
Print title of ellipse on ostream.
RooEllipse()=default
void printClassName(std::ostream &os) const override
Print class name of ellipse on ostream.
void printName(std::ostream &os) const override
Print name of ellipse on ostream.
void printMultiline(std::ostream &os, Int_t contents, bool verbose=false, TString indent="") const override
Print detailed multi line information on ellipse on ostreamx.
void printMultiline(std::ostream &os, Int_t contents, bool verbose=false, TString indent="") const override
Print detailed information.
void updateYAxisLimits(double y)
Definition RooPlotable.h:30
void setYAxisLimits(double ymin, double ymax)
Definition RooPlotable.h:34
Int_t fNpoints
Number of points <= fMaxSize.
Definition TGraph.h:46
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition TGraph.cxx:2291
Double_t * fY
[fNpoints] array of Y points
Definition TGraph.h:48
void SetName(const char *name="") override
Set graph name.
Definition TGraph.cxx:2330
Double_t * fX
[fNpoints] array of X points
Definition TGraph.h:47
void SetTitle(const char *title="") override
Change (i.e.
Definition TGraph.cxx:2346
const char * GetName() const override
Returns name of object.
Definition TNamed.h:49
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition TObject.cxx:226
Basic string class.
Definition TString.h:139
constexpr Double_t Pi()
Definition TMath.h:38