MathTools/src/Funcs.cc

// $Id: Funcs.cc,v 1.3 2009/07/20 04:55:44 loizides Exp $

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

ClassImp(mithep::Funcs)

using namespace mithep;

//--------------------------------------------------------------------------------------------------
Double_t Funcs::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 Funcs::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 Funcs::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 Funcs::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 Funcs::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 Funcs::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 Funcs::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 Funcs::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 Funcs::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 *Funcs::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::Funcs::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.1,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,5);
  f->SetParLimits(5,0.1,5);
  f->SetParLimits(6,0,1);
  f->SetParLimits(7,0,10);
  f->FixParameter(8,xmin);
  f->FixParameter(9,xmax);

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