Fixed warnings [-Wunused-but-set-variable] from GCC 4.6 -
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardFlowTaskQC.cxx
CommitLineData
d2bea14e 1//
d226802c 2// Calculate flow in the forward and central regions using the Q cumulants method.
d2bea14e 3//
4// Inputs:
5// - AliAODEvent
6//
7// Outputs:
8// - AnalysisResults.root
9//
d2bea14e 10#include <TROOT.h>
11#include <TSystem.h>
12#include <TInterpreter.h>
13#include <TChain.h>
14#include <TFile.h>
15#include <TList.h>
16#include <iostream>
17#include <TMath.h>
2f9be372 18#include "THnSparse.h"
58f5fae2 19#include "TH3D.h"
2f9be372 20#include "TProfile2D.h"
d2bea14e 21#include "AliLog.h"
22#include "AliForwardFlowTaskQC.h"
23#include "AliAnalysisManager.h"
24#include "AliAODHandler.h"
25#include "AliAODInputHandler.h"
26#include "AliAODMCParticle.h"
d2bea14e 27#include "AliAODForwardMult.h"
58f5fae2 28#include "AliAODEvent.h"
29
30//
31// Enumeration for adding and retrieving stuff from the histogram
32//
d226802c 33enum { kW2Two = 1, kW2, kW4Four, kW4, kQnRe, kQnIm, kM,
34 kCosphi1phi2, kSinphi1phi2, kCosphi1phi2phi3m, kSinphi1phi2phi3m, kMm1m2,
35 kw2two, kw2, kw4four, kw4, kpnRe, kpnIm, kmp,
36 kCospsi1phi2, kSinpsi1phi2, kCospsi1phi2phi3m, kSinpsi1phi2phi3m,
37 kmpmq, kCospsi1phi2phi3p, kSinpsi1phi2phi3p };
d2bea14e 38
39ClassImp(AliForwardFlowTaskQC)
9d05ffeb 40#if 0
41; // For emacs
42#endif
d2bea14e 43
44AliForwardFlowTaskQC::AliForwardFlowTaskQC()
58f5fae2 45 : fOutputList(0), // Output list
46 fFlowUtil(0), // AliForwardFlowUtil
9d05ffeb 47 fAOD(0), // AOD input event
58f5fae2 48 fMC(kFALSE), // MC flag
d226802c 49 fEtaBins(48), // # of etaBin bins in histograms
50 fEtaRef(12), // # of Eta bins for reference flow
58f5fae2 51 fAddFlow(0), // Add flow string
52 fAddType(0), // Add flow type #
53 fAddOrder(0), // Add flow order
d226802c 54 fZvertex(1111), // Z vertex range
55 fCent(-1) // Centrality
d2bea14e 56{
57 //
58 // Default constructor
59 //
60}
61//_____________________________________________________________________
62AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name) :
63 AliAnalysisTaskSE(name),
d2bea14e 64 fOutputList(0), // Output list
58f5fae2 65 fFlowUtil(0), // AliForwardFlowUtil
d2bea14e 66 fAOD(0), // AOD input event
67 fMC(kFALSE), // MC flag
d226802c 68 fEtaBins(48), // # of Eta bins
69 fEtaRef(12), // # of Eta bins for reference flow
58f5fae2 70 fAddFlow(0), // Add flow string
71 fAddType(0), // Add flow type #
72 fAddOrder(0), // Add flow order
d226802c 73 fZvertex(1111), // Z vertex range
74 fCent(-1) // Centrality
d2bea14e 75{
76 //
77 // Constructor
78 //
79 // Parameters:
80 // name: Name of task
81 //
d226802c 82 for (Int_t n = 0; n <= 6; n++) fv[n] = kTRUE;
d2bea14e 83 DefineOutput(1, TList::Class());
84}
85//_____________________________________________________________________
86AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o) :
87 AliAnalysisTaskSE(o),
d2bea14e 88 fOutputList(o.fOutputList), // Output list
58f5fae2 89 fFlowUtil(o.fFlowUtil), // AliForwardFlowUtil
d2bea14e 90 fAOD(o.fAOD), // AOD input event
91 fMC(o.fMC), // MC flag
58f5fae2 92 fEtaBins(o.fEtaBins), // # of Eta bins
d226802c 93 fEtaRef(o.fEtaRef), // # of Eta bins for reference flow
58f5fae2 94 fAddFlow(o.fAddFlow), // Add flow string
95 fAddType(o.fAddType), // Add flow type #
96 fAddOrder(o.fAddOrder), // Add flow order
97 fZvertex(o.fZvertex), // Z vertex range
2f9be372 98 fCent(o.fCent) // Centrality
d2bea14e 99{
100 //
101 // Copy constructor
102 //
103 // Parameters:
104 // o Object to copy from
105 //
d226802c 106 for (Int_t n = 0; n <= 6; n++) fv[n] = o.fv[n];
d2bea14e 107 DefineOutput(1, TList::Class());
108}
109//_____________________________________________________________________
9453b19e 110void AliForwardFlowTaskQC::UserCreateOutputObjects()
d2bea14e 111{
112 //
113 // Create output objects
114 //
115 if (!fOutputList)
116 fOutputList = new TList();
117 fOutputList->SetName("QCumulants");
58f5fae2 118 fOutputList->SetOwner();
119 fFlowUtil = new AliForwardFlowUtil(fOutputList);
d226802c 120
121 Double_t x[101] = { 0. };
b4db5d3b 122 // Double_t etaMin = -6;
123 // Double_t etaMax = 6;
d226802c 124
b4db5d3b 125 // First we have a number of options for eta binning, also if it is
126 // not really eta, but centrality or pt we want to do flow as a
127 // function of, then this is possible:
d226802c 128 if (fEtaBins == 5) {
129 x[0] = 0.;
130 x[1] = 1.;
131 x[2] = 2.;
132 x[3] = 3.;
133 x[4] = 4.5;
134 x[5] = 6.0;
b4db5d3b 135 // etaMin = 0;
136 // etaMax = 6;
2f9be372 137 }
2f9be372 138
d226802c 139 else if (fEtaBins == 100) {
140 for (Int_t e = 0; e<= 100; e++) {
141 x[e] = e;
142 }
b4db5d3b 143 // etaMin = 0;
144 // etaMax = 100;
d226802c 145 }
146
147 else {
148 if (fEtaBins % 6) fEtaBins = 48;
149 for (Int_t e = 0; e <=fEtaBins; e++) {
150 x[e] = -6. + e*(12./(Double_t)fEtaBins);
151 }
152 }
153
154 // Reference flow binning
155 Double_t xR[101] = { 0. };
156 for (Int_t r = 0; r <= fEtaRef; r++) {
157 xR[r] = -6 + r*(12./(Double_t)fEtaRef);
158 }
159
160 // Phi binning
161 Double_t phi[21] = { 0. };
162 for (Int_t p = 0; p <= 20; p++) {
163 phi[p] = p*2.*TMath::Pi() / 20.;
164 }
165 Double_t phiMC[201] = { 0. };
166 for (Int_t p = 0; p <= 200; p++) {
167 phiMC[p] = p*2.*TMath::Pi() / 200.;
168 }
d2bea14e 169
170 // Histograms for cumulants analysis
171
172 // We loop over flow histograms here to add different orders of harmonics
58f5fae2 173 TString type = "";
2f9be372 174 for (Int_t loop = 1; loop <= 5; loop++) {
58f5fae2 175
2f9be372 176 if (loop == 1) type = "FMD";
58f5fae2 177 if (loop == 2) type = "SPD";
178 if (loop == 3) type = "MC";
2f9be372 179 if (loop == 4) type = "FMDTR";
180 if (loop == 5) type = "SPDTR";
181
182 TList* tList = new TList();
183 tList->SetName(type.Data());
184 fOutputList->Add(tList);
d2bea14e 185
d226802c 186 for (Int_t n = 1; n <= 6; n++) {
58f5fae2 187 if (!fv[n]) continue;
2f9be372 188
189 TList* vList = new TList();
190 vList->SetName(Form("v%d", n));
191 tList->Add(vList);
192
58f5fae2 193 TString tag = TString();
d226802c 194 for (Int_t c = 0; c < 100; c++) {
58f5fae2 195 // Output histograms
d226802c 196 tag = Form("%d_%d", c, c+1);
197
198 TH3D* hFlowHist = new TH3D(Form("hQ%dCumuHist%s_%s", n, type.Data(), tag.Data()),
199 Form("hQ%dCumuHist%s_%s", n, type.Data(), tag.Data()), fEtaBins, -6, 6, 10, -5, 5, 26, 0.5, 26.5);
200 TAxis* xAxis = hFlowHist->GetXaxis();
201 xAxis->Set(fEtaBins, x);
202 hFlowHist->Sumw2();
203 vList->Add(hFlowHist);
204
2f9be372 205 TProfile* hCumulant2DiffFlow = new TProfile(Form("hQ%dCumulant2DiffFlow%s_%s", n, type.Data(), tag.Data()),
58f5fae2 206 Form("hQ%dCumulant2DiffFlow%s_%s", n, type.Data(), tag.Data()), fEtaBins, x);
207 hCumulant2DiffFlow->Sumw2();
2f9be372 208 vList->Add(hCumulant2DiffFlow);
58f5fae2 209
2f9be372 210 TProfile* hCumulant4DiffFlow = new TProfile(Form("hQ%dCumulant4DiffFlow%s_%s", n, type.Data(), tag.Data()),
58f5fae2 211 Form("hQ%dCumulant4DiffFlow%s_%s", n, type.Data(), tag.Data()), fEtaBins, x);
212 hCumulant4DiffFlow->Sumw2();
2f9be372 213 vList->Add(hCumulant4DiffFlow);
58f5fae2 214 } // end of centrality loop
215
216 } // end of v_{n} loop
217
218 // Single Event histograms
d226802c 219 TH2D* hdNdphiSE = new TH2D(Form("hdNdphiSE%s", type.Data()),
220 Form("hdNdphiSE%s", type.Data()), fEtaRef, xR, (type.Contains("MC") ? 200 : 20), (type.Contains("MC") ? phiMC : phi));
58f5fae2 221 hdNdphiSE->Sumw2();
222 fOutputList->Add(hdNdphiSE);
d2bea14e 223
2f9be372 224 TH2D* hdNdetadphiSE = new TH2D(Form("hdNdetadphiSE%s", type.Data()),
d226802c 225 Form("hdNdetadphiSE%s", type.Data()), fEtaBins, x, (type.Contains("MC") ? 200 : 20), (type.Contains("MC") ? phiMC : phi));
2f9be372 226 hdNdetadphiSE->Sumw2();
227 fOutputList->Add(hdNdetadphiSE);
58f5fae2 228 } // end of type loop
d2bea14e 229
d226802c 230 TProfile2D* pMCTruth = new TProfile2D("pMCTruth", "pMCTruth", 48, -6, 6, 100, 0, 100);
231 TAxis* xAxis = pMCTruth->GetXaxis();
232 xAxis->Set(fEtaBins, x);
2f9be372 233 pMCTruth->Sumw2();
234 fOutputList->Add(pMCTruth);
235
d226802c 236 // Monitoring plots
2f9be372 237 TH1D* cent = new TH1D("Centralities", "Centralities", 100, 0, 100);
58f5fae2 238 fOutputList->Add(cent);
d2bea14e 239
d226802c 240 TH1D* vertexSel = new TH1D("VertexSelected", "VertexSelected", 50, -10, 10);
241 fOutputList->Add(vertexSel);
242 TH1D* vertexAll = new TH1D("VertexAll", "VertexAll", 50, -10, 10);
243 fOutputList->Add(vertexAll);
244
245 TH2D* vertex2D = new TH2D("CoverageVsVertex", "CoverageVsVertex", fEtaBins, -6, 6, 20, -10, 10);
246 fOutputList->Add(vertex2D);
2f9be372 247
d2bea14e 248 PostData(1, fOutputList);
249}
250//_____________________________________________________________________
251void AliForwardFlowTaskQC::UserExec(Option_t */*option*/)
252{
253 //
254 // Process each event
255 //
256 // Parameters:
257 // option: Not used
258 //
259
d226802c 260 // Reset data members
261 fCent = -1;
262 fZvertex = 1111;
263
d2bea14e 264 // Get input event
265 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
266 if (!fAOD) return;
267
2f9be372 268 // Fill histograms
58f5fae2 269 if (!fFlowUtil->LoopAODFMD(fAOD)) return;
d226802c 270 if (!fFlowUtil->LoopAODSPD(fAOD)) return;
d2bea14e 271
d226802c 272 // Get centrality and vertex from flow utility and fill monitoring histograms
2f9be372 273 fCent = fFlowUtil->GetCentrality();
274 fZvertex = fFlowUtil->GetVertex();
d226802c 275
2f9be372 276 TH1D* cent = (TH1D*)fOutputList->FindObject("Centralities");
d226802c 277
2f9be372 278 cent->Fill(fCent);
d226802c 279 TH1D* vertex = (TH1D*)fOutputList->FindObject("VertexSelected");
280 vertex->Fill(fZvertex);
2f9be372 281
d226802c 282 // Run analysis on FMD and SPD
283 for (Int_t n = 1; n <= 6; n++) {
284 if (fv[n]) {
285 CumulantsMethod("FMD", n);
2f9be372 286 CumulantsMethod("SPD", n);
d226802c 287 }
2f9be372 288 }
d226802c 289
290 // And do analysis if there is.
291 if (fMC) {
d2bea14e 292 ProcessPrimary();
d226802c 293 }
d2bea14e 294
d226802c 295 PostData(1, fOutputList);
d2bea14e 296}
297//_____________________________________________________________________
9d05ffeb 298void AliForwardFlowTaskQC::CumulantsMethod(TString type = "",
299 Int_t harmonic = 2)
d2bea14e 300{
301 //
302 // Calculate the Q cumulant of order n
303 //
304 // Parameters:
305 // type: Determines which histograms should be used
d226802c 306 // - "FMD" or "SPD" = data histograms
307 // - "FMDTR" or "SPDTR" = track reference histograms
d2bea14e 308 // - "MC" = MC truth histograms
309 // harmonic: Which harmonic to calculate
310 //
311 Double_t n = harmonic;
312
2f9be372 313 TList* tList = (TList*)fOutputList->FindObject(type.Data());
314 TList* vList = (TList*)tList->FindObject(Form("v%d", harmonic));
d226802c 315 // We get the histograms
316 TH3D* flowHist = (TH3D*)vList->FindObject(Form("hQ%dCumuHist%s_%d_%d", harmonic, type.Data(), (Int_t)fCent, (Int_t)fCent+1));
2f9be372 317 TH2D* dNdetadphi = (TH2D*)fOutputList->FindObject(Form("hdNdetadphiSE%s", type.Data()));
d226802c 318 TH2D* dNdphi = (TH2D*)fOutputList->FindObject(Form("hdNdphiSE%s", /*"FMD"*/type.Data()));
2f9be372 319
d226802c 320 // We create the coordinate array used to fill the THnSparse, centrality and vertex is set from the beginning.
321 Double_t coord[4] = { 0., fCent, fZvertex, 0. };
d2bea14e 322 // We create the objects needed for the analysis
d226802c 323 Double_t dQnRe = 0, dQ2nRe = 0, dQnIm = 0, dQ2nIm = 0, mult = 0;
58f5fae2 324 Double_t pnRe = 0, p2nRe = 0, qnRe = 0, q2nRe = 0, pnIm = 0, p2nIm = 0, qnIm = 0, q2nIm = 0;
2f9be372 325 Double_t two = 0, four = 0, twoPrime = 0, fourPrime = 0;
d226802c 326 Double_t cosPhi1Phi2 = 0, cosPhi1Phi2Phi3m = 0;
327 Double_t sinPhi1Phi2 = 0, sinPhi1Phi2Phi3m = 0;
328 Double_t cosPsi1Phi2 = 0, cosPsi1Phi2Phi3m = 0, cosPsi1Phi2Phi3p = 0;
329 Double_t sinPsi1Phi2 = 0, sinPsi1Phi2Phi3m = 0, sinPsi1Phi2Phi3p = 0;
330 Double_t phi = 0, eta = 0;
d2bea14e 331 Double_t multi = 0, multp = 0, mp = 0, mq = 0;
58f5fae2 332 Double_t w2 = 0, w4 = 0, w2p = 0, w4p = 0;
d226802c 333 Int_t refEtaBin = 0;
334 Bool_t kEventCount = kFALSE;
d2bea14e 335
336 // We loop over the data 1 time!
2f9be372 337 for (Int_t etaBin = 1; etaBin <= dNdetadphi->GetNbinsX(); etaBin++) {
d226802c 338 eta = dNdetadphi->GetXaxis()->GetBinCenter(etaBin);
339 refEtaBin = dNdphi->GetXaxis()->FindBin(eta);
58f5fae2 340 // The values for each individual etaBin bins are reset
d2bea14e 341 mp = 0;
342 pnRe = 0;
343 p2nRe = 0;
344 pnIm = 0;
345 p2nIm = 0;
346
d226802c 347 mult = 0;
348 dQnRe = 0;
349 dQnIm = 0;
350 dQ2nRe = 0;
351 dQ2nIm = 0;
d2bea14e 352
d226802c 353 for (Int_t phiBin = 1; phiBin <= dNdphi->GetNbinsY(); phiBin++) {
354 phi = dNdphi->GetYaxis()->GetBinCenter(phiBin);
355
356 // Reference flow
357 multi = dNdphi->GetBinContent(refEtaBin, phiBin);
358 dQnRe += multi*TMath::Cos(n*phi);
359 dQnIm += multi*TMath::Sin(n*phi);
360 dQ2nRe += multi*TMath::Cos(2.*n*phi);
361 dQ2nIm += multi*TMath::Sin(2.*n*phi);
362 mult += multi;
d2bea14e 363
d226802c 364 // For each etaBin bin the necessary values for differential flow
d2bea14e 365 // is calculated. Here is the loop over the phi's.
d226802c 366 multp = dNdetadphi->GetBinContent(etaBin, phiBin);
d2bea14e 367 mp += multp;
58f5fae2 368 pnRe += multp*TMath::Cos(n*phi);
369 pnIm += multp*TMath::Sin(n*phi);
370 p2nRe += multp*TMath::Cos(2.*n*phi);
d226802c 371 p2nIm += multp*TMath::Sin(2.*n*phi);
d2bea14e 372 }
d226802c 373 if (mult == 0) continue;
374
375 if (!kEventCount) {
376 // Count number of events
377 coord[0] = -7;
378 coord[3] = -0.5;
379 flowHist->Fill(coord[0], coord[2], coord[3], 1.);
380 kEventCount = kTRUE;
381 }
d2bea14e 382
d226802c 383 // The reference flow is calculated
384 coord[0] = eta;
385
386 // 2-particle
387 w2 = mult * (mult - 1.);
388 two = dQnRe*dQnRe + dQnIm*dQnIm - mult;
389 coord[3] = kW2Two;
390 flowHist->Fill(coord[0], coord[2], coord[3], two);
391 coord[3] = kW2;
392 flowHist->Fill(coord[0], coord[2], coord[3], w2);
d2bea14e 393
d226802c 394 coord[3] = kQnRe;
395 flowHist->Fill(coord[0], coord[2], coord[3], dQnRe);
396 coord[3] = kQnIm;
397 flowHist->Fill(coord[0], coord[2], coord[3], dQnIm);
398 coord[3] = kM;
399 flowHist->Fill(coord[0], coord[2], coord[3], mult);
400
401 // 4-particle
402 w4 = mult * (mult - 1.) * (mult - 2.) * (mult - 3.);
403
404 four = 2.*mult*(mult-3.) + TMath::Power((TMath::Power(dQnRe,2.)+TMath::Power(dQnIm,2.)),2.)
405 -4.*(mult-2.)*(TMath::Power(dQnRe,2.) + TMath::Power(dQnIm,2.))
406 -2.*(TMath::Power(dQnRe,2.)*dQ2nRe+2.*dQnRe*dQnIm*dQ2nIm-TMath::Power(dQnIm,2.)*dQ2nRe)
407 +(TMath::Power(dQ2nRe,2.)+TMath::Power(dQ2nIm,2.));
2f9be372 408
d226802c 409 coord[3] = kW4Four;
410 flowHist->Fill(coord[0], coord[2], coord[3], four);
411 coord[3] = kW4;
412 flowHist->Fill(coord[0], coord[2], coord[3], w4);
2f9be372 413
d226802c 414 cosPhi1Phi2 = dQnRe*dQnRe - dQnIm*dQnIm - dQ2nRe;
415 sinPhi1Phi2 = 2.*dQnRe*dQnIm - dQ2nIm;
d2bea14e 416
d226802c 417 cosPhi1Phi2Phi3m = dQnRe*(TMath::Power(dQnRe,2)+TMath::Power(dQnIm,2))-dQnRe*dQ2nRe-dQnIm*dQ2nIm-2.*(mult-1)*dQnRe;
d2bea14e 418
d226802c 419 sinPhi1Phi2Phi3m = -dQnIm*(TMath::Power(dQnRe,2)+TMath::Power(dQnIm,2))+dQnRe*dQ2nIm-dQnIm*dQ2nRe+2.*(mult-1)*dQnIm;
d2bea14e 420
d226802c 421 coord[3] = kCosphi1phi2;
422 flowHist->Fill(coord[0], coord[2], coord[3], cosPhi1Phi2);
423 coord[3] = kSinphi1phi2;
424 flowHist->Fill(coord[0], coord[2], coord[3], sinPhi1Phi2);
425 coord[3] = kCosphi1phi2phi3m;
426 flowHist->Fill(coord[0], coord[2], coord[3], cosPhi1Phi2Phi3m);
427 coord[3] = kSinphi1phi2phi3m;
428 flowHist->Fill(coord[0], coord[2], coord[3], sinPhi1Phi2Phi3m);
429 coord[3] = kMm1m2;
430 flowHist->Fill(coord[0], coord[2], coord[3], mult*(mult-1.)*(mult-2.));
d2bea14e 431
58f5fae2 432 // Differential flow calculations for each etaBin bin is done:
d2bea14e 433 mq = mp;
434 qnRe = pnRe;
435 qnIm = pnIm;
58f5fae2 436 q2nRe = p2nRe;
437 q2nIm = p2nIm;
d2bea14e 438
439 // 2-particle differential flow
58f5fae2 440 w2p = mp * mult - mq;
2f9be372 441 twoPrime = pnRe*dQnRe + pnIm*dQnIm - mq;
d2bea14e 442
d226802c 443 coord[3] = kw2two;
444 flowHist->Fill(coord[0], coord[2], coord[3], twoPrime);
445 coord[3] = kw2;
446 flowHist->Fill(coord[0], coord[2], coord[3], w2p);
447
448 coord[3] = kpnRe;
449 flowHist->Fill(coord[0], coord[2], coord[3], pnRe);
450 coord[3] = kpnIm;
451 flowHist->Fill(coord[0], coord[2], coord[3], pnIm);
452 coord[3] = kmp;
453 flowHist->Fill(coord[0], coord[2], coord[3], mp);
d2bea14e 454
455 // 4-particle differential flow
58f5fae2 456 w4p = (mp * mult - 3.*mq)*(mult - 1.)*(mult - 2.);
d226802c 457
458 fourPrime = (TMath::Power(dQnRe,2.)+TMath::Power(dQnIm,2.))*(pnRe*dQnRe+pnIm*dQnIm)
459 - q2nRe*(TMath::Power(dQnRe,2.)-TMath::Power(dQnIm,2.))
460 - 2.*q2nIm*dQnRe*dQnIm
461 - pnRe*(dQnRe*dQ2nRe+dQnIm*dQ2nIm)
462 + pnIm*(dQnIm*dQ2nRe-dQnRe*dQ2nIm)
463 - 2.*mult*(pnRe*dQnRe+pnIm*dQnIm)
464 - 2.*(TMath::Power(dQnRe,2.)+TMath::Power(dQnIm,2.))*mq
465 + 6.*(qnRe*dQnRe+qnIm*dQnIm)
466 + 1.*(q2nRe*dQ2nRe+q2nIm*dQ2nIm)
467 + 2.*(pnRe*dQnRe+pnIm*dQnIm)
468 + 2.*mq*mult
469 - 6.*mq;
470
471 coord[3] = kw4four;
472 flowHist->Fill(coord[0], coord[2], coord[3], fourPrime);
473 coord[3] = kw4;
474 flowHist->Fill(coord[0], coord[2], coord[3], w4p);
58f5fae2 475
d226802c 476 cosPsi1Phi2 = pnRe*dQnRe - pnIm*dQnIm - q2nRe;
477 sinPsi1Phi2 = pnRe*dQnIm + pnIm*dQnRe - q2nIm;
58f5fae2 478
d226802c 479 cosPsi1Phi2Phi3p = pnRe*(TMath::Power(dQnIm,2.)+TMath::Power(dQnRe,2.)-mult)
480 - 1.*(q2nRe*dQnRe+q2nIm*dQnIm)
481 - mq*dQnRe+2.*qnRe;
482
483 sinPsi1Phi2Phi3p = pnIm*(TMath::Power(dQnIm,2.)+TMath::Power(dQnRe,2.)-mult)
484 - 1.*(q2nIm*dQnRe-q2nRe*dQnIm)
485 - mq*dQnIm+2.*qnIm;
58f5fae2 486
d226802c 487 cosPsi1Phi2Phi3m = pnRe*(TMath::Power(dQnRe,2.)-TMath::Power(dQnIm,2.))+2.*pnIm*dQnRe*dQnIm
488 - 1.*(pnRe*dQ2nRe+pnIm*dQ2nIm)
489 - 2.*mq*dQnRe+2.*qnRe;
490
491 sinPsi1Phi2Phi3m = pnIm*(TMath::Power(dQnRe,2.)-TMath::Power(dQnIm,2.))-2.*pnRe*dQnRe*dQnIm
492 - 1.*(pnIm*dQ2nRe-pnRe*dQ2nIm)
493 + 2.*mq*dQnIm-2.*qnIm;
494
495 coord[3] = kCospsi1phi2;
496 flowHist->Fill(coord[0], coord[2], coord[3], cosPsi1Phi2);
497 coord[3] = kSinpsi1phi2;
498 flowHist->Fill(coord[0], coord[2], coord[3], sinPsi1Phi2);
499 coord[3] = kCospsi1phi2phi3m;
500 flowHist->Fill(coord[0], coord[2], coord[3], cosPsi1Phi2Phi3m);
501 coord[3] = kSinpsi1phi2phi3m;
502 flowHist->Fill(coord[0], coord[2], coord[3], sinPsi1Phi2Phi3m);
503 coord[3] = kmpmq;
504 flowHist->Fill(coord[0], coord[2], coord[3], (mp*mult-2.*mq)*(mult-1.));
505 coord[3] = kCospsi1phi2phi3p;
506 flowHist->Fill(coord[0], coord[2], coord[3], cosPsi1Phi2Phi3p);
507 coord[3] = kSinpsi1phi2phi3p;
508 flowHist->Fill(coord[0], coord[2], coord[3], sinPsi1Phi2Phi3p);
d2bea14e 509
510 }
511
512}
513//_____________________________________________________________________
514void AliForwardFlowTaskQC::Terminate(Option_t */*option*/)
515{
516 //
517 // End of job
518 // Finalizes the Q cumulant calculations
519 //
520 // Parameters:
521 // option Not used
522 //
2f9be372 523
d226802c 524 TH3D* cumulantsHist = 0;
2f9be372 525 TProfile* cumulant2diffHist = 0;
526 TProfile* cumulant4diffHist = 0;
527 TList* tList = 0;
528 TList* vList = 0;
58f5fae2 529
530 // For flow calculations
b4db5d3b 531 Double_t two = 0, qc2 = 0, /* vnTwo = 0, */ four = 0, qc4 = 0 /*, vnFour = 0*/;
2f9be372 532 Double_t twoPrime = 0, qc2Prime = 0, vnTwoDiff = 0, fourPrime = 0, qc4Prime = 0, vnFourDiff = 0;
533 Double_t w2 = 0, w4 = 0, w2p = 0, w4p = 0;
2f9be372 534 Double_t w2Two = 0, w2pTwoPrime = 0, w4Four = 0, w4pFourPrime = 0;
58f5fae2 535 Double_t cosP1nPhi = 0, sinP1nPhi = 0, mult = 0, cosP1nPhi1P1nPhi2 = 0, sinP1nPhi1P1nPhi2 = 0;
536 Double_t cosP1nPhi1M1nPhi2M1nPhi3 = 0, sinP1nPhi1M1nPhi2M1nPhi3 = 0, multm1m2 = 0;
537 Double_t cosP1nPsi = 0, sinP1nPsi = 0, mp = 0, cosP1nPsi1P1nPhi2 = 0, sinP1nPsi1P1nPhi2 = 0;
538 Double_t cosP1nPsi1M1nPhi2M1nPhi3 = 0, sinP1nPsi1M1nPhi2M1nPhi3 = 0, mpqMult = 0;
539 Double_t cosP1nPsi1P1nPhi2M1nPhi3 = 0, sinP1nPsi1P1nPhi2M1nPhi3 = 0;
2f9be372 540
541 Int_t nLoops = (fMC ? 5 : 2);
58f5fae2 542 TString type = "";
d2bea14e 543
544 // Do a loop over the difference analysis types, calculating flow
58f5fae2 545 // 2 loops for real data, 3 for MC data
d2bea14e 546 // inside each is a nested loop over each harmonic (1, 2, 3 and 4 at the moment)
547 for (Int_t loop = 1; loop <= nLoops; loop++) {
548
2f9be372 549 if (loop == 1) type = "FMD";
58f5fae2 550 if (loop == 2) type = "SPD";
551 if (loop == 3) type = "MC";
2f9be372 552 if (loop == 4) type = "FMDTR";
553 if (loop == 5) type = "SPDTR";
d2bea14e 554
d226802c 555 for (Int_t n = 1; n <= 6; n++) {
d2bea14e 556 if (!fv[n]) continue;
2f9be372 557
558 tList = (TList*)fOutputList->FindObject(type.Data());
559 vList = (TList*)tList->FindObject(Form("v%d", n));
d226802c 560 TString tag = "";
d2bea14e 561
58f5fae2 562 // Centrality loop
d226802c 563 for (Int_t c = 0; c < 100; c++) {
2f9be372 564
565 Double_t nEv = 0;
d226802c 566 tag = Form("%d_%d", c, c+1);
567 cumulantsHist = (TH3D*)vList->FindObject(Form("hQ%dCumuHist%s_%s", n, type.Data(), tag.Data()));
568
569 cumulant2diffHist = (TProfile*)vList->FindObject(Form("hQ%dCumulant2DiffFlow%s_%s", n, type.Data(), tag.Data()));
570 cumulant4diffHist = (TProfile*)vList->FindObject(Form("hQ%dCumulant4DiffFlow%s_%s", n, type.Data(), tag.Data()));
2f9be372 571
572 for (Int_t vertexBin = 1; vertexBin <= cumulantsHist->GetNbinsY(); vertexBin++) {
573
2f9be372 574 for (Int_t etaBin = 1; etaBin <= cumulantsHist->GetNbinsX(); etaBin++) {
d226802c 575 Double_t eta = cumulantsHist->GetXaxis()->GetBinCenter(etaBin);
576 // 2-particle reference flow
577 w2Two = cumulantsHist->GetBinContent(etaBin, vertexBin, kW2Two);
578 if (!w2Two) continue;
579 w2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kW2);
580 cosP1nPhi = cumulantsHist->GetBinContent(etaBin, vertexBin, kQnRe);
581 sinP1nPhi = cumulantsHist->GetBinContent(etaBin, vertexBin, kQnIm);
582 mult = cumulantsHist->GetBinContent(etaBin, vertexBin, kM);
2f9be372 583
d226802c 584 cosP1nPhi /= mult;
585 sinP1nPhi /= mult;
586 two = w2Two / w2;
587 qc2 = two - TMath::Power(cosP1nPhi, 2) - TMath::Power(sinP1nPhi, 2);
588 if (qc2 <= 0) continue;
b4db5d3b 589 // vnTwo = TMath::Sqrt(qc2);
d226802c 590 // if (!TMath::IsNaN(vnTwo*mult))
591 // cumulant2diffHist->Fill(eta, vnTwo, cumulantsHist->GetBinContent(0,vertexBin,0));
592
593 // 4-particle reference flow
594 w4Four = cumulantsHist->GetBinContent(etaBin, vertexBin, kW4Four);
595 w4 = cumulantsHist->GetBinContent(etaBin, vertexBin, kW4);
596 cosP1nPhi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kCosphi1phi2);
597 sinP1nPhi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kSinphi1phi2);
598 cosP1nPhi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kCosphi1phi2phi3m);
599 sinP1nPhi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kSinphi1phi2phi3m);
600 multm1m2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kMm1m2);
601
602 cosP1nPhi1P1nPhi2 /= w2;
603 sinP1nPhi1P1nPhi2 /= w2;
604 cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
605 sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
606 four = w4Four / w4;
607 qc4 = four-2.*TMath::Power(two,2.)
608 - 4.*cosP1nPhi*cosP1nPhi1M1nPhi2M1nPhi3
609 + 4.*sinP1nPhi*sinP1nPhi1M1nPhi2M1nPhi3-TMath::Power(cosP1nPhi1P1nPhi2,2.)-TMath::Power(sinP1nPhi1P1nPhi2,2.)
610 + 4.*cosP1nPhi1P1nPhi2*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.))
611 + 8.*sinP1nPhi1P1nPhi2*sinP1nPhi*cosP1nPhi
612 + 8.*two*(TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.))
613 - 6.*TMath::Power((TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.)),2.);
2f9be372 614
d226802c 615 if (qc4 >= 0) continue;
b4db5d3b 616 // vnFour = TMath::Power(-qc4, 0.25);
d226802c 617 // if (!TMath::IsNaN(vnFour*mult))
618 // cumulant4diffHist->Fill(eta, vnFour, cumulantsHist->GetBinContent(0,vertexBin,0));
2f9be372 619
620 // 2-particle differential flow
d226802c 621 w2pTwoPrime = cumulantsHist->GetBinContent(etaBin, vertexBin, kw2two);
2f9be372 622 if (!w2pTwoPrime) continue;
d226802c 623 w2p = cumulantsHist->GetBinContent(etaBin, vertexBin, kw2);
624 cosP1nPsi = cumulantsHist->GetBinContent(etaBin, vertexBin, kpnRe);
625 sinP1nPsi = cumulantsHist->GetBinContent(etaBin, vertexBin, kpnIm);
626 mp = cumulantsHist->GetBinContent(etaBin, vertexBin, kmp);
627
2f9be372 628 cosP1nPsi /= mp;
629 sinP1nPsi /= mp;
630 twoPrime = w2pTwoPrime / w2p;
631 qc2Prime = twoPrime - sinP1nPsi*sinP1nPhi - cosP1nPsi*cosP1nPhi;
632
633 vnTwoDiff = qc2Prime / TMath::Sqrt(qc2);
d226802c 634 if (!TMath::IsNaN(vnTwoDiff*mp) && vnTwoDiff > 0) cumulant2diffHist->Fill(eta, vnTwoDiff, cumulantsHist->GetBinContent(0,vertexBin,0));
2f9be372 635
636 // 4-particle differential flow
d226802c 637 w4pFourPrime = cumulantsHist->GetBinContent(etaBin, vertexBin, kw4four);
638 w4p = cumulantsHist->GetBinContent(etaBin, vertexBin, kw4);
639 cosP1nPsi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kCospsi1phi2);
640 sinP1nPsi1P1nPhi2 = cumulantsHist->GetBinContent(etaBin, vertexBin, kSinpsi1phi2);
641 cosP1nPsi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kCospsi1phi2phi3m);
642 sinP1nPsi1M1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kSinpsi1phi2phi3m);
643 cosP1nPsi1P1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kCospsi1phi2phi3p);
644 sinP1nPsi1P1nPhi2M1nPhi3 = cumulantsHist->GetBinContent(etaBin, vertexBin, kSinpsi1phi2phi3p);
645 mpqMult = cumulantsHist->GetBinContent(etaBin, vertexBin, kmpmq);
646
2f9be372 647 cosP1nPsi1P1nPhi2 /= w2p;
648 sinP1nPsi1P1nPhi2 /= w2p;
649 cosP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
650 sinP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
651 cosP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
652 sinP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
653 fourPrime = w4pFourPrime / w4p;
d226802c 654
2f9be372 655 qc4Prime = fourPrime-2.*twoPrime*two
656 - cosP1nPsi*cosP1nPhi1M1nPhi2M1nPhi3
657 + sinP1nPsi*sinP1nPhi1M1nPhi2M1nPhi3
658 - cosP1nPhi*cosP1nPsi1M1nPhi2M1nPhi3
659 + sinP1nPhi*sinP1nPsi1M1nPhi2M1nPhi3
660 - 2.*cosP1nPhi*cosP1nPsi1P1nPhi2M1nPhi3
661 - 2.*sinP1nPhi*sinP1nPsi1P1nPhi2M1nPhi3
662 - cosP1nPsi1P1nPhi2*cosP1nPhi1P1nPhi2
663 - sinP1nPsi1P1nPhi2*sinP1nPhi1P1nPhi2
664 + 2.*cosP1nPhi1P1nPhi2*(cosP1nPsi*cosP1nPhi-sinP1nPsi*sinP1nPhi)
665 + 2.*sinP1nPhi1P1nPhi2*(cosP1nPsi*sinP1nPhi+sinP1nPsi*cosP1nPhi)
666 + 4.*two*(cosP1nPsi*cosP1nPhi+sinP1nPsi*sinP1nPhi)
667 + 2.*cosP1nPsi1P1nPhi2*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.))
668 + 4.*sinP1nPsi1P1nPhi2*cosP1nPhi*sinP1nPhi
669 + 4.*twoPrime*(TMath::Power(cosP1nPhi,2.)+TMath::Power(sinP1nPhi,2.))
670 - 6.*(TMath::Power(cosP1nPhi,2.)-TMath::Power(sinP1nPhi,2.))
671 * (cosP1nPsi*cosP1nPhi-sinP1nPsi*sinP1nPhi)
672 - 12.*cosP1nPhi*sinP1nPhi
673 * (sinP1nPsi*cosP1nPhi+cosP1nPsi*sinP1nPhi);
d2bea14e 674
2f9be372 675 vnFourDiff = - qc4Prime / TMath::Power(-qc4, 0.75);
d226802c 676 if (!TMath::IsNaN(vnFourDiff*mp) && vnFourDiff > 0) cumulant4diffHist->Fill(eta, vnFourDiff, cumulantsHist->GetBinContent(0,vertexBin,0));
2f9be372 677 } // End of etaBin loop
678 nEv += cumulantsHist->GetBinContent(0,vertexBin,0);
679 } // End of vertexBin loop
58f5fae2 680 // Number of events:
2f9be372 681 cumulant2diffHist->Fill(7., nEv);
682 cumulant4diffHist->Fill(7., nEv);
683
58f5fae2 684 } // End of centrality loop
d2bea14e 685 } // End of harmonics loop
686 } // End of type loop
58f5fae2 687
d226802c 688 PostData(1, fOutputList);
689
690 return;
d2bea14e 691}
692// _____________________________________________________________________
693void AliForwardFlowTaskQC::ProcessPrimary()
694{
695 //
696 // If fMC == kTRUE this function takes care of organizing the input
697 // Monte Carlo data and histograms so AliForwardFlowTaskQC::QCumulants
698 // can be run on it.
699 //
700
2f9be372 701 if (fFlowUtil->LoopAODFMDtrrefHits(fAOD)) {
58f5fae2 702 // Run analysis on TrackRefs
d226802c 703 for (Int_t n = 1; n <= 6; n++) {
58f5fae2 704 if (fv[n])
2f9be372 705 CumulantsMethod("FMDTR", n);
706 }
707 }
708
709 if (fFlowUtil->LoopAODSPDtrrefHits(fAOD)) {
710 // Run analysis on SPD TrackRefs
d226802c 711 for (Int_t n = 1; n <= 6; n++) {
2f9be372 712 if (fv[n])
713 CumulantsMethod("SPDTR", n);
58f5fae2 714 }
d2bea14e 715 }
716
58f5fae2 717 if (fFlowUtil->LoopAODmc(fAOD, fAddFlow, fAddType, fAddOrder)) {
718 // Run analysis on MC truth
d226802c 719 for (Int_t n = 1; n <= 6; n++) {
58f5fae2 720 if (fv[n])
721 CumulantsMethod("MC", n);
722 }
d2bea14e 723 }
724
d2bea14e 725}
726//_____________________________________________________________________
727//
728//
729// EOF