Fixed warnings [-Wunused-but-set-variable] from GCC 4.6 -
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardFlowTaskQC.cxx
1 //
2 // Calculate flow in the forward and central regions using the Q cumulants method.
3 //
4 // Inputs:
5 //  - AliAODEvent
6 //
7 // Outputs:
8 //  - AnalysisResults.root
9 //
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>
18 #include "THnSparse.h"
19 #include "TH3D.h"
20 #include "TProfile2D.h"
21 #include "AliLog.h"
22 #include "AliForwardFlowTaskQC.h"
23 #include "AliAnalysisManager.h"
24 #include "AliAODHandler.h"
25 #include "AliAODInputHandler.h"
26 #include "AliAODMCParticle.h"
27 #include "AliAODForwardMult.h"
28 #include "AliAODEvent.h"
29
30 //
31 // Enumeration for adding and retrieving stuff from the histogram
32 //
33 enum { 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 };
38
39 ClassImp(AliForwardFlowTaskQC)
40 #if 0
41 ; // For emacs 
42 #endif
43
44 AliForwardFlowTaskQC::AliForwardFlowTaskQC()
45   : fOutputList(0),     // Output list
46     fFlowUtil(0),       // AliForwardFlowUtil
47     fAOD(0),            // AOD input event
48     fMC(kFALSE),        // MC flag
49     fEtaBins(48),       // # of etaBin bins in histograms
50     fEtaRef(12),        // # of Eta bins for reference flow
51     fAddFlow(0),        // Add flow string
52     fAddType(0),        // Add flow type #
53     fAddOrder(0),       // Add flow order
54     fZvertex(1111),     // Z vertex range
55     fCent(-1)           // Centrality
56 {
57   // 
58   // Default constructor
59   //
60 }
61 //_____________________________________________________________________
62 AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name) :
63   AliAnalysisTaskSE(name),
64   fOutputList(0),       // Output list
65   fFlowUtil(0),         // AliForwardFlowUtil
66   fAOD(0),              // AOD input event
67   fMC(kFALSE),          // MC flag
68   fEtaBins(48),         // # of Eta bins
69   fEtaRef(12),          // # of Eta bins for reference flow
70   fAddFlow(0),          // Add flow string
71   fAddType(0),          // Add flow type #
72   fAddOrder(0),         // Add flow order
73   fZvertex(1111),       // Z vertex range
74   fCent(-1)             // Centrality
75 {
76   // 
77   // Constructor
78   //
79   // Parameters:
80   //  name: Name of task
81   //
82   for (Int_t n = 0; n <= 6; n++) fv[n] = kTRUE;
83   DefineOutput(1, TList::Class());
84 }
85 //_____________________________________________________________________
86 AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o) :
87   AliAnalysisTaskSE(o),
88   fOutputList(o.fOutputList),   // Output list
89   fFlowUtil(o.fFlowUtil),       // AliForwardFlowUtil
90   fAOD(o.fAOD),                 // AOD input event
91   fMC(o.fMC),                   // MC flag
92   fEtaBins(o.fEtaBins),         // # of Eta bins
93   fEtaRef(o.fEtaRef),           // # of Eta bins for reference flow
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
98   fCent(o.fCent)                // Centrality
99 {
100   // 
101   // Copy constructor 
102   // 
103   // Parameters:
104   //    o Object to copy from 
105   //
106   for (Int_t n = 0; n <= 6; n++) fv[n] = o.fv[n];
107   DefineOutput(1, TList::Class());
108 }
109 //_____________________________________________________________________
110 void AliForwardFlowTaskQC::UserCreateOutputObjects()
111 {
112   //
113   // Create output objects
114   // 
115   if (!fOutputList)
116     fOutputList = new TList();
117   fOutputList->SetName("QCumulants");
118   fOutputList->SetOwner();
119   fFlowUtil = new AliForwardFlowUtil(fOutputList);
120   
121   Double_t x[101] = { 0. };
122   // Double_t etaMin = -6;
123   // Double_t etaMax = 6;
124
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:
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;
135     // etaMin = 0;
136     // etaMax = 6;
137   }
138
139   else if (fEtaBins == 100) { 
140     for (Int_t e = 0; e<= 100; e++) {
141       x[e] = e;
142     }
143     // etaMin = 0;
144     // etaMax = 100;
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   }
169
170   // Histograms for cumulants analysis
171
172   // We loop over flow histograms here to add different orders of harmonics
173   TString type = "";
174   for (Int_t loop = 1; loop <= 5; loop++) {
175     
176     if (loop == 1) type = "FMD";
177     if (loop == 2) type = "SPD";
178     if (loop == 3) type = "MC";
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);
185
186     for (Int_t n = 1; n <= 6; n++) {
187       if (!fv[n]) continue;
188       
189       TList* vList = new TList();
190       vList->SetName(Form("v%d", n));
191       tList->Add(vList);
192       
193       TString tag = TString();
194       for (Int_t c = 0; c < 100; c++) {
195         // Output histograms
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
205         TProfile* hCumulant2DiffFlow = new TProfile(Form("hQ%dCumulant2DiffFlow%s_%s", n, type.Data(), tag.Data()), 
206                      Form("hQ%dCumulant2DiffFlow%s_%s", n, type.Data(), tag.Data()), fEtaBins, x);
207         hCumulant2DiffFlow->Sumw2();
208         vList->Add(hCumulant2DiffFlow);
209    
210         TProfile* hCumulant4DiffFlow = new TProfile(Form("hQ%dCumulant4DiffFlow%s_%s", n, type.Data(), tag.Data()), 
211                      Form("hQ%dCumulant4DiffFlow%s_%s", n, type.Data(), tag.Data()), fEtaBins, x);
212         hCumulant4DiffFlow->Sumw2();
213         vList->Add(hCumulant4DiffFlow);
214       } // end of centrality loop
215
216     } // end of v_{n} loop
217  
218     // Single Event histograms
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));
221     hdNdphiSE->Sumw2();
222     fOutputList->Add(hdNdphiSE);
223
224     TH2D* hdNdetadphiSE = new TH2D(Form("hdNdetadphiSE%s", type.Data()), 
225                  Form("hdNdetadphiSE%s", type.Data()), fEtaBins, x, (type.Contains("MC") ? 200 : 20), (type.Contains("MC") ? phiMC : phi));
226     hdNdetadphiSE->Sumw2();
227     fOutputList->Add(hdNdetadphiSE);
228   } // end of type loop
229
230   TProfile2D* pMCTruth = new TProfile2D("pMCTruth", "pMCTruth", 48, -6, 6, 100, 0, 100);
231   TAxis* xAxis = pMCTruth->GetXaxis();
232   xAxis->Set(fEtaBins, x);
233   pMCTruth->Sumw2();
234   fOutputList->Add(pMCTruth);
235
236   // Monitoring plots
237   TH1D* cent = new TH1D("Centralities", "Centralities", 100, 0, 100);
238   fOutputList->Add(cent);
239
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);
247  
248   PostData(1, fOutputList);
249 }
250 //_____________________________________________________________________
251 void AliForwardFlowTaskQC::UserExec(Option_t */*option*/)
252 {
253   // 
254   // Process each event
255   //
256   // Parameters:
257   //  option: Not used
258   //
259
260   // Reset data members
261   fCent = -1;
262   fZvertex = 1111;
263
264   // Get input event
265   fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
266   if (!fAOD) return;
267
268   // Fill histograms
269   if (!fFlowUtil->LoopAODFMD(fAOD)) return;
270   if (!fFlowUtil->LoopAODSPD(fAOD)) return;
271
272   // Get centrality and vertex from flow utility and fill monitoring histograms
273   fCent = fFlowUtil->GetCentrality();
274   fZvertex = fFlowUtil->GetVertex();
275
276   TH1D* cent = (TH1D*)fOutputList->FindObject("Centralities");
277
278   cent->Fill(fCent);
279   TH1D* vertex = (TH1D*)fOutputList->FindObject("VertexSelected");
280   vertex->Fill(fZvertex);
281
282   // Run analysis on FMD and SPD
283   for (Int_t n = 1; n <= 6; n++) {
284     if (fv[n]) {
285       CumulantsMethod("FMD", n);
286       CumulantsMethod("SPD", n);
287     }
288   }
289
290   // And do analysis if there is.
291   if (fMC) { 
292     ProcessPrimary();
293   }
294
295   PostData(1, fOutputList);
296 }
297 //_____________________________________________________________________
298 void AliForwardFlowTaskQC::CumulantsMethod(TString type = "", 
299                                            Int_t harmonic = 2)
300 {
301   // 
302   // Calculate the Q cumulant of order n
303   //
304   // Parameters:
305   //  type: Determines which histograms should be used
306   //        - "FMD" or "SPD" = data histograms
307   //        - "FMDTR" or "SPDTR" = track reference histograms
308   //        - "MC" = MC truth histograms
309   //  harmonic: Which harmonic to calculate
310   //
311   Double_t n = harmonic;
312   
313   TList* tList = (TList*)fOutputList->FindObject(type.Data());
314   TList* vList = (TList*)tList->FindObject(Form("v%d", harmonic));
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));
317   TH2D* dNdetadphi = (TH2D*)fOutputList->FindObject(Form("hdNdetadphiSE%s", type.Data()));
318   TH2D* dNdphi = (TH2D*)fOutputList->FindObject(Form("hdNdphiSE%s", /*"FMD"*/type.Data()));
319   
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. };
322   // We create the objects needed for the analysis
323   Double_t dQnRe = 0, dQ2nRe = 0, dQnIm = 0, dQ2nIm = 0, mult = 0;
324   Double_t pnRe = 0, p2nRe = 0, qnRe = 0, q2nRe = 0, pnIm = 0, p2nIm = 0, qnIm = 0, q2nIm = 0;
325   Double_t two = 0, four = 0, twoPrime = 0, fourPrime = 0;
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;
331   Double_t multi = 0, multp = 0, mp = 0, mq = 0;
332   Double_t w2 = 0, w4 = 0, w2p = 0, w4p = 0;
333   Int_t refEtaBin = 0;
334   Bool_t kEventCount = kFALSE;
335
336   // We loop over the data 1 time!
337   for (Int_t etaBin = 1; etaBin <= dNdetadphi->GetNbinsX(); etaBin++) {
338     eta = dNdetadphi->GetXaxis()->GetBinCenter(etaBin);
339     refEtaBin = dNdphi->GetXaxis()->FindBin(eta);
340     // The values for each individual etaBin bins are reset
341     mp = 0;
342     pnRe = 0;
343     p2nRe = 0;
344     pnIm = 0;
345     p2nIm = 0;
346
347     mult = 0;
348     dQnRe = 0;
349     dQnIm = 0;
350     dQ2nRe = 0;
351     dQ2nIm = 0;
352
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;
363       
364       // For each etaBin bin the necessary values for differential flow
365       // is calculated. Here is the loop over the phi's.
366       multp = dNdetadphi->GetBinContent(etaBin, phiBin);
367       mp += multp;
368       pnRe += multp*TMath::Cos(n*phi);
369       pnIm += multp*TMath::Sin(n*phi);
370       p2nRe += multp*TMath::Cos(2.*n*phi);
371       p2nIm += multp*TMath::Sin(2.*n*phi);
372     }
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     } 
382
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);
393
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.));
408
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);
413
414     cosPhi1Phi2 = dQnRe*dQnRe - dQnIm*dQnIm - dQ2nRe;
415     sinPhi1Phi2 = 2.*dQnRe*dQnIm - dQ2nIm;
416       
417     cosPhi1Phi2Phi3m = dQnRe*(TMath::Power(dQnRe,2)+TMath::Power(dQnIm,2))-dQnRe*dQ2nRe-dQnIm*dQ2nIm-2.*(mult-1)*dQnRe;
418
419     sinPhi1Phi2Phi3m = -dQnIm*(TMath::Power(dQnRe,2)+TMath::Power(dQnIm,2))+dQnRe*dQ2nIm-dQnIm*dQ2nRe+2.*(mult-1)*dQnIm; 
420
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.));
431
432     // Differential flow calculations for each etaBin bin is done:
433     mq = mp;
434     qnRe = pnRe;
435     qnIm = pnIm;
436     q2nRe = p2nRe;
437     q2nIm = p2nIm;
438
439     // 2-particle differential flow
440     w2p = mp * mult - mq;
441     twoPrime = pnRe*dQnRe + pnIm*dQnIm - mq;
442     
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);
454
455     // 4-particle differential flow
456     w4p = (mp * mult - 3.*mq)*(mult - 1.)*(mult - 2.);
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);
475
476     cosPsi1Phi2 = pnRe*dQnRe - pnIm*dQnIm - q2nRe;
477     sinPsi1Phi2 = pnRe*dQnIm + pnIm*dQnRe - q2nIm;
478
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; 
486
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); 
509
510   }
511
512 }
513 //_____________________________________________________________________
514 void 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   //
523
524   TH3D* cumulantsHist = 0;
525   TProfile* cumulant2diffHist = 0;
526   TProfile* cumulant4diffHist = 0;
527   TList* tList = 0;
528   TList* vList = 0;
529
530   // For flow calculations
531   Double_t two = 0, qc2 = 0, /* vnTwo = 0, */ four = 0, qc4 = 0 /*, vnFour = 0*/; 
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;
534   Double_t w2Two = 0, w2pTwoPrime = 0, w4Four = 0, w4pFourPrime = 0;
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;
540
541   Int_t nLoops = (fMC ? 5 : 2);
542   TString type = "";
543
544   // Do a loop over the difference analysis types, calculating flow
545   // 2 loops for real data, 3 for MC data
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
549     if (loop == 1) type = "FMD";
550     if (loop == 2) type = "SPD";
551     if (loop == 3) type = "MC";
552     if (loop == 4) type = "FMDTR";
553     if (loop == 5) type = "SPDTR";
554     
555     for (Int_t n = 1; n <= 6; n++) {
556       if (!fv[n]) continue;
557      
558       tList = (TList*)fOutputList->FindObject(type.Data());
559       vList = (TList*)tList->FindObject(Form("v%d", n));
560       TString tag = "";
561  
562       // Centrality loop
563       for (Int_t c = 0; c < 100; c++) {
564
565         Double_t nEv = 0;
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()));
571         
572         for (Int_t vertexBin = 1; vertexBin <= cumulantsHist->GetNbinsY(); vertexBin++) {
573         
574           for (Int_t etaBin = 1; etaBin <= cumulantsHist->GetNbinsX(); etaBin++) {
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);
583               
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;
589             // vnTwo = TMath::Sqrt(qc2);
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.);
614             
615             if (qc4 >= 0) continue;
616             // vnFour = TMath::Power(-qc4, 0.25);
617        //     if (!TMath::IsNaN(vnFour*mult)) 
618        //         cumulant4diffHist->Fill(eta, vnFour, cumulantsHist->GetBinContent(0,vertexBin,0));
619   
620             // 2-particle differential flow
621             w2pTwoPrime = cumulantsHist->GetBinContent(etaBin, vertexBin, kw2two);
622             if (!w2pTwoPrime) continue;
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
628             cosP1nPsi /= mp;
629             sinP1nPsi /= mp;
630             twoPrime = w2pTwoPrime / w2p;
631             qc2Prime = twoPrime - sinP1nPsi*sinP1nPhi - cosP1nPsi*cosP1nPhi;
632
633             vnTwoDiff = qc2Prime / TMath::Sqrt(qc2);
634             if (!TMath::IsNaN(vnTwoDiff*mp) && vnTwoDiff > 0) cumulant2diffHist->Fill(eta, vnTwoDiff, cumulantsHist->GetBinContent(0,vertexBin,0));
635
636             // 4-particle differential flow
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
647             cosP1nPsi1P1nPhi2 /= w2p;
648             sinP1nPsi1P1nPhi2 /= w2p;
649             cosP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
650             sinP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
651             cosP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
652             sinP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
653             fourPrime = w4pFourPrime / w4p;
654
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);
674   
675             vnFourDiff = - qc4Prime / TMath::Power(-qc4, 0.75);
676             if (!TMath::IsNaN(vnFourDiff*mp) && vnFourDiff > 0) cumulant4diffHist->Fill(eta, vnFourDiff, cumulantsHist->GetBinContent(0,vertexBin,0));
677           } // End of etaBin loop
678           nEv += cumulantsHist->GetBinContent(0,vertexBin,0);
679         } // End of vertexBin loop
680         // Number of events:
681         cumulant2diffHist->Fill(7., nEv);
682         cumulant4diffHist->Fill(7., nEv);
683   
684       } // End of centrality loop
685     } // End of harmonics loop
686   }  // End of type loop
687
688   PostData(1, fOutputList);
689   
690   return;
691 }
692 // _____________________________________________________________________
693 void 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
701   if (fFlowUtil->LoopAODFMDtrrefHits(fAOD)) {
702    // Run analysis on TrackRefs
703     for (Int_t n = 1; n <= 6; n++) {
704       if (fv[n])
705         CumulantsMethod("FMDTR", n);
706     }
707   }
708
709   if (fFlowUtil->LoopAODSPDtrrefHits(fAOD)) {
710     // Run analysis on SPD TrackRefs
711     for (Int_t n = 1; n <= 6; n++) {
712       if (fv[n])
713         CumulantsMethod("SPDTR", n);
714     }
715   }
716
717   if (fFlowUtil->LoopAODmc(fAOD, fAddFlow, fAddType, fAddOrder)) {
718     // Run analysis on MC truth
719     for (Int_t n = 1; n <= 6; n++) {
720       if (fv[n])
721         CumulantsMethod("MC", n);
722     }
723   }
724
725 }
726 //_____________________________________________________________________
727 //
728 //
729 // EOF