Skip to content

Commit 31d2605

Browse files
JackHLindonlmoneta
JackHLindon
authored and
lmoneta
committed
Add Relativistic Breit Wigner Function to TMath
Implemented relativistic version of breit wigner (non-relativistic case already exists in TMath). Define BreitWignerRelativistic(Double_t x, Double_t median=0, Double_t gamma=1) in TMath.h and then calculate in TMath.cxx, a similar non relativistic function BreitWigner already exists in the same location which was used as a template. A tutorial BreitWigner.C has been added in tutorials/math which produces plots comparing the non relativistic and relativistic case
1 parent 60ec158 commit 31d2605

File tree

4 files changed

+210
-0
lines changed

4 files changed

+210
-0
lines changed

‎math/mathcore/inc/TMath.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -482,6 +482,7 @@ struct Limits {
482482
Double_t Binomial(Int_t n,Int_t k); // Calculate the binomial coefficient n over k
483483
Double_t BinomialI(Double_t p, Int_t n, Int_t k);
484484
Double_t BreitWigner(Double_t x, Double_t mean=0, Double_t gamma=1);
485+
Double_t BreitWignerRelativistic(Double_t x, Double_t median=0, Double_t gamma=1);
485486
Double_t CauchyDist(Double_t x, Double_t t=0, Double_t s=1);
486487
Double_t ChisquareQuantile(Double_t p, Double_t ndf);
487488
Double_t FDist(Double_t F, Double_t N, Double_t M);

‎math/mathcore/src/TMath.cxx‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -445,6 +445,24 @@ Double_t TMath::BreitWigner(Double_t x, Double_t mean, Double_t gamma)
445445
return bw/(2*Pi());
446446
}
447447

448+
////////////////////////////////////////////////////////////////////////////////
449+
/// Calculates a Relativistic Breit Wigner function with median and gamma.
450+
// \f$ BW(E) = \frac{2\sqrt{2}}{\pi}\frac{M^{2}\gamma\sqrt{M^{2} + \gamma^{2}}}{\left(\sqrt{M^{2}+M\sqrt{M^{2} + \gamma^{2}}}\right)\left(\left(E^{2} - M^{2}\right)^{2} + M^{2}\gamma^{2}\right)} \f$
451+
452+
Double_t TMath::BreitWignerRelativistic(Double_t x, Double_t median, Double_t gamma)
453+
{
454+
Double_t mm = median*median;
455+
Double_t gg = gamma*gamma;
456+
Double_t mg = median*gamma;
457+
Double_t xxMinusmm = x*x - mm;
458+
459+
Double_t y = sqrt(mm * (mm + gg));
460+
Double_t k = (0.90031631615710606*mg*y)/(sqrt(mm+y)); //2*sqrt(2)/pi = 0.90031631615710606
461+
462+
Double_t bw = k/(xxMinusmm*xxMinusmm + mg*mg);
463+
return bw;
464+
}
465+
448466
////////////////////////////////////////////////////////////////////////////////
449467
/// Calculates a gaussian function with mean and sigma.
450468
/// If norm=kTRUE (default is kFALSE) the result is divided

‎math/mathcore/test/testTMath.cxx‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -141,6 +141,20 @@ void testPlane()
141141
<< endl;
142142
}
143143

144+
void testBreitWignerRelativistic()
145+
{
146+
Double_t median = 5000;
147+
Double_t gamma = 100;
148+
Int_t nPoints = 10;
149+
Double_t xMinimum = 0; Double_t xMaximum = 10000;
150+
Double_t xStepSize = (xMaximum-xMinimum)/nPoints;
151+
152+
for (Int_t i=0;i<=nPoints;i++) {
153+
Double_t currentX = xMinimum+i*xStepSize;
154+
cout << "BreitWignerRelativistic(" << currentX << "," << median << "," << gamma << ") = " << BreitWignerRelativistic(currentX,median,gamma) << endl;
155+
}
156+
}
157+
144158
void testTMath()
145159
{
146160
cout << "Starting tests on TMath..." << endl;
@@ -179,6 +193,10 @@ void testTMath()
179193

180194
testPlane<Double_t>();
181195
testPlane<Float_t>();
196+
197+
cout << "\nBreitWignerRelativistic tests: " << endl;
198+
199+
testBreitWignerRelativistic();
182200
}
183201

184202
int main()

‎tutorials/math/BreitWigner.C‎

Lines changed: 173 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,173 @@
1+
/// \file
2+
/// \ingroup tutorial_math
3+
/// \notebook
4+
/// Tutorial illustrating how to create a plot comparing a Breit Wigner to a Relativistic Breit Wigner
5+
///
6+
/// can be run with:
7+
///
8+
/// ~~~{.cpp}
9+
/// root[0] .x BreitWigner.C
10+
/// ~~~
11+
///
12+
/// \macro_image
13+
/// \macro_code
14+
///
15+
/// \author Jack Lindon
16+
17+
#include "TMath.h"
18+
19+
#include <limits>
20+
#include <string>
21+
#include "TAxis.h"
22+
#include "TGraph.h"
23+
#include "TCanvas.h"
24+
#include "TLatex.h"
25+
#include "TLegend.h"
26+
#include "TStyle.h" //For gStyle to remove stat box.
27+
28+
29+
void plotTwoTGraphs(Double_t x[], Double_t y1[], Double_t y2[], const Int_t nPoints
30+
, Double_t lowerXLimit, Double_t upperXLimit
31+
, Double_t lowerYLimit, Double_t upperYLimit
32+
, std::string legend1, std::string legend2
33+
, std::string plotTitle1, std::string plotTitle2, std::string plotTitle3
34+
, std::string pdfTitle){
35+
36+
///////////////////////////////////////////////////////
37+
//Define variables for plot aesthetics and positioning
38+
Double_t legendXPos = 0.63; Double_t legendYPos = 0.85;
39+
Double_t legendXWidth = 0.29; Double_t legendYHeight = 0.1;
40+
Double_t plotTitleXPos = 0.23; Double_t plotTitleYPos = 0.25;
41+
Double_t fontSize = 0.04;
42+
Double_t lineWidth = 2;
43+
bool setLimitPlotLogScale = true;
44+
std::string xAxisTitle = "E [GeV]"; std::string yAxisTitle = "Events";
45+
Double_t xAxisTitleOffset = 1; Double_t yAxisTitleOffset = 1.3;
46+
gStyle->SetOptStat(0);
47+
48+
///////////////////////////////////////////////////////
49+
// Initialize TGraphs
50+
TGraph* gr1 = new TGraph(nPoints,x,y1);
51+
TGraph* gr2 = new TGraph(nPoints,x,y2);
52+
gr1->SetLineWidth(lineWidth);
53+
gr2->SetLineWidth(lineWidth);
54+
gr1->SetLineColor(kBlack);
55+
gr2->SetLineColor(kBlue);
56+
57+
/////////////////////////////////////////////////////////
58+
// Initialize canvas
59+
TCanvas* c1 = new TCanvas("c1","transparent pad",200,10,600,600);
60+
c1->SetLogy(setLimitPlotLogScale);
61+
c1->SetTicks(1,1);
62+
c1->SetRightMargin(0.02);
63+
c1->SetTopMargin(0.02);
64+
65+
66+
///////////////////////////////////////////////////////
67+
//Make just a basic invisible TGraph just for the axes
68+
const Double_t axis_x[2] = {lowerXLimit,upperXLimit};
69+
const Double_t axis_y[2] = {lowerYLimit,upperYLimit};
70+
TGraph* grAxis = new TGraph(2, axis_x, axis_y);
71+
grAxis->SetTitle("");
72+
grAxis->GetYaxis()->SetTitle(yAxisTitle.c_str());
73+
grAxis->GetXaxis()->SetTitle(xAxisTitle.c_str());
74+
grAxis->GetXaxis()->SetRangeUser(lowerXLimit,upperXLimit);
75+
grAxis->GetYaxis()->SetRangeUser(lowerYLimit,upperYLimit);
76+
grAxis->GetXaxis()->SetLabelSize(fontSize);
77+
grAxis->GetYaxis()->SetLabelSize(fontSize);
78+
grAxis->GetXaxis()->SetTitleSize(fontSize);
79+
grAxis->GetYaxis()->SetTitleSize(fontSize);
80+
grAxis->GetXaxis()->SetTitleOffset(xAxisTitleOffset);
81+
grAxis->GetYaxis()->SetTitleOffset(yAxisTitleOffset);
82+
grAxis->SetLineWidth(0);//So invisible
83+
84+
///////////////////////////////////////////////////////////
85+
// Make legend and set aesthetics
86+
auto legend = new TLegend(legendXPos,legendYPos,legendXPos+legendXWidth,legendYPos+legendYHeight);
87+
legend->SetFillStyle(0);
88+
legend->SetBorderSize(0);
89+
legend->SetTextSize(fontSize);
90+
legend->AddEntry(gr1,legend1.c_str(),"L");
91+
legend->AddEntry(gr2,legend2.c_str(),"L");
92+
93+
94+
/////////////////////////////////////////////////////////////
95+
// Add plot title to plot. Make in three lines so not crowded.
96+
// Shift each line down by shiftY
97+
float shiftY{0.037};
98+
TLatex* tex_Title = new TLatex(plotTitleXPos,plotTitleYPos-0*shiftY,plotTitle1.c_str());
99+
tex_Title->SetNDC();
100+
tex_Title->SetTextFont(42);
101+
tex_Title->SetTextSize(fontSize);
102+
tex_Title->SetLineWidth(lineWidth);
103+
TLatex* tex_Title2 = new TLatex(plotTitleXPos,plotTitleYPos-1*shiftY,plotTitle2.c_str());
104+
tex_Title2->SetNDC();
105+
tex_Title2->SetTextFont(42);
106+
tex_Title2->SetTextSize(fontSize);
107+
tex_Title2->SetLineWidth(lineWidth);
108+
TLatex* tex_Title3 = new TLatex(plotTitleXPos,plotTitleYPos-2*shiftY,plotTitle3.c_str());
109+
tex_Title3->SetNDC();
110+
tex_Title3->SetTextFont(42);
111+
tex_Title3->SetTextSize(fontSize);
112+
tex_Title3->SetLineWidth(lineWidth);
113+
114+
115+
/////////////////////////////////////
116+
// Draw everything
117+
grAxis->Draw("AL");
118+
gr1->Draw("L same");
119+
gr2->Draw("L same");
120+
legend->Draw();
121+
tex_Title->Draw();
122+
tex_Title2->Draw();
123+
tex_Title3->Draw();
124+
c1->RedrawAxis(); //Be sure to redraw axis AFTER plotting TGraphs otherwise TGraphs will be on top of tick marks and axis borders.
125+
126+
gPad->Print(pdfTitle.c_str());
127+
128+
}
129+
130+
131+
void BreitWigner(){
132+
133+
/////////////////////////////////////////////////////////
134+
// Define x axis limits and steps for each plotted point
135+
const Int_t nPoints = 1000;
136+
Double_t xMinimum = 0; Double_t xMaximum = 13000;
137+
Double_t xStepSize = (xMaximum-xMinimum)/nPoints;
138+
139+
///////////////////////////////////////////////////////
140+
// Define arrays of (x,y) points.
141+
Double_t x[nPoints];
142+
Double_t y_nonRelBW[nPoints], y_relBW[nPoints];
143+
144+
145+
//////////////////////////////////
146+
// Define Breit-Wigner parameters
147+
Double_t width = 1350;
148+
Double_t sigma = 269.7899;
149+
Double_t median = 9000;
150+
151+
///////////////////////////////////////////////////
152+
// Loop over x axis range, filling in (x,y) points,
153+
// and finding y minimums and maximums for axis limit.
154+
Double_t yMinimum = std::numeric_limits<Double_t>::max();
155+
Double_t yMaximum = TMath::BreitWignerRelativistic(median,median,width); //y maximum is at x=median (and non relativistic = relativistic at median so choice of function does not matter).
156+
for (Int_t i=0;i<nPoints;i++) {
157+
Double_t currentX = xMinimum+i*xStepSize;
158+
x[i] = currentX;
159+
y_nonRelBW[i] = TMath::BreitWigner(currentX,median,width);
160+
y_relBW[i] = TMath::BreitWignerRelativistic(currentX,median,width);
161+
162+
if (y_nonRelBW[i]<yMinimum){yMinimum = y_nonRelBW[i];}
163+
if (y_relBW[i]<yMinimum){yMinimum = y_relBW[i];}
164+
}
165+
166+
plotTwoTGraphs(x, y_nonRelBW, y_relBW, nPoints
167+
, xMinimum, xMaximum //xAxis limits
168+
, yMinimum/4, yMaximum*4 //yAxis limits, expand for aesthetics.
169+
,"NonRel BW", "Rel BW" //Legend entries
170+
, "Comparing BW", "M = " + std::to_string(int(round(median))) + " GeV", "#Gamma = " + std::to_string(int(round(width))) + " GeV" //Plot Title entry (three lines)
171+
, "BW_M"+std::to_string(int(round(median)))+"_Gamma" + std::to_string(int(round(width))) +".pdf)" //PDF file title.
172+
);
173+
}

0 commit comments

Comments
 (0)