Fixes for reading zero-suppressed data. These should be propagated to
[u/mrichter/AliRoot.git] / FMD / scripts / AliFMDAnaFlow.C
1 #include <TH1D.h>
2 #include <TH2D.h>
3 #include <TBrowser.h>
4 #include <TArrayF.h>
5 #include <TCanvas.h>
6 #include <TStyle.h>
7 #include <TLegend.h>
8 #include <TProfile.h>
9 #include <TProfile2D.h>
10 #include <TROOT.h>
11 #include "AliFMDAnaESD.h"
12 #include "AliFMDAnaFlowRing.h"
13 #include <TFile.h>
14 #if 1
15 # include "background.C"
16 // # include "background2.C"
17 #else
18 TH2* bgFMD1i() { return 0; }
19 TH2* bgFMD2i() { return 0; }
20 TH2* bgFMD2o() { return 0; }
21 TH2* bgFMD3o() { return 0; }
22 TH2* bgFMD3o() { return 0; }
23 #endif
24 #if 1
25 # include "reactionplane.C"
26 #else
27 Float_t reactionplane(int ev) { return 0; }
28 #endif
29 #include <iostream>
30
31
32 //====================================================================
33 class AliFMDAnaFlow : public AliFMDAnaESD
34 {
35 public:
36   //__________________________________________________________________
37   /** Constructor */
38   AliFMDAnaFlow(Int_t n=2, Bool_t bg=false, Float_t cut0=.68) 
39     : fAxis(n, -3.7, 5.3),
40       fV2("all", "Full FMD", 2, 1, fAxis)
41   {
42     for (int i = 0; i < 5; i++) { 
43       UShort_t det = (i == 0 ? 1 : (i <= 2 ? 2 : 3));
44       Char_t   rng = (i == 0 ? 'I' :  (i % 2 == 1 ? 'I' : 'O'));
45       AddRing(new AliFMDAnaFlowRing(det,rng,n,EtaMin(i)+.1,EtaMax(i)-.1,
46                            Bg(bg,i),cut0,Cut1(i)));
47     }
48   }
49   //__________________________________________________________________
50   /** Constructor */
51   AliFMDAnaFlow(Float_t deta, Bool_t bg=false, Float_t cut0=.68) 
52     : fAxis(Int_t((5.3- -3.7)/deta), -3.7, 5.3),
53       fV2("all", "Full FMD", 2, 1, fAxis)
54   {
55     for (int i = 0; i < 5; i++) { 
56       UShort_t det = (i == 0 ? 1 : i <= 2 ? 2 : 3);
57       Char_t   rng = (i % 2 == 0 ? 'I' : 'O');
58       AddRing(new AliFMDAnaFlowRing(det,rng,fAxis,Bg(bg,i),cut0,Cut1(i)));
59     }
60   }
61   //__________________________________________________________________
62   /** Initialize */
63   Bool_t Init()
64   {
65     fPhis.Set(51200);
66     fEtas.Set(51200);
67     fWeights.Set(51200);
68     fNObserv = 0;
69
70     Float_t pi2 = 2 * TMath::Pi();
71
72     fPsi  = new TH2D("all_psi", "Event planes", 
73                      fAxis.N(), fAxis.Bins(), 40, 0, pi2);
74     fPsi->SetXTitle("#eta");
75     fPsi->SetYTitle("#Psi");
76     fPsi->SetZTitle("Frequency");
77     fPsi->SetMarkerStyle(20);
78
79     fDPsi = new TProfile("all_dpsi", "Distance to real #Psi_{R}", 
80                          fAxis.N(), fAxis.Bins(), -pi2, pi2);
81     fDPsi->SetXTitle("#eta");
82     fDPsi->SetYTitle("#Delta#Psi=|#Psi-#Psi_{R}|");
83
84     fDPhi = new TProfile2D("all_dphi","Distribution of phi signals", 
85                           fAxis.N(), -3.7, 5.3, 40, 0, pi2, 0, 20);
86     fDPhi->SetXTitle("#eta");
87     fDPhi->SetYTitle("#varphi");
88     fDPhi->SetZTitle("signal");
89     fDPhi->SetMarkerStyle(20);
90
91     return AliFMDAnaESD::Init();
92   }
93   //__________________________________________________________________
94   /** Called at the beginning of an event */
95   Bool_t Begin(Int_t ev) 
96   {
97     fNObserv = 0;
98     Float_t psir = reactionplane(fNEvents);
99     AliFMDAnaFlowRing* fring = 0;
100     for (size_t i = 0; i < 5; i++) {
101       if (!fRing[i]) continue;
102       fring = static_cast<AliFMDAnaFlowRing*>(fRing[i]);
103       fring->fTrue.SetPsi(psir);
104     }
105     return AliFMDAnaESD::Begin(ev);
106   }
107   //__________________________________________________________________
108   /** Fill in one strip */
109   void Fill(Float_t phi, Float_t eta, Float_t mult)
110   {
111     Int_t imult = Int_t(mult+.33);
112     if (fNObserv + imult >= fEtas.fN) { 
113       ULong_t n = Int_t(fEtas.fN * 1.5)+imult;
114       fEtas.Set(n);
115       fPhis.Set(n);
116       fWeights.Set(n);
117     }
118     for (Int_t i = 0; i < imult; i++) { 
119       fEtas[fNObserv]    = eta;
120       fPhis[fNObserv]    = phi;
121       fWeights[fNObserv] = 1;
122       fNObserv++;
123     }
124     fDPhi->Fill(eta, phi);
125   }
126   //__________________________________________________________________
127   /** Called at the end of an event */
128   Bool_t End() 
129   {
130     float       psir    = reactionplane(fNEvents-1);
131     // std::cout << "PsiR=" << psir << std::endl;
132     fV2.Event(fNObserv, fPhis.fArray, fEtas.fArray, 
133               fWeights.fArray, fWeights.fArray);
134     for (UShort_t i = 0; i < fAxis.N(); i++) { 
135       AliFMDFlowBin* bin = fV2.GetBin(i);
136       if (bin && bin->Counts() > 0) {
137         fPsi->Fill(fAxis.BinCenter(i),  bin->Psi() - psir);
138         fDPsi->Fill(fAxis.BinCenter(i), bin->Psi() - psir);
139       }
140     }
141     return AliFMDAnaESD::End();
142   }
143   //__________________________________________________________________
144   /** Make a frame */
145   TH1* MakeFrame(const char* xname, const char* yname, 
146                  Float_t ymin, Float_t ymax)
147   {
148     static int id = 0;
149     id++;
150     TH1* frame = new TH1D(Form("frame%d", id), "Frame", 
151                           fAxis.N(), fAxis.Bins()[0], fAxis.Bins()[fAxis.N()]);
152     frame->SetMinimum(ymin);
153     frame->SetMaximum(ymax);
154     frame->SetXTitle(xname);
155     frame->SetYTitle(yname);
156     frame->SetLabelFont(132, "XYZ");
157     frame->SetTitleFont(132, "XYZ");
158     frame->SetLabelSize(.08, "XYZ");
159     frame->SetTitleSize(.08, "XYZ");
160     frame->SetTitleOffset(.5, "XYZ");
161     frame->SetNdivisions(8, "XYZ");
162     frame->Draw();
163     return frame;
164   }
165   //__________________________________________________________________
166   /** called at the end of a run */
167   Bool_t Finish()
168   {
169     gStyle->SetPalette(1);
170     gStyle->SetOptTitle(0);
171     gStyle->SetOptStat(0);
172     AliFMDAnaFlowRing* fring[5];
173     for (size_t i = 0; i < 5; i++) 
174       fring[i] = static_cast<AliFMDAnaFlowRing*>(fRing[i]);
175     
176     TCanvas* c = new TCanvas("c1", "Result", 600, 750);
177     c->cd();
178     c->SetGridy();
179     c->SetTopMargin(0.02);
180     c->SetRightMargin(0.02);
181     c->SetFillColor(0);
182     c->SetBorderSize(0);
183     c->SetBorderMode(0);
184     c->SetHighLightColor(0);
185     c->Divide(1,3,0,0);
186
187     TVirtualPad* p = 0;
188     TH1*         f = 0;
189
190     p = c->cd(1);
191     p->SetGridx();
192     p->SetGridy();
193     p->SetFillColor(0);
194     f = MakeFrame("#eta", "|#Psi_{2}-#Psi_{R}|", -TMath::Pi(), TMath::Pi());
195     for (int i = 0; i < 5; i++) fring[i]->fDPsi.Draw("same");
196     fDPsi->Draw("same");
197
198     p = c->cd(2);
199     p->SetGridx();
200     p->SetGridy();
201     p->SetFillColor(0);
202     f = MakeFrame("#eta", "R_{k}", .7, 1);
203     for (int i = 0; i < 5; i++) fring[i]->fV2.Draw("tr:same");
204     fV2.Draw("tr:same");
205
206     p = c->cd(3);
207     p->SetGridx();
208     p->SetGridy();
209     p->SetFillColor(0);
210     f = MakeFrame("#eta", "v_{2}", 0, .1);
211     for (int i = 0; i < 5; i++) { 
212       fring[i]->fV2.Draw("th:same");
213       fring[i]->fTrue.Draw("th:same");
214     }
215     fV2.Draw("th:same");
216
217
218     for (int i = 0; i < 5; i++) { 
219       fring[i]->fV2.Print("s");
220       fring[i]->fTrue.Print("s");
221     }
222     fV2.Print("s");
223     gROOT->GetListOfSpecials()->Add(this);
224     return kTRUE;
225   }
226   void ToFile(const char* file) { 
227     TFile* out = TFile::Open(file, "RECREATE");
228     for (int i = 0; i < 5; i++) fRing[i]->Write(fRing[i]->Name());
229     fV2.Write();
230     if (fPsi) fPsi->Write();
231     if (fDPsi) fDPsi->Write();
232     if (fDPhi) fDPhi->Write();
233     out->Close();
234   }
235     
236     
237
238   //__________________________________________________________________
239   /** Get min pseudo rapidity of ring idx */
240   Float_t EtaMin(Int_t idx) const { 
241     switch (idx) { 
242     case 0: return  3.61728; break; // FMD1I
243     case 1: return  2.18122; break; // FMD2I
244     case 2: return  1.79821; break; // FMD2O
245     case 3: return -3.39913; break; // FMD3I
246     case 4: return -2.28925; break; // FMD3O
247     }
248     return 0;
249   }
250   //__________________________________________________________________
251   /** Get max pseudo rapidity of ring idx */
252   Float_t EtaMax(Int_t idx) const { 
253     switch (idx) { 
254     case 0: return  5.02643; break; // FMD1I
255     case 1: return  3.57899; break; // FMD2I
256     case 2: return  2.39085; break; // FMD2O
257     case 3: return -2.00644; break; // FMD3I
258     case 4: return -1.70080; break; // FMD3O
259     }
260     return 0;
261   }
262   //__________________________________________________________________
263   /** Get cut 1 of ring idx */
264   Float_t Cut1(Int_t idx) const { 
265     switch (idx) { 
266     case 0: return 1.33857; break; // FMD1I
267     case 1: return 1.30464; break; // FMD2I
268     case 2: return 1.27263; break; // FMD2O
269     case 3: return 1.28052; break; // FMD3I
270     case 4: return 1.23041; break; // FMD3O
271     }
272     return 0;
273   }
274   //__________________________________________________________________
275   /** Get background correction of ring idx */
276   TH2* Bg(Bool_t use, Int_t idx) const { 
277     if (!use) return 0;
278     switch (idx) { 
279     case 0: return bgFMD1i(); break;
280     case 1: return bgFMD2i(); break;
281     case 2: return bgFMD2o(); break;
282     case 3: return bgFMD3i(); break;
283     case 4: return bgFMD3o(); break;
284     }
285     return 0;
286   }
287   /** Flow axis */
288   AliFMDFlowAxis fAxis;
289   /** Flow object */
290   AliFMDFlowBinned1D fV2;
291   /** Cache of phi */
292   TArrayD            fPhis;
293   /** Cache of eta */
294   TArrayD            fEtas;
295   /** Cache of weights */
296   TArrayD            fWeights;
297   /** Number of observables */
298   Int_t              fNObserv;
299   /** Histogram of Psi */
300   TH2D*              fPsi;       // Event planes
301   /** Histogram of dPsi */
302   TProfile*          fDPsi;   // Mean distance to Psi
303   /** Histogram of dPhi */
304   TProfile2D*        fDPhi;
305   ClassDef(AliFMDAnaFlow,1) 
306 };
307
308 //____________________________________________________________________
309 //
310 // EOF
311 //