Introducing a new AOD class: AliAODcascade (A.Maire)
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / macros / AliRsnView.C
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
9 Double_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
19 Double_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
60 Double_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
66 void 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