bug fix
[u/mrichter/AliRoot.git] / FMD / scripts / AliFMDAnaFlowKine.C
CommitLineData
9b98d361 1#include <AliFMDHit.h>
2#include <TH1D.h>
3#include <TH2D.h>
4#include <TParticle.h>
5#include <TBrowser.h>
6#include <TArrayF.h>
7#include <TCanvas.h>
8#include <TStyle.h>
9#include <TLegend.h>
10#include <TProfile.h>
11#include <TProfile2D.h>
12#include <TROOT.h>
13#include <AliFMDAnaESD.h>
14#include "AliFMDAnaFlowRing.h"
15#if 1
16# include "background.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 AliFMDAnaFlowKine : public AliFMDInput
34{
35public:
36 //__________________________________________________________________
37 /** Constructor */
38 AliFMDAnaFlowKine(Int_t n=1, Bool_t segments=true,
39 Bool_t primary_only=true, Float_t cut0=.68)
40 : fAxis(n, -3.7, 5.3),
41 fV2("all", "Full FMD", 2, fAxis),
42 fPrimaryOnly(primary_only),
43 fSegmented(segments)
44 {
45 AddLoad(kHits);
46 AddLoad(kKinematics);
47 for (int i = 0; i < 5; i++) {
48 UShort_t det = (i == 0 ? 1 : (i <= 2 ? 2 : 3));
49 Char_t rng = (i == 0 ? 'I' : (i % 2 == 1 ? 'I' : 'O'));
50 fRing[i] = new AliFMDAnaFlowRing(det,rng,n,EtaMin(i)-.1,EtaMax(i)+.1,
51 0,cut0,Cut1(i));
52 }
53 }
54 //__________________________________________________________________
55 /** Constructor */
56 AliFMDAnaFlowKine(Float_t deta, Bool_t segments=true,
57 Bool_t primary_only=true, Float_t cut0=.68)
58 : fAxis(Int_t((5.3- -3.7)/deta), -3.7, 5.3),
59 fV2("all", "Full FMD", 2, fAxis),
60 fPrimaryOnly(primary_only),
61 fSegmented(segments)
62 {
63 for (int i = 0; i < 5; i++) {
64 UShort_t det = (i == 0 ? 1 : i <= 2 ? 2 : 3);
65 Char_t rng = (i % 2 == 0 ? 'I' : 'O');
66 fRing[i] = new AliFMDAnaFlowRing(det,rng,fAxis,0,cut0,Cut1(i));
67 }
68 }
69 //__________________________________________________________________
70 Int_t FindRing(UShort_t det, Char_t ring) const
71 {
72 Int_t idx = -1;
73 switch (det) {
74 case 1: idx = 0; break;
75 case 2: idx = 1 + (ring == 'O' || ring == 'o' ? 1 : 0); break;
76 case 3: idx = 3 + (ring == 'O' || ring == 'o' ? 1 : 0); break;
77 }
78 return idx;
79 }
80 //__________________________________________________________________
81 /** Initialize */
82 Bool_t Init()
83 {
84 fPhis.Set(51200);
85 fEtas.Set(51200);
86 fWeights.Set(51200);
87 fNObserv = 0;
88 fNEvents = 0;
89
90 Float_t pi2 = 2 * TMath::Pi();
91
92 fPsi = new TH2D("all_psi", "Event planes",
93 fAxis.N(), fAxis.Bins(), 40, 0, pi2);
94 fPsi->SetXTitle("#eta");
95 fPsi->SetYTitle("#Psi");
96 fPsi->SetZTitle("Frequency");
97 fPsi->SetMarkerStyle(20);
98
99 fDPsi = new TProfile("all_dpsi", "Distance to real #Psi_{R}",
100 fAxis.N(), fAxis.Bins(), -pi2, pi2);
101 fDPsi->SetXTitle("#eta");
102 fDPsi->SetYTitle("#Delta#Psi=|#Psi-#Psi_{R}|");
103
104 fDPhi = new TProfile2D("all_dphi","Distribution of phi signals",
105 fAxis.N(), -3.7, 5.3, 40, 0, pi2, 0, 20);
106 fDPhi->SetXTitle("#eta");
107 fDPhi->SetYTitle("#varphi");
108 fDPhi->SetZTitle("signal");
109 fDPhi->SetMarkerStyle(20);
110
111 for (size_t i = 0; i < 5; i++) fRing[i]->Init();
112 return AliFMDInput::Init();
113 }
114 //__________________________________________________________________
115 /** Called at the beginning of an event */
116 Bool_t Begin(Int_t ev)
117 {
118 fNObserv = 0;
119 fNEvents++;
120 Float_t psir = reactionplane(fNEvents-1);
121 for (size_t i = 0; i < 5; i++) {
122 fRing[i]->Begin();
123 fRing[i]->fTrue.SetPsi(psir);
124 }
125 return AliFMDInput::Begin(ev);
126 }
127 //__________________________________________________________________
128 /** Called at the beginning of an event */
129 Bool_t ProcessHit(AliFMDHit* hit, TParticle* p)
130 {
131 if (fPrimaryOnly && !p->IsPrimary()) return kTRUE;
132 if (hit->IsStop()) return kTRUE;
133 UShort_t det = hit->Detector();
134 Char_t rng = hit->Ring();
135 UShort_t sec = hit->Sector();
136 AliFMDAnaFlowRing* frn = fRing[FindRing(det, rng)];
137 Double_t the = p->Theta();
138 Double_t phi = p->Phi();
139 Double_t eta = -TMath::Log(TMath::Tan(the/2));
140 // Phi
141 if (fSegmented) {
142 Double_t dPhi = 2*TMath::Pi() / frn->NSeq();
143 if (det == 3) phi = TMath::Pi() - (sec + .5) * dPhi;
144 else phi = (sec + .5) * dPhi;
145 if (phi < 0) phi += 2 * TMath::Pi();
146 }
147
148 Double_t mult = 1; // hit->Edep();
149 frn->Fill(phi, eta, mult);
150 Fill(phi, eta, mult);
151 return kTRUE;
152 }
153
154 //__________________________________________________________________
155 /** Fill in one strip */
156 void Fill(Float_t phi, Float_t eta, Float_t mult)
157 {
158 if (fNObserv >= fEtas.fN) {
159 ULong_t n = Int_t(fEtas.fN * 1.5)+1;
160 fEtas.Set(n);
161 fPhis.Set(n);
162 fWeights.Set(n);
163 }
164 fEtas[fNObserv] = eta;
165 fPhis[fNObserv] = phi;
166 fWeights[fNObserv] = mult;
167 fNObserv++;
168 fDPhi->Fill(eta, phi);
169 }
170 //__________________________________________________________________
171 /** Called at the end of an event */
172 Bool_t End()
173 {
174 for (size_t i = 0; i < 5; i++) fRing[i]->End();
175 float psir = reactionplane(fNEvents-1)/2;
176 // std::cout << "PsiR=" << psir << std::endl;
177 fV2.Event(fPhis.fArray, fEtas.fArray, fWeights.fArray, fNObserv);
178 for (UShort_t i = 0; i < fAxis.N(); i++) {
179 AliFMDFlowBin* bin = fV2.GetBin(i);
180 if (bin && bin->Counts() > 0) {
181 fPsi->Fill(fAxis.BinCenter(i), bin->Psi() - psir);
182 fDPsi->Fill(fAxis.BinCenter(i), bin->Psi() - psir);
183 }
184 }
185 return AliFMDInput::End();
186 }
187 //__________________________________________________________________
188 /** Make a frame */
189 TH1* MakeFrame(const char* xname, const char* yname,
190 Float_t ymin, Float_t ymax)
191 {
192 static int id = 0;
193 id++;
194 TH1* frame = new TH1D(Form("frame%d", id), "Frame",
195 fAxis.N(), fAxis.Bins()[0], fAxis.Bins()[fAxis.N()]);
196 frame->SetMinimum(ymin);
197 frame->SetMaximum(ymax);
198 frame->SetXTitle(xname);
199 frame->SetYTitle(yname);
200 frame->SetLabelFont(132, "XYZ");
201 frame->SetTitleFont(132, "XYZ");
202 frame->SetLabelSize(.08, "XYZ");
203 frame->SetTitleSize(.08, "XYZ");
204 frame->SetTitleOffset(.5, "XYZ");
205 frame->SetNdivisions(8, "XYZ");
206 frame->Draw();
207 return frame;
208 }
209 //__________________________________________________________________
210 /** called at the end of a run */
211 Bool_t Finish()
212 {
213 gStyle->SetPalette(1);
214 gStyle->SetOptTitle(0);
215 gStyle->SetOptStat(0);
216 AliFMDAnaFlowRing* fring[5];
217 for (size_t i = 0; i < 5; i++)
218 fring[i] = static_cast<AliFMDAnaFlowRing*>(fRing[i]);
219
220 TCanvas* c = new TCanvas("c1", "Result", 600, 750);
221 c->cd();
222 c->SetGridy();
223 c->SetTopMargin(0.02);
224 c->SetRightMargin(0.02);
225 c->SetFillColor(0);
226 c->SetBorderSize(0);
227 c->SetBorderMode(0);
228 c->SetHighLightColor(0);
229 c->Divide(1,3,0,0);
230
231 TVirtualPad* p = 0;
232 TH1* f = 0;
233
234 p = c->cd(1);
235 p->SetGridx();
236 p->SetGridy();
237 p->SetFillColor(0);
238 f = MakeFrame("#eta", "v_{2}", 0, .1);
239 for (int i = 0; i < 5; i++) {
240 fring[i]->fV2.Draw("th:same");
241 fring[i]->fTrue.Draw("th:same");
242 }
243 fV2.Draw("th:same");
244
245 p = c->cd(2);
246 p->SetGridx();
247 p->SetGridy();
248 p->SetFillColor(0);
249 f = MakeFrame("#eta", "R_{k}", .7, 1);
250 for (int i = 0; i < 5; i++) fring[i]->fV2.Draw("tr:same");
251 fV2.Draw("tr:same");
252
253 p = c->cd(3);
254 p->SetGridx();
255 p->SetGridy();
256 p->SetFillColor(0);
257 f = MakeFrame("#eta", "Counts", 0, 2e5);
258 f->GetYaxis()->SetNoExponent();
259 for (int i = 0; i < 5; i++) fring[i]->fV2.Draw("tc:same");
260 // fV2.Draw("tc:same");
261
262 c = new TCanvas("c2", "Event plane");
263 c->cd();
264 c->SetGridy();
265 c->SetTopMargin(0.02);
266 c->SetRightMargin(0.02);
267 c->SetFillColor(0);
268 c->SetBorderSize(0);
269 c->SetBorderMode(0);
270 c->SetHighLightColor(0);
271 float pi2 = 2*TMath::Pi();
272 TH2* h = new TH2D("h", "h", fAxis.N(), fAxis.Bins(), 80, -pi2, pi2);
273 h->Draw();
274 for (int i = 0; i < 5; i++) fring[i]->fDPsi.Draw("same col");
275 // fDPhi->Draw("same");
276
277 for (int i = 0; i < 5; i++) {
278 fring[i]->fV2.Print("s");
279 fring[i]->fTrue.Print("s");
280 }
281 fV2.Print("s");
282 // gROOT->GetListOfSpecials()->Add(this);
283 return kTRUE;
284 }
285 //__________________________________________________________________
286 Bool_t IsFolder() const { return kTRUE; }
287 //__________________________________________________________________
288 void Browse(TBrowser* b)
289 {
290 for (size_t i = 0; i < 5; i++)
291 b->Add(fRing[i], fRing[i]->Name());
292 b->Add(&fV2);
293 if (fPsi) b->Add(fPsi);
294 if (fDPsi) b->Add(fDPsi);
295 if (fDPhi) b->Add(fDPhi);
296 }
297 //__________________________________________________________________
298 /** Get min pseudo rapidity of ring idx */
299 Float_t EtaMin(Int_t idx) const {
300 switch (idx) {
301 case 0: return 3.61728; break; // FMD1I
302 case 1: return 2.18122; break; // FMD2I
303 case 2: return 1.79821; break; // FMD2O
304 case 3: return -3.39913; break; // FMD3I
305 case 4: return -2.28925; break; // FMD3O
306 }
307 return 0;
308 }
309 //__________________________________________________________________
310 /** Get max pseudo rapidity of ring idx */
311 Float_t EtaMax(Int_t idx) const {
312 switch (idx) {
313 case 0: return 5.02643; break; // FMD1I
314 case 1: return 3.57899; break; // FMD2I
315 case 2: return 2.39085; break; // FMD2O
316 case 3: return -2.00644; break; // FMD3I
317 case 4: return -1.70080; break; // FMD3O
318 }
319 return 0;
320 }
321 //__________________________________________________________________
322 /** Get cut 1 of ring idx */
323 Float_t Cut1(Int_t idx) const {
324 switch (idx) {
325 case 0: return 1.33857; break; // FMD1I
326 case 1: return 1.30464; break; // FMD2I
327 case 2: return 1.27263; break; // FMD2O
328 case 3: return 1.28052; break; // FMD3I
329 case 4: return 1.23041; break; // FMD3O
330 }
331 return 0;
332 }
333 /** Number of events */
334 Int_t fNEvents;
335 /** Flow rings */
336 AliFMDAnaFlowRing* fRing[5];
337 /** Flow axis */
338 AliFMDFlowAxis fAxis;
339 /** Flow object */
340 AliFMDFlowBinned1D fV2;
341 /** Cache of phi */
342 TArrayD fPhis;
343 /** Cache of eta */
344 TArrayD fEtas;
345 /** Cache of weights */
346 TArrayD fWeights;
347 /** Number of observables */
348 Int_t fNObserv;
349 /** Histogram of Psi */
350 TH2D* fPsi; // Event planes
351 /** Histogram of dPsi */
352 TProfile* fDPsi; // Mean distance to Psi
353 /** Histogram of dPhi */
354 TProfile2D* fDPhi;
355 /** Wether to include primaries only */
356 Bool_t fPrimaryOnly;
357 /** Whether to segment the phi space according to detector */
358 Bool_t fSegmented;
359 ClassDef(AliFMDAnaFlowKine,1)
360};
361
362//____________________________________________________________________
363//
364// EOF
365//