MathTools/src/MathUtils.cc

// $Id: MathUtils.cc,v 1.10 2009/11/03 10:28:08 sixie Exp $

#include "MitCommon/MathTools/interface/MathUtils.h"
#include <TError.h>
#include <TH1D.h>
#include <TGraphAsymmErrors.h>
#include "PhysicsTools/RooStatsCms/interface/ClopperPearsonBinomialInterval.h"
#include "PhysicsTools/RooStatsCms/interface/FeldmanCousinsBinomialInterval.h"


ClassImp(mithep::MathUtils)

using namespace mithep;

//--------------------------------------------------------------------------------------------------
Double_t MathUtils::AddInQuadrature(Double_t a, Double_t b)
{
  // Add quantities in quadrature.

  return(TMath::Sqrt(a*a + b*b));
}

//--------------------------------------------------------------------------------------------------
void MathUtils::CalcRatio(Double_t n1, Double_t n2, Double_t &r, Double_t &rlow, Double_t &rup, Int_t type = 0)
{
  // Calculate ratio and lower/upper errors from given values using various methods, dependent on type:
  //  0: Bayes
  //  1: Feldman-Cousins
  //  2: Clopper-Pearson
  
  if (n1>n2) {
    Error("CalcRatio", "First value should be smaller than second: %f > %f", n1, n2);
    r = n1/n2;
    rlow = 0;
    rup  = 0;
    return;
  }
 
  if (n2 >= 1) {
    if (type == 1) {
      //compute errors using Feldman-Cousins
      r = double(n1) / double(n2);
      FeldmanCousinsBinomialInterval fc;
      const double alpha = (1-0.682);
      fc.init(alpha);
      fc.calculate(n1, n2);
      rlow = r - fc.lower();
      rup = fc.upper() - r;    
      
    } else if (type == 2) {
      //compute errors using Clopper-Pearson
      r = double(n1) / double(n2);
      ClopperPearsonBinomialInterval cp;
      const double alpha = (1-0.682);
      cp.init(alpha);
      cp.calculate(n1, n2);
      rlow = r - cp.lower();
      rup = cp.upper() - r;
      
    } else {
      //compute using Bayes
      TH1D h1("dummy1","",1,1,2);
      h1.SetBinContent(1,n1);
      
      TH1D h2("dummy2","",1,1,2);
      h2.SetBinContent(1,n2);
      
      TGraphAsymmErrors g;
      g.BayesDivide(&h1,&h2);
      r = g.GetY()[0];
      rup = g.GetErrorYhigh(0);
      rlow = g.GetErrorYlow(0);
    } 
  } else {
    r = 0;
    rlow = 0;
    rup  = 0;
  }
}       


//--------------------------------------------------------------------------------------------------
Double_t MathUtils::DeltaPhi(Double_t phi1, Double_t phi2)
{
  // Compute DeltaPhi between two given angles. Results is in [-pi/2,pi/2].

  Double_t dphi = TMath::Abs(phi1-phi2);
  while (dphi>TMath::Pi())
    dphi = TMath::Abs(dphi - TMath::TwoPi());
  return(dphi);
}

//--------------------------------------------------------------------------------------------------
Double_t MathUtils::DeltaR(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
{
  // Compute DeltaR between two given points in the eta/phi plane.

  Double_t dR = TMath::Sqrt(DeltaR2(phi1,eta1,phi2,eta2));
  return(dR);
}

//--------------------------------------------------------------------------------------------------
Double_t MathUtils::DeltaR2(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
{
  // Compute DeltaR between two given points in the eta/phi plane.

  Double_t dphi = DeltaPhi(phi1, phi2);
  Double_t deta = eta1-eta2;
  Double_t dR = dphi*dphi + deta*deta;
  return(dR);
}

//--------------------------------------------------------------------------------------------------
Double_t MathUtils::Eta2Theta(Double_t eta) 
{ 
  // Compute theta from given eta value.

  return 2.*TMath::ATan(exp(-eta)); 
}

//--------------------------------------------------------------------------------------------------
Double_t MathUtils::Theta2Eta(Double_t theta) 
{ 
  // Compute eta from given theta value.

  return -TMath::Log(TMath::Tan(theta/2.)); 
}
Revision:	1.11
Committed:	Wed Nov 4 12:44:33 2009 UTC (15 years, 6 months ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_025c_branch2, Mit_025c_branch1, Mit_025c_branch0, Mit_025c, Mit_025b, Mit_025a, Mit_025, Mit_025pre2, Mit_024b, Mit_025pre1, Mit_024a, Mit_024, Mit_023, Mit_022a, Mit_022, Mit_020d, TMit_020d, Mit_020c, Mit_021, Mit_021pre2, Mit_021pre1, Mit_020b, Mit_020a, Mit_020, Mit_020pre1, Mit_018, Mit_017, Mit_017pre3, Mit_017pre2, Mit_017pre1, V07-05-00, Mit_016, Mit_015b, Mit_015a, Mit_015, Mit_014e, Mit_014d, Mit_014c, Mit_014b, ConvRejection-10-06-09, Mit_014a, Mit_014, Mit_014pre3, Mit_014pre2, Mit_014pre1, Mit_013d, Mit_013c, Mit_013b, Mit_013a, Mit_013, Mit_013pre1, Mit_012i, Mit_012g, Mit_012f, Mit_012e, Mit_012d, Mit_012c, Mit_012b, Mit_012a, Mit_012
Branch point for:	Mit_025c_branch
Changes since 1.10:	+8 -6 lines
Log Message:	Improved docu
#	Content
1	// $Id: MathUtils.cc,v 1.10 2009/11/03 10:28:08 sixie Exp $
2
3	#include "MitCommon/MathTools/interface/MathUtils.h"
4	#include <TError.h>
5	#include <TH1D.h>
6	#include <TGraphAsymmErrors.h>
7	#include "PhysicsTools/RooStatsCms/interface/ClopperPearsonBinomialInterval.h"
8	#include "PhysicsTools/RooStatsCms/interface/FeldmanCousinsBinomialInterval.h"
9
10
11	ClassImp(mithep::MathUtils)
12
13	using namespace mithep;
14
15	//--------------------------------------------------------------------------------------------------
16	Double_t MathUtils::AddInQuadrature(Double_t a, Double_t b)
17	{
18	// Add quantities in quadrature.
19
20	return(TMath::Sqrt(aa + bb));
21	}
22
23	//--------------------------------------------------------------------------------------------------
24	void MathUtils::CalcRatio(Double_t n1, Double_t n2, Double_t &r, Double_t &rlow, Double_t &rup, Int_t type = 0)
25	{
26	// Calculate ratio and lower/upper errors from given values using various methods, dependent on type:
27	// 0: Bayes
28	// 1: Feldman-Cousins
29	// 2: Clopper-Pearson
30
31	if (n1>n2) {
32	Error("CalcRatio", "First value should be smaller than second: %f > %f", n1, n2);
33	r = n1/n2;
34	rlow = 0;
35	rup = 0;
36	return;
37	}
38
39	if (n2 >= 1) {
40	if (type == 1) {
41	//compute errors using Feldman-Cousins
42	r = double(n1) / double(n2);
43	FeldmanCousinsBinomialInterval fc;
44	const double alpha = (1-0.682);
45	fc.init(alpha);
46	fc.calculate(n1, n2);
47	rlow = r - fc.lower();
48	rup = fc.upper() - r;
49
50	} else if (type == 2) {
51	//compute errors using Clopper-Pearson
52	r = double(n1) / double(n2);
53	ClopperPearsonBinomialInterval cp;
54	const double alpha = (1-0.682);
55	cp.init(alpha);
56	cp.calculate(n1, n2);
57	rlow = r - cp.lower();
58	rup = cp.upper() - r;
59
60	} else {
61	//compute using Bayes
62	TH1D h1("dummy1","",1,1,2);
63	h1.SetBinContent(1,n1);
64
65	TH1D h2("dummy2","",1,1,2);
66	h2.SetBinContent(1,n2);
67
68	TGraphAsymmErrors g;
69	g.BayesDivide(&h1,&h2);
70	r = g.GetY()[0];
71	rup = g.GetErrorYhigh(0);
72	rlow = g.GetErrorYlow(0);
73	}
74	} else {
75	r = 0;
76	rlow = 0;
77	rup = 0;
78	}
79	}
80
81
82	//--------------------------------------------------------------------------------------------------
83	Double_t MathUtils::DeltaPhi(Double_t phi1, Double_t phi2)
84	{
85	// Compute DeltaPhi between two given angles. Results is in [-pi/2,pi/2].
86
87	Double_t dphi = TMath::Abs(phi1-phi2);
88	while (dphi>TMath::Pi())
89	dphi = TMath::Abs(dphi - TMath::TwoPi());
90	return(dphi);
91	}
92
93	//--------------------------------------------------------------------------------------------------
94	Double_t MathUtils::DeltaR(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
95	{
96	// Compute DeltaR between two given points in the eta/phi plane.
97
98	Double_t dR = TMath::Sqrt(DeltaR2(phi1,eta1,phi2,eta2));
99	return(dR);
100	}
101
102	//--------------------------------------------------------------------------------------------------
103	Double_t MathUtils::DeltaR2(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
104	{
105	// Compute DeltaR between two given points in the eta/phi plane.
106
107	Double_t dphi = DeltaPhi(phi1, phi2);
108	Double_t deta = eta1-eta2;
109	Double_t dR = dphidphi + detadeta;
110	return(dR);
111	}
112
113	//--------------------------------------------------------------------------------------------------
114	Double_t MathUtils::Eta2Theta(Double_t eta)
115	{
116	// Compute theta from given eta value.
117
118	return 2.*TMath::ATan(exp(-eta));
119	}
120
121	//--------------------------------------------------------------------------------------------------
122	Double_t MathUtils::Theta2Eta(Double_t theta)
123	{
124	// Compute eta from given theta value.
125
126	return -TMath::Log(TMath::Tan(theta/2.));
127	}