A lot of changes after detector review:
[u/mrichter/AliRoot.git] / FMD / scripts / AliFMDAnaFlowKine.C
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
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 AliFMDAnaFlowKine : public AliFMDInput
34 {
35 public:
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 //