Introducing a new AOD class: AliAODcascade (A.Maire)
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / macros / AliRsnView.C
CommitLineData
ec908cae 1#include <Riostream.h>
2#include <TMath.h>
3#include <TH1.h>
4#include <TF1.h>
5#include <TFile.h>
6#include <TCanvas.h>
7#include <TStyle.h>
8
9Double_t bw(Double_t *x, Double_t *par)
10{
11 // Fit parameters:
12 // par[0] = normalization factor
13 // par[1] = peak position
14 // par[2] = FWHM
15
16 return par[0] * TMath::BreitWigner(x[0], par[1], par[2]);
17}
18
19Double_t bwgaus(Double_t *x, Double_t *par)
20{
21 // Fit parameters:
22 // par[0] = global normalization constant
23 // par[1] = BW peak position
24 // par[2] = FWHM
25 // par[3] = sigma of convoluted gaussian
26
27 // Numeric constants
28 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
29
30 // Control constants
31 Double_t nsteps = 100.0; // number of convolution steps
32 Double_t nsigma = 3.0; // convolution extends to +-sc Gaussian sigmas
33
34 // Range of convolution integral
35 Double_t x1 = x[0] - nsigma * par[3];
36 Double_t x2 = x[0] + nsigma * par[3];
37 Double_t dx = (x2 - x1) / nsteps;
38
39 // Variables
40 Double_t i, xx, fbw, fgaus, sum = 0.0;
41
42 // Convolution integral of Breit-Wigner and Gaussian by sum
43 for(i = 1.0; i <= 0.5 * nsteps; i++) {
44
45 xx = x1 + (i - 0.5) * dx;
46 fbw = TMath::BreitWigner(xx, par[1], par[2]);
47 fgaus = TMath::Gaus(x[0], xx, par[3]);
48 sum += fbw * fgaus;
49
50 xx = x2 - (i - 0.5) * dx;
51 fbw = TMath::BreitWigner(xx, par[1], par[2]);
52 fgaus = TMath::Gaus(x[0], xx, par[3]);
53 sum += fbw * fgaus;
54
55 }
56
57 return (par[0] * dx * sum * invsq2pi / par[3]);
58}
59
60Double_t bwgaus_line(Double_t *x, Double_t *par)
61{
62 Double_t bwg = bwgaus(x, par);
63 return bwg + par[4] + par[5] * x[0];
64}
65
66void view
67(Int_t rebin = 8, Double_t mult = 4.0, Double_t fitmin = 0.7, Double_t fitmax = 0.8,
68 const char *filename = "kstar.invmass.root")
69{
70 TFile *f = TFile::Open(filename);
71
72 TH1D * hsign = (TH1D*)f->Get("h_Pi(+)_K(-)");
73 TH1D * hsign_2 = (TH1D*)f->Get("h_Pi(-)_K(+)");
74
75 TH1D * hmix = (TH1D*)f->Get("hmix_Pi(+)_K(-)");
76 TH1D * hmix_2 = (TH1D*)f->Get("hmix_Pi(-)_K(+)");
77 TH1D * hmix_3 = (TH1D*)f->Get("hmix_K(+)_Pi(-)");
78 TH1D * hmix_4 = (TH1D*)f->Get("hmix_K(-)_Pi(+)");
79
80 TH1D * htrue = (TH1D*)f->Get("h_Pi(+)_K(-)_true");
81 TH1D * htrue_2 = (TH1D*)f->Get("h_Pi(-)_K(+)_true");
82
83 hsign->Rebin(rebin);
84 hsign_2->Rebin(rebin);
85
86 htrue->Rebin(rebin);
87 htrue_2->Rebin(rebin);
88
89 hmix->Rebin(rebin);
90 hmix_2->Rebin(rebin);
91 hmix_3->Rebin(rebin);
92 hmix_4->Rebin(rebin);
93
94 hsign->Add(hsign_2);
95 htrue->Add(htrue_2);
96
97 hmix->Add(hmix_2);
98 hmix->Add(hmix_3);
99 hmix->Add(hmix_4);
100
101 hsign->Sumw2();
102 hmix->Sumw2();
103 htrue->Sumw2();
104
105 cout << "# signal entries: " << hsign->Integral() << endl;
106 cout << "# backgr. entries: " << hmix->Integral() << endl;
107 cout << "# true pairs: " << htrue->Integral() << endl;
108
109 // mean and width of lambda*
110 Double_t mean = 0.896;
111 Double_t gamma = 0.05;
112
113 // find bin limits for fit
114 Int_t b1 = hsign->GetXaxis()->FindFixBin(fitmin);
115 Int_t b2 = hsign->GetXaxis()->FindFixBin(fitmax);
116 hmix->Scale( hsign->Integral(b1, b2) / hmix->Integral(b1, b2) );
117
118 // drawing options
119 gStyle->SetOptStat(0);
120 gROOT->SetStyle("Plain");
121
122 // draw signal and background together
123 TCanvas *c1 = new TCanvas("c1", "", 0, 0, 640, 480);
124 TH1D *hsign1 = (TH1D*)hsign->Clone();
125 TH1D *hmix1 = (TH1D*)hmix->Clone();
126 hsign1->GetXaxis()->SetTitle("Inv. mass (GeV/c^{2})");
127 hsign1->GetYaxis()->SetTitle("counts");
128 hsign1->SetTitle("");
129 hsign1->SetMarkerStyle(8);
130 hsign1->SetMarkerSize(0.8);
131 hsign1->SetStats(0);
132 hsign1->Draw("PE");
133 //hmix1->SetLineColor(kRed);
134 hmix1->SetMarkerStyle(4);
135 hmix1->SetStats(0);
136 hmix1->GetXaxis()->SetTitle("Inv. mass (GeV/c^{2})");
137 hmix1->GetYaxis()->SetTitle("counts");
138 hmix1->Draw("CESAME");
139 c1->Update();
140
141 // draw subtraction and fit with BW+gaus
142 TCanvas *c2 = new TCanvas("c2", "", 50, 50, 640, 480);
143 TH1D *hdiff = (TH1D*)hsign1->Clone();
144 hdiff->Add(hmix, -1.0);
145 hdiff->GetXaxis()->SetRangeUser(mean - mult*gamma, mean + mult*gamma);
146// TF1 *fcn = new TF1("fcn", bwgaus_line, mean - mult*gamma, mean + mult*gamma, 6);
147 TF1 *fcn = new TF1("fcn", bwgaus, mean - mult*gamma, mean + mult*gamma, 4);
148 fcn->SetParameter(0, hdiff->GetMaximum() * 0.6);
149 fcn->SetParameter(1, mean);
150 fcn->SetParameter(2, gamma);
151 fcn->SetParameter(3, 0.001);
152// fcn->SetParameter(4, 20);
153// fcn->SetParameter(5, -50);
154// fcn->SetParNames("Constant", "BW_peak", "BW_gamma", "Gaus_sigma", "Line_A", "Line_B");
155 fcn->SetParNames("Constant", "BW_peak", "BW_gamma", "Gaus_sigma");
156 fcn->SetLineColor(kGreen);
157// hdiff->Fit(fcn, "RE");
158 hdiff->SetStats(0);
159 hdiff->GetXaxis()->SetTitle("Inv. mass (GeV/c^{2})");
160 hdiff->GetYaxis()->SetTitle("counts");
161 hdiff->Draw("PE");
162 cout << "Fit results: " << endl;
163 cout << "Peak center: " << fcn->GetParameter(1) << " +/- " << fcn->GetParError(1) << endl;
164 cout << "Peak width : " << fcn->GetParameter(2) << " +/- " << fcn->GetParError(2) << endl;
165 cout << "Gaus sigma : " << fcn->GetParameter(3) << " +/- " << fcn->GetParError(3) << endl;
166 c2->Update();
167
168 // draw true pairs
169 TCanvas *c3 = new TCanvas("c3", "", 100, 100, 640, 480);
170 TH1D *htrue1 = (TH1D*)htrue->Clone();
171 htrue1->SetMarkerStyle(21);
172 htrue1->SetMarkerColor(kRed);
173 htrue1->GetXaxis()->SetRangeUser(mean - mult*gamma, mean + mult*gamma);
174 htrue1->SetStats(0);
175 htrue1->GetXaxis()->SetTitle("Inv. mass (GeV/c^{2})");
176 htrue1->GetYaxis()->SetTitle("counts");
177 htrue1->Draw("PE");
178 c3->Update();
179
180 // draw true pairs with diff
181 TCanvas *c4 = new TCanvas("c4", "", 150, 150, 640, 480);
182 hdiff->Draw("PE");
183 htrue1->SetMarkerStyle(4);
184 htrue1->SetMarkerColor(1);
185 htrue1->Draw("PEsame");
186 c4->Update();
187
188 // print true and reconstructed pairs
189 b1 = htrue->GetXaxis()->FindFixBin(mean - 1.5*gamma);
190 b2 = htrue->GetXaxis()->FindFixBin(mean + 1.5*gamma);
191 TF1 *fbw = new TF1("fbw", bw, mean - mult*gamma, mean + mult*gamma, 3);
192 fbw->SetParameter(0, fcn->GetParameter(0));
193 fbw->SetParameter(1, fcn->GetParameter(1));
194 fbw->SetParameter(2, fcn->GetParameter(2));
195 Double_t S = htrue->Integral(b1, b2);
196 Double_t B = hsign->Integral(b1, b2);
197 cout << "True pairs: " << S << endl;
198 cout << "Reconstructed: " << fbw->Integral(mean - 1.5*gamma, mean + 1.5*gamma) / htrue->GetBinWidth(1) << endl;
199 cout << "Background pairs: " << B << endl;
200 cout << "S/B (+/- 1.5 G): " << S / B << endl;
201 cout << "sign. (+/- 1.5 G): " << S / TMath::Sqrt(S + B) << endl;
202
203 // save pictures
204 Text_t answer;
205 cout << "Save canvases as eps files (y/n)? ";
206 cin >> answer;
207 if (answer == 'y' || answer == 'Y') {
208 c1->SaveAs("lambda_sign_bg.eps");
209 c2->SaveAs("lambda_diff.eps");
210 c3->SaveAs("lambda_true.eps");
211 c4->SaveAs("lambda_comparison.eps");
212 }
213}
214