A lot of changes after detector review:
[u/mrichter/AliRoot.git] / FMD / scripts / AliFMDAnaFlow.C
CommitLineData
9b98d361 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
18TH2* bgFMD1i() { return 0; }
19TH2* bgFMD2i() { return 0; }
20TH2* bgFMD2o() { return 0; }
21TH2* bgFMD3o() { return 0; }
22TH2* bgFMD3o() { return 0; }
23#endif
24#if 1
25# include "reactionplane.C"
26#else
27Float_t reactionplane(int ev) { return 0; }
28#endif
29#include <iostream>
30
31
32//====================================================================
33class AliFMDAnaFlow : public AliFMDAnaESD
34{
35public:
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//