PWG2/FORWARD/analysis2/AliForwardFlowTaskQC.cxx

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