bwnorm.C
Go to the documentation of this file.
1 //
2 // Test Breit-Wigner normalization
3 //
4 // To run:
5 // root[0] .L bwnorm.C++
6 // root[1] bwnorm()
7 //
8 // Costas Andreopoulos
9 //
10 
11 #include <iostream>
12 
13 #include <TH1F.h>
14 #include <TMath.h>
15 #include <TCanvas.h>
16 #include <TGraph.h>
17 #include <TLegend.h>
18 
19 //
20 // constants
21 //
22 
23 const int kNRes = 18;
24 
25 const char * kResName[kNRes] = {
26  "P33(1232)", "S11(1535)", "D13(1520)", "S11(1650)", "D13(1700)",
27  "D15(1675)", "S31(1620)", "D33(1700)", "P11(1440)", "P33(1600)",
28  "P13(1720)", "F15(1680)", "P31(1910)", "P33(1920)", "F35(1905)",
29  "F37(1950)", "P11(1710)", "F17(1970)"
30 };
31 
32 const int kResN[kNRes] = {
33  0, 1, 1, 1, 1,
34  1, 1, 1, 1, 9,
35  2, 2, 2, 2, 2,
36  2, 2, 0
37 };
38 
39 const int kResL[kNRes] = {
40  1, 0, 2, 0, 2,
41  2, 0, 2, 1, 1,
42  1, 3, 1, 1, 3,
43  3, 1, 0
44 };
45 
46 const double kResMass[kNRes] = {
47  1.232, 1.535, 1.520, 1.650, 1.700,
48  1.675, 1.620, 1.700, 1.440, 1.600,
49  1.720, 1.680, 1.910, 1.920, 1.905,
50  1.950, 1.710, 1.970
51 };
52 
53 const double kResWidth[kNRes] = {
54  0.120, 0.150, 0.120, 0.150, 0.100,
55  0.150, 0.150, 0.300, 0.350, 0.350,
56  0.150, 0.130, 0.250, 0.200, 0.350,
57  0.300, 0.100, 0.325
58 };
59 
60 const int kResColor[kNRes] = {
61  1, 8, 2, 8, 2,
62  2, 8, 2, 1, 1,
63  1, 5, 1, 1, 5,
64  5, 1, 8
65 };
66 
67 const int kResStyle[kNRes] = {
68  1, 1, 1, 2, 2,
69  3, 3, 4, 2, 3,
70  4, 1, 2, 2, 2,
71  3, 3, 4
72 };
73 
74 //
75 // func prototypes
76 //
77 
78 void bwnorm (void);
79 void print_norm (int method, double Wmax = 3.0);
80 double bwfunc (double mR, double gRo, int L, double W);
81 double bwintegrate (double mR, double gRo, int L, double Wmin, double Wmax);
82 double bwintegrate_neugen (double mR, double gRo, int L, int n);
83 
84 
85 void bwnorm(void)
86 {
87  TCanvas * c = new TCanvas("c","",20,20,500,500);
88  c->SetFillColor(0);
89  c->SetBorderMode(0);
90 
91  TLegend * lg = new TLegend(0.15, 0.30, 0.35, 0.85);
92  lg->SetFillColor(0);
93 
94  const int kNW = 150;
95 
96  double norm[kNRes][kNW];
97  double wmax[kNRes][kNW];
98 
99  TGraph * grnorm[kNRes];
100 
101  for(int ires=0; ires<kNRes; ires++) {
102 
103  double m = kResMass [ires];
104  double w = kResWidth [ires];
105  int l = kResL [ires];
106 
107  double Wmin = 0.01;
108  double Wmax = m+w;
109  for(int iw=0; iw<kNW; iw++) {
110  Wmax+=w;
111  wmax[ires][iw] = Wmax;
112  norm[ires][iw] = bwintegrate(m,w,l,Wmin,Wmax);
113  }//w
114 
115  grnorm[ires] = new TGraph(kNW, wmax[ires], norm[ires]);
116 
117  grnorm[ires]->SetLineWidth(2);
118  grnorm[ires]->SetLineStyle(kResStyle[ires]);
119  grnorm[ires]->SetLineColor(kResColor[ires]);
120 
121  lg->AddEntry(grnorm[ires], kResName[ires], "L");
122 
123  }//res
124 
125  c->cd();
126  TH1F * hframe = (TH1F*)c->DrawFrame(1.0, 0.6, 3.0, 1.2);
127  hframe->GetXaxis()->SetTitle("W_{max} (GeV) in integration");
128  hframe->GetYaxis()->SetTitle("normalization");
129  hframe->Draw();
130  for(int ires=0; ires<kNRes; ires++) {
131  grnorm[ires]->Draw("L");
132  }
133  lg->Draw();
134  c->Update();
135 }
136 //__________________________________________________________________________
137 void print_norm(int method, double Wmax)
138 {
139  for(int ires=0; ires<kNRes; ires++) {
140 
141  double m = kResMass [ires];
142  double w = kResWidth [ires];
143  int l = kResL [ires];
144  int n = kResN [ires];
145 
146  double norm=0;
147  if (method==0) norm = bwintegrate_neugen(m,w,l,n);
148  else if (method==1) norm = bwintegrate(m,w,l,0.01,Wmax);
149 
150  cout << kResName[ires] << " --> Norm = " << norm
151  << " (multiplied by pi -> "
152  << norm * TMath::Pi() << ")" << endl;
153  }
154 }
155 //__________________________________________________________________________
156 double bwintegrate(
157  double mR, double gRo, int L, double Wmin, double Wmax)
158 {
159 // integrate within input W range
160 //
161  int N = 1000* TMath::Nint( (Wmax-Wmin)/gRo );
162  if(N%2==0) N++;
163 
164  double dW = (Wmax-Wmin)/(N-1);
165 
166  double sum = 0.5 * (bwfunc(mR,gRo,L,Wmin) + bwfunc(mR,gRo,L,Wmax));
167 
168  for(int i=1; i<N-1; i++) {
169  double W = Wmin + i*dW;
170  sum += ( bwfunc(mR,gRo,L,W) * (i%2+1) );
171  }
172  sum *= (2.*dW/3.);
173 
174  return sum;
175 }
176 //__________________________________________________________________________
177 double bwintegrate_neugen(double mR, double gRo, int L, int n)
178 {
179 // integrate using std neugen W range
180 //
181  int NW = 4;
182  if(n==2) NW=2;
183  if(n==0) NW=6;
184 
185  double Wmin = 0.001;
186  double Wmax = mR + NW*gRo;
187 
188  return bwintegrate(mR,gRo,L,Wmin,Wmax);
189 }
190 //__________________________________________________________________________
191 double bwfunc(double mR, double gRo, int L, double W)
192 {
193 // breit-wigner function
194 //
195  double pi = 3.141;
196  double mN = 0.938;
197  double mPi = 0.140;
198  double mR2 = TMath::Power(mR, 2);
199  double mN2 = TMath::Power(mN, 2);
200  double mPi2 = TMath::Power(mPi,2);
201  double W2 = TMath::Power(W, 2);
202  double qpW2 = TMath::Power(W2 - mN2 - mPi2, 2) - 4*mN2*mPi2;
203  double qpM2 = TMath::Power(mR2 - mN2 - mPi2, 2) - 4*mN2*mPi2;
204  double qpW = TMath::Sqrt(TMath::Max(0.,qpW2)) / (2*W);
205  double qpM = TMath::Sqrt(TMath::Max(0.,qpM2)) / (2*mR);
206  double gR = gRo * TMath::Power( qpW/qpM, 2*L+1 );
207  double gR2 = TMath::Power(gR, 2);
208  double brwg = (0.5/pi)*gR / (TMath::Power(W-mR, 2) + 0.25*gR2);
209  return brwg;
210 }
211 //__________________________________________________________________________
double bwintegrate_neugen(double mR, double gRo, int L, int n)
Definition: bwnorm.C:177
double Wmax[kNWBins]
const int kResN[kNRes]
Definition: bwnorm.C:32
double bwfunc(double mR, double gRo, int L, double W)
Definition: bwnorm.C:191
const char * kResName[kNRes]
Definition: bwnorm.C:25
static constexpr double L
const int kResL[kNRes]
Definition: bwnorm.C:39
const double kResMass[kNRes]
Definition: bwnorm.C:46
void print_norm(int method, double Wmax=3.0)
Definition: bwnorm.C:137
OStream cout
Definition: OStream.cxx:6
const int kResColor[kNRes]
Definition: bwnorm.C:60
void bwnorm(void)
Definition: bwnorm.C:85
Float_t norm
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
double Wmin[kNWBins]
const int kNRes
Definition: bwnorm.C:23
Double_t sum
Definition: plot.C:31
double bwintegrate(double mR, double gRo, int L, double Wmin, double Wmax)
Definition: bwnorm.C:156
const double kResWidth[kNRes]
Definition: bwnorm.C:53
Float_t w
Definition: plot.C:20
#define W(x)
const int kResStyle[kNRes]
Definition: bwnorm.C:67