MathTools/src/TPFuncs.cc

// $Id: TPFuncs.cc,v 1.2 2009/06/11 20:26:57 loizides Exp $

#include "MitCommon/MathTools/interface/TPFuncs.h"
#include <TF1.h>

using namespace mithep;

//--------------------------------------------------------------------------------------------------
Double_t TPFuncs::BreitGaus(Double_t x, Double_t m, Double_t mwidth, Double_t msig, 
                            Double_t fintf, Double_t xmin, Double_t xmax)
{
  // Breit-Wigner convoluted with Gaussian.
  // Fit parameters:
  //  m: Most probable location) Breit mean
  //  mwidth: Width (scale) Breit-Wigner
  //  msig: Width (sigma) of convoluted Gaussian function
  //  fintf: interference fraction between Z and gamma
  //  xmin: minimum x (for normalization)
  //  xmax: maximum x (for normalization)

  // Numeric constants
  const Double_t invsq2pi = 0.3989422804014; //(2pi)^(-1/2)

  // Control constants
  const Double_t np = 300; // number of convolution steps
  const Double_t sc = 3.0; // convolution extends to +-sc Gaussian sigmas

  // Range of convolution integral
  const Double_t xlow = x - sc * msig;
  const Double_t xupp = x + sc * msig;
  const Double_t step = (xupp-xlow) / np;

  // Convolution integral of Breit and Gaussian by sum
  Double_t sum0 = 0.0;
  if (fintf<1) {
    for(Double_t i=1.0; i<=np/2; i++) {
      Double_t xx = xlow + (i-0.5) * step;
      Double_t fland = BreitWignerZ(xx,m,mwidth);
      sum0 += fland * TMath::Gaus(x,xx,msig);
      xx = xupp - (i-0.5) * step;
      fland = BreitWignerZ(xx,m,mwidth);
      sum0 += fland * TMath::Gaus(x,xx,msig);
    }
    sum0 = (step * sum0 * invsq2pi / msig);
  }

  Double_t sum1 = 0.0;
  if (fintf>0) {
    for(Double_t i=1.0; i<=np/2; i++) {
      Double_t xx = xlow + (i-0.5) * step;
      Double_t fland = BreitWignerGamma(xx,m,mwidth);
      sum1 += fland * TMath::Gaus(x,xx,msig);
      xx = xupp - (i-0.5) * step;
      fland = BreitWignerGamma(xx,m,mwidth);
      sum1 += fland * TMath::Gaus(x,xx,msig);
    }
    sum1 = (step * sum1 * invsq2pi / msig / (xmax-xmin));
  }
  return ((1.-fintf)*sum0+fintf*sum1);
}

//--------------------------------------------------------------------------------------------------
Double_t TPFuncs::BreitGaus(Double_t *x, Double_t *par)
{
  // Breit-Wigner convoluted with Gaussian.
  // Fit parameters:
  //  norm: Normalization factor
  //  m: Most probable location) Breit mean
  //  mwidth: Width (scale) Breit-Wigner
  //  msig: Width (sigma) of convoluted Gaussian function
  //  fintf: interference fraction between Z and gamma
  //  xmin: minimum x (for normalization)
  //  xmax: maximum x (for normalization)

  Double_t xx     = x[0];
  Double_t norm   = par[0];
  Double_t m      = par[1];
  Double_t mwidth = par[2];
  Double_t msig   = par[3];
  Double_t fintf  = par[4];
  Double_t xmin   = par[5];
  Double_t xmax   = par[6];
  Double_t ret = norm * BreitGaus(xx, m, mwidth, msig, fintf, xmin, xmax);
  return ret;
}

//--------------------------------------------------------------------------------------------------
Double_t TPFuncs::BreitWignerZ(Double_t x, Double_t mean, Double_t gamma)
{
  // Breit-Wigner for Z peak.

  Double_t bw =  (gamma*mean*mean)/
                 ((x*x-mean*mean)*(x*x-mean*mean) + x*x*x*x*(gamma*gamma)/(mean*mean));
  return bw / (TMath::Pi()/2.0);
}

//--------------------------------------------------------------------------------------------------
Double_t TPFuncs::BreitWignerZ(Double_t *x, Double_t *par)
{
  // Breit-Wigner for Z peak.

  Double_t xx    = x[0];
  Double_t norm  = par[0];
  Double_t mean  = par[1];
  Double_t gamma = par[2];

  Double_t ret = norm * BreitWignerZ(xx,mean,gamma);
  return ret;
}

//--------------------------------------------------------------------------------------------------
Double_t TPFuncs::BreitWignerGamma(Double_t x, Double_t mean, Double_t gamma)
{
  // Breit-Wigner for gamma interference.

  Double_t bw =  ((x*x-mean*mean)*(x*x-mean*mean))/
                 ((x*x-mean*mean)*(x*x-mean*mean) + x*x*x*x*(gamma*gamma)/(mean*mean));
  return bw;
}

//--------------------------------------------------------------------------------------------------
Double_t TPFuncs::BreitWignerGamma(Double_t *x, Double_t *par)
{
  // Breit-Wigner for gamma interference.

  Double_t xx    = x[0];
  Double_t norm  = par[0];
  Double_t mean  = par[1];
  Double_t gamma = par[2];

  Double_t ret = norm * BreitWignerGamma(xx,mean,gamma);
  return ret;
}

//--------------------------------------------------------------------------------------------------
Double_t TPFuncs::ExpRange(Double_t mass, Double_t lambda, Double_t xmin, Double_t xmax)
{
  // Probability of an exponential in a given interval.

  if (mass < xmin || mass > xmax)
    return 0.0;

  if (lambda == 0.0)
    return 1.0/(xmax-xmin);

  Double_t xmed  = 0.5*(xmin+xmax);
  Double_t integ = lambda * TMath::Exp(-lambda*xmed) /
    (TMath::Exp(-lambda*xmin)-TMath::Exp(-lambda*xmax));

  return integ*TMath::Exp(-lambda*(mass-xmed));
}

//--------------------------------------------------------------------------------------------------
Double_t TPFuncs::ExpRange(Double_t *x, Double_t *par)
{
  // Probability of an exponential in a given interval.

  Double_t xx     = x[0];
  Double_t norm   = par[0];
  Double_t lambda = par[1];
  Double_t xmin   = par[2];
  Double_t xmax   = par[3];

  Double_t ret = norm * ExpRange(xx,lambda,xmin,xmax);
  return ret;
}

//--------------------------------------------------------------------------------------------------
Double_t TPFuncs::ZLineShapePlusBkg(Double_t *x, Double_t *par)
{ 
  // Z line shape (BreitGaus) plus background (ExpRange).

  // Parameters:
  //  0: overall norm
  //  1: relative fraction of signal and background
  //  2: relative fraction of exponential and polynomial
  //  3: most probable location) Breit mean
  //  4: width (scale) Breit-Wigner
  //  5: width (sigma) of convoluted Gaussian function
  //  6: interference fraction between Z and gamma
  //  7: lambda of exponential
  //  8: minimum x (for normalization)
  //  9: maximum x (for normalization)

  Double_t xmin=par[8];
  Double_t xmax=par[9];
  Double_t s = (1.0 - par[1]) * BreitGaus(x[0], par[3], par[4], par[5], par[6], xmin, xmax);
  Double_t b = par[1] * ((1.0 - par[2]) * ExpRange(x[0], par[7], xmin, xmax) + par[2]/(xmax-xmin));

  return (s+b)*par[0];
}

//--------------------------------------------------------------------------------------------------
TF1 *TPFuncs::CreateZLineShapePlusBkg(Double_t norm, Double_t xmin, Double_t xmax, const char *n)
{
  // Create TF1 for Z line shape and background.

  TF1 *f = new TF1(n,mithep::TPFuncs::ZLineShapePlusBkg,xmin,xmax,10);
  f->SetParName(0,"norm");
  f->SetParName(1,"f_sb");
  f->SetParName(2,"f_eb");
  f->SetParName(3,"mass");
  f->SetParName(4,"width");
  f->SetParName(5,"gsig");
  f->SetParName(6,"f_int");
  f->SetParName(7,"lambda");
  f->SetParName(8,"min");
  f->SetParName(9,"max");

  f->SetParameters(norm,0.5,0.5,91.,1.,1.,0.,0.1,xmin,xmax);
  f->SetParLimits(0,0,1e6);
  f->SetParLimits(1,0,1);
  f->SetParLimits(2,0,1);
  f->SetParLimits(3,xmin,xmax);
  f->SetParLimits(4,0.1,3);
  f->SetParLimits(5,0.1,3);
  f->SetParLimits(6,0,0);
  f->SetParLimits(7,0,10);
  f->SetParLimits(8,xmin,xmin);
  f->SetParLimits(9,xmax,xmax);

  return f;
}
Revision:	1.3
Committed:	Thu Jun 11 20:31:19 2009 UTC (15 years, 10 months ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.2:	+1 -1 lines
State:	*FILE REMOVED*
Log Message:	Renamed since neither class nor function names was really precise.
#	User	Rev	Content
1	loizides	1.3	// $Id: TPFuncs.cc,v 1.2 2009/06/11 20:26:57 loizides Exp $
2	loizides	1.1
3			#include "MitCommon/MathTools/interface/TPFuncs.h"
4			#include <TF1.h>
5
6			using namespace mithep;
7
8			//--------------------------------------------------------------------------------------------------
9			Double_t TPFuncs::BreitGaus(Double_t x, Double_t m, Double_t mwidth, Double_t msig,
10			Double_t fintf, Double_t xmin, Double_t xmax)
11			{
12			// Breit-Wigner convoluted with Gaussian.
13			// Fit parameters:
14			// m: Most probable location) Breit mean
15			// mwidth: Width (scale) Breit-Wigner
16			// msig: Width (sigma) of convoluted Gaussian function
17			// fintf: interference fraction between Z and gamma
18			// xmin: minimum x (for normalization)
19			// xmax: maximum x (for normalization)
20
21			// Numeric constants
22			const Double_t invsq2pi = 0.3989422804014; //(2pi)^(-1/2)
23
24			// Control constants
25			const Double_t np = 300; // number of convolution steps
26			const Double_t sc = 3.0; // convolution extends to +-sc Gaussian sigmas
27
28			// Range of convolution integral
29			const Double_t xlow = x - sc * msig;
30			const Double_t xupp = x + sc * msig;
31			const Double_t step = (xupp-xlow) / np;
32
33			// Convolution integral of Breit and Gaussian by sum
34			Double_t sum0 = 0.0;
35			if (fintf<1) {
36			for(Double_t i=1.0; i<=np/2; i++) {
37			Double_t xx = xlow + (i-0.5) * step;
38			Double_t fland = BreitWignerZ(xx,m,mwidth);
39			sum0 += fland * TMath::Gaus(x,xx,msig);
40			xx = xupp - (i-0.5) * step;
41			fland = BreitWignerZ(xx,m,mwidth);
42			sum0 += fland * TMath::Gaus(x,xx,msig);
43			}
44			sum0 = (step * sum0 * invsq2pi / msig);
45			}
46
47			Double_t sum1 = 0.0;
48			if (fintf>0) {
49			for(Double_t i=1.0; i<=np/2; i++) {
50			Double_t xx = xlow + (i-0.5) * step;
51			Double_t fland = BreitWignerGamma(xx,m,mwidth);
52			sum1 += fland * TMath::Gaus(x,xx,msig);
53			xx = xupp - (i-0.5) * step;
54			fland = BreitWignerGamma(xx,m,mwidth);
55			sum1 += fland * TMath::Gaus(x,xx,msig);
56			}
57			sum1 = (step * sum1 * invsq2pi / msig / (xmax-xmin));
58			}
59			return ((1.-fintf)sum0+fintfsum1);
60			}
61
62			//--------------------------------------------------------------------------------------------------
63			Double_t TPFuncs::BreitGaus(Double_t x, Double_t par)
64			{
65			// Breit-Wigner convoluted with Gaussian.
66			// Fit parameters:
67			// norm: Normalization factor
68			// m: Most probable location) Breit mean
69			// mwidth: Width (scale) Breit-Wigner
70			// msig: Width (sigma) of convoluted Gaussian function
71			// fintf: interference fraction between Z and gamma
72			// xmin: minimum x (for normalization)
73			// xmax: maximum x (for normalization)
74
75			Double_t xx = x[0];
76			Double_t norm = par[0];
77			Double_t m = par[1];
78			Double_t mwidth = par[2];
79			Double_t msig = par[3];
80			Double_t fintf = par[4];
81			Double_t xmin = par[5];
82			Double_t xmax = par[6];
83			Double_t ret = norm * BreitGaus(xx, m, mwidth, msig, fintf, xmin, xmax);
84			return ret;
85			}
86
87			//--------------------------------------------------------------------------------------------------
88			Double_t TPFuncs::BreitWignerZ(Double_t x, Double_t mean, Double_t gamma)
89			{
90			// Breit-Wigner for Z peak.
91
92			Double_t bw = (gammameanmean)/
93			((xx-meanmean)(xx-meanmean) + xxxx(gammagamma)/(mean*mean));
94			return bw / (TMath::Pi()/2.0);
95			}
96
97			//--------------------------------------------------------------------------------------------------
98			Double_t TPFuncs::BreitWignerZ(Double_t x, Double_t par)
99			{
100			// Breit-Wigner for Z peak.
101
102			Double_t xx = x[0];
103			Double_t norm = par[0];
104			Double_t mean = par[1];
105			Double_t gamma = par[2];
106
107			Double_t ret = norm * BreitWignerZ(xx,mean,gamma);
108			return ret;
109			}
110
111			//--------------------------------------------------------------------------------------------------
112			Double_t TPFuncs::BreitWignerGamma(Double_t x, Double_t mean, Double_t gamma)
113			{
114			// Breit-Wigner for gamma interference.
115
116			Double_t bw = ((xx-meanmean)(xx-mean*mean))/
117			((xx-meanmean)(xx-meanmean) + xxxx(gammagamma)/(mean*mean));
118			return bw;
119			}
120
121			//--------------------------------------------------------------------------------------------------
122			Double_t TPFuncs::BreitWignerGamma(Double_t x, Double_t par)
123			{
124			// Breit-Wigner for gamma interference.
125
126			Double_t xx = x[0];
127			Double_t norm = par[0];
128			Double_t mean = par[1];
129			Double_t gamma = par[2];
130
131			Double_t ret = norm * BreitWignerGamma(xx,mean,gamma);
132			return ret;
133			}
134
135			//--------------------------------------------------------------------------------------------------
136			Double_t TPFuncs::ExpRange(Double_t mass, Double_t lambda, Double_t xmin, Double_t xmax)
137			{
138			// Probability of an exponential in a given interval.
139
140			if (mass < xmin \|\| mass > xmax)
141			return 0.0;
142
143			if (lambda == 0.0)
144			return 1.0/(xmax-xmin);
145
146			Double_t xmed = 0.5*(xmin+xmax);
147			Double_t integ = lambda * TMath::Exp(-lambda*xmed) /
148			(TMath::Exp(-lambdaxmin)-TMath::Exp(-lambdaxmax));
149
150			return integTMath::Exp(-lambda(mass-xmed));
151			}
152
153			//--------------------------------------------------------------------------------------------------
154			Double_t TPFuncs::ExpRange(Double_t x, Double_t par)
155			{
156			// Probability of an exponential in a given interval.
157
158			Double_t xx = x[0];
159			Double_t norm = par[0];
160			Double_t lambda = par[1];
161			Double_t xmin = par[2];
162			Double_t xmax = par[3];
163
164			Double_t ret = norm * ExpRange(xx,lambda,xmin,xmax);
165			return ret;
166			}
167
168			//--------------------------------------------------------------------------------------------------
169	loizides	1.2	Double_t TPFuncs::ZLineShapePlusBkg(Double_t x, Double_t par)
170	loizides	1.1	{
171	loizides	1.2	// Z line shape (BreitGaus) plus background (ExpRange).
172	loizides	1.1
173			// Parameters:
174			// 0: overall norm
175			// 1: relative fraction of signal and background
176			// 2: relative fraction of exponential and polynomial
177			// 3: most probable location) Breit mean
178			// 4: width (scale) Breit-Wigner
179			// 5: width (sigma) of convoluted Gaussian function
180			// 6: interference fraction between Z and gamma
181			// 7: lambda of exponential
182			// 8: minimum x (for normalization)
183			// 9: maximum x (for normalization)
184
185			Double_t xmin=par[8];
186			Double_t xmax=par[9];
187			Double_t s = (1.0 - par[1]) * BreitGaus(x[0], par[3], par[4], par[5], par[6], xmin, xmax);
188			Double_t b = par[1] * ((1.0 - par[2]) * ExpRange(x[0], par[7], xmin, xmax) + par[2]/(xmax-xmin));
189
190			return (s+b)*par[0];
191			}
192
193			//--------------------------------------------------------------------------------------------------
194	loizides	1.2	TF1 TPFuncs::CreateZLineShapePlusBkg(Double_t norm, Double_t xmin, Double_t xmax, const char n)
195	loizides	1.1	{
196	loizides	1.2	// Create TF1 for Z line shape and background.
197	loizides	1.1
198	loizides	1.2	TF1 *f = new TF1(n,mithep::TPFuncs::ZLineShapePlusBkg,xmin,xmax,10);
199	loizides	1.1	f->SetParName(0,"norm");
200			f->SetParName(1,"f_sb");
201			f->SetParName(2,"f_eb");
202			f->SetParName(3,"mass");
203			f->SetParName(4,"width");
204			f->SetParName(5,"gsig");
205			f->SetParName(6,"f_int");
206			f->SetParName(7,"lambda");
207			f->SetParName(8,"min");
208			f->SetParName(9,"max");
209
210			f->SetParameters(norm,0.5,0.5,91.,1.,1.,0.,0.1,xmin,xmax);
211			f->SetParLimits(0,0,1e6);
212			f->SetParLimits(1,0,1);
213			f->SetParLimits(2,0,1);
214			f->SetParLimits(3,xmin,xmax);
215			f->SetParLimits(4,0.1,3);
216			f->SetParLimits(5,0.1,3);
217			f->SetParLimits(6,0,0);
218			f->SetParLimits(7,0,10);
219			f->SetParLimits(8,xmin,xmin);
220			f->SetParLimits(9,xmax,xmax);
221
222			return f;
223			}