Bug fix for reading raw data - thanks Frederic YERMIA, Diego, and Gines MARTINEZ
[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   TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max) 
51   {
52     TArrayF bins(n+1);
53     Float_t dp   = n / TMath::Log10(max / min);
54     Float_t pmin = TMath::Log10(min);
55     bins[0]      = min;
56     for (Int_t i = 1; i < n+1; i++) {
57       Float_t p = pmin + i / dp;
58       bins[i]   = TMath::Power(10, p);
59     }
60     return bins;
61   }
62   //__________________________________________________________________
63   DrawESD(Int_t n=420, Double_t mmin=-0.5, Double_t mmax=20.5) 
64     : fCorr(1) // 0.68377 / 1.1)
65   { 
66     AddLoad(kESD);
67     fMult = new TH1D("mult", " Multiplicity (strip)", n, mmin, mmax);
68     fMult->SetXTitle("Strip Multiplicity");
69   }
70   //__________________________________________________________________
71   /** Begining of event
72       @param ev Event number
73       @return @c false on error */
74   Bool_t Begin(Int_t ev) 
75   {
76     return AliFMDInput::Begin(ev);
77   }
78   //__________________________________________________________________
79   Bool_t ProcessESD(UShort_t /* det */, 
80                     Char_t   /* ring */, 
81                     UShort_t /* sec */, 
82                     UShort_t /* strip */, 
83                     Float_t  /* eta */, 
84                     Float_t  mult)
85   {
86     // Cache the energy loss 
87     fMult->Fill(mult/fCorr);
88     return kTRUE;
89   }
90   //__________________________________________________________________
91   Bool_t Finish()
92   {
93     gStyle->SetPalette(1);
94     gStyle->SetOptTitle(0);
95     gStyle->SetOptFit(1111111);
96     gStyle->SetCanvasColor(0);
97     gStyle->SetCanvasBorderSize(0);
98     gStyle->SetPadColor(0);
99     gStyle->SetPadBorderSize(0);
100     
101     TCanvas* c = new TCanvas("c", "C");
102     c->cd();
103     c->SetLogy();
104     // fMult->GetXaxis()->SetRangeUser(0.2,4);
105     // fMult->Sumw2();
106     // fMult->Scale(1. / fMult->GetMaximum());
107     fMult->SetStats(kFALSE);
108     fMult->SetFillColor(2);
109     fMult->SetFillStyle(3001);
110     fMult->Draw();
111     Double_t x1 = .75; // .8;  // .65 / fCorr;
112     Double_t x2 = 1.3; // 1.7; // fCorr;
113     Double_t x3 = 2.5; // 2.7; // fCorr;
114     Double_t x4 = 3.7; // fCorr;
115     TF1* l1 = new TF1("landau1", "landau", x1, x2);
116     TF1* l2 = new TF1("landau2", "landau", x2, x3);
117     // TF1* l3 = new TF1("landau3", "landau", x3, x4);
118     TF1* f  = new TF1("user", "landau(0)+landau(3)", x1, x3);
119
120     fMult->SetStats(kTRUE);
121     fMult->GetXaxis()->SetRangeUser(0, 4);
122     fMult->Fit(l1, "E0L", "", x1, x2);
123     fMult->Fit(l2, "E0L", "", x2, x3);
124     // fMult->Fit(l3, "E0L", "", x3, x4);
125     f->SetParNames("A_{1}", "Mpv_{1}", "#sigma_{1}",
126                    "A_{2}", "Mpv_{2}", "#sigma_{2}",
127                    "A_{3}", "Mpv_{3}", "#sigma_{3}");
128     f->SetParameters(l1->GetParameter(0), 
129                      l1->GetParameter(1), 
130                      l1->GetParameter(2), 
131                      l2->GetParameter(0), 
132                      l2->GetParameter(1), 
133                      l2->GetParameter(2),
134                      l1->GetParameter(0), 
135                      l1->GetParameter(1) * 3, 
136                      l1->GetParameter(2) * 3);
137     fMult->Fit(f, "E0L", "", x1, x3);
138     fMult->Fit(f, "MEL", "E1", x1, x3);
139
140     TH1* h = static_cast<TH1*>(fMult->DrawClone("same hist"));
141
142     // l1->SetLineStyle(2);
143     l1->SetLineColor(3);
144     l1->SetLineWidth(2);
145     l1->SetRange(0, 4);
146     l1->Draw("same");
147     // l2->SetLineStyle(3);
148     l2->SetLineColor(4);
149     l2->SetLineWidth(2);
150     l2->SetRange(0, 4);
151     l2->Draw("same");
152     // l3->SetLineStyle(3);
153     // l3->SetLineWidth(2);
154     // l3->SetRange(0, 5);
155     // l3->Draw("same");
156     f->SetLineWidth(2);
157     f->SetRange(0, 4);
158     f->Draw("same");
159
160     TLegend* l = new TLegend(0.2, 0.6, .6, .89);
161     l->AddEntry(l1, "1 particle Landau", "l");
162     l->AddEntry(l2, "2 particle Landau", "l");
163     l->AddEntry(f,  "1+2 particle Landau", "l");
164     l->SetFillColor(0);
165     l->Draw("same");
166
167
168     c = new TCanvas("c2", "Landaus");
169     c->SetLogy();
170     fMult->DrawClone("axis");
171     f->Draw("same");
172     TF1* ll1 = new TF1("ll1", "landau", 0, 4);
173     ll1->SetParameters(f->GetParameter(0), 
174                        f->GetParameter(1), 
175                        f->GetParameter(2));
176     ll1->SetLineColor(3);
177     ll1->Draw("same");
178     TF1* ll2 = new TF1("ll2", "landau", 0, 4);
179     ll2->SetParameters(f->GetParameter(3), 
180                        f->GetParameter(4), 
181                        f->GetParameter(5));
182     ll2->SetLineColor(4);
183     ll2->Draw("same");
184
185     Double_t y1  = fMult->GetMinimum() * 1.1;
186     Double_t y2  = fMult->GetMaximum() * .9;
187     Double_t xc1 = f->GetParameter(1)-3*f->GetParameter(2);
188     Double_t xc2 = f->GetParameter(4)-2*f->GetParameter(5);
189     TLine* c1 = new TLine(xc1, y1, xc1, y2);
190     c1->Draw("same");
191     TLine* c2 = new TLine(xc2, y1, xc2, y2);
192     c2->Draw("same");
193
194     l = new TLegend(0.2, 0.6, .6, .89);
195     l->AddEntry(ll1, "1 particle Landau", "l");
196     l->AddEntry(ll2, "2 particle Landau", "l");
197     l->AddEntry(f,  "1+2 particle Landau", "l");
198     l->SetFillColor(0);
199     l->Draw("same");
200
201 #if 0
202     c = new TCanvas("s", "Spectrum");
203     TSpectrum* s = new TSpectrum(16);
204     h->GetXaxis()->SetRangeUser(0.3, 20);
205     s->Search(h, .15, "", 0.01);
206     c->Update();
207     TH1* b = s->Background(h);
208     b->SetFillColor(4);
209     b->SetFillStyle(3001);
210     b->Draw("same");
211 #endif
212
213     return kTRUE;
214   }
215
216   ClassDef(DrawESD,0);
217   
218 };
219
220 //____________________________________________________________________
221 //
222 // EOF
223 //