94924fca07b6443f78e4bc28963aa52940aee0d5
[u/mrichter/AliRoot.git] / FMD / scripts / DrawESD.C
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  */
43 class DrawESD : public AliFMDInput
44 {
45 private:
46   TH1D* fMult; // Histogram 
47   const Double_t fCorr;
48 public:
49   //__________________________________________________________________
50   DrawESD(Int_t n=1000, Double_t mmin=-0.5, Double_t mmax=20.5) 
51     : fCorr(1) // 0.68377 / 1.1)
52   { 
53     AddLoad(kESD);
54     fMult = new TH1D("mult", " Multiplicity (strip)", n, mmin, mmax);
55     fMult->Sumw2();
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   //__________________________________________________________________
67   Bool_t ProcessESD(UShort_t det, 
68                     Char_t   rng, 
69                     UShort_t sec, 
70                     UShort_t str, 
71                     Float_t  /* eta */, 
72                     Float_t  mult)
73   {
74     // Cache the energy loss 
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);
78     return kTRUE;
79   }
80   //__________________________________________________________________
81   TF1* FitPeak(Int_t n, TH1D* hist, Double_t min, Double_t& max)
82   {
83     if (TMath::Abs(max-min) < .25) return 0;
84     std::cout << "Fit peack in range " << min << " to " << max << std::endl;
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   //__________________________________________________________________
117   Bool_t Finish()
118   {
119     Info("Finish", "Will draw results");
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     
128     if (fMult->GetEntries() <= 0) return kFALSE;
129     
130     TCanvas* c = new TCanvas("c", "C");
131     c->cd();
132     // c->SetLogy();
133     fMult->GetXaxis()->SetRangeUser(0,4);
134     fMult->Scale(1. / fMult->GetEntries());
135     fMult->SetStats(kFALSE);
136     fMult->SetFillColor(2);
137     fMult->SetFillStyle(3001);
138     fMult->Draw("hist e");
139     
140     return kTRUE;
141     
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);
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     
176     fMult->GetXaxis()->SetRangeUser(0, 4);
177     fMult->Fit(f, "0QR", "", x1, x3);
178     fMult->Fit(f, "ME0R", "E1", x1, x3);
179     fMult->DrawClone("same hist");
180
181     l1->SetLineColor(3);
182     l1->SetLineWidth(2);
183     l1->SetRange(0, 4);
184     l1->Draw("same");
185     if (l2) {
186       l2->SetLineColor(4);
187       l2->SetLineWidth(2);
188       l2->SetRange(0, 4);
189       l2->Draw("same");
190     }
191     f->SetLineWidth(2);
192     f->SetRange(0, 4);
193     f->Draw("same");
194
195     TLegend* l = new TLegend(0.6, 0.6, .89, .89);
196     l->AddEntry(l1, "1 particle Landau", "l");
197     if (l2) l->AddEntry(l2, "2 particle Landau", "l");
198     l->AddEntry(f,  "1+2 particle Landau", "l");
199     l->SetFillColor(0);
200     l->Draw("same");
201
202
203 #if 0
204     c = new TCanvas("c2", "Landaus");
205     // c->SetLogy();
206     fMult->DrawClone("axis");
207     f->Draw("same");
208     Double_t* p1 = f->GetParameters();
209     Double_t* p2 = &(p1[3]);
210     TF1* ll1 = new TF1("ll1", "landau", 0, 4);
211     ll1->SetParameters(p1);
212     ll1->SetLineColor(3);
213     ll1->Draw("same");
214     TF1* ll2 = new TF1("ll2", "landau", 0, 4);
215     ll2->SetParameters(p2);
216     ll2->SetLineColor(4);
217     ll2->Draw("same");
218
219     Double_t y1  = fMult->GetMinimum() * 1.1;
220     Double_t y2  = fMult->GetMaximum() * .9;
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;
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");
228     TLine* c3 = new TLine(xc3, y1, xc3, y2);
229     c3->Draw("same");
230
231     l = new TLegend(0.6, 0.6, .89, .89);
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");
237 #endif
238
239
240     return kTRUE;
241   }
242
243   ClassDef(DrawESD,0);
244   
245 };
246
247 //____________________________________________________________________
248 //
249 // EOF
250 //