Try to fix Coverity defect
[u/mrichter/AliRoot.git] / FMD / scripts / DrawESD.C
CommitLineData
2b893216 1//____________________________________________________________________
2//
3// $Id$
4//
5// Script that contains a class to draw eloss from hits, versus ADC
6// counts from digits, using the AliFMDInputHits class in the util library.
7//
8// It draws the energy loss versus the p/(mq^2). It can be overlayed
9// with the Bethe-Bloc curve to show how the simulation behaves
10// relative to the expected.
11//
12// Use the script `Compile.C' to compile this class using ACLic.
13//
14#include <TH1D.h>
15#include <AliFMDHit.h>
16#include <AliFMDDigit.h>
17#include <AliFMDInput.h>
18#include <AliFMDUShortMap.h>
19#include <AliFMDFloatMap.h>
20#include <AliFMDRecPoint.h>
21#include <AliESDFMD.h>
22#include <AliLog.h>
23#include <iostream>
24#include <TStyle.h>
25#include <TArrayF.h>
26#include <TCanvas.h>
27#include <TMath.h>
28#include <TF1.h>
29#include <TSpectrum.h>
30#include <TLegend.h>
31#include <TLine.h>
32
33/** @class DrawESD
34 @brief Draw digit ADC versus Rec point mult
35 @code
36 Root> .L Compile.C
37 Root> Compile("DrawESD.C")
38 Root> DrawESD c
39 Root> c.Run();
40 @endcode
41 @ingroup FMD_script
42 */
43class DrawESD : public AliFMDInput
44{
45private:
46 TH1D* fMult; // Histogram
47 const Double_t fCorr;
48public:
49 //__________________________________________________________________
f48d9b11 50 DrawESD(Int_t n=1000, Double_t mmin=-0.5, Double_t mmax=20.5)
2b893216 51 : fCorr(1) // 0.68377 / 1.1)
52 {
53 AddLoad(kESD);
54 fMult = new TH1D("mult", " Multiplicity (strip)", n, mmin, mmax);
69893a66 55 fMult->Sumw2();
2b893216 56 fMult->SetXTitle("Strip Multiplicity");
57 }
58 //__________________________________________________________________
59 /** Begining of event
60 @param ev Event number
61 @return @c false on error */
62 Bool_t Begin(Int_t ev)
63 {
64 return AliFMDInput::Begin(ev);
65 }
66 //__________________________________________________________________
faf80567 67 Bool_t ProcessESD(UShort_t det,
68 Char_t rng,
69 UShort_t sec,
70 UShort_t str,
2b893216 71 Float_t /* eta */,
72 Float_t mult)
73 {
74 // Cache the energy loss
faf80567 75 if (mult > 0)
76 Info("ProcessESD", "FMD%d%c[%2d,%3d]=%f", det, rng, sec, str, mult);
77 if (mult/fCorr > 0.001) fMult->Fill(mult/fCorr);
2b893216 78 return kTRUE;
79 }
80 //__________________________________________________________________
f48d9b11 81 TF1* FitPeak(Int_t n, TH1D* hist, Double_t min, Double_t& max)
82 {
faf80567 83 if (TMath::Abs(max-min) < .25) return 0;
84 std::cout << "Fit peack in range " << min << " to " << max << std::endl;
f48d9b11 85 TF1* l = new TF1(Form("l%02d", n), "landau", min, max);
86 hist->GetXaxis()->SetRangeUser(0, 4);
87 hist->Fit(l, "0Q", "", min, max);
88 Double_t mpv = l->GetParameter(1);
89 Double_t empv = l->GetParError(1);
90 Double_t sigma = l->GetParameter(2);
91 l->SetRange(mpv-empv, mpv+3*sigma);
92 hist->Fit(l, "EMQ0", "", mpv-3*empv, mpv+3*sigma);
93 std::cout << "Peak # " << n << " [" << min << "," << max << "]\n"
94 << " MPV: " << l->GetParameter(1)
95 << " +/- " << l->GetParError(1)
96 << " Var: " << l->GetParameter(2)
97 << " +/- " << l->GetParError(2)
98 << " Chi^2/NDF: " << l->GetChisquare() / l->GetNDF()
99 << std::endl;
100 mpv = l->GetParameter(1);
101 sigma = l->GetParameter(2);
102 min = mpv - sigma * 2; // (n==1 ? 3 : 2);
103 max = mpv + sigma * 3;
104 // l->SetRange(min, max);
105 l->Draw("same");
106 return l;
107 }
108 //__________________________________________________________________
109 void MaxInRange(TH1D* hist, Double_t min, Double_t& mean, Double_t& var)
110 {
111 hist->GetXaxis()->SetRangeUser(min, 4);
112 mean = hist->GetMean();
113 var = hist->GetRMS();
114 }
115
116 //__________________________________________________________________
2b893216 117 Bool_t Finish()
118 {
faf80567 119 Info("Finish", "Will draw results");
2b893216 120 gStyle->SetPalette(1);
121 gStyle->SetOptTitle(0);
122 gStyle->SetOptFit(1111111);
123 gStyle->SetCanvasColor(0);
124 gStyle->SetCanvasBorderSize(0);
125 gStyle->SetPadColor(0);
126 gStyle->SetPadBorderSize(0);
127
faf80567 128 if (fMult->GetEntries() <= 0) return kFALSE;
129
2b893216 130 TCanvas* c = new TCanvas("c", "C");
131 c->cd();
faf80567 132 // c->SetLogy();
f48d9b11 133 fMult->GetXaxis()->SetRangeUser(0,4);
69893a66 134 fMult->Scale(1. / fMult->GetEntries());
2b893216 135 fMult->SetStats(kFALSE);
136 fMult->SetFillColor(2);
137 fMult->SetFillStyle(3001);
f48d9b11 138 fMult->Draw("hist e");
faf80567 139
140 return kTRUE;
141
f48d9b11 142 Double_t mean, rms;
143 MaxInRange(fMult, 0.2, mean, rms);
144 Double_t x1 = mean-rms/2; // .75; // .8; // .65 / fCorr;
145 Double_t x2 = mean+rms/2; // 1.3; // 1.7; // fCorr;
146 TF1* l1 = FitPeak(1, fMult, x1, x2);
147 x2 = TMath::Max(mean+rms/2, x2);
148 MaxInRange(fMult, x2, mean, rms);
149 Double_t x3 = mean + rms;
150 TF1* l2 = FitPeak(2, fMult, x2, x3);
faf80567 151 TF1* f = 0;
152 Double_t diff = 0;
153 if (l2) {
154 diff = l2->GetParameter(1)-l1->GetParameter(1);
155 f = new TF1("user", "landau(0)+landau(3)", x1, x3);
156 f->SetParNames("A_{1}", "Mpv_{1}", "#sigma_{1}",
157 "A_{2}", "Mpv_{2}", "#sigma_{2}",
158 "A_{3}", "Mpv_{3}", "#sigma_{3}");
159 f->SetParameters(l1->GetParameter(0),
160 l1->GetParameter(1),
161 l1->GetParameter(2),
162 l2 ? l2->GetParameter(0) : 0,
163 l2 ? l2->GetParameter(1) : 0,
164 l2 ? l2->GetParameter(2) : 0,
165 l2->GetParameter(0)/10,
166 l2->GetParameter(1) + diff,
167 l2->GetParameter(2));
168 }
169 else {
170 x3 = x2;
171 f = new TF1("usr", "landau", x1, x3);
172 }
173
174 std::cout << "Range: " << x1 << "-" << x3 << std::endl;
175
2b893216 176 fMult->GetXaxis()->SetRangeUser(0, 4);
faf80567 177 fMult->Fit(f, "0QR", "", x1, x3);
178 fMult->Fit(f, "ME0R", "E1", x1, x3);
f48d9b11 179 fMult->DrawClone("same hist");
2b893216 180
2b893216 181 l1->SetLineColor(3);
182 l1->SetLineWidth(2);
183 l1->SetRange(0, 4);
184 l1->Draw("same");
faf80567 185 if (l2) {
186 l2->SetLineColor(4);
187 l2->SetLineWidth(2);
188 l2->SetRange(0, 4);
189 l2->Draw("same");
190 }
2b893216 191 f->SetLineWidth(2);
192 f->SetRange(0, 4);
193 f->Draw("same");
194
f48d9b11 195 TLegend* l = new TLegend(0.6, 0.6, .89, .89);
2b893216 196 l->AddEntry(l1, "1 particle Landau", "l");
faf80567 197 if (l2) l->AddEntry(l2, "2 particle Landau", "l");
2b893216 198 l->AddEntry(f, "1+2 particle Landau", "l");
199 l->SetFillColor(0);
200 l->Draw("same");
201
202
faf80567 203#if 0
2b893216 204 c = new TCanvas("c2", "Landaus");
faf80567 205 // c->SetLogy();
2b893216 206 fMult->DrawClone("axis");
207 f->Draw("same");
f48d9b11 208 Double_t* p1 = f->GetParameters();
209 Double_t* p2 = &(p1[3]);
2b893216 210 TF1* ll1 = new TF1("ll1", "landau", 0, 4);
f48d9b11 211 ll1->SetParameters(p1);
2b893216 212 ll1->SetLineColor(3);
213 ll1->Draw("same");
214 TF1* ll2 = new TF1("ll2", "landau", 0, 4);
f48d9b11 215 ll2->SetParameters(p2);
2b893216 216 ll2->SetLineColor(4);
217 ll2->Draw("same");
218
219 Double_t y1 = fMult->GetMinimum() * 1.1;
220 Double_t y2 = fMult->GetMaximum() * .9;
f48d9b11 221 Double_t xc1 = p1[1]-3*p1[2];
222 Double_t xc2 = p2[1]-2*p2[2];
223 Double_t xc3 = p2[1]-2*p2[2]+diff;
2b893216 224 TLine* c1 = new TLine(xc1, y1, xc1, y2);
225 c1->Draw("same");
226 TLine* c2 = new TLine(xc2, y1, xc2, y2);
227 c2->Draw("same");
f48d9b11 228 TLine* c3 = new TLine(xc3, y1, xc3, y2);
229 c3->Draw("same");
2b893216 230
f48d9b11 231 l = new TLegend(0.6, 0.6, .89, .89);
2b893216 232 l->AddEntry(ll1, "1 particle Landau", "l");
233 l->AddEntry(ll2, "2 particle Landau", "l");
234 l->AddEntry(f, "1+2 particle Landau", "l");
235 l->SetFillColor(0);
236 l->Draw("same");
2b893216 237#endif
238
f48d9b11 239
2b893216 240 return kTRUE;
241 }
242
243 ClassDef(DrawESD,0);
244
245};
246
247//____________________________________________________________________
248//
249// EOF
250//