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