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