1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
|
// Novosibirsk function
// See H. Ikeda et al. / Nuclear Instruments and Methods in Physics Research A 441 (2000) 401-426
// Also implemented as RooNovosibirsk, with tail --> -sigma0
// Adapted from code by Chris Hearty by Jean-François Caron.
#include <cmath>
#include <TMath.h>
#include <TF1.h>
Double_t novosibirsk(Double_t * x, Double_t * par)
{ // parameters are:
// [0] -> amplitude/normalization
// [1] -> peak location
// [2] -> width parameter (FWMH/2.36)
// [3] -> tail parameter epsilon
// In the Belle paper, the parameters are
// [0] = N
// [1] = x_p
// [2] = sigma_E
// [3] = eta
Double_t qa=0,qb=0,qc=0,qx=0,qy=0;
const Double_t peak = par[1];
const Double_t width = par[2];
static const Double_t sln4 = TMath::Sqrt(TMath::Log(4));
const Double_t y = par[3]*sln4;
const Double_t tail = -TMath::Log(y + TMath::Sqrt(1 + y*y))/sln4;
if(TMath::Abs(tail) < 1e-7)
{ // If the tail variable is small enough, this is just a Gaussian.
//qc = 0.5*TMath::Power(((x[0]-peak)/width),2);
qc = 0.5*((x[0]-peak)*(x[0]-peak)/width/width);
}
else
{
qa = tail*sln4;
qb = sinh(qa)/qa;
qx = -1*(x[0]-peak)/width*qb; // The -1 here swaps the X-axis from Chris' original.
qy = 1.0+tail*qx;
if( qy > 1e-7 )
{
//qc = 0.5*(TMath::Power((TMath::Log(qy)/tail),2) + tail*tail);
const Double_t lqt = TMath::Log(qy)/tail;
qc = 0.5*(lqt*lqt + tail*tail);
}
else
{
qc = 15.0;
}
}
return par[0]*TMath::Exp(-qc);
}
double belle_novosibirsk(const double * const x, const double * const c)
{
const double N = c[0];
const double x_p = c[1];
const double sigma_E = c[2];
const double eta = c[3];
static const double xi = 2*std::sqrt(std::log(4));
const double sigma_0 = (2.0/xi)*std::asinh(eta*xi/2.0);
const double sigma_0_squared = sigma_0*sigma_0;
const double logterm = std::log1p((x[0]-x_p)*eta/sigma_E);//log1p is
//log(1+arg);
// log1p( (x-x_p)*eta/sigma_E ) limits the values of x.
// x > x_p - sigma_E/eta
const double expo_arg = -0.5*(logterm*logterm/sigma_0_squared + sigma_0_squared);
return N*std::exp(expo_arg);
}
TF1 f_novosibirsk()
{
TF1 f("f_novo",novosibirsk,0,10,4);
f.SetParameters(1,1,1,1);
f.SetParNames("N","x_p","omega (FWHM/2.36)","epsilon");
return f;
}
TF1 f_belle()
{
TF1 f("f_belle",belle_novosibirsk,0,10,4);
f.SetParameters(1,1,1,1);
f.SetParNames("N","x_0","sigma_E (FWHM/2.36)","eta");
return f;
}
// Double_t simple_novosibirsk(Double_t * x, Double_t * par)
// {
// const double N = par[0];
// const Double_t x_p = par[1];
// const Double_t omega = par[2];
// const double epsilon = par[3];
// static const Double_t xi = 2*TMath::Sqrt(TMath::Log(4));
// const Double_t tail = -2.0/xi*TMath::ASinH(epsilon*xi/2.0);
// const Double_t lqt = TMath::Log(1+(x[0]-x_p)*epsilon/omega)/tail;
// return N*TMath::Exp(-0.5*(lqt*lqt + tail*tail));
// }
// void test()
// {
// TF1 * f_chris = new TF1("f_chris",novosibirsk,0,10,4);
// TF1 * f_belle = new TF1("f_belle",Real_Novosibirsk,0,10,4);
// TF1 * f_simple = new TF1("f_simple",simple_novosibirsk,0,10,4);
// f_chris->SetParameters(1,1,1,1);
// f_belle->SetParameters(1,1,1,1);
// f_simple->SetParameters(1,1,1,1);
// f_chris->SetParNames("N","x_o","omega","epsilon");
// f_belle->SetParNames("N","x_p","sigma_E","eta");
// f_simple->SetParNames("N","x_o","omega","epsilon");
// f_chris->SetLineColor(kRed);
// f_belle->SetLineColor(kBlue);
// f_simple->SetLineColor(kGreen);
// TCanvas * c1 = new TCanvas();
// c1->cd();
// f_belle->Draw();
// f_chris->Draw("same");
// //f_simple->Draw("same");
// }
|