MathTools/src/MathUtils.cc

// $Id: MathUtils.cc,v 1.11 2009/11/04 12:44:33 loizides 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 = phi1-phi2;
  while (dphi > TMath::Pi())
    dphi -= TMath::TwoPi();
  while (dphi <= -TMath::Pi())
    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.12
Committed:	Sat Dec 10 23:59:51 2011 UTC (13 years, 4 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_032, Mit_031, Mit_030, Mit_029c, Mit_030_pre1, Mit_029a, Mit_029, Mit_029_pre1, Mit_028a, Mit_028, Mit_027a, Mit_027, Mit_026, Mit_025e, Mit_025d, HEAD
Changes since 1.11:	+10 -4 lines
Log Message:	deltaphi on proper interval
#	User	Rev	Content
1	bendavid	1.12	// $Id: MathUtils.cc,v 1.11 2009/11/04 12:44:33 loizides Exp $
2	sixie	1.1
3			#include "MitCommon/MathTools/interface/MathUtils.h"
4	loizides	1.8	#include <TError.h>
5			#include <TH1D.h>
6			#include <TGraphAsymmErrors.h>
7	sixie	1.10	#include "PhysicsTools/RooStatsCms/interface/ClopperPearsonBinomialInterval.h"
8			#include "PhysicsTools/RooStatsCms/interface/FeldmanCousinsBinomialInterval.h"
9
10	sixie	1.1
11	loizides	1.9	ClassImp(mithep::MathUtils)
12
13	loizides	1.3	using namespace mithep;
14	sixie	1.1
15	loizides	1.2	//--------------------------------------------------------------------------------------------------
16	loizides	1.4	Double_t MathUtils::AddInQuadrature(Double_t a, Double_t b)
17	loizides	1.2	{
18	loizides	1.5	// Add quantities in quadrature.
19
20	loizides	1.2	return(TMath::Sqrt(aa + bb));
21			}
22	sixie	1.1
23	loizides	1.2	//--------------------------------------------------------------------------------------------------
24	sixie	1.10	void MathUtils::CalcRatio(Double_t n1, Double_t n2, Double_t &r, Double_t &rlow, Double_t &rup, Int_t type = 0)
25	loizides	1.8	{
26	loizides	1.11	// 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	loizides	1.8
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	sixie	1.10
39			if (n2 >= 1) {
40			if (type == 1) {
41	loizides	1.11	//compute errors using Feldman-Cousins
42	sixie	1.10	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	loizides	1.11	//compute errors using Clopper-Pearson
52	sixie	1.10	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	loizides	1.11	//compute using Bayes
62	sixie	1.10	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	loizides	1.8
81
82			//--------------------------------------------------------------------------------------------------
83	loizides	1.4	Double_t MathUtils::DeltaPhi(Double_t phi1, Double_t phi2)
84	sixie	1.1	{
85	loizides	1.5	// Compute DeltaPhi between two given angles. Results is in [-pi/2,pi/2].
86
87	bendavid	1.12	Double_t dphi = phi1-phi2;
88			while (dphi > TMath::Pi())
89			dphi -= TMath::TwoPi();
90			while (dphi <= -TMath::Pi())
91			dphi += TMath::TwoPi();
92
93
94	loizides	1.4	return(dphi);
95	bendavid	1.12
96
97	sixie	1.1	}
98
99	loizides	1.2	//--------------------------------------------------------------------------------------------------
100	loizides	1.4	Double_t MathUtils::DeltaR(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
101			{
102	loizides	1.5	// Compute DeltaR between two given points in the eta/phi plane.
103
104	loizides	1.7	Double_t dR = TMath::Sqrt(DeltaR2(phi1,eta1,phi2,eta2));
105			return(dR);
106			}
107
108			//--------------------------------------------------------------------------------------------------
109			Double_t MathUtils::DeltaR2(Double_t phi1, Double_t eta1, Double_t phi2, Double_t eta2)
110			{
111			// Compute DeltaR between two given points in the eta/phi plane.
112
113	loizides	1.4	Double_t dphi = DeltaPhi(phi1, phi2);
114			Double_t deta = eta1-eta2;
115	loizides	1.7	Double_t dR = dphidphi + detadeta;
116	sixie	1.1	return(dR);
117			}
118
119	loizides	1.2	//--------------------------------------------------------------------------------------------------
120	loizides	1.4	Double_t MathUtils::Eta2Theta(Double_t eta)
121	loizides	1.2	{
122	loizides	1.5	// Compute theta from given eta value.
123
124	loizides	1.2	return 2.*TMath::ATan(exp(-eta));
125	sixie	1.1	}
126
127	loizides	1.2	//--------------------------------------------------------------------------------------------------
128	loizides	1.4	Double_t MathUtils::Theta2Eta(Double_t theta)
129	loizides	1.2	{
130	loizides	1.5	// Compute eta from given theta value.
131
132	loizides	1.2	return -TMath::Log(TMath::Tan(theta/2.));
133	sixie	1.1	}