git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fixed a warning
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardFlowTaskQC.cxx
CommitLineData
d2bea14e 1//
2// Calculate flow in the forward regions using the Q cumulants method
3//
4// Inputs:
5// - AliAODEvent
6//
7// Outputs:
8// - AnalysisResults.root
9//
10// TODO!
11// - Add centrality stuff
12// END OF TODO!
13//
14#include <TROOT.h>
15#include <TSystem.h>
16#include <TInterpreter.h>
17#include <TChain.h>
18#include <TFile.h>
19#include <TList.h>
20#include <iostream>
21#include <TMath.h>
22#include "AliLog.h"
23#include "AliForwardFlowTaskQC.h"
24#include "AliAnalysisManager.h"
25#include "AliAODHandler.h"
26#include "AliAODInputHandler.h"
27#include "AliAODMCParticle.h"
28#include "AliForwardFlowBase.h"
29#include "AliAODForwardMult.h"
30
31ClassImp(AliForwardFlowTaskQC)
32
33AliForwardFlowTaskQC::AliForwardFlowTaskQC()
34: fDebug(0), // Debug flag
35 fOutputList(0), // Output list
36 fAOD(0), // AOD input event
37 fMC(kFALSE), // MC flag
38 fEtaBins(20) // # of eta bins in histograms
39{
40 //
41 // Default constructor
42 //
43}
44//_____________________________________________________________________
45AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name) :
46 AliAnalysisTaskSE(name),
47 fDebug(0), // Debug flag
48 fOutputList(0), // Output list
49 fAOD(0), // AOD input event
50 fMC(kFALSE), // MC flag
51 fEtaBins(20) // # of Eta bins
52{
53 //
54 // Constructor
55 //
56 // Parameters:
57 // name: Name of task
58 //
59 for (Int_t n = 0; n <= 4; n++) fv[n] = kTRUE;
60 DefineOutput(1, TList::Class());
61}
62//_____________________________________________________________________
63AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o) :
64 AliAnalysisTaskSE(o),
65 fDebug(o.fDebug), // Debug flag
66 fOutputList(o.fOutputList), // Output list
67 fAOD(o.fAOD), // AOD input event
68 fMC(o.fMC), // MC flag
69 fEtaBins(o.fEtaBins) // # of Eta bins
70{
71 //
72 // Copy constructor
73 //
74 // Parameters:
75 // o Object to copy from
76 //
77 for (Int_t n = 0; n <= 4; n++) fv[n] = o.fv[n];
78 DefineOutput(1, TList::Class());
79}
80//_____________________________________________________________________
81void AliForwardFlowTaskQC::CreateOutputObjects()
82{
83 //
84 // Create output objects
85 //
86 if (!fOutputList)
87 fOutputList = new TList();
88 fOutputList->SetName("QCumulants");
89
90 if (fEtaBins % 20) fEtaBins = 20;
91
92 // Histograms for cumulants analysis
93
94 // We loop over flow histograms here to add different orders of harmonics
95 for (Int_t n = 1; n <= 4; n++) {
96 if (!fv[n]) continue;
97 // Only one flow histogram is needed for each type of data;
98 // x-axis is eta-bins with differential flow, integrated is in underflowbin
99 // y-axis bin 1: (w_<2> * <2>).Re()
100 // y-axis bin 2: (w_<2> * <2>).Im()
101 // y-axis bin 3: w_<2> = M(M-1)
102 // y-axis bin 4: (w_<2> * <2> * <2>).Re()
103 // y-axis bin 5: (w_<2> * <2> * <2>).Im()
104 // y-axis bin 6: (w_<2> * w_<2'> * <2> * <2'>).Re()
105 // y-axis bin 7: (w_<2> * w_<2'> * <2> * <2'>).Im()
106 // y-axis bin 8: w_<2> * w_<2'>
107 // y-axis bin 9: (w_<4> * <4>).Re()
108 // y-axis bin 10: (w_<4> * <4>).Im()
109 // y-axis bin 11: w_<4>
110 // y-axis bin 12: (w_<4> * <4> * <4>).Re()
111 // y-axis bin 13: (w_<4> * <4> * <4>).Im()
112 // y-axis bin 14: w_<2> * w_<4>
113 // y-axis bin 15: (w_<2> * w_<4> * <2> * <4>).Re()
114 // y-axis bin 16: (w_<2> * w_<4> * <2> * <4>).Im()
115 // y-axis bin 17: w_<2> * w_<4'>
116 // y-axis bin 18: (w_<2> * w_<4'> * <2> * <4'>).Re()
117 // y-axis bin 19: (w_<2> * w_<4'> * <2> * <4'>).Im()
118 // y-axis bin 20: w_<4> * w_<2'>
119 // y-axis bin 21: (w_<4> * w_<2'> * <4> * <2'>).Re()
120 // y-axis bin 22: (w_<4> * w_<2'> * <4> * <2'>).Im()
121 // y-axis bin 23: w_<4> * w_<4'>
122 // y-axis bin 24: (w_<4> * w_<4'> * <4> * <4'>).Re()
123 // y-axis bin 25: (w_<4> * w_<4'> * <4> * <4'>).Im()
124 // y-axis bin 26: Qn or pn.Re() = <<cos2phi or psi>>
125 // y-axis bin 27: Qn or pn.Im() = <<sin2phi or psi>>
126 // y-axis bin 28: M or mp
127 // y-axis bin 29: (Qn*Qn-Q2n).Re() = <<cos(2(phi1 or psi1+phi2))>>
128 // y-axis bin 30: (Qn*Qn-Q2n).Im() = <<sin(2(phi1 or psi1+phi2))>>
129 // y-axis bin 31: <<cos(2(phi1 or psi1-phi2-phi3))>>
130 // y-axis bin 32: <<sin(2(phi1 or psi1-phi2-phi3))>>
131 // y-axis bin 33: M*(M-1)*(M-2) or similar for diff
132 // y-axis bin 34: <<cos(2(psi1+phi2-phi3>>
133 // y-axis bin 35: <<sin(2(psi1+phi2-phi3>>
134 TH2D* hFlowHist = new TH2D(Form("hQ%dCumuHist", n), Form("hQ%dCumuHist", n), fEtaBins, -4, 6, 35, 0.5, 35.5);
135 hFlowHist->Sumw2();
136 fOutputList->Add(hFlowHist);
137
138 TH2D* hFlowHistMC = new TH2D(Form("hQ%dCumuHistMC", n), Form("hQ%dCumuHistMC", n), fEtaBins, -4, 6, 35, 0.5, 35.5);
139 hFlowHistMC->Sumw2();
140 fOutputList->Add(hFlowHistMC);
141
142 TH2D* hFlowHistTrRef = new TH2D(Form("hQ%dCumuHistTrRef", n), Form("hQ%dCumuHistTrRef", n), fEtaBins, -4, 6, 35, 0.5, 35.5);
143 hFlowHistTrRef->Sumw2();
144 fOutputList->Add(hFlowHistTrRef);
145
146 // Output histograms
147 TH1D* hCumulant2Flow = new TH1D(Form("hQ%dCumulant2Flow", n), Form("hQ%dCumulant2Flow", n), fEtaBins, -4, 6);
148 hCumulant2Flow->Sumw2();
149 fOutputList->Add(hCumulant2Flow);
150
151 TH1D* hCumulant2FlowMC = new TH1D(Form("hQ%dCumulant2FlowMC", n),Form("hQ%dCumulant2FlowMC", n), fEtaBins, -4, 6);
152 hCumulant2FlowMC->Sumw2();
153 fOutputList->Add(hCumulant2FlowMC);
154
155 TH1D* hCumulant2FlowTrRef = new TH1D(Form("hQ%dCumulant2FlowTrRef", n), Form("hQ%dCumulant2FlowTrRef", n), fEtaBins, -4, 6);
156 hCumulant2FlowTrRef->Sumw2();
157 fOutputList->Add(hCumulant2FlowTrRef);
158
159
160 TH1D* hCumulant4Flow = new TH1D(Form("hQ%dCumulant4Flow", n), Form("hQ%dCumulant4Flow", n), fEtaBins, -4, 6);
161 hCumulant4Flow->Sumw2();
162 fOutputList->Add(hCumulant4Flow);
163
164 TH1D* hCumulant4FlowMC = new TH1D(Form("hQ%dCumulant4FlowMC", n), Form("hQ%dCumulant4FlowMC", n), fEtaBins, -4, 6);
165 hCumulant4FlowMC->Sumw2();
166 fOutputList->Add(hCumulant4FlowMC);
167
168 TH1D* hCumulant4FlowTrRef = new TH1D(Form("hQ%dCumulant4FlowTrRef", n), Form("hQ%dCumulant4FlowTrRef", n), fEtaBins, -4, 6);
169 hCumulant4FlowTrRef->Sumw2();
170 fOutputList->Add(hCumulant4FlowTrRef);
171 }
172
173 // Single Event histograms
174 TH1D* hdNdphiSE = new TH1D("hdNdphiSE","hdNdphiSE", 20, 0, 2*TMath::Pi());
175 hdNdphiSE->Sumw2();
176 fOutputList->Add(hdNdphiSE);
177
178 TH2D* hdNdetadphiSE = new TH2D("hdNdetadphiSE", "hdNdetadphiSE", fEtaBins, -4, 6, 20, 0, 2*TMath::Pi());
179 hdNdetadphiSE->Sumw2();
180 fOutputList->Add(hdNdetadphiSE);
181
182 TH1D* hdNdphiSEMC = new TH1D("hdNdphiSEMC","hdNdphiSEMC", 20, 0, 2*TMath::Pi());
183 hdNdphiSEMC->Sumw2();
184 fOutputList->Add(hdNdphiSEMC);
185
186 TH2D* hdNdetadphiSEMC = new TH2D("hdNdetadphiSEMC", "hdNdetadphiSEMC", fEtaBins, -4, 6, 20, 0, 2*TMath::Pi());
187 hdNdetadphiSEMC->Sumw2();
188 fOutputList->Add(hdNdetadphiSEMC);
189
190 TH1D* hdNdphiSETrRef = new TH1D("hdNdphiSETrRef","hdNdphiSETrRef", 20, 0, 2*TMath::Pi());
191 hdNdphiSETrRef->Sumw2();
192 fOutputList->Add(hdNdphiSETrRef);
193
194 TH2D* hdNdetadphiSETrRef = new TH2D("hdNdetadphiSETrRef", "hdNdetadphiSETrRef", fEtaBins, -4, 6, 20, 0, 2*TMath::Pi());
195 hdNdetadphiSETrRef->Sumw2();
196 fOutputList->Add(hdNdetadphiSETrRef);
197
198 PostData(1, fOutputList);
199}
200//_____________________________________________________________________
201void AliForwardFlowTaskQC::UserExec(Option_t */*option*/)
202{
203 //
204 // Process each event
205 //
206 // Parameters:
207 // option: Not used
208 //
209
210 // Get input event
211 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
212 if (!fAOD) return;
213
214 // Load histograms and reset from last event
215 TH1D* dNdphi = (TH1D*)fOutputList->FindObject("hdNdphiSE");
216 TH2D* dNdetadphi = (TH2D*)fOutputList->FindObject("hdNdetadphiSE");
217
218 dNdphi->Reset();
219 dNdetadphi->Reset();
220
221 // Initiate FlowCommon and fill histograms
222 AliForwardFlowBase* common = new AliForwardFlowBase(fOutputList);
223
224 if (!common->LoopAODFMD(fAOD)) return;
225// else if (!common->LoopAODSPD(fAOD)) return;
226// if (!common->LoopAODFMDandSPD(fAOD)) return;
227
228 // Run analysis
229 for (Int_t n = 1; n <= 4; n++) {
230 if (fv[n])
231 CumulantsMethod("", n);
232 }
233
234 // Find out if there's any MC data present
235 if (!fMC) {
236 TClonesArray* mcArray = 0;
237 mcArray = (TClonesArray*)fAOD->FindListObject(AliAODMCParticle::StdBranchName());
238 if (mcArray) fMC = kTRUE;
239 }
240 if (fMC)
241 ProcessPrimary();
242
243}
244//_____________________________________________________________________
245void AliForwardFlowTaskQC::CumulantsMethod(TString type = "", Int_t harmonic = 2)
246{
247 //
248 // Calculate the Q cumulant of order n
249 //
250 // Parameters:
251 // type: Determines which histograms should be used
252 // - "" = data histograms
253 // - "TrRef" = track reference histograms
254 // - "MC" = MC truth histograms
255 // harmonic: Which harmonic to calculate
256 //
257 Double_t n = harmonic;
258
259 // We get histograms depending on if it's real data or MC truth data
260 TH2D* flowHist = (TH2D*)fOutputList->FindObject(Form("hQ%dCumuHist%s", harmonic, type.Data()));
261 TH1D* dNdphi = (TH1D*)fOutputList->FindObject(Form("hdNdphiSE%s", type.Data()));
262 TH2D* dNdetadphi = (TH2D*)fOutputList->FindObject(Form("hdNdetadphiSE%s", type.Data()));
263
264 // We create the objects needed for the analysis
265 Double_t Mult = dNdphi->GetBinContent(0);
266
267 Double_t QnRe = 0, Q2nRe = 0, QnIm = 0, Q2nIm = 0;
268 Double_t pnRe = 0, p2nRe = 0, qnRe = 0, qnnRe = 0, pnIm = 0, p2nIm = 0, qnIm = 0, qnnIm = 0;
269 Double_t avg2 = 0, avg4 = 0, avg2p = 0, avg4p = 0;
270 Double_t w2avg2sq = 0, w2pavg2psq = 0, w4avg4sq = 0, w4pavg4psq = 0;
271 Double_t w2w2pavg2avg2p = 0, w2w4avg2avg4 = 0, w2pw4pavg2pavg4p = 0;
272 Double_t w2w4pavg2avg4p = 0, w4w2pavg4avg2p = 0, w4w4pavg4avg4p = 0;
273 Double_t CosPhi1Phi2 = 0, CosPhi1Phi2Phi3m = 0, CosPhi1Phi2Phi3p = 0;
274 Double_t SinPhi1Phi2 = 0, SinPhi1Phi2Phi3m = 0, SinPhi1Phi2Phi3p = 0;
275 Double_t Phii = 0;
276 Double_t multi = 0, multp = 0, mp = 0, mq = 0;
277 Double_t W2 = 0, W4 = 0, W2p = 0, W4p = 0;
278
279 // We loop over the data 1 time!
280 for (Int_t eta = 1; eta <= dNdetadphi->GetNbinsX(); eta++) {
281 // The values for each individual eta bins are reset
282 mp = 0;
283 pnRe = 0;
284 p2nRe = 0;
285 pnIm = 0;
286 p2nIm = 0;
287
288 for (Int_t phii = 1; phii <= dNdphi->GetNbinsX()+1; phii++) {
289 Phii = dNdphi->GetXaxis()->GetBinCenter(phii);
290 multi = dNdphi->GetBinContent(phii);
291
292 // In the phi loop on the first eta loop the integrated flow
293 // is calculated from the dNdphi histogram
294 if(eta == 1) {
295 QnRe += multi*TMath::Cos(n*Phii);
296 QnIm += multi*TMath::Sin(n*Phii);
297 Q2nRe += multi*TMath::Cos(2.*n*Phii);
298 Q2nIm += multi*TMath::Sin(2.*n*Phii);
299 }
300
301 // For each eta bin the necesarry values for differential flow
302 // is calculated. Here is the loop over the phi's.
303 multp = dNdetadphi->GetBinContent(eta, phii);
304 mp += multp;
305 pnRe += multp*TMath::Cos(n*Phii);
306 pnIm += multp*TMath::Sin(n*Phii);
307 p2nRe += multp*TMath::Cos(2.*n*Phii);
308 p2nIm += multp*TMath::Sin(2.*n*Phii);
309 }
310
311 // The integrated flow is calculated
312 if (eta == 1) {
313 Double_t Eta = flowHist->GetXaxis()->GetBinCenter(0);
314
315 // 2-particle
316 W2 = Mult * (Mult - 1.);
317 avg2 = QnRe*QnRe + QnIm*QnIm - Mult;
318 avg2 /= W2;
319 w2avg2sq = W2 * avg2 * avg2;
320
321 flowHist->Fill(Eta, 1, W2 * avg2);
322 flowHist->Fill(Eta, 3, W2);
323 flowHist->Fill(Eta, 4, w2avg2sq);
324
325 flowHist->Fill(Eta, 26, QnRe);
326 flowHist->Fill(Eta, 27, QnIm);
327 flowHist->Fill(Eta, 28, Mult);
328
329 // 4-particle
330 W4 = Mult * (Mult - 1.) * (Mult - 2.) * (Mult - 3.);
331 Double_t real = Q2nRe*QnRe*QnRe - Q2nRe*QnIm*QnIm + 2.*Q2nIm*QnRe*QnIm;
332
333 avg4 = TMath::Power(QnRe*QnRe + QnIm*QnIm, 2);
334 avg4 += Q2nRe*Q2nRe + Q2nIm*Q2nIm - 2.*real;
335 avg4 -= 4.*(Mult - 2.)*(QnRe*QnRe + QnIm*QnIm) - 2.*Mult*(Mult - 3.);
336
337 avg4 /= W4;
338 w4avg4sq = W4 * avg4 * avg4;
339 w2w4avg2avg4 = W2 * W4 * avg2 * avg4;
340
341 flowHist->Fill(Eta, 9, W4 * avg4);
342 flowHist->Fill(Eta, 11, W4);
343 flowHist->Fill(Eta, 12, w4avg4sq);
344 flowHist->Fill(Eta, 14, W2 * W4);
345 flowHist->Fill(Eta, 15, w2w4avg2avg4);
346
347 CosPhi1Phi2 = QnRe*QnRe - QnIm*QnIm - Q2nRe;
348 SinPhi1Phi2 = 2.*QnRe*QnIm - Q2nIm;
349
350 CosPhi1Phi2Phi3m = TMath::Power(QnRe, 3) + QnRe*QnIm*QnIm;
351 CosPhi1Phi2Phi3m -= QnRe*Q2nRe + QnIm*Q2nIm + 2.*(Mult - 1.)*QnRe;
352 SinPhi1Phi2Phi3m = -TMath::Power(QnIm, 3) - QnRe*QnRe*QnIm;
353 SinPhi1Phi2Phi3m -= QnIm*Q2nRe - QnRe*Q2nIm + 2.*(Mult - 1.)*QnIm;
354
355 flowHist->Fill(Eta, 29, CosPhi1Phi2);
356 flowHist->Fill(Eta, 30, SinPhi1Phi2);
357 flowHist->Fill(Eta, 31, CosPhi1Phi2Phi3m);
358 flowHist->Fill(Eta, 32, SinPhi1Phi2Phi3m);
359 flowHist->Fill(Eta, 33, Mult*(Mult-1.)*(Mult-2.));
360
361 // Count number of events
362 flowHist->Fill(Eta, 0., 1.);
363 } // end of harmonics loop
364
365 // Differential flow calculations for each eta bin is done:
366 if (mp == 0) continue;
367 Double_t Eta = dNdetadphi->GetXaxis()->GetBinCenter(eta);
368
369 mq = mp;
370 qnRe = pnRe;
371 qnIm = pnIm;
372 qnnRe = p2nRe;
373 qnnIm = p2nIm;
374
375 // 2-particle differential flow
376 W2p = mp * Mult - mq;
377 avg2p = pnRe*QnRe + pnIm*QnIm - mq;
378 avg2p /= W2p;
379 w2pavg2psq = W2p * avg2p * avg2p;
380 w2w2pavg2avg2p = W2 * W2p * avg2 * avg2p;
381
382 flowHist->Fill(Eta, 1, W2p * avg2p);
383 flowHist->Fill(Eta, 3, W2p);
384 flowHist->Fill(Eta, 4, w2pavg2psq);
385 flowHist->Fill(Eta, 6, w2w2pavg2avg2p);
386 flowHist->Fill(Eta, 8, W2 * W2p);
387
388 flowHist->Fill(Eta, 26, pnRe);
389 flowHist->Fill(Eta, 27, pnIm);
390 flowHist->Fill(Eta, 28, mp);
391
392 // 4-particle differential flow
393 W4p = (mp * Mult - 3.*mq)*(Mult - 1.)*(Mult - 2.);
394
395 avg4p = pnRe*QnRe*(QnRe*QnRe + QnIm*QnIm) + pnIm*QnIm*(QnRe*QnRe + QnIm*QnIm);
396 avg4p -= qnnRe*QnRe*QnRe - qnnRe*QnIm*QnIm + 2.*qnnIm*QnRe*QnIm;
397 avg4p -= pnRe*QnRe*Q2nRe - pnRe*QnIm*Q2nIm + pnIm*QnRe*Q2nIm + pnIm*QnIm*Q2nRe;
398 avg4p -= 2.*Mult*(pnRe*QnRe + pnIm*QnIm);
399
400 avg4p += - 2.*mq*(QnRe*QnRe + QnIm*QnIm) + 7.*(qnRe*QnRe + qnIm*QnIm);
401 avg4p += - (QnRe*qnRe + QnIm*qnIm) + (qnnRe*Q2nRe + qnnIm*Q2nIm);
402 avg4p += 2.*(pnRe*QnRe + pnIm*QnIm) + 2.*mq*Mult - 6.*mq;
403 avg4p /= W4p;
404
405 w4pavg4psq = W4p * avg4p * avg4p;
406 w2w4pavg2avg4p = W2 * W4p * avg2 * avg4p;
407 w4w2pavg4avg2p = W4 * W2p * avg4 * avg2p;
408 w4w4pavg4avg4p = W4 * W4p * avg4 * avg4p;
409 w2pw4pavg2pavg4p = W2p * W4p * avg2p * avg4p;
410
411 flowHist->Fill(Eta, 9, W4p * avg4p);
412 flowHist->Fill(Eta, 11, W4p);
413 flowHist->Fill(Eta, 12, w4pavg4psq);
414 flowHist->Fill(Eta, 14, W2p * W4p);
415 flowHist->Fill(Eta, 15, w2pw4pavg2pavg4p);
416 flowHist->Fill(Eta, 17, W2 * W4p);
417 flowHist->Fill(Eta, 18, w2w4pavg2avg4p);
418 flowHist->Fill(Eta, 20, W4 * W2p);
419 flowHist->Fill(Eta, 21, w4w2pavg4avg2p);
420 flowHist->Fill(Eta, 23, W4 * W4p);
421 flowHist->Fill(Eta, 24, w4w4pavg4avg4p);
422
423 CosPhi1Phi2 = pnRe*QnRe - pnIm*QnIm - qnnRe;
424 SinPhi1Phi2 = pnRe*QnIm + pnIm*QnRe - qnnIm;
425
426 CosPhi1Phi2Phi3p = pnRe*(QnRe*QnRe + QnIm*QnIm - Mult);
427 CosPhi1Phi2Phi3p -= qnnRe*QnRe - qnnIm*QnIm + mq*QnRe - 2.*qnRe;
428 SinPhi1Phi2Phi3p = pnIm*(QnRe*QnRe + QnIm*QnIm - Mult);
429 SinPhi1Phi2Phi3p -= qnnIm*QnRe - qnnRe*QnIm + mq*QnIm - 2.*qnIm;
430
431 CosPhi1Phi2Phi3m = pnRe*(QnRe*QnRe - QnIm*QnIm) + 2.*pnIm*QnRe*QnIm;
432 CosPhi1Phi2Phi3m -= pnRe*Q2nRe + pnIm*Q2nIm + 2.*mq*QnRe - 2.*qnRe;
433 SinPhi1Phi2Phi3m = pnIm*(QnRe*QnRe - QnIm*QnIm) - 2.*pnRe*QnRe*QnIm;
434 SinPhi1Phi2Phi3m += - pnIm*Q2nRe + pnRe*Q2nIm + 2.*mq*QnIm - 2.*qnIm;
435
436 flowHist->Fill(Eta, 29, CosPhi1Phi2);
437 flowHist->Fill(Eta, 30, SinPhi1Phi2);
438 flowHist->Fill(Eta, 31, CosPhi1Phi2Phi3m);
439 flowHist->Fill(Eta, 32, SinPhi1Phi2Phi3m);
440 flowHist->Fill(Eta, 33, (mp*Mult-2.*mq)*(Mult-1.));
441 flowHist->Fill(Eta, 34, CosPhi1Phi2Phi3p);
442 flowHist->Fill(Eta, 35, SinPhi1Phi2Phi3p);
443
444 }
445
446}
447//_____________________________________________________________________
448void AliForwardFlowTaskQC::Terminate(Option_t */*option*/)
449{
450 //
451 // End of job
452 // Finalizes the Q cumulant calculations
453 //
454 // Parameters:
455 // option Not used
456 //
457
458 TH2D* cumulantsHist;
459 TH1D* cumulant2Hist;
460 TH1D* cumulant4Hist;
461
462 Int_t nLoops = (fMC ? 3 : 1);
463
464 // Do a loop over the difference analysis types, calculating flow
465 // 1 loop for real data, 3 for MC data
466 // inside each is a nested loop over each harmonic (1, 2, 3 and 4 at the moment)
467 for (Int_t loop = 1; loop <= nLoops; loop++) {
468
469 TString type;
470 if (loop == 1) type = "";
471 if (loop == 2) type = "MC";
472 if (loop == 3) type = "TrRef";
473
474 for (Int_t n = 1; n <= 4; n++) {
475 if (!fv[n]) continue;
476
477 cumulantsHist = (TH2D*)fOutputList->FindObject(Form("hQ%dCumuHist%s", n, type.Data()));
478 cumulant2Hist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant2Flow%s", n, type.Data()));
479 cumulant4Hist = (TH1D*)fOutputList->FindObject(Form("hQ%dCumulant4Flow%s", n, type.Data()));
480
481 // For flow calculations
482 Double_t Avg2 = 0, c2 = 0, vTwo2 = 0, Avg4 = 0, c4 = 0, vTwo4 = 0;
483 Double_t Avg2p = 0, d2 = 0, vTwo2diff = 0, Avg4p = 0, d4 = 0, vTwo4diff = 0;
484 Double_t W2 = 0, W4 = 0, W2p = 0, W4p = 0, sqrtW2sq = 0, sqrtW2psq = 0, W2W2p = 0;
485 Double_t W2W4 = 0, W2W4p = 0, W4W2p = 0, W4W4p = 0, W2pW4p = 0;
486 Double_t sqrtW4sq = 0, sqrtW4psq = 0;
487 Double_t W2avg2 = 0, W2pavg2p = 0, W4avg4 = 0, W4pavg4p = 0;
488 Double_t AvgCos2Phi = 0, AvgSin2Phi = 0, Mult = 0, AvgCos2Phi1Phi2 = 0, AvgSin2Phi1Phi2 = 0;
489 Double_t AvgCos2Phi1Phi2Phi3m = 0, AvgSin2Phi1Phi2Phi3m = 0, Multm1m2 = 0;
490 Double_t AvgCos2Psi = 0, AvgSin2Psi = 0, mp = 0, AvgCos2Psi1Phi2 = 0, AvgSin2Psi1Phi2 = 0;
491 Double_t AvgCos2Psi1Phi2Phi3m = 0, AvgSin2Psi1Phi2Phi3m = 0, mpqMult = 0;
492 Double_t AvgCos2Psi1Phi2Phi3p = 0, AvgSin2Psi1Phi2Phi3p = 0;
493
494 // For error calculations
495 Double_t W2avg2sq = 0, W2W2pavg2avg2p = 0, W2pavg2psq = 0;
496 Double_t W4avg4sq = 0, W2W4avg2avg4 = 0, W4W4pavg4avg4p = 0, W4pavg4psq = 0;
497 Double_t W2W4pavg2avg4p = 0, W4W2pavg4avg2p = 0;
498 Double_t sAvg2sq = 0, sAvg2psq = 0, sAvg4sq = 0, sAvg4psq = 0;
499 Double_t vTwo2err = 0, vTwo2diffErr = 0, vTwo4err = 0, vTwo4diffErr = 0;
500 Double_t Cov22p = 0, Cov24 = 0, Cov24p = 0, Cov42p = 0, Cov44p = 0, Cov2p2np = 0;
501 Double_t W2pW4pavg2pavg4p = 0;
502
503 for (Int_t eta = 0; eta <= cumulantsHist->GetNbinsX(); eta++) {
504 if (eta == 0) {
505 // 2-particle reference flow
506 W2avg2 = cumulantsHist->GetBinContent(eta, 1);
507 if (!W2avg2) continue;
508 W2 = cumulantsHist->GetBinContent(eta, 3);
509 AvgCos2Phi = cumulantsHist->GetBinContent(eta, 26);
510 AvgSin2Phi = cumulantsHist->GetBinContent(eta, 27);
511 Mult = cumulantsHist->GetBinContent(eta, 28);
512 AvgCos2Phi /= Mult;
513 AvgSin2Phi /= Mult;
514 Avg2 = W2avg2 / W2;
515 c2 = Avg2 - TMath::Power(AvgCos2Phi, 2) - TMath::Power(AvgSin2Phi, 2);
516 vTwo2 = TMath::Sqrt(c2);
517 cumulant2Hist->SetBinContent(eta, vTwo2);
518
519 // 4-particle reference flow
520 W4avg4 = cumulantsHist->GetBinContent(eta, 9);
521 W4 = cumulantsHist->GetBinContent(eta, 11);
522 AvgCos2Phi1Phi2 = cumulantsHist->GetBinContent(eta, 29);
523 AvgSin2Phi1Phi2 = cumulantsHist->GetBinContent(eta, 30);
524 AvgCos2Phi1Phi2 /= W2;
525 AvgSin2Phi1Phi2 /= W2;
526 AvgCos2Phi1Phi2Phi3m = cumulantsHist->GetBinContent(eta, 31);
527 AvgSin2Phi1Phi2Phi3m = cumulantsHist->GetBinContent(eta, 32);
528 Multm1m2 = cumulantsHist->GetBinContent(eta, 33);
529 AvgCos2Phi1Phi2Phi3m /= Multm1m2;
530 AvgSin2Phi1Phi2Phi3m /= Multm1m2;
531 Avg4 = W4avg4 / W4;
532 c4 = Avg4 - 2. * Avg2 * Avg2;
533 c4 -= 4.*AvgCos2Phi*AvgCos2Phi1Phi2Phi3m;
534 c4 += 4.*AvgSin2Phi*AvgSin2Phi1Phi2Phi3m;
535 c4 -= TMath::Power(AvgCos2Phi1Phi2, 2) + TMath::Power(AvgSin2Phi1Phi2 , 2);
536 c4 += 4.*AvgCos2Phi1Phi2*(TMath::Power(AvgCos2Phi, 2) - TMath::Power(AvgSin2Phi, 2));
537 c4 += 8.*AvgSin2Phi1Phi2*AvgSin2Phi*AvgCos2Phi;
538 c4 += 8.*Avg2*(TMath::Power(AvgCos2Phi, 2) + TMath::Power(AvgSin2Phi, 2));
539 c4 -= 6.*TMath::Power(TMath::Power(AvgCos2Phi, 2)+TMath::Power(AvgSin2Phi, 2), 2);
540
541 vTwo4 = TMath::Power(-c4, 0.25);
542 cumulant4Hist->SetBinContent(eta, vTwo4);
543
544 // 2-particle reference flow error calculations
545 W2avg2sq = cumulantsHist->GetBinContent(eta, 4);
546 sqrtW2sq = cumulantsHist->GetBinError(eta, 3);
547
548 sAvg2sq = VarSQ(W2avg2sq, Avg2, W2, W2avg2, sqrtW2sq);
549 vTwo2err = sqrtW2sq * TMath::Sqrt(sAvg2sq) / (2. * TMath::Sqrt(Avg2) * W2);
550 cumulant2Hist->SetBinError(eta, vTwo2err);
551
552 // 4-particle reference flow error calculations
553 W4avg4sq = cumulantsHist->GetBinContent(eta, 12);
554 sqrtW4sq = cumulantsHist->GetBinError(eta, 11);
555 W2W4 = cumulantsHist->GetBinContent(eta, 14);
556 W2W4avg2avg4 = cumulantsHist->GetBinContent(eta, 15);
557
558 sAvg4sq = VarSQ(W4avg4sq, Avg4, W4, W4avg4, sqrtW4sq);
559 Cov24 = CovXY(W2W4avg2avg4, W2W4, Avg2*Avg4, W2, W4);
560
561 vTwo4err = Avg2*Avg2 * TMath::Power(sqrtW2sq, 2) * sAvg2sq / (W2*W2);
562 vTwo4err += TMath::Power(sqrtW4sq, 2) * sAvg4sq / (16. * W4*W4);
563 vTwo4err -= Avg2 * W2W4 * Cov24 / (2. * W2 * W4);
564 vTwo4err /= TMath::Power(2. * Avg2*Avg2 - Avg4, 1.5);
565 vTwo4err = TMath::Sqrt(vTwo4err);
566 cumulant4Hist->SetBinError(eta, vTwo4err);
567
568 continue;
569 }
570
571 // 2-particle differential flow
572 W2pavg2p = cumulantsHist->GetBinContent(eta, 1);
573 if (!W2pavg2p) continue;
574 W2p = cumulantsHist->GetBinContent(eta, 3);
575 AvgCos2Psi = cumulantsHist->GetBinContent(eta, 26);
576 AvgSin2Psi = cumulantsHist->GetBinContent(eta, 27);
577 mp = cumulantsHist->GetBinContent(eta, 28);
578 AvgCos2Psi /= mp;
579 AvgSin2Psi /= mp;
580 Avg2p = W2pavg2p / W2p;
581 d2 = Avg2p - AvgCos2Psi*AvgCos2Phi - AvgSin2Psi*AvgSin2Phi;
582 vTwo2diff = d2 / TMath::Sqrt(c2);
583 cumulant2Hist->SetBinContent(eta, vTwo2diff);
584
585 // 4-particle differential flow
586 W4pavg4p = cumulantsHist->GetBinContent(eta, 9);
587 W4p = cumulantsHist->GetBinContent(eta, 11);
588 AvgCos2Psi1Phi2 = cumulantsHist->GetBinContent(eta, 29);
589 AvgSin2Psi1Phi2 = cumulantsHist->GetBinContent(eta, 30);
590 AvgCos2Psi1Phi2 /= W2p;
591 AvgSin2Psi1Phi2 /= W2p;
592 AvgCos2Psi1Phi2Phi3m = cumulantsHist->GetBinContent(eta, 31);
593 AvgSin2Psi1Phi2Phi3m = cumulantsHist->GetBinContent(eta, 32);
594 mpqMult = cumulantsHist->GetBinContent(eta, 33);
595 AvgCos2Psi1Phi2Phi3m /= mpqMult;
596 AvgSin2Psi1Phi2Phi3m /= mpqMult;
597 AvgCos2Psi1Phi2Phi3p = cumulantsHist->GetBinContent(eta, 34);
598 AvgSin2Psi1Phi2Phi3p = cumulantsHist->GetBinContent(eta, 35);
599 AvgCos2Psi1Phi2Phi3p /= mpqMult;
600 AvgSin2Psi1Phi2Phi3p /= mpqMult;
601
602 Avg4p = W4pavg4p / W4p;
603 d4 = Avg4p - 2. * Avg2p * Avg2;
604 d4 -= AvgCos2Psi*AvgCos2Phi1Phi2Phi3m;
605 d4 += AvgSin2Psi*AvgSin2Phi1Phi2Phi3m;
606 d4 -= AvgCos2Phi*AvgCos2Psi1Phi2Phi3m;
607 d4 += AvgSin2Phi*AvgSin2Psi1Phi2Phi3m;
608 d4 -= 2.*AvgCos2Phi*AvgCos2Psi1Phi2Phi3p;
609 d4 -= 2.*AvgSin2Phi*AvgSin2Psi1Phi2Phi3p;
610 d4 -= AvgCos2Psi1Phi2*AvgCos2Phi1Phi2;
611 d4 -= AvgSin2Psi1Phi2*AvgSin2Phi1Phi2;
612 d4 += 2.*AvgCos2Phi1Phi2*(AvgCos2Psi*AvgCos2Phi - AvgSin2Psi*AvgSin2Phi);
613 d4 += 2.*AvgSin2Phi1Phi2*(AvgCos2Psi*AvgSin2Phi + AvgSin2Psi*AvgCos2Phi);
614 d4 += 4.*Avg2*(AvgCos2Psi*AvgCos2Phi + AvgSin2Psi*AvgSin2Phi);
615 d4 += 2.*AvgCos2Psi1Phi2*(TMath::Power(AvgCos2Phi, 2) - TMath::Power(AvgSin2Phi, 2));
616 d4 += 4.*AvgSin2Psi1Phi2*AvgCos2Phi*AvgSin2Phi;
617 d4 += 4.*Avg2p*(TMath::Power(AvgCos2Phi, 2) + TMath::Power(AvgSin2Phi, 2));
618 d4 -= 6.*(TMath::Power(AvgCos2Phi, 2) - TMath::Power(AvgSin2Phi, 2))
619 *(AvgCos2Psi*AvgCos2Phi-AvgSin2Psi*AvgSin2Phi);
620 d4 -= 12.*AvgCos2Phi*AvgSin2Phi*(AvgSin2Psi*AvgCos2Phi+AvgCos2Psi*AvgSin2Phi);
621
622 vTwo4diff = - d4 / TMath::Power(-c4, 0.75);
623 cumulant4Hist->SetBinContent(eta, vTwo4diff);
624
625 // 2-particle differential flow error calculations
626 W2pavg2psq = cumulantsHist->GetBinContent(eta, 4);
627 sqrtW2psq = cumulantsHist->GetBinError(eta, 3);
628 W2W2pavg2avg2p = cumulantsHist->GetBinContent(eta, 6);
629 W2W2p = cumulantsHist->GetBinContent(eta, 8);
630
631 Cov22p = CovXY(W2W2pavg2avg2p, W2W2p, Avg2*Avg2p, W2, W2p);
632 sAvg2psq = VarSQ(W2pavg2psq, Avg2p, W2p, W2pavg2p, sqrtW2psq);
633
634 vTwo2diffErr = Avg2p*Avg2p*TMath::Power(sqrtW2psq, 2)*sAvg2sq/(W2*W2);
635 vTwo2diffErr += 4.*Avg2*Avg2*TMath::Power(sqrtW2psq, 2)*sAvg2psq/(W2p*W2p);
636 vTwo2diffErr -= 4.*Avg2*Avg2p*W2W2p*Cov22p/(W2*W2p);
637 vTwo2diffErr /= (4. * TMath::Power(Avg2, 3));
638 vTwo2diffErr = TMath::Sqrt(vTwo2diffErr);
639 cumulant2Hist->SetBinError(eta, vTwo2diffErr);
640
641 // 4-particle differential flow error calculations
642 sqrtW4psq = cumulantsHist->GetBinError(eta, 11);
643 W4pavg4psq = cumulantsHist->GetBinContent(eta, 12);
644 W2pW4p = cumulantsHist->GetBinContent(eta, 14);
645 W2pW4pavg2pavg4p = cumulantsHist->GetBinContent(eta, 15);
646 W2W4p = cumulantsHist->GetBinContent(eta, 17);
647 W2W4pavg2avg4p = cumulantsHist->GetBinContent(eta, 18);
648 W4W2p = cumulantsHist->GetBinContent(eta, 20);
649 W4W2pavg4avg2p = cumulantsHist->GetBinContent(eta, 21);
650 W4W4p = cumulantsHist->GetBinContent(eta, 23);
651 W4W4pavg4avg4p = cumulantsHist->GetBinContent(eta, 24);
652
653 sAvg4psq = VarSQ(W4pavg4psq, Avg4p, W4p, W4pavg4p, sqrtW4psq);
654 Cov24p = CovXY(W2W4pavg2avg4p, W2W4p, Avg2*Avg4p, W2, W4p);
655 Cov42p = CovXY(W4W2pavg4avg2p, W4W2p, Avg4*Avg2p, W4, W2p);
656 Cov44p = CovXY(W4W4pavg4avg4p, W4W4p, Avg4*Avg4p, W4, W4p);
657 Cov2p2np = CovXY(W2pW4pavg2pavg4p, W2pW4p, Avg2p*Avg4p, W2p, W4p);
658
659 // Numbers on the side reference term number in paper (cite needed) loosely
660/*1*/ vTwo4diffErr = TMath::Power(2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p, 2)
661 * TMath::Power(sqrtW2sq, 2) * sAvg2sq / (W2*W2);
662/*2*/ vTwo4diffErr += 9. * TMath::Power(2.*Avg2*Avg2p - Avg4p, 2) * TMath::Power(sqrtW4sq, 2)
663 * sAvg4sq / (16. * W4*W4);
664/*3*/ vTwo4diffErr += 4. * Avg2*Avg2 * TMath::Power(2.*Avg2*Avg2 - Avg4, 2) * TMath::Power(sqrtW2psq, 2)
665 * sAvg2psq / (W2p*W2p);
666/*4*/ vTwo4diffErr += TMath::Power(2.*Avg2*Avg2 - Avg4, 2) * TMath::Power(sqrtW4psq, 2) * sAvg4psq
667 / (W4p*W4p);
668/*5*/ vTwo4diffErr -= 1.5 * (2.*Avg2*Avg2p - Avg4p) * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p)
669 * W2W4 * Cov24 / (W2*W4);
670/*6*/ vTwo4diffErr -= 4. * Avg2 * (2.*Avg2*Avg2 - Avg4)
671 * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p)
672 * W2W2p * Cov22p / (W2 * W2p);
673/*7*/ vTwo4diffErr += 2. * (2.*Avg2*Avg2 - Avg4)
674 * (2.*Avg2*Avg2*Avg2p - 3.*Avg2*Avg4p + 2.*Avg4*Avg2p)
675 * W2W4p * Cov24p / (W2 * W4p);
676/*8*/ vTwo4diffErr += 3.*Avg2*(2.*Avg2*Avg2 - Avg4)*(2.*Avg2*Avg2p - Avg4p)
677 * W4W2p * Cov42p / (W4*W2p);
678/*9*/ vTwo4diffErr -= 1.5 * (2.*Avg2*Avg2 - Avg4)*(2.*Avg2*Avg2p - Avg4p)
679 * W4W4p * Cov44p / (W4 * W4p);
680/*10*/vTwo4diffErr -= 4.*Avg2*TMath::Power(2.*Avg2*Avg2 - Avg4, 2)
681 * W2pW4p * Cov2p2np / (W2p * W4p);
682/*11*/vTwo4diffErr /= TMath::Power(2.*Avg2*Avg2 - Avg4, 3.5);
683 vTwo4diffErr = TMath::Sqrt(vTwo4diffErr);
684
685 cumulant4Hist->SetBinError(eta, vTwo4diffErr);
686 } // End of eta loop
687
688 // Number of events:
a0a945d3 689 Int_t nEv = (Int_t)cumulantsHist->GetBinContent(0,0);
d2bea14e 690 cumulant2Hist->SetBinContent(cumulant2Hist->GetNbinsX() + 1, nEv);
691 cumulant4Hist->SetBinContent(cumulant4Hist->GetNbinsX() + 1, nEv);
692 } // End of harmonics loop
693 } // End of type loop
694}
695// _____________________________________________________________________
696void AliForwardFlowTaskQC::ProcessPrimary()
697{
698 //
699 // If fMC == kTRUE this function takes care of organizing the input
700 // Monte Carlo data and histograms so AliForwardFlowTaskQC::QCumulants
701 // can be run on it.
702 //
703
704 // Histograms are loaded and reset
705 TH1D* dNdphi = (TH1D*)fOutputList->FindObject("hdNdphiSEMC");
706 TH2D* dNdetadphi = (TH2D*)fOutputList->FindObject("hdNdetadphiSEMC");
707 TH1D* dNdphiTrRef = (TH1D*)fOutputList->FindObject("hdNdphiSETrRef");
708 TH2D* dNdetadphiTrRef = (TH2D*)fOutputList->FindObject("hdNdetadphiSETrRef");
709
710 dNdphi->Reset();
711 dNdetadphi->Reset();
712 dNdphiTrRef->Reset();
713 dNdetadphiTrRef->Reset();
714
715 // Loads AliFMDFlowCommon and fills histograms and runs analysis.
716 // AOD events also get a TrackRef histogram
717 AliForwardFlowBase* common = new AliForwardFlowBase(fOutputList);
718
719 if (fAOD) {
720 if (!common->LoopAODmc(fAOD)) return;
721 if (!common->LoopAODtrrefHits(fAOD)) return;
722// if (!common->LoopMCaddptFlow(fAOD)) return;
723// if (!common->LoopMCaddpdgFlow(fAOD)) return;
724// if (!common->LoopMCaddetaFlow(fAOD)) return;
725 }
726
727 // Run analysis on MC truth
728 for (Int_t n = 1; n <= 4; n++) {
729 if (fv[n])
730 CumulantsMethod("MC", n);
731 }
732
733 // Run analysis on TrackRefs
734 for (Int_t n = 1; n <= 4; n++) {
735 if (fv[n])
736 CumulantsMethod("TrRef", n);
737 }
738
739}
740//_____________________________________________________________________
741Double_t AliForwardFlowTaskQC::VarSQ(Double_t wxx2, Double_t x, Double_t wx, Double_t wxx, Double_t sqrtwx2)
742{
743 //
744 // Small function to compute the variance squared - used by Terminte()
745 //
746 Double_t sx;
747
748 sx = wxx2 + x*x*wx - (Double_t)2*x*wxx;
749 sx *= (Double_t)1 / wx;
750 sx *= (Double_t)1 / (1 - TMath::Power(sqrtwx2, 2) / (wx*wx));
751
752 return sx;
753}
754//_____________________________________________________________________
755Double_t AliForwardFlowTaskQC::CovXY(Double_t wxwyxy, Double_t wxwy, Double_t xy, Double_t wx, Double_t wy)
756{
757 //
758 // Small function to compute the covariance between two numbers
759 // - used by Terminate()
760 //
761 Double_t Cov, denominator, numerator;
762
763 denominator = (wxwyxy / wxwy) - xy;
764 numerator = 1 - (wxwy / (wx * wy));
765
766 Cov = denominator / numerator;
767 return Cov;
768}
769//_____________________________________________________________________
770//
771//
772// EOF