PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.cxx

   1
   2 //
   3 // Calculate flow in the forward and central regions using the Q cumulants method.
   4 //
   5 // Inputs:
   6 //  - AliAODEvent
   7 //
   8 // Outputs:
   9 //  - AnalysisResults.root or forward_flow.root
  10 //
  11 #include <TROOT.h>
  12 #include <TSystem.h>
  13 #include <TInterpreter.h>
  14 #include <TChain.h>
  15 #include <TFile.h>
  16 #include <TList.h>
  17 #include <TMath.h>
  18 #include <TH3D.h>
  19 #include <TProfile2D.h>
  20 #include <TParameter.h>
  21 #include <TMatrixD.h>
  22 #include <TVectorD.h>
  23 #include <TGraph.h>
  24 #include "AliLog.h"
  25 #include "AliForwardFlowTaskQC.h"
  26 #include "AliAnalysisManager.h"
  27 #include "AliAODHandler.h"
  28 #include "AliAODInputHandler.h"
  29 #include "AliAODForwardMult.h"
  30 #include "AliAODCentralMult.h"
  31 #include "AliAODEvent.h"
  32 #include "AliForwardUtil.h"
  33 #include "AliVVZERO.h"
  34 #include "AliAODVertex.h"
  35 #include "AliCentrality.h"
  36 #include "AliESDEvent.h"
  37 #include "AliVTrack.h"
  38 #include "AliESDtrack.h"
  39 #include "AliAODTrack.h"
  40 #include "AliAnalysisFilter.h"
  41 #include "AliAODMCHeader.h"
  42 #include "AliForwardFlowWeights.h"
  43
  44 ClassImp(AliForwardFlowTaskQC)
  45 #if 0
  46 ; // For emacs
  47 #endif
  48
  49 AliForwardFlowTaskQC::AliForwardFlowTaskQC()
  50   : AliAnalysisTaskSE(),
  51     fVtxAxis(),          // Axis to control vertex binning
  52     fCentAxis(),         // Axis to control centrality/multiplicity binning
  53     fFMDCut(-1),         // FMD sigma cut
  54     fSPDCut(-1),         // SPD sigma cut
  55     fFlowFlags(0),       // Flow flags
  56     fEtaGap(0.),         // Eta gap value
  57     fBinsForward(),      // List with forward flow hists
  58     fBinsCentral(),      // List with central flow hists
  59     fSumList(0),         // Event sum list
  60     fOutputList(0),      // Result output list
  61     fAOD(0),             // AOD input event
  62     fTrackCuts(0),    // ESD track cuts
  63     fMaxMoment(0),       // Max flow moment
  64     fVtx(1111),          // Z vertex coordinate
  65     fCent(-1),           // Centrality
  66     fHistdNdedpV0(),     // Hist for v0
  67     fHistdNdedp3Cor(),   // Hist for combining detectors
  68     fHistFMDSPDCorr(),   // FMD SPD correlation
  69     fHistCent(),         // Hist for centrality
  70     fHistVertexSel(),    // Hist for selected vertices
  71     fHistEventSel()      // Hist for event selection
  72 {
  73   //
  74   //  Default constructor
  75   //
  76 }
  77 //_____________________________________________________________________
  78 AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name)
  79   : AliAnalysisTaskSE(name),
  80     fVtxAxis(),         // Axis to control vertex binning
  81     fCentAxis(),        // Axis to control centrality/multiplicity binning
  82     fFMDCut(-1),        // FMD sigma cut
  83     fSPDCut(-1),        // SPD sigma cut
  84     fFlowFlags(0x0),    // Flow flags
  85     fEtaGap(0.),        // Eta gap value
  86     fBinsForward(),     // List with forward flow hists
  87     fBinsCentral(),     // List with central flow hists
  88     fSumList(0),        // Event sum list
  89     fOutputList(0),     // Result output list
  90     fAOD(0),            // AOD input event
  91     fTrackCuts(0),      // ESD track cuts
  92     fMaxMoment(4),      // Max flow moment
  93     fVtx(1111),         // Z vertex coordinate
  94     fCent(-1),          // Centrality
  95     fHistdNdedpV0(),    // Histo for v0
  96     fHistdNdedp3Cor(),  // Histo for combining detectors
  97     fHistFMDSPDCorr(),  // FMD SPD correlation
  98     fHistCent(),        // Hist for centrality
  99     fHistVertexSel(),   // Hist for selected vertices
 100     fHistEventSel()     // Hist for event selection
 101 {
 102   //
 103   //  Constructor
 104   //
 105   //  Parameters:
 106   //   name: Name of task
 107   //
 108   DefineOutput(1, TList::Class());
 109   DefineOutput(2, TList::Class());
 110 }
 111 //_____________________________________________________________________
 112 AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o)
 113   : AliAnalysisTaskSE(o),
 114     fVtxAxis(o.fVtxAxis),              // Axis to control vertex binning
 115     fCentAxis(o.fCentAxis),            // Array to control centrality/multiplicity binning
 116     fFMDCut(o.fFMDCut),                // FMD sigma cut
 117     fSPDCut(o.fSPDCut),                // SPD sigma cut
 118     fFlowFlags(o.fFlowFlags),          // Flow flags
 119     fEtaGap(o.fEtaGap),                // Eta gap value
 120     fBinsForward(),                    // List with forward flow hists
 121     fBinsCentral(),                    // List with central flow hists
 122     fSumList(o.fSumList),              // Event sum list
 123     fOutputList(o.fOutputList),        // Result output list
 124     fAOD(o.fAOD),                      // AOD input event
 125     fTrackCuts(o.fTrackCuts),          // ESD track cuts
 126     fMaxMoment(o.fMaxMoment),          // Flow moments
 127     fVtx(o.fVtx),                      // Z vertex coordinate
 128     fCent(o.fCent),                    // Centrality
 129     fHistdNdedpV0(o.fHistdNdedpV0),    // Histo for v0
 130     fHistdNdedp3Cor(o.fHistdNdedp3Cor),// Histo for combining detectors
 131     fHistFMDSPDCorr(o.fHistFMDSPDCorr),// FMD SPD correlation
 132     fHistCent(o.fHistCent),            // Hist for centrality
 133     fHistVertexSel(o.fHistVertexSel),  // Hist for selected vertices
 134     fHistEventSel(o.fHistEventSel)     // Hist for event selection
 135 {
 136   //
 137   //  Copy constructor
 138   //
 139   //  Parameters:
 140   //   o: Object to copy from
 141   //
 142 }
 143 //_____________________________________________________________________
 144 AliForwardFlowTaskQC&
 145 AliForwardFlowTaskQC::operator=(const AliForwardFlowTaskQC& o)
 146 {
 147   //
 148   //  Assignment operator
 149   //
 150   if (&o == this) return *this;
 151   fVtxAxis        = o.fVtxAxis;
 152   fCentAxis       = o.fCentAxis;
 153   fFMDCut         = o.fFMDCut;
 154   fSPDCut         = o.fSPDCut;
 155   fFlowFlags      = o.fFlowFlags;
 156   fEtaGap         = o.fEtaGap;
 157   fSumList        = o.fSumList;
 158   fOutputList     = o.fOutputList;
 159   fAOD            = o.fAOD;
 160   fTrackCuts      = o.fTrackCuts;
 161   fMaxMoment      = o.fMaxMoment;
 162   fVtx            = o.fVtx;
 163   fCent           = o.fCent;
 164   fHistdNdedpV0   = o.fHistdNdedpV0;
 165   fHistdNdedp3Cor = o.fHistdNdedp3Cor;
 166   fHistFMDSPDCorr = o.fHistFMDSPDCorr;
 167   fHistCent       = o.fHistCent;
 168   fHistVertexSel  = o.fHistVertexSel;
 169   fHistEventSel   = o.fHistEventSel;
 170
 171   return *this;
 172 }
 173 //_____________________________________________________________________
 174 void AliForwardFlowTaskQC::SetFlowFlags(UShort_t flags)
 175 {
 176   //
 177   //  Set flow flags, making sure the detector setup is right
 178   //
 179   //  Parameters:
 180   //   flags: Flow flags
 181   //
 182   if ((flags & kFMD) && (flags & kVZERO))
 183     AliFatal("Cannot do analysis on more than one forward detector!");
 184   else if (!(flags & kFMD) && !(flags & kVZERO))
 185     AliFatal("You need to add a forward detector!");
 186   else fFlowFlags = flags;
 187 }
 188 //_____________________________________________________________________
 189 void AliForwardFlowTaskQC::UserCreateOutputObjects()
 190 {
 191   //
 192   //  Create output objects
 193   //
 194   InitVertexBins();
 195   InitHists();
 196   if ((fFlowFlags & kTracks) && !fTrackCuts) AliFatal("No track cuts set!");
 197   PrintFlowSetup();
 198
 199   PostData(1, fSumList);
 200 }
 201 //_____________________________________________________________________
 202 void AliForwardFlowTaskQC::InitVertexBins()
 203 {
 204   //
 205   //  Init vertexbin objects for forward and central detectors, and add them to the lists
 206   //
 207   for (Int_t v = 1; v <= fVtxAxis->GetNbins(); v++) {
 208     Int_t vL = Int_t(fVtxAxis->GetBinLowEdge(v));
 209     Int_t vH = Int_t(fVtxAxis->GetBinUpEdge(v));
 210     if ((fFlowFlags & kFMD)) {
 211       fBinsForward.Add(new VertexBin(vL, vH, fMaxMoment, "FMD", fFlowFlags, fFMDCut, fEtaGap));
 212       if (!(fFlowFlags & k3Cor))
 213         fBinsCentral.Add(new VertexBin(vL, vH, fMaxMoment, "SPD-FMD", fFlowFlags|kNUAcorr|kSPD, fSPDCut, fEtaGap));
 214     }
 215     else if ((fFlowFlags & kVZERO)) {
 216       fBinsForward.Add(new VertexBin(vL, vH, fMaxMoment, "VZERO", fFlowFlags, 0, fEtaGap));
 217       if ((fFlowFlags & kEtaGap) && !(fFlowFlags & kTracks))
 218         fBinsCentral.Add(new VertexBin(vL, vH, fMaxMoment, "SPD-VZERO", fFlowFlags|kNUAcorr|kSPD, fSPDCut, fEtaGap));
 219     }
 220   }
 221 }
 222 //_____________________________________________________________________
 223 void AliForwardFlowTaskQC::InitHists()
 224 {
 225   //
 226   //  Init histograms and add vertex bin histograms to the sum list
 227   //
 228   if (!fSumList)
 229     fSumList = new TList();
 230   fSumList->SetName("Sums");
 231   fSumList->SetOwner();
 232
 233   if (!fVtxAxis) fVtxAxis = new TAxis(20, -10, 10);
 234   fVtxAxis->SetName("VtxAxis");
 235   if (!fCentAxis) fCentAxis = new TAxis(20, 0, 100);
 236   fVtxAxis->SetName("CentAxis");
 237
 238   fHistCent         = new TH1D("hCent", "Centralities", 100, 0, 100);
 239   fHistVertexSel    = new TH1D("hVertexSel", "Selected vertices", fVtxAxis->GetNbins(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax());
 240   fHistEventSel     = new TH1I("hEventSel", "Event Selection", kOK, 0.5, kOK+0.5);
 241   fHistEventSel->GetXaxis()->SetBinLabel(kNoEvent, "No AOD event");
 242   fHistEventSel->GetXaxis()->SetBinLabel(kNoForward, "No forward det");
 243   fHistEventSel->GetXaxis()->SetBinLabel(kNoCentral, "No central det");
 244   fHistEventSel->GetXaxis()->SetBinLabel(kNoTrigger, "Not triggered");
 245   fHistEventSel->GetXaxis()->SetBinLabel(kNoCent, "No centrality");
 246   fHistEventSel->GetXaxis()->SetBinLabel(kInvCent, "Centrality outside range");
 247   fHistEventSel->GetXaxis()->SetBinLabel(kNoVtx, "No vertex");
 248   fHistEventSel->GetXaxis()->SetBinLabel(kInvVtx, "Vtx outside range");
 249   fHistEventSel->GetXaxis()->SetBinLabel(kOK, "OK!");
 250
 251   fHistFMDSPDCorr = new TH2D("hFMDSPDCorr", "hFMDSPCCorr", 200, 0., 20000., 200, 0, 7500);
 252
 253   TList* dList = new TList();
 254   dList->SetName("Diagnostics");
 255   dList->Add(fHistCent);
 256   dList->Add(fHistVertexSel);
 257   dList->Add(fHistEventSel);
 258   dList->Add(fHistFMDSPDCorr);
 259   fSumList->Add(dList);
 260
 261   fHistdNdedp3Cor = TH2D(Form("hdNdedpCombined_%s", GetQCType(fFlowFlags)), Form("hdNdedpCombined_%s", GetQCType(fFlowFlags)),
 262                            200, -4., 6., 20, 0., TMath::TwoPi());
 263   if ((fFlowFlags & kVZERO)) {
 264     Double_t bins[12] = { -6, -3.7, -3.2, -2.7, -2.2, -1.7,
 265                           2.8, 3.4,  3.9,  4.5,  5.1, 6 };
 266     fHistdNdedpV0 = TH2D(Form("hdNdedpv0%s", GetQCType(fFlowFlags)), Form("hdNdedpv0%s", GetQCType(fFlowFlags)),
 267                     11, -6, 6, 8, 0, TMath::TwoPi());
 268     fHistdNdedpV0.GetXaxis()->Set(11, bins);
 269     if ((fFlowFlags & k3Cor)) {
 270       Double_t bins2[20] = { -6, -3.7, -3.2, -2.7, -2.2, // VZERO
 271                              -2.0, -1.5, -1.0, -0.5 , 0., 0.5, 1.0, 1.5, 2.0, // SPD
 272                             2.8, 3.4,  3.9,  4.5,  5.1, 6 }; // VZERO
 273       fHistdNdedp3Cor.GetXaxis()->Set(19, bins2);
 274       fHistdNdedp3Cor.GetYaxis()->Set(8, 0., TMath::TwoPi());
 275     }
 276   }
 277
 278   TIter nextForward(&fBinsForward);
 279   VertexBin* bin = 0;
 280   while ((bin = static_cast<VertexBin*>(nextForward()))) {
 281     bin->AddOutput(fSumList, fCentAxis);
 282   }
 283   TIter nextCentral(&fBinsCentral);
 284   while ((bin = static_cast<VertexBin*>(nextCentral()))) {
 285     bin->AddOutput(fSumList, fCentAxis);
 286   }
 287 }
 288 //_____________________________________________________________________
 289 void AliForwardFlowTaskQC::UserExec(Option_t */*option*/)
 290 {
 291   //
 292   //  Calls the analyze function - called every event
 293   //
 294   //  Parameters:
 295   //   option: Not used
 296   //
 297
 298   // Reset data members
 299   fCent = -1;
 300   fVtx = 1111;
 301
 302   Analyze();
 303
 304   PostData(1, fSumList);
 305
 306   return;
 307 }
 308 //_____________________________________________________________________
 309 Bool_t AliForwardFlowTaskQC::Analyze()
 310 {
 311   //
 312   //  Load forward and central detector objects from aod tree and call the
 313   //  cumulants calculation for the correct vertex bin
 314   //
 315   //  Return: true on success
 316   //
 317
 318   // Get input event
 319   fAOD = dynamic_cast<AliAODEvent*>(AliForwardUtil::GetAODEvent(this));
 320   if (!fAOD) {
 321     fHistEventSel->Fill(kNoEvent);
 322     return kFALSE;
 323   }
 324
 325   // Get detector objects
 326   AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("Forward"));
 327   AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClusters"));
 328   AliVVZERO* vzero = GetVZERO();
 329   if ((fFlowFlags & kVZERO)) {
 330     if (vzero) {
 331       fHistdNdedpV0.Reset();
 332       FillVZEROHist(vzero);
 333     }
 334   }
 335
 336   // We make sure that the necessary forward object is there
 337   if ((fFlowFlags & kFMD) && !aodfmult) {
 338     fHistEventSel->Fill(kNoForward);
 339     return kFALSE;
 340   }
 341   else if ((fFlowFlags & kVZERO) && !vzero) {
 342     fHistEventSel->Fill(kNoForward);
 343 //    return kFALSE;
 344   }
 345   if (!aodcmult) fHistEventSel->Fill(kNoCentral);
 346
 347    // Check event for triggers, get centrality, vtx etc.
 348   if (!CheckEvent(aodfmult)) return kFALSE;
 349   Int_t vtx = fVtxAxis->FindBin(fVtx)-1;
 350
 351   // Then we assign a reference to the forward histogram of interest
 352   TH2D& forwarddNdedp = ((fFlowFlags & kFMD) ? aodfmult->GetHistogram() : fHistdNdedpV0);
 353   TH2D& spddNdedp = aodcmult->GetHistogram();
 354   if ((fFlowFlags & kStdQC)) {
 355     FillVtxBinList(fBinsForward, forwarddNdedp, vtx);
 356   } else if ((fFlowFlags & kEtaGap)) {
 357     FillVtxBinListEtaGap(fBinsForward, forwarddNdedp, forwarddNdedp, vtx);
 358   }
 359   // At the moment only clusters are supported for the central region (some day add tracks?)
 360   // So no extra checks necessary
 361   if (aodcmult) {
 362     if ((fFlowFlags & kStdQC)) {
 363       FillVtxBinList(fBinsCentral, spddNdedp, vtx);
 364     } else if ((fFlowFlags & kEtaGap)) {
 365       FillVtxBinListEtaGap(fBinsCentral, forwarddNdedp, spddNdedp, vtx);
 366     } else if ((fFlowFlags & k3Cor)) {
 367       FillVtxBinList3Cor(fBinsForward, spddNdedp, forwarddNdedp, vtx);
 368     }
 369     // Diagnostics
 370     if (aodfmult) {
 371       Double_t totForward = forwarddNdedp.Integral(1, forwarddNdedp.GetNbinsX(), 1, forwarddNdedp.GetNbinsY());
 372       Double_t totSPD = spddNdedp.Integral(1, spddNdedp.GetNbinsX(), 1, spddNdedp.GetNbinsY());
 373       fHistFMDSPDCorr->Fill(totForward, totSPD);
 374     }
 375   }
 376
 377   return kTRUE;
 378 }
 379 //_____________________________________________________________________
 380 void AliForwardFlowTaskQC::FillVtxBinList(const TList& list, TH2D& h, Int_t vtx, UShort_t flags) const
 381 {
 382   //
 383   //  Loops over list of VtxBins, fills hists of bins for current vertex
 384   //  and runs analysis on those bins
 385   //
 386   //  Parameters:
 387   //   list:  list of VtxBins
 388   //   h:     dN/detadphi histogram
 389   //   vtx:   current vertex bin
 390   //   flags: extra flags to handle calculations.
 391   //
 392   //   Note: The while loop is used in this function and the next 2 for historical reasons,
 393   //         as originally each moment had it's own VertexBin object.
 394   VertexBin* bin = 0;
 395   Int_t i = 0;
 396   Int_t nVtxBins = fVtxAxis->GetNbins();
 397
 398   while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
 399     i++;
 400     // If no tracks do things normally
 401     if (!(fFlowFlags & kTracks) || (flags & kMC)) {
 402       if (!bin->FillHists(h, fCent, kFillBoth|flags|kReset)) continue;
 403     }
 404     // if tracks things are more complicated
 405     else if ((fFlowFlags & kTracks)) {
 406       if (!FillTracks(bin, kFillRef|kReset|flags)) continue;
 407       if (!bin->FillHists(h, fCent, kFillDiff|kReset|flags)) continue;
 408     }
 409     bin->CumulantsAccumulate(fCent);
 410   }
 411
 412   return;
 413 }
 414 //_____________________________________________________________________
 415 void AliForwardFlowTaskQC::FillVtxBinListEtaGap(const TList& list, TH2D& href,
 416                                                   TH2D& hdiff, Int_t vtx, UShort_t flags) const
 417 {
 418   //
 419   //  Loops over list of VtxBins, fills hists of bins for current vertex
 420   //  and runs analysis on those bins
 421   //
 422   //  Parameters:
 423   //   list:  list of VtxBins
 424   //   href:  dN/detadphi histogram for ref. flow
 425   //   hdiff: dN/detadphi histogram for diff. flow
 426   //   vBin:  current vertex bin
 427   //   flags: extra flags to handle calculations.
 428   //
 429   VertexBin* bin = 0;
 430   Int_t i = 0;
 431   Int_t nVtxBins = fVtxAxis->GetNbins();
 432
 433   while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
 434     i++;
 435     if (!(fFlowFlags & kTracks) || (flags & kMC)) {
 436       if(!bin->FillHists(href, fCent, kFillRef|flags|kReset)) continue;
 437     }
 438     else if ((fFlowFlags & kTracks)) {
 439       if (!FillTracks(bin, kFillRef|kReset|flags)) continue;
 440     }
 441     if (!bin->FillHists(hdiff, fCent, kFillDiff|kReset|flags)) continue;
 442     bin->CumulantsAccumulate(fCent);
 443   }
 444
 445   return;
 446 }
 447 //_____________________________________________________________________
 448 void AliForwardFlowTaskQC::FillVtxBinList3Cor(const TList& list, TH2D& hcent,
 449                                               TH2D& hfwd, Int_t vtx, UShort_t flags)
 450 {
 451   //
 452   //  Loops over list of VtxBins, fills hists of bins for current vertex
 453   //  and runs analysis on those bins
 454   //
 455   //  Parameters:
 456   //   list:  list of VtxBins
 457   //   hcent: dN/detadphi histogram for central coverage
 458   //   hfwd:  dN/detadphi histogram for forward coverage
 459   //   vBin:  current vertex bin
 460   //   flags: extra flags to handle calculations.
 461   //
 462   VertexBin* bin = 0;
 463   Int_t i = 0;
 464   Int_t nVtxBins = fVtxAxis->GetNbins();
 465
 466   TH2D& h = CombineHists(hcent, hfwd);
 467
 468   while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
 469     i++;
 470     if (!bin->FillHists(h, fCent, kFillBoth|flags|kReset)) continue;
 471     bin->CumulantsAccumulate3Cor(fCent);
 472   }
 473
 474   return;
 475 }
 476 //_____________________________________________________________________
 477 TH2D& AliForwardFlowTaskQC::CombineHists(TH2D& hcent, TH2D& hfwd)
 478 {
 479   //
 480   //  Combines a forward and central d^2N/detadphi histogram.
 481   //  At some point it might need a flag to choose which histogram gets
 482   //  priority when there is an overlap, at the moment the average is chosen
 483   //
 484   //  Parameters:
 485   //    hcent: Central barrel detector
 486   //    hfwd: Forward detector
 487   //
 488   //  Return: reference to combined hist
 489   //
 490
 491   // If hists are the same (MC input) don't do anything
 492   if (&hcent == &hfwd) return hcent;
 493
 494   fHistdNdedp3Cor.Reset();
 495   // FMD, SPD input
 496   if ((fFlowFlags & kFMD)) {
 497     for (Int_t e = 1; e <= fHistdNdedp3Cor.GetNbinsX(); e++)  {
 498       Double_t eta = fHistdNdedp3Cor.GetXaxis()->GetBinCenter(e);
 499       Bool_t fwdCov = (hfwd.GetBinContent(e, 0) != 0);
 500       Bool_t centCov = (hcent.GetBinContent(e, 0) != 0);
 501       if (!fwdCov && !centCov) continue;
 502       else fHistdNdedp3Cor.SetBinContent(e, 0, 1);
 503       for (Int_t p = 1; p <= fHistdNdedp3Cor.GetNbinsY(); p++) {
 504         Double_t phi = fHistdNdedp3Cor.GetYaxis()->GetBinCenter(p);
 505         Int_t n = 0;
 506         Double_t cont = 0.;
 507         if (fwdCov) {
 508           cont += hfwd.GetBinContent(e, p);
 509           n++;
 510         }
 511         if (centCov) {
 512           cont += hcent.GetBinContent(e, p);
 513           n++;
 514         }
 515         if (cont == 0 || n == 0) continue;
 516         cont /= n;
 517         fHistdNdedp3Cor.Fill(eta, phi, cont);
 518       }
 519     }
 520   // VZERO, SPD input, here we do not average but cut to avoid
 521   // (too much) overlap.
 522   } else if ((fFlowFlags & kVZERO)) {
 523     // VZERO loop
 524     for (Int_t eV = 1; eV <= hfwd.GetNbinsX(); eV++) {
 525       Double_t eta = hfwd.GetXaxis()->GetBinLowEdge(eV)+0.1;
 526       if (hfwd.GetBinContent(eV, 0) == 0) continue;
 527       else {
 528         Int_t he = fHistdNdedp3Cor.GetXaxis()->FindBin(eta);
 529         fHistdNdedp3Cor.SetBinContent(he, 0, 1);
 530       }
 531       for (Int_t p = 1; p <= hfwd.GetNbinsY(); p++) {
 532         Double_t phi = hfwd.GetYaxis()->GetBinCenter(p);
 533         Double_t cont = hfwd.GetBinContent(eV, p);
 534         fHistdNdedp3Cor.Fill(eta, phi, cont);
 535       }
 536     }
 537     // SPD loop
 538     Int_t eSs = hcent.GetXaxis()->FindBin(-1.99);
 539     Int_t eSe = hcent.GetXaxis()->FindBin(1.99);
 540     for (Int_t eS = eSs; eS <= eSe; eS++) {
 541       Double_t eta = hcent.GetXaxis()->GetBinCenter(eS);
 542       if (hcent.GetBinContent(eS, 0) == 0) continue;
 543       else {
 544         Int_t he = fHistdNdedp3Cor.GetXaxis()->FindBin(eta);
 545         fHistdNdedp3Cor.SetBinContent(he, 0, 1);
 546       }
 547       for (Int_t p = 1; p <= hcent.GetNbinsY(); p++) {
 548         Double_t phi = hcent.GetYaxis()->GetBinCenter(p);
 549         Double_t cont = hcent.GetBinContent(eS, p);
 550         fHistdNdedp3Cor.Fill(eta, phi, cont);
 551       }
 552     }
 553   }
 554   return fHistdNdedp3Cor;
 555 }
 556 //_____________________________________________________________________
 557 Bool_t AliForwardFlowTaskQC::FillTracks(VertexBin* bin, UShort_t mode) const
 558 {
 559   //
 560   //  Get TPC tracks to use for reference flow.
 561   //
 562   //  Return: TObjArray with tracks
 563   //
 564   TObjArray* trList = 0;
 565   AliESDEvent* esdEv = 0;
 566   if (AliForwardUtil::CheckForAOD() == 1) // AOD tracks
 567     trList = static_cast<TObjArray*>(fAOD->GetTracks());
 568   else
 569     esdEv = dynamic_cast<AliESDEvent*>(InputEvent());
 570
 571   Bool_t useEvent = bin->FillTracks(trList, esdEv, fTrackCuts, mode);
 572   return useEvent;
 573 }
 574 //_____________________________________________________________________
 575 void AliForwardFlowTaskQC::Terminate(Option_t */*option*/)
 576 {
 577   //
 578   //  Calls the finalize function, done at the end of the analysis
 579   //
 580   //  Parameters:
 581   //   option: Not used
 582   //
 583
 584   // Make sure pointers are set to the correct lists
 585   fSumList = dynamic_cast<TList*> (GetOutputData(1));
 586   if(!fSumList) {
 587     AliError("Could not retrieve TList fSumList");
 588     return;
 589   }
 590   if (!fOutputList)
 591     fOutputList = new TList();
 592   fOutputList->SetName("Results");
 593   fOutputList->SetOwner();
 594
 595   if ((fFlowFlags & kEtaGap) || (fFlowFlags & kTracks)) {
 596     TParameter<Double_t>* etaGap = new TParameter<Double_t>("EtaGap", fEtaGap);
 597     fOutputList->Add(etaGap);
 598   }
 599   // We only add axes in terminate, as TAxis object do not merge well,
 600   // and so we get a mess when running on the grid if we put them in the sum list...
 601   fVtxAxis->SetName("VtxAxis");
 602   fOutputList->Add(fVtxAxis);
 603   fCentAxis->SetName("CentAxis");
 604   fOutputList->Add(fCentAxis);
 605
 606   // Run finalize on VertexBins
 607   Finalize();
 608
 609   // Loop over output to get dN/deta hists - used for diagnostics
 610   TIter next(fOutputList);
 611   TObject* o = 0;
 612   TString name;
 613   TH2D* dNdeta = 0;
 614   TH1D* cent = 0;
 615   while ((o = next())) {
 616     name = o->GetName();
 617     if (name.Contains("dNdeta")) {
 618       dNdeta = dynamic_cast<TH2D*>(o);
 619       name.ReplaceAll("dNdeta", "cent");
 620       name.ReplaceAll("Ref", "");
 621       name.ReplaceAll("Diff", "");
 622       cent = dynamic_cast<TH1D*>(fOutputList->FindObject(name.Data()));
 623       if (!dNdeta || !cent) continue;
 624       for (Int_t cBin = 1; cBin <= dNdeta->GetNbinsY(); cBin++) {
 625         Double_t nEvents = cent->GetBinContent(cBin);
 626         if (nEvents == 0) continue;
 627         for (Int_t eBin = 1; eBin <= dNdeta->GetNbinsX(); eBin++) {
 628           dNdeta->SetBinContent(eBin, cBin, dNdeta->GetBinContent(eBin, cBin)/nEvents);
 629           dNdeta->SetBinError(eBin, cBin, dNdeta->GetBinError(eBin, cBin)/nEvents);
 630         }
 631       }
 632     }
 633   }
 634
 635   // Loop over output and make 1D projections for fast look at results
 636   MakeCentralityHists(fOutputList);
 637   TList* vtxList = (TList*)fOutputList->FindObject("vtxList");
 638   if (vtxList) MakeCentralityHists(vtxList);
 639   TList* nuaList = (TList*)fOutputList->FindObject("NUATerms");
 640   TIter nextNua(nuaList);
 641   o = 0;
 642   TH2D* h = 0;
 643   while ((o = nextNua())) {
 644     if (!(h = dynamic_cast<TH2D*>(o))) continue;
 645     Double_t evts = h->GetBinContent(0, 0);
 646     if (evts != 0) h->Scale(1./evts);
 647   }
 648   if (nuaList) MakeCentralityHists(nuaList);
 649
 650   PostData(2, fOutputList);
 651
 652   return;
 653 }
 654 //_____________________________________________________________________
 655 void AliForwardFlowTaskQC::Finalize()
 656 {
 657   //
 658   //  Finalize command, called by Terminate()
 659   //
 660
 661   // Reinitiate vertex bins if Terminate is called separately!
 662   if (fBinsForward.GetEntries() == 0) InitVertexBins();
 663
 664   // Iterate over all vertex bins objects and finalize cumulants
 665   // calculations
 666   EndVtxBinList(fBinsForward);
 667   EndVtxBinList(fBinsCentral);
 668
 669   return;
 670 }
 671 //_____________________________________________________________________
 672 void AliForwardFlowTaskQC::EndVtxBinList(const TList& list) const
 673 {
 674   //
 675   //  Loop over VertexBin list and call terminate on each
 676   //
 677   //  Parameters:
 678   //   list: VertexBin list
 679   //
 680   TIter next(&list);
 681   VertexBin* bin = 0;
 682   while ((bin = static_cast<VertexBin*>(next()))) {
 683     bin->CumulantsTerminate(fSumList, fOutputList);
 684   }
 685   return;
 686 }
 687 // _____________________________________________________________________
 688 void AliForwardFlowTaskQC::MakeCentralityHists(TList* list) const
 689 {
 690   //
 691   // Loop over a list containing a TH2D with flow results
 692   // and project to TH1's in specific centrality bins
 693   //
 694   // Parameters:
 695   //  list: Flow results list
 696   //
 697   TH2D* hist2D = 0;
 698   TList* centList = 0;
 699   TH1D* hist1D = 0;
 700   TObject* helper = 0;
 701   TIter nextHist(list);
 702   while ((helper = dynamic_cast<TObject*>(nextHist()))) {
 703     if (!(hist2D = dynamic_cast<TH2D*>(helper))) continue;
 704     for (Int_t cBin = 1; cBin <= hist2D->GetNbinsY(); cBin++) {
 705       Int_t cMin = Int_t(hist2D->GetYaxis()->GetBinLowEdge(cBin));
 706       Int_t cMax = Int_t(hist2D->GetYaxis()->GetBinUpEdge(cBin));
 707       TString name = Form("cent_%d-%d", cMin, cMax);
 708       centList = (TList*)list->FindObject(name.Data());
 709       if (!centList) {
 710         centList = new TList();
 711         centList->SetName(name.Data());
 712         list->Add(centList);
 713       }
 714       hist1D = hist2D->ProjectionX(Form("%s_%s", hist2D->GetName(), name.Data()),
 715                                          cBin, cBin, "E");
 716       hist1D->SetTitle(hist1D->GetName());
 717       centList->Add(hist1D);
 718     }
 719   }
 720 }
 721 // _____________________________________________________________________
 722 Bool_t AliForwardFlowTaskQC::CheckEvent(const AliAODForwardMult* aodfm)
 723 {
 724   //
 725   //  Function to check that an AOD event meets the cuts
 726   //
 727   //  Parameters:
 728   //   AliAODForwardMult: forward mult object with trigger and vertex info
 729   //
 730   //  Return: false if there is no trigger or if the centrality or vertex
 731   //  is out of range. Otherwise true.
 732   //
 733
 734   // First check for trigger
 735   if (!CheckTrigger(aodfm)) {
 736     fHistEventSel->Fill(kNoTrigger);
 737     return kFALSE;
 738   }
 739
 740   // Then check for centrality
 741   if (!GetCentrality(aodfm)) {
 742     return kFALSE;
 743   }
 744
 745   // And finally check for vertex
 746   if (!GetVertex(aodfm)) {
 747     return kFALSE;
 748   }
 749
 750   // Ev. accepted - filling diag. hists
 751   fHistCent->Fill(fCent);
 752   fHistVertexSel->Fill(fVtx);
 753   fHistEventSel->Fill(kOK);
 754
 755   return kTRUE;
 756 }
 757 // _____________________________________________________________________
 758 Bool_t AliForwardFlowTaskQC::CheckTrigger(const AliAODForwardMult* aodfm) const
 759 {
 760   //
 761   //  Function to look for a trigger string in the event.
 762   //  First check for info in forward mult object, if not there, use the AOD header
 763   //
 764   //  Parameters:
 765   //   AliAODForwardMult: forward mult object with trigger and vertex info
 766   //
 767   //  Return: true if offline trigger is present
 768   //
 769   if (aodfm) return aodfm->IsTriggerBits(AliAODForwardMult::kOffline);
 770   // this may need to be changed for 2011 data to handle kCentral and so on...
 771   else return (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
 772                  ->IsEventSelected() & AliVEvent::kMB);
 773 }
 774 // _____________________________________________________________________
 775 Bool_t AliForwardFlowTaskQC::GetCentrality(const AliAODForwardMult* aodfm)
 776 {
 777   //
 778   //  Function to look get centrality of the event.
 779   //  First check for info in forward mult object, if not there, use the AOD header
 780   //
 781   //  Parameters:
 782   //   AliAODForwardMult: forward mult object with trigger and vertex info
 783   //
 784   //  Return: true if centrality determination is present
 785   //
 786   if (aodfm) {
 787     if (aodfm->HasCentrality()) {
 788       fCent = (Double_t)aodfm->GetCentrality();
 789       if (fCentAxis->GetXmin() > fCent || fCent >= fCentAxis->GetXmax()) {
 790         fHistEventSel->Fill(kInvCent);
 791         return kFALSE;
 792       }
 793     }
 794     else {
 795       fCent = 97.5;
 796       fHistEventSel->Fill(kNoCent);
 797     }
 798     return kTRUE;
 799   } else {
 800     AliCentrality* aodCent = fAOD->GetCentrality();
 801     if (aodCent) {
 802       fCent = (Double_t)aodCent->GetCentralityPercentile("V0M");
 803       if (fCentAxis->GetXmin() > fCent || fCent >= fCentAxis->GetXmax()) {
 804         fHistEventSel->Fill(kInvCent);
 805         return kFALSE;
 806       }
 807     }
 808     else {
 809       fCent = 97.5;
 810       fHistEventSel->Fill(kNoCent);
 811     }
 812     return kTRUE;
 813   }
 814 }
 815 //_____________________________________________________________________
 816 Bool_t AliForwardFlowTaskQC::GetVertex(const AliAODForwardMult* aodfm)
 817 {
 818   //
 819   //  Function to look for vertex determination in the event.
 820   //  First check for info in forward mult object, if not there, use the AOD header
 821   //
 822   //  Parameters:
 823   //   AliAODForwardMult: forward mult object with trigger and vertex info
 824   //
 825   //  Return: true if vertex is determined
 826   //
 827   if (aodfm) {
 828     if (aodfm->HasIpZ()) {
 829       fVtx = aodfm->GetIpZ();
 830       if (fVtx < fVtxAxis->GetXmin() || fVtx >= fVtxAxis->GetXmax()) {
 831         fHistEventSel->Fill(kInvVtx);
 832         return kFALSE;
 833       }
 834     } else {
 835       fVtx = 9999;
 836       fHistEventSel->Fill(kNoVtx);
 837       return kFALSE;
 838     }
 839     return kTRUE;
 840   } else {
 841     AliAODVertex* aodVtx = fAOD->GetPrimaryVertex();
 842     if (aodVtx) {
 843       fVtx = aodVtx->GetZ();
 844       if (fVtx < fVtxAxis->GetXmin() || fVtx >= fVtxAxis->GetXmax()) {
 845         fHistEventSel->Fill(kInvVtx);
 846         return kFALSE;
 847       }
 848     } else {
 849       fVtx = 9999;
 850       fHistEventSel->Fill(kNoVtx);
 851       return kFALSE;
 852     }
 853     return kTRUE;
 854   }
 855 }
 856 // _____________________________________________________________________
 857 AliVVZERO* AliForwardFlowTaskQC::GetVZERO() const
 858 {
 859   //
 860   //  Get VZERO object from ESD or AOD
 861   //
 862   //  Return: VZERO data object
 863   //
 864   AliVVZERO* vzero = 0;
 865   // Get input type
 866   UShort_t input = AliForwardUtil::CheckForAOD();
 867   switch (input) {
 868     // If AOD input, simply get the track array from the event
 869     case 1: vzero = (AliVVZERO*)fAOD->GetVZEROData();
 870             break;
 871     case 2: {
 872     // If ESD input get event, apply track cuts
 873               AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 874               if (!esd) return 0;
 875               vzero = (AliVVZERO*)esd->GetVZEROData();
 876               break;
 877             }
 878     default: AliFatal("Neither ESD or AOD input. This should never happen");
 879             break;
 880   }
 881   return vzero;
 882 }
 883 // _____________________________________________________________________
 884 void AliForwardFlowTaskQC::FillVZEROHist(AliVVZERO* vzero)
 885 {
 886   //
 887   //  Loops over VZERO data object and fill up d^2N/detadphi histogram for flow analysis
 888   //
 889   //  Parameters:
 890   //   vzero: VZERO AOD data object
 891   //
 892   Int_t ring = 0;
 893   Int_t bin = 0;
 894   Double_t eta = 0;
 895   // Sort of right for 2010 data, do not use for 2011!
 896   Double_t eq[64] = { 1.43536, 1.45727, 1.44993, 1.30051, 1.17425, 1.2335, 1.22247, 1.14362,
 897                       1.14647, 1.25208, 1.17681, 1.21642, 1.16604, 1.05532, 1.03212, 1.1032,
 898                       1.22941, 1.36986, 1.14652, 1.20056, 0.927086, 1.10809, 1.03343, 1.29472,
 899                       1.21204, 1.29217, 1.2003, 2.10382, 1.28513, 1.40558, 1.25784, 1.21848,
 900                       0.475162, 0.50421, 0.503617, 0.512471, 0.515276, 0.39831, 0.415199, 0.444664,
 901                       0.521922, 0.785915, 0.703658, 0.832479, 0.77461, 0.73129, 0.778697, 0.710265,
 902                       0.89686, 0.967688, 0.974225, 0.873445, 0.811096, 0.828493, 0.889609, 0.586056,
 903                       1.15877, 0.954656, 0.914557, 0.979028, 1.04907, 0.748518, 0.928043, 0.98175 };
 904   for (Int_t i = 0; i < 64; i++) {
 905     if (i % 8 == 0) {
 906       ring++;
 907       bin = (ring < 5 ? ring+1 : 15-ring);
 908       eta = fHistdNdedpV0.GetXaxis()->GetBinCenter(bin);
 909       fHistdNdedpV0.SetBinContent(bin, 0, 1);
 910     }
 911     Float_t amp = vzero->GetMultiplicity(i);
 912     amp /= eq[i];
 913     Double_t phi = TMath::Pi()/8.+TMath::TwoPi()*i/8.;
 914     while (phi > TMath::TwoPi()) phi -= TMath::TwoPi();
 915     fHistdNdedpV0.Fill(eta, phi, amp);
 916   }
 917 }
 918 //_____________________________________________________________________
 919 AliForwardFlowTaskQC::VertexBin::VertexBin()
 920   : TNamed(),
 921     fMaxMoment(0),   // Max flow moment for this vertexbin
 922     fVzMin(0),       // Vertex z-coordinate min [cm]
 923     fVzMax(0),       // Vertex z-coordinate max [cm]
 924     fType(),         // Data type name e.g., FMD/SPD/FMDTR/SPDTR/MC
 925     fFlags(0),       // Flow flags
 926     fSigmaCut(-1),   // Sigma cut to remove outlier events
 927     fEtaGap(-1),     // Eta gap value
 928     fEtaLims(),      // Limits for binning in 3Cor method
 929     fCumuRef(),      // Histogram for reference flow
 930     fCumuDiff(),     // Histogram for differential flow
 931     fCumuHists(0,0), // CumuHists object for keeping track of results
 932     fCumuNUARef(),   // Histogram for ref NUA terms
 933     fCumuNUADiff(),  // Histogram for diff NUA terms
 934     fdNdedpRefAcc(), // Diagnostics histogram for acc. maps
 935     fdNdedpDiffAcc(),// Diagnostics histogram for acc. maps
 936     fOutliers(),     // Histogram for sigma distribution
 937     fDebug()         // Debug level
 938 {
 939   //
 940   //  Default constructor
 941   //
 942 }
 943 //_____________________________________________________________________
 944 AliForwardFlowTaskQC::VertexBin::VertexBin(Int_t vLow, Int_t vHigh,
 945                                            UShort_t moment, TString name,
 946                                            UShort_t flags, Double_t cut,
 947                                            Double_t etaGap)
 948   : TNamed("", ""),
 949     fMaxMoment(moment),  // Max flow moment for this vertexbin
 950     fVzMin(vLow),        // Vertex z-coordinate min [cm]
 951     fVzMax(vHigh),       // Vertex z-coordinate max [cm]
 952     fType(name),         // Data type name e.g., FMD/SPD/FMDTR/SPDTR/MC
 953     fFlags(flags),       // Flow flags
 954     fSigmaCut(cut),      // Sigma cut to remove outlier events
 955     fEtaGap(etaGap),     // Eta gap value
 956     fEtaLims(),          // Limits for binning in 3Cor method
 957     fCumuRef(),          // Histogram for reference flow
 958     fCumuDiff(),         // Histogram for differential flow
 959     fCumuHists(moment,0),// CumuHists object for keeping track of results
 960     fCumuNUARef(),       // Histogram for ref NUA terms
 961     fCumuNUADiff(),      // Histogram for diff NUA terms
 962     fdNdedpRefAcc(),     // Diagnostics histogram for acc. maps
 963     fdNdedpDiffAcc(),    // Diagnostics histogram for acc. maps
 964     fOutliers(),         // Histogram for sigma distribution
 965     fDebug(0)            // Debug level
 966 {
 967   //
 968   //  Constructor
 969   //
 970   //  Parameters
 971   //   vLow: min z-coordinate
 972   //   vHigh: max z-coordinate
 973   //   moment: max flow moment
 974   //   name: data type name (FMD/SPD/FMDTR/SPDTR/MC)
 975   //   flags: flow flags
 976   //   cut: sigma cut
 977   //   etaGap: eta-gap value
 978   //
 979   fType.ToUpper();
 980
 981   SetName(Form("%svertexBin%d_%d_%d%s", fType.Data(), moment, vLow, vHigh, GetQCType(fFlags)));
 982   SetTitle(Form("%svertexBin%d_%d_%d%s", fType.Data(), moment, vLow, vHigh, GetQCType(fFlags)));
 983
 984   fDebug = AliAnalysisManager::GetAnalysisManager()->GetDebugLevel();
 985   if (fDebug > 0) Printf("AliForwardFlowTaskQC::VertexBin()\tDebugMode: %d", fDebug);
 986
 987   // Set limits for 3 correlator method
 988   if ((fFlags & kFMD)) {
 989     fEtaLims[0] = -6.;
 990     fEtaLims[1] = -1.5;
 991     fEtaLims[2] = -0.5;
 992     fEtaLims[3] = 2.;
 993     fEtaLims[4] = 3.;
 994     fEtaLims[5] = 6.;
 995   } else if ((fFlags & kVZERO)) {
 996     fEtaLims[0] = -6;
 997     fEtaLims[1] = -2.7;
 998     fEtaLims[2] = -2.0;
 999     fEtaLims[3] = 2.0;
1000     fEtaLims[4] = 3.9;
1001     fEtaLims[5] = 6;
1002   }
1003 }
1004 //_____________________________________________________________________
1005 AliForwardFlowTaskQC::VertexBin&
1006 AliForwardFlowTaskQC::VertexBin::operator=(const AliForwardFlowTaskQC::VertexBin& o)
1007 {
1008   //
1009   //  Assignment operator
1010   //
1011   //  Parameters
1012   //   o: AliForwardFlowTaskQC::VertexBin
1013   //
1014   if (&o == this) return *this;
1015   fMaxMoment     = o.fMaxMoment;
1016   fVzMin         = o.fVzMin;
1017   fVzMax         = o.fVzMax;
1018   fType          = o.fType;
1019   fFlags         = o.fFlags;
1020   fSigmaCut      = o.fSigmaCut;
1021   fEtaGap        = o.fEtaGap;
1022   fCumuRef       = o.fCumuRef;
1023   fCumuDiff      = o.fCumuDiff;
1024   fCumuHists     = o.fCumuHists;
1025   fCumuNUARef    = o.fCumuNUARef;
1026   fCumuNUADiff   = o.fCumuNUADiff;
1027   fdNdedpRefAcc  = o.fdNdedpRefAcc;
1028   fdNdedpDiffAcc = o.fdNdedpDiffAcc;
1029   fOutliers      = o.fOutliers;
1030   fDebug         = o.fDebug;
1031   for (UInt_t i = 0; i < sizeof(fEtaLims)/sizeof(Double_t); i++) fEtaLims[i] = o.fEtaLims[i];
1032
1033   return *this;
1034 }
1035 //_____________________________________________________________________
1036 void AliForwardFlowTaskQC::VertexBin::AddOutput(TList* outputlist, TAxis* centAxis)
1037 {
1038   //
1039   //  Add histograms to outputlist
1040   //
1041   //  Parameters
1042   //   outputlist: list of histograms
1043   //   centAxis: centrality axis
1044   //
1045
1046   // First we try to find an outputlist for this vertexbin
1047   TList* list = (TList*)outputlist->FindObject(Form("%svertex_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1048   // If it doesn't exist we make one
1049   if (!list) {
1050     list = new TList();
1051     list->SetName(Form("%svertex_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1052     outputlist->Add(list);
1053   }
1054
1055   // Get bin numbers and binning defined
1056   Int_t nHBins = GetBinNumberSin();
1057   Int_t nEtaBins = 48;
1058   if ((fFlags & k3Cor)) {
1059     if ((fFlags & kFMD)) nEtaBins = 24;
1060     else if ((fFlags & kVZERO)) nEtaBins = 19;
1061   }
1062   else if ((fFlags & kVZERO) && (fFlags & kEtaGap)) nEtaBins = 19;
1063   else if ((fFlags & kVZERO)) nEtaBins = 11;
1064
1065   Double_t vzeroBins[12] = { -6, -3.7, -3.2, -2.7, -2.2, -1.7,
1066                              2.8, 3.4,  3.9,  4.5,  5.1, 6 };
1067   Double_t vzeroBins2[20] = { -6, -3.7, -3.2, -2.7, -2.2, // VZERO
1068                              -2.0, -1.5, -1.0, -0.5 , 0., 0.5, 1.0, 1.5, 2.0, // SPD
1069                             2.8, 3.4,  3.9,  4.5,  5.1, 6 }; // VZERO
1070
1071   Int_t nRefBins = nEtaBins; // needs to be something as default
1072   if ((fFlags & kStdQC)) {
1073     if ((fFlags & kSymEta) && !((fFlags & kTracks) && (fFlags & kSPD))) nRefBins = 1;
1074     else nRefBins = 2;
1075   } else if ((fFlags & kEtaGap )) {
1076     nRefBins = 2;
1077   } else if ((fFlags & k3Cor)) {
1078     nRefBins = 24;
1079   }
1080
1081   // We initiate the reference histogram
1082   fCumuRef = new TH2D(Form("%s_%d_%d%s_ref", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1083                       Form("%s_%d_%d%s_ref", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1084                         nRefBins, -6., 6.,
1085                         nHBins, 0.5, nHBins+0.5);
1086   if ((fFlags & kVZERO) && (fFlags & k3Cor)) fCumuRef->GetXaxis()->Set(nEtaBins, vzeroBins2);
1087   SetupNUALabels(fCumuRef->GetYaxis());
1088   list->Add(fCumuRef);
1089   // We initiate the differential histogram
1090   fCumuDiff = new TH2D(Form("%s_%d_%d%s_diff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1091                        Form("%s_%d_%d%s_diff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1092                        nEtaBins, -6., 6., nHBins, 0.5, nHBins+0.5);
1093   if ((fFlags & kVZERO)) {
1094     if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1095       fCumuDiff->GetXaxis()->Set(nEtaBins, vzeroBins2);
1096     else fCumuDiff->GetXaxis()->Set(nEtaBins, vzeroBins);
1097   }
1098   SetupNUALabels(fCumuDiff->GetYaxis());
1099   list->Add(fCumuDiff);
1100
1101   // Cumulants sum hists
1102   Int_t cBins = centAxis->GetNbins();
1103   fCumuHists.ConnectList(Form("%sCumu_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)), list);
1104   TH3D* cumuHist = 0;
1105   Int_t nC2Bins = ((fFlags & kEtaGap) || (fFlags & k3Cor) ? kW2 : k3pWeight);
1106   Int_t nC4Bins = ((fFlags & kEtaGap) ? kW2 : ((fFlags & k3Cor) ? kW4 : kSinphi1phi2phi3p));
1107   for (Int_t n = 2; n <= fMaxMoment; n++) {
1108     // Initiate the ref cumulant sum histogram
1109     cumuHist = new TH3D(Form("%sv%d_vertex_%d_%d%s_cumuRef", fType.Data(), n, fVzMin, fVzMax, GetQCType(fFlags)),
1110                         Form("%sv%d_vertex_%d_%d%s_cumuRef", fType.Data(), n, fVzMin, fVzMax, GetQCType(fFlags)),
1111                         nRefBins, -6., 6.,
1112                         cBins, 0., 100., nC2Bins, 0.5, nC2Bins+0.5);
1113     if ((fFlags & kVZERO) && (fFlags & k3Cor)) cumuHist->GetXaxis()->Set(nEtaBins, vzeroBins2);
1114     cumuHist->GetYaxis()->Set(cBins, centAxis->GetXbins()->GetArray());
1115     fCumuHists.Add(cumuHist);
1116     // Initiate the diff cumulant sum histogram
1117     cumuHist = new TH3D(Form("%sv%d_vertex_%d_%d%s_cumuDiff", fType.Data(), n, fVzMin, fVzMax, GetQCType(fFlags)),
1118                         Form("%sv%d_vertex_%d_%d%s_cumuDiff", fType.Data(), n, fVzMin, fVzMax, GetQCType(fFlags)),
1119                         nEtaBins, -6., 6., cBins, 0., 100., nC4Bins, 0.5, nC4Bins+0.5);
1120     if ((fFlags & kVZERO)) {
1121       if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1122         cumuHist->GetXaxis()->Set(nEtaBins, vzeroBins2);
1123       else cumuHist->GetXaxis()->Set(nEtaBins, vzeroBins);
1124     }
1125     cumuHist->GetYaxis()->Set(cBins, centAxis->GetXbins()->GetArray());
1126     fCumuHists.Add(cumuHist);
1127   }
1128
1129   // Common NUA histograms
1130   fCumuNUARef = new TH3D(Form("%s_vertex_%d_%d%s_cumuNUARef", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1131                          Form("%s_vertex_%d_%d%s_cumuNUARef", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1132                          nRefBins, -6., 6.,
1133                           cBins, 0., 100., nHBins, 0.5, nHBins+0.5);
1134   if ((fFlags & kVZERO) && (fFlags & k3Cor)) fCumuNUARef->GetXaxis()->Set(nEtaBins, vzeroBins2);
1135   fCumuNUARef->GetYaxis()->Set(cBins, centAxis->GetXbins()->GetArray());
1136   SetupNUALabels(fCumuNUARef->GetZaxis());
1137   fCumuNUARef->Sumw2();
1138   list->Add(fCumuNUARef);
1139
1140   fCumuNUADiff = new TH3D(Form("%s_vertex_%d_%d%s_cumuNUADiff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1141                           Form("%s_vertex_%d_%d%s_cumuNUADiff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1142                           nEtaBins, -6., 6., cBins, 0., 100., nHBins, 0.5, nHBins+0.5);
1143   if ((fFlags & kVZERO)) {
1144     if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1145       fCumuNUADiff->GetXaxis()->Set(nEtaBins, vzeroBins2);
1146     else fCumuNUADiff->GetXaxis()->Set(nEtaBins, vzeroBins);
1147   }
1148   fCumuNUADiff->GetYaxis()->Set(cBins, centAxis->GetXbins()->GetArray());
1149   SetupNUALabels(fCumuNUADiff->GetZaxis());
1150   fCumuNUADiff->Sumw2();
1151   list->Add(fCumuNUADiff);
1152
1153   // We create diagnostic histograms.
1154   TList* dList = (TList*)outputlist->FindObject("Diagnostics");
1155   if (!dList) AliFatal("No diagnostics list found");
1156
1157   // Acceptance hist
1158   Int_t nPhiBins = ((fFlags & kFMD) ? 20 : 8);
1159   fdNdedpRefAcc = new TH2F(Form("h%sdNdedpRefAcc_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1160     Form("%s reference flow acceptance map for %d cm < v_{z} < %d cm", fType.Data(), fVzMin, fVzMax),
1161     nEtaBins, -6., 6., nPhiBins, 0, TMath::TwoPi());
1162   if ((fFlags & kVZERO)) {
1163     if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1164       fdNdedpRefAcc->GetXaxis()->Set(nEtaBins, vzeroBins2);
1165     else fdNdedpRefAcc->GetXaxis()->Set(nEtaBins, vzeroBins);
1166   }
1167   dList->Add(fdNdedpRefAcc);
1168
1169   fdNdedpDiffAcc = new TH2F(Form("h%sdNdedpDiffAcc_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1170     Form("%s differential flow acceptance map for %d cm < v_{z} < %d cm", fType.Data(), fVzMin, fVzMax),
1171     nEtaBins, -6., 6., nPhiBins, 0, TMath::TwoPi());
1172   if ((fFlags & kVZERO)) {
1173     if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1174       fdNdedpDiffAcc->GetXaxis()->Set(nEtaBins, vzeroBins2);
1175     else fdNdedpDiffAcc->GetXaxis()->Set(nEtaBins, vzeroBins);
1176   }
1177   dList->Add(fdNdedpDiffAcc);
1178
1179   fOutliers = new TH2F(Form("hOutliers_%s_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1180                        Form("Maximum #sigma from mean N_{ch} pr. bin - %s, %d < v_{z} < %d",
1181                        fType.Data(), fVzMin, fVzMax),
1182                        20, 0., 100., 500, 0., ((fFlags & kMC) ? 15. : 5.));
1183   dList->Add(fOutliers);
1184
1185   return;
1186 }
1187 //_____________________________________________________________________
1188 Bool_t AliForwardFlowTaskQC::VertexBin::FillHists(TH2D& dNdetadphi, Double_t cent, UShort_t mode)
1189 {
1190   //
1191   //  Fill reference and differential eta-histograms
1192   //
1193   //  Parameters:
1194   //   dNdetadphi: 2D histogram with input data
1195   //   cent: centrality
1196   //   mode: filling mode: kFillRef/kFillDiff/kFillBoth
1197   //
1198   if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1199   Bool_t useEvent = kTRUE;
1200
1201   // Fist we reset histograms
1202   if ((mode & kReset)) {
1203     if ((mode & kFillRef))  fCumuRef->Reset();
1204     if ((mode & kFillDiff)) fCumuDiff->Reset();
1205   }
1206   // Then we loop over the input and calculate sum cos(k*n*phi)
1207   // and fill it in the reference and differential histograms
1208   Int_t nBadBins = 0;
1209   Double_t limit = 9999.;
1210   for (Int_t etaBin = 1; etaBin <= dNdetadphi.GetNbinsX(); etaBin++) {
1211     Double_t eta = dNdetadphi.GetXaxis()->GetBinCenter(etaBin);
1212     // Numbers to cut away bad events
1213     Double_t runAvg = 0;
1214     Double_t max = 0;
1215     Int_t nInAvg = 0;
1216     Double_t avgSqr = 0;
1217     for (Int_t phiBin = 0; phiBin <= dNdetadphi.GetNbinsY(); phiBin++) {
1218       // Check for acceptance
1219       if (phiBin == 0) {
1220         if (dNdetadphi.GetBinContent(etaBin, 0) == 0) break;
1221         // Central limit for eta gap break for reference flow
1222         if ((fFlags & kEtaGap) && (mode & kFillRef) &&
1223              TMath::Abs(eta) < fEtaGap) break;
1224         // Backward and forward eta gap break for reference flow
1225         if ((fFlags & kEtaGap) && (mode & kFillRef) && TMath::Abs(eta) > TMath::Abs(limit)) break;
1226         if ((fFlags & kStdQC) && (fFlags & kMC) && !(fFlags & kTracks)) {
1227           if (!(fFlags & kSPD) && TMath::Abs(eta) < 1.75) break;
1228           if ((fFlags & kSPD) && TMath::Abs(eta) > 2.00) break;
1229         }
1230         if (limit > 1e3) limit = dNdetadphi.GetXaxis()->GetBinLowEdge(etaBin);
1231         continue;
1232       } // End of phiBin == 0
1233       Double_t phi = dNdetadphi.GetYaxis()->GetBinCenter(phiBin);
1234       Double_t weight = dNdetadphi.GetBinContent(etaBin, phiBin);
1235
1236       // We calculate the average Nch per. bin
1237       avgSqr += weight*weight;
1238       runAvg += weight;
1239       nInAvg++;
1240       if (weight == 0) continue;
1241       if (weight > max) max = weight;
1242       // Fill into Cos() and Sin() hists
1243       if ((mode & kFillRef) && !((fFlags & kTracks) && (fFlags & kMC) && TMath::Abs(eta) > 0.75)) {
1244         fCumuRef->Fill(eta, 0., weight);// mult goes in underflowbin - no visual, but not needed?
1245         fdNdedpRefAcc->Fill(eta, phi, weight);
1246       }
1247       if ((mode & kFillDiff)) {
1248         fCumuDiff->Fill(eta, 0., weight);
1249         fdNdedpDiffAcc->Fill(eta, phi, weight);
1250       }
1251       for (Int_t n = 1; n <= 2*fMaxMoment; n++) {
1252         Double_t cosBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberCos(n));
1253         Double_t sinBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberSin(n));
1254         Double_t cosnPhi = weight*TMath::Cos(n*phi);
1255         Double_t sinnPhi = weight*TMath::Sin(n*phi);
1256         // fill ref
1257         if ((mode & kFillRef) && !((fFlags & kTracks) && (fFlags & kMC) && TMath::Abs(eta) > 0.75)) {
1258           fCumuRef->Fill(eta, cosBin, cosnPhi);
1259           fCumuRef->Fill(eta, sinBin, sinnPhi);
1260         }
1261         // fill diff
1262         if ((mode & kFillDiff)) {
1263           fCumuDiff->Fill(eta, cosBin, cosnPhi);
1264           fCumuDiff->Fill(eta, sinBin, sinnPhi);
1265         }
1266       } // End of NUA loop
1267     } // End of phi loop
1268     // Outlier cut calculations
1269     if (nInAvg > 0) {
1270       runAvg /= nInAvg;
1271       avgSqr /= nInAvg;
1272       Double_t stdev = (nInAvg > 1 ? TMath::Sqrt(nInAvg/(nInAvg-1))*TMath::Sqrt(avgSqr - runAvg*runAvg) : 0);
1273       Double_t nSigma = (stdev == 0 ? 0 : (max-runAvg)/stdev);
1274       if (fSigmaCut > 0. && nSigma >= fSigmaCut && cent < 60) nBadBins++;
1275       else nBadBins = 0;
1276       fOutliers->Fill(cent, nSigma);
1277       // We still finish the loop, for fOutliers to make sense,
1278       // but we do no keep the event for analysis
1279       if (nBadBins > 3) useEvent = kFALSE;
1280     }
1281   } // End of eta bin
1282
1283   return useEvent;
1284 }
1285 //_____________________________________________________________________
1286 Bool_t AliForwardFlowTaskQC::VertexBin::FillTracks(TObjArray* trList, AliESDEvent* esd,
1287                                                    AliAnalysisFilter* trFilter, UShort_t mode)
1288 {
1289   //
1290   //  Fill reference and differential eta-histograms
1291   //
1292   //  Parameters:
1293   //   trList: list of tracks
1294   //   mode: filling mode: kFillRef/kFillDiff/kFillBoth
1295   //
1296   if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1297   if (!trList && !esd) {
1298     AliError("FillTracks: No AOD track list or ESD event - something might be wrong!");
1299     return kFALSE;
1300   }
1301
1302   // Fist we reset histograms
1303   if ((mode & kReset)) {
1304     if ((mode & kFillRef))  fCumuRef->Reset();
1305     if ((mode & kFillDiff)) fCumuDiff->Reset();
1306   }
1307
1308   // Then we loop over the input and calculate sum cos(k*n*phi)
1309   // and fill it in the reference and differential histograms
1310   Int_t nTr = 0;
1311   if (trList) nTr = trList->GetEntries();
1312   if (esd) nTr = esd->GetNumberOfTracks();
1313   if (nTr == 0) return kFALSE;
1314   AliVTrack* tr = 0;
1315   // Cuts for AOD tracks (have already been applied to ESD tracks) - except dEdx
1316 //  const tpcdEdx = 10;
1317   for (Int_t i = 0; i < nTr; i++) { // track loop
1318     tr = (trList ? (AliVTrack*)trList->At(i) : (AliVTrack*)esd->GetTrack(i));
1319     if (!tr) continue;
1320     if (esd) {
1321       AliESDtrack* esdTr = (AliESDtrack*)tr;
1322       if (!trFilter->IsSelected(esdTr)) continue;
1323     }
1324     else if (trList) { // If AOD input
1325       Double_t pTMin = 0, pTMax = 0, etaMin = 0, etaMax = 0, minNCl = 0;
1326       UInt_t bit = 0;
1327       if ((fFlags & kTPC) == kTPC)    pTMin = 0.2, pTMax = 5., etaMin = -0.8, etaMax = 0.8, minNCl = 70, bit = 128;
1328       if ((fFlags & kHybrid) == kHybrid) pTMin = 0.2, pTMax = 5., etaMin = -0.8, etaMax = 0.8, minNCl = 70, bit = 272;
1329
1330       AliAODTrack* aodTr = (AliAODTrack*)tr;
1331       if (aodTr->GetID() > -1) continue;
1332       if (!aodTr->TestFilterBit(bit) || aodTr->Pt() > pTMax || aodTr->Pt() < pTMin ||
1333         aodTr->Eta() > etaMax || aodTr->Eta() < etaMin || aodTr->GetTPCNcls() < minNCl) continue;
1334     }
1335
1336 //    if (tr->GetTPCsignal() < tpcdEdx) continue;
1337     // Track accepted
1338     Double_t eta = tr->Eta();
1339     if (((fFlags & kSPD) || (fFlags & kEtaGap)) && TMath::Abs(eta) < fEtaGap) continue;
1340     Double_t phi = tr->Phi();
1341 //    Double_t pT = tr->Pt();
1342 //    AliAODMCHeader* head = static_cast<AliAODMCHeader*>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
1343 //    Double_t rp = head->GetReactionPlaneAngle();
1344 //    Double_t b = head->GetImpactParameter();
1345     Double_t weight = 1.;//AliForwardFlowWeights::Instance().CalcWeight(eta, pT, phi, 0, rp, b);
1346
1347     if ((mode & kFillRef)) {
1348       fCumuRef->Fill(eta, 0., weight);// mult goes in underflowbin - no visual, but not needed?
1349       fdNdedpRefAcc->Fill(eta, phi, weight);
1350     }
1351     if ((mode & kFillDiff)) {
1352       fCumuDiff->Fill(eta, 0., weight);
1353       fdNdedpDiffAcc->Fill(eta, phi, weight);
1354     }
1355     for (Int_t n = 1; n <= 2*fMaxMoment; n++) {
1356       Double_t cosBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberCos(n));
1357       Double_t sinBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberSin(n));
1358       Double_t cosnPhi = weight*TMath::Cos(n*phi);
1359       Double_t sinnPhi = weight*TMath::Sin(n*phi);
1360       // fill ref
1361       if ((mode & kFillRef)) {
1362         fCumuRef->Fill(eta, cosBin, cosnPhi);
1363         fCumuRef->Fill(eta, sinBin, sinnPhi);
1364       }
1365       // fill diff
1366       if ((mode & kFillDiff)) {
1367         fCumuDiff->Fill(eta, cosBin, cosnPhi);
1368         fCumuDiff->Fill(eta, sinBin, sinnPhi);
1369       }
1370     } // End of NUA loop
1371   } // End of track loop
1372   return kTRUE;
1373 }
1374 //_____________________________________________________________________
1375 void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate(Double_t cent)
1376 {
1377   //
1378   //  Calculate the Q cumulant up to order fMaxMoment
1379   //
1380   //  Parameters:
1381   //   cent: centrality of event
1382   //
1383   if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1384
1385   // Fill out NUA hists
1386   for (Int_t etaBin = 1; etaBin <= fCumuRef->GetNbinsX(); etaBin++) {
1387     Double_t eta = fCumuRef->GetXaxis()->GetBinCenter(etaBin);
1388     if (fCumuRef->GetBinContent(etaBin, 0) == 0) continue;
1389     if ((fFlags & kTracks) && (fFlags && kSPD) && !(fFlags & kEtaGap)) eta = -eta;
1390     for (Int_t qBin = 0; qBin <= fCumuRef->GetNbinsY(); qBin++) {
1391       fCumuNUARef->Fill(eta, cent, Double_t(qBin), fCumuRef->GetBinContent(etaBin, qBin));
1392     }
1393   }
1394   for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1395     Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1396     if (fCumuDiff->GetBinContent(etaBin, 0) == 0) continue;
1397     for (Int_t qBin = 0; qBin <= fCumuDiff->GetNbinsY(); qBin++) {
1398       fCumuNUADiff->Fill(eta, cent, Double_t(qBin), fCumuDiff->GetBinContent(etaBin, qBin));
1399     }
1400   }
1401
1402   // We create the objects needed for the analysis
1403   TH3D* cumuRef = 0;
1404   TH3D* cumuDiff = 0;
1405   // For each n we loop over the hists
1406   for (Int_t n = 2; n <= fMaxMoment; n++) {
1407     cumuRef  = (TH3D*)fCumuHists.Get('r',n);
1408     cumuDiff = (TH3D*)fCumuHists.Get('d',n);
1409     Int_t prevRefEtaBin = 0;
1410
1411     // Per mom. quantities
1412     Double_t dQnReA = 0, dQnImA = 0, multA = 0;
1413     Double_t dQnReB = 0, dQnImB = 0, multB = 0;
1414     Double_t dQ2nReA = 0, dQ2nImA = 0;
1415     Double_t two = 0, w2 = 0, four = 0, w4 = 0;
1416     Double_t cosPhi1Phi2 = 0, cosPhi1Phi2Phi3m = 0;
1417     Double_t sinPhi1Phi2 = 0, sinPhi1Phi2Phi3m = 0;
1418     for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1419       Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1420       Double_t refEta = eta;
1421       if ((fFlags & kTracks) && (fFlags && kSPD) && !(fFlags & kEtaGap)) refEta = -eta;
1422       Int_t refEtaBinA = fCumuRef->GetXaxis()->FindBin(refEta);
1423       if ((fFlags & kEtaGap)) refEta = -eta;
1424       Int_t refEtaBinB = fCumuRef->GetXaxis()->FindBin(refEta);
1425       if (refEtaBinA != prevRefEtaBin) {
1426         prevRefEtaBin = refEtaBinA;
1427         // Reference flow
1428         multA  = fCumuRef->GetBinContent(refEtaBinA, 0);
1429         dQnReA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberCos(n));
1430         dQnImA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberSin(n));
1431         dQ2nReA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberCos(n*2));
1432         dQ2nImA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberSin(n*2));
1433
1434         multB  = fCumuRef->GetBinContent(refEtaBinB, 0);
1435         dQnReB = fCumuRef->GetBinContent(refEtaBinB, GetBinNumberCos(n));
1436         dQnImB = fCumuRef->GetBinContent(refEtaBinB, GetBinNumberSin(n));
1437
1438         if (multA <= 3 || multB <= 3) return;
1439         // The reference flow is calculated
1440         // 2-particle
1441         if ((fFlags & kStdQC)) {
1442           w2 = multA * (multA - 1.);
1443           two = dQnReA*dQnReA + dQnImA*dQnImA - multA;
1444         } else {
1445           w2 = multA * multB;
1446           two = dQnReA*dQnReB + dQnImA*dQnImB;
1447         }
1448         cumuRef->Fill(eta, cent, kW2Two, two);
1449         cumuRef->Fill(eta, cent, kW2, w2);
1450
1451         // The reference flow is calculated
1452         // 4-particle
1453         if ((fFlags & kStdQC)) {
1454           w4 = multA * (multA - 1.) * (multA - 2.) * (multA - 3.);
1455
1456           four = 2.*multA*(multA-3.) + TMath::Power((TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.)),2.)
1457                    -4.*(multA-2.)*(TMath::Power(dQnReA,2.) + TMath::Power(dQnImA,2.))
1458                    -2.*(TMath::Power(dQnReA,2.)*dQ2nReA+2.*dQnReA*dQnImA*dQ2nImA-TMath::Power(dQnImA,2.)*dQ2nReA)
1459                    +(TMath::Power(dQ2nReA,2.)+TMath::Power(dQ2nImA,2.));
1460
1461           cumuRef->Fill(eta, cent, kW4Four, four);
1462           cumuRef->Fill(eta, cent, kW4, w4);
1463
1464           // NUA
1465           cosPhi1Phi2 = dQnReA*dQnReA - dQnImA*dQnImA - dQ2nReA;
1466           sinPhi1Phi2 = 2.*dQnReA*dQnImA - dQ2nImA;
1467
1468           cosPhi1Phi2Phi3m = dQnReA*(TMath::Power(dQnReA,2)+TMath::Power(dQnImA,2))
1469                               -dQnReA*dQ2nReA-dQnImA*dQ2nImA-2.*(multA-1)*dQnReA;
1470
1471           sinPhi1Phi2Phi3m = -dQnImA*(TMath::Power(dQnReA,2)+TMath::Power(dQnImA,2))
1472                               +dQnReA*dQ2nImA-dQnImA*dQ2nReA+2.*(multA-1)*dQnImA;
1473
1474           cumuRef->Fill(eta, cent, kCosphi1phi2, cosPhi1Phi2);
1475           cumuRef->Fill(eta, cent, kSinphi1phi2, sinPhi1Phi2);
1476           cumuRef->Fill(eta, cent, kCosphi1phi2phi3m, cosPhi1Phi2Phi3m);
1477           cumuRef->Fill(eta, cent, kSinphi1phi2phi3m, sinPhi1Phi2Phi3m);
1478           cumuRef->Fill(eta, cent, k3pWeight, multA*(multA-1.)*(multA-2.));
1479         } // End of QC{4}
1480       } // End of reference flow
1481       // For each etaBin bin the necessary values for differential flow is calculated
1482       Double_t mp = fCumuDiff->GetBinContent(etaBin, 0);
1483       Double_t pnRe = fCumuDiff->GetBinContent(etaBin, GetBinNumberCos(n));
1484       Double_t pnIm = fCumuDiff->GetBinContent(etaBin, GetBinNumberSin(n));
1485       Double_t p2nRe = fCumuDiff->GetBinContent(etaBin, GetBinNumberCos(n*2));
1486       Double_t p2nIm = fCumuDiff->GetBinContent(etaBin, GetBinNumberSin(n*2));
1487       if (mp == 0) continue;
1488       Double_t mq = 0;
1489       Double_t qnRe = 0;
1490       Double_t qnIm = 0;
1491       Double_t q2nRe = 0;
1492       Double_t q2nIm = 0;
1493
1494       // Differential flow calculations for each eta bin is done:
1495       // 2-particle differential flow
1496       if ((fFlags & kStdQC) && (!(fFlags & kTracks) || ((fFlags & kTracks) && (fFlags & kMC) && !(fFlags & kSPD) && TMath::Abs(eta) < 0.75))) {
1497         mq = mp;
1498         qnRe = pnRe;
1499         qnIm = pnIm;
1500         q2nRe = p2nRe;
1501         q2nIm = p2nIm;
1502       }
1503
1504       Double_t w2p = mp * multB - mq;
1505       Double_t twoPrime = pnRe*dQnReB + pnIm*dQnImB - mq;
1506
1507       cumuDiff->Fill(eta, cent, kW2Two, twoPrime);
1508       cumuDiff->Fill(eta, cent, kW2, w2p);
1509
1510       if ((fFlags & kEtaGap)) continue;
1511       // Differential flow calculations for each eta bin bin is done:
1512       // 4-particle differential flow
1513       Double_t w4p = (mp * multA - 3.*mq)*(multA - 1.)*(multA - 2.);
1514
1515       Double_t fourPrime = (TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.))*(pnRe*dQnReA+pnIm*dQnImA)
1516                           - q2nRe*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))
1517                           - 2.*q2nIm*dQnReA*dQnImA
1518                           - pnRe*(dQnReA*dQ2nReA+dQnImA*dQ2nImA)
1519                           + pnIm*(dQnImA*dQ2nReA-dQnReA*dQ2nImA)
1520                           - 2.*multA*(pnRe*dQnReA+pnIm*dQnImA)
1521                           - 2.*(TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.))*mq
1522                           + 6.*(qnRe*dQnReA+qnIm*dQnImA)
1523                           + 1.*(q2nRe*dQ2nReA+q2nIm*dQ2nImA)
1524                           + 2.*(pnRe*dQnReA+pnIm*dQnImA)
1525                           + 2.*mq*multA
1526                           - 6.*mq;
1527
1528       cumuDiff->Fill(eta, cent, kW4Four, fourPrime);
1529       cumuDiff->Fill(eta, cent, kW4, w4p);
1530
1531       // NUA
1532       Double_t cosPsi1Phi2 = pnRe*dQnReA - pnIm*dQnImA - q2nRe;
1533       Double_t sinPsi1Phi2 = pnRe*dQnImA + pnIm*dQnReA - q2nIm;
1534
1535       Double_t cosPsi1Phi2Phi3p = pnRe*(TMath::Power(dQnImA,2.)+TMath::Power(dQnReA,2.)-multA)
1536                             - 1.*(q2nRe*dQnReA+q2nIm*dQnImA)
1537                             - mq*dQnReA+2.*qnRe;
1538
1539       Double_t sinPsi1Phi2Phi3p = pnIm*(TMath::Power(dQnImA,2.)+TMath::Power(dQnReA,2.)-multA)
1540                             - 1.*(q2nIm*dQnReA-q2nRe*dQnImA)
1541                             - mq*dQnImA+2.*qnIm;
1542
1543       Double_t cosPsi1Phi2Phi3m = pnRe*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))+2.*pnIm*dQnReA*dQnImA
1544                             - 1.*(pnRe*dQ2nReA+pnIm*dQ2nImA)
1545                             - 2.*mq*dQnReA+2.*qnRe;
1546
1547       Double_t sinPsi1Phi2Phi3m = pnIm*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))-2.*pnRe*dQnReA*dQnImA
1548                             - 1.*(pnIm*dQ2nReA-pnRe*dQ2nImA)
1549                             + 2.*mq*dQnImA-2.*qnIm;
1550
1551       cumuDiff->Fill(eta, cent, kCosphi1phi2, cosPsi1Phi2);
1552       cumuDiff->Fill(eta, cent, kSinphi1phi2, sinPsi1Phi2);
1553       cumuDiff->Fill(eta, cent, kCosphi1phi2phi3m, cosPsi1Phi2Phi3m);
1554       cumuDiff->Fill(eta, cent, kSinphi1phi2phi3m, sinPsi1Phi2Phi3m);
1555       cumuDiff->Fill(eta, cent, k3pWeight, (mp*multA-2.*mq)*(multA-1.));
1556       cumuDiff->Fill(eta, cent, kCosphi1phi2phi3p, cosPsi1Phi2Phi3p);
1557       cumuDiff->Fill(eta, cent, kSinphi1phi2phi3p, sinPsi1Phi2Phi3p);
1558     } // End of eta loop
1559     // Event count
1560     cumuRef->Fill(-7., cent, -0.5, 1.);
1561   } // End of moment loop
1562   return;
1563 }
1564 //_____________________________________________________________________
1565 void AliForwardFlowTaskQC::VertexBin::GetLimits(Int_t bin, Int_t& aLow, Int_t& aHigh,
1566                                                 Int_t& bLow, Int_t& bHigh) const
1567 {
1568   //
1569   //  Get the limits for the 3 correlator method
1570   //
1571   //  Parameters:
1572   //   bin  : reference bin #
1573   //   aLow : Lowest bin to be included in v_A calculations
1574   //   aHigh: Highest bin to be included in v_A calculations
1575   //   bLow : Lowest bin to be included in v_B calculations
1576   //   bHigh: Highest bin to be included in v_B calculations
1577   //
1578   if ((fFlags & kFMD)) {
1579     switch(bin) {
1580       case 0:
1581         aLow = 14; aHigh = 15;
1582         bLow = 20; bHigh = 22;
1583         break;
1584       case 1:
1585         aLow = 16; aHigh = 16;
1586         bLow = 21; bHigh = 22;
1587         break;
1588       case 2:
1589         aLow =  6; aHigh =  7;
1590         bLow = 21; bHigh = 22;
1591         break;
1592       case 3:
1593         aLow =  6; aHigh =  7;
1594         bLow = 12; bHigh = 12;
1595         break;
1596       case 4:
1597         aLow =  6; aHigh =  8;
1598         bLow = 13; bHigh = 14;
1599         break;
1600       default:
1601         AliFatal(Form("No limits for this eta region! (%d)", bin));
1602     }
1603   }
1604   else if ((fFlags & kVZERO)) {
1605     switch(bin) {
1606       case 0:
1607         aLow =  6; aHigh = 13;
1608         bLow = 17; bHigh = 18;
1609         break;
1610       case 1:
1611         aLow =  6; aHigh =  9;
1612         bLow = 17; bHigh = 18;
1613         break;
1614       case 2:
1615         aLow =  2; aHigh =  3;
1616         bLow = 17; bHigh = 18;
1617         break;
1618       case 3:
1619         aLow =  2; aHigh =  3;
1620         bLow =  6; bHigh =  9;
1621         break;
1622       case 4:
1623         aLow =  2; aHigh =  3;
1624         bLow =  6; bHigh = 13;
1625         break;
1626       default:
1627         AliFatal(Form("No limits for this eta region! (%d)", bin));
1628     }
1629   }
1630   // Try to catch cases where fEtaLimits and these values do not correspond to each other
1631   if (aHigh > fCumuNUARef->GetNbinsX() || bHigh > fCumuNUARef->GetNbinsX())
1632     AliFatal(Form("Limits outside vtx range! (%d) - aHigh = %d, bHigh = %d, Nbins = %d",
1633       bin, aHigh, bHigh, fCumuNUARef->GetNbinsX()));
1634 }
1635 //_____________________________________________________________________
1636 void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate3Cor(Double_t cent)
1637 {
1638   //
1639   //  Calculate the Q cumulant up to order fMaxMoment
1640   //
1641   //  Parameters:
1642   //   cent: centrality of event
1643   //
1644   if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1645
1646   // Fill out NUA hists
1647   for (Int_t etaBin = 1; etaBin <= fCumuRef->GetNbinsX(); etaBin++) {
1648     Double_t eta = fCumuRef->GetXaxis()->GetBinCenter(etaBin);
1649     if (fCumuRef->GetBinContent(etaBin, 0) == 0) continue;
1650     for (Int_t qBin = 0; qBin <= fCumuRef->GetNbinsY(); qBin++) {
1651       fCumuNUARef->Fill(eta, cent, Double_t(qBin), fCumuRef->GetBinContent(etaBin, qBin));
1652     }
1653   }
1654   for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1655     Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1656     if (fCumuDiff->GetBinContent(etaBin, 0) == 0) continue;
1657     for (Int_t qBin = 0; qBin <= fCumuDiff->GetNbinsY(); qBin++) {
1658       fCumuNUADiff->Fill(eta, cent, Double_t(qBin), fCumuDiff->GetBinContent(etaBin, qBin));
1659     }
1660   }
1661
1662   // We create the objects needed for the analysis
1663   TH3D* cumuRef = 0;
1664   TH3D* cumuDiff = 0;
1665   // For each n we loop over the hists
1666   for (Int_t n = 2; n <= fMaxMoment; n++) {
1667     cumuRef  = (TH3D*)fCumuHists.Get('r',n);
1668     cumuDiff = (TH3D*)fCumuHists.Get('d',n);
1669
1670     // Per mom. quantities
1671     Int_t prevLim = 0;
1672     Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
1673     Double_t dQnReA = 0, dQnImA = 0, multA = 0;
1674     Double_t dQnReB = 0, dQnImB = 0, multB = 0;
1675     Double_t two = 0, w2 = 0;
1676     for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1677       Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1678       if (fEtaLims[prevLim] < eta) {
1679         GetLimits(prevLim, aLow, aHigh, bLow, bHigh);
1680         prevLim++;
1681         multA = 0; dQnReA = 0; dQnImA = 0;
1682         multB = 0; dQnReB = 0; dQnImB = 0;
1683         // Reference flow
1684         for (Int_t a = aLow; a <= aHigh; a++) {
1685           multA  += fCumuRef->GetBinContent(a, 0);
1686           dQnReA += fCumuRef->GetBinContent(a, GetBinNumberCos(n));
1687           dQnImA += fCumuRef->GetBinContent(a, GetBinNumberSin(n));
1688         }
1689         for (Int_t b = bLow; b <= bHigh; b++) {
1690           multB  += fCumuRef->GetBinContent(b, 0);
1691           dQnReB += fCumuRef->GetBinContent(b, GetBinNumberCos(n));
1692           dQnImB += fCumuRef->GetBinContent(b, GetBinNumberSin(n));
1693         }
1694         // The reference flow is calculated
1695         // 2-particle
1696         w2 = multA * multB;
1697         two = dQnReA*dQnReB + dQnImA*dQnImB;
1698       } // End of reference flow
1699       cumuRef->Fill(eta, cent, kW2Two, two);
1700       cumuRef->Fill(eta, cent, kW2, w2);
1701
1702       // For each etaBin bin the necessary values for differential flow is calculated
1703       Double_t mp = fCumuDiff->GetBinContent(etaBin, 0);
1704       Double_t pnRe = fCumuDiff->GetBinContent(etaBin, GetBinNumberCos(n));
1705       Double_t pnIm = fCumuDiff->GetBinContent(etaBin, GetBinNumberSin(n));
1706       if (mp == 0) continue;
1707
1708       // Differential flow calculations for each eta bin is done:
1709       // 2-particle differential flow
1710       Double_t w2pA = mp * multA;
1711       Double_t twoPrimeA = pnRe*dQnReA + pnIm*dQnImA;
1712       cumuDiff->Fill(eta, cent, kW2Two, twoPrimeA);
1713       cumuDiff->Fill(eta, cent, kW2, w2pA);
1714
1715       Double_t w2pB = mp * multB;
1716       Double_t twoPrimeB = pnRe*dQnReB + pnIm*dQnImB;
1717       cumuDiff->Fill(eta, cent, kW4Four, twoPrimeB);
1718       cumuDiff->Fill(eta, cent, kW4, w2pB);
1719      } // End of eta loop
1720     // Event count
1721     cumuRef->Fill(-7., cent, -0.5, 1.);
1722   } // End of moment loop
1723   return;
1724
1725 }
1726 //_____________________________________________________________________
1727 void AliForwardFlowTaskQC::VertexBin::CumulantsTerminate(TList* inlist, TList* outlist)
1728 {
1729   //
1730   //  Finalizes the Q cumulant calculations
1731   //
1732   //  Parameters:
1733   //   inlist: input sumlist
1734   //   outlist: output result list
1735   //
1736
1737   // Re-find cumulants hist if Terminate is called separately
1738   if (!fCumuHists.IsConnected()) {
1739     TList* list = (TList*)inlist->FindObject(Form("%svertex_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1740     fCumuHists.ConnectList(Form("%sCumu_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)), list);
1741
1742     if (!fCumuNUARef)
1743       fCumuNUARef = (TH3D*)list->FindObject(Form("%s_vertex_%d_%d%s_cumuNUARef", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1744     if (!fCumuNUADiff)
1745       fCumuNUADiff = (TH3D*)list->FindObject(Form("%s_vertex_%d_%d%s_cumuNUADiff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1746   }
1747   // Clone to avoid normalization problems when redoing terminate locally
1748   fCumuNUARef = (TH3D*)fCumuNUARef->Clone(Form("%s_vertex_%d_%d%s_cumuNUARefNorm", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1749   fCumuNUADiff = (TH3D*)fCumuNUADiff->Clone(Form("%s_vertex_%d_%d%s_cumuNUADiffNorm", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1750
1751   // Diagnostics histograms
1752   TH2I* quality = (TH2I*)outlist->FindObject(Form("hQCQuality%s%s", fType.Data(), GetQCType(fFlags)));
1753   if (!quality) {
1754     quality = MakeQualityHist(Form("hQCQuality%s%s", fType.Data(), GetQCType(fFlags)));
1755     outlist->Add(quality);
1756   }
1757   TH1D* cent = (TH1D*)outlist->FindObject(Form("%s%s_cent", fType.Data(), GetQCType(fFlags)));
1758   if (!cent) {
1759     cent = new TH1D(Form("%s%s_cent", fType.Data(), GetQCType(fFlags)),
1760                     Form("%s%s_cent", fType.Data(), GetQCType(fFlags)),
1761                 fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXmin(), fCumuNUARef->GetYaxis()->GetXmax());
1762     cent->GetXaxis()->Set(fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXbins()->GetArray());
1763     outlist->Add(cent);
1764   }
1765   TH2D* dNdetaRef = (TH2D*)outlist->FindObject(Form("%s%s_dNdetaRef", fType.Data(), GetQCType(fFlags)));
1766   if (!dNdetaRef) {
1767     dNdetaRef = new TH2D(Form("%s%s_dNdetaRef", fType.Data(), GetQCType(fFlags)),
1768                                Form("%s%s_dNdetaRef", fType.Data(), GetQCType(fFlags)),
1769                 fCumuNUARef->GetNbinsX(), fCumuNUARef->GetXaxis()->GetXmin(), fCumuNUARef->GetXaxis()->GetXmax(),
1770                 fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXmin(), fCumuNUARef->GetYaxis()->GetXmax());
1771     dNdetaRef->GetYaxis()->Set(fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXbins()->GetArray());
1772     dNdetaRef->Sumw2();
1773     outlist->Add(dNdetaRef);
1774   }
1775   TH2D* dNdetaDiff = (TH2D*)outlist->FindObject(Form("%s%s_dNdetaDiff", fType.Data(), GetQCType(fFlags)));
1776   if (!dNdetaDiff) {
1777     dNdetaDiff = new TH2D(Form("%s%s_dNdetaDiff", fType.Data(), GetQCType(fFlags)),
1778                                 Form("%s%s_dNdetaDiff", fType.Data(), GetQCType(fFlags)),
1779                 fCumuNUADiff->GetNbinsX(), fCumuNUADiff->GetXaxis()->GetXmin(), fCumuNUADiff->GetXaxis()->GetXmax(),
1780                 fCumuNUADiff->GetNbinsY(), fCumuNUADiff->GetYaxis()->GetXmin(), fCumuNUADiff->GetYaxis()->GetXmax());
1781     dNdetaDiff->GetYaxis()->Set(fCumuNUADiff->GetNbinsY(), fCumuNUADiff->GetYaxis()->GetXbins()->GetArray());
1782     dNdetaDiff->Sumw2();
1783     outlist->Add(dNdetaDiff);
1784   }
1785
1786   // Setting up outputs
1787   // Create output lists and diagnostics
1788   TList* vtxList = (TList*)outlist->FindObject("vtxList");
1789   if (!vtxList) {
1790     vtxList = new TList();
1791     vtxList->SetName("vtxList");
1792     outlist->Add(vtxList);
1793   }
1794   vtxList->Add(fCumuNUARef);
1795   vtxList->Add(fCumuNUADiff);
1796
1797   // Setup output profiles
1798   CumuHistos cumu2(fMaxMoment, ((fFlags & kNUAcorr) ? 2 : 0));
1799   CumuHistos cumu4(fMaxMoment, ((fFlags & kNUAcorr) ? 1 : 0));
1800
1801   cumu2.ConnectList(Form("%sQC2_Cumu%s_vtx_%d_%d", fType.Data(), GetQCType(fFlags), fVzMin, fVzMax), vtxList);
1802   if ((fFlags & kStdQC))
1803     cumu4.ConnectList(Form("%sQC4_Cumu%s_vtx_%d_%d", fType.Data(), GetQCType(fFlags), fVzMin, fVzMax), vtxList);
1804
1805   for (Int_t n = 2; n <= fMaxMoment; n++) {
1806     // 2-particle
1807     cumu2.Add(MakeOutputHist(2, n, "Ref", CumuHistos::kNoNUA));
1808     if ((fFlags & k3Cor)){
1809       cumu2.Add(MakeOutputHist(2, n, "DiffA", CumuHistos::kNoNUA));
1810       cumu2.Add(MakeOutputHist(2, n, "DiffB", CumuHistos::kNoNUA));
1811     } else {
1812       cumu2.Add(MakeOutputHist(2, n, "Diff", CumuHistos::kNoNUA));
1813     }
1814     // 4-particle
1815     if ((fFlags & kStdQC)) {
1816       cumu4.Add(MakeOutputHist(4, n, "Ref", CumuHistos::kNoNUA));
1817       cumu4.Add(MakeOutputHist(4, n, "Diff", CumuHistos::kNoNUA));
1818     }
1819   } // End of v_n result loop
1820   // NUA corrected
1821   if ((fFlags & kNUAcorr)) {
1822     for (Int_t n = 2; n <= fMaxMoment; n++) {
1823       // 2-particle
1824       cumu2.Add(MakeOutputHist(2, n, "Ref", CumuHistos::kNUAOld));
1825       if ((fFlags & k3Cor)) {
1826         cumu2.Add(MakeOutputHist(2, n, "DiffA", CumuHistos::kNUAOld));
1827         cumu2.Add(MakeOutputHist(2, n, "DiffB", CumuHistos::kNUAOld));
1828       } else {
1829         cumu2.Add(MakeOutputHist(2, n, "Diff", CumuHistos::kNUAOld));
1830       }
1831       // 4-particle
1832       if ((fFlags & kStdQC)) {
1833         cumu4.Add(MakeOutputHist(4, n, "Ref", CumuHistos::kNUAOld));
1834         cumu4.Add(MakeOutputHist(4, n, "Diff", CumuHistos::kNUAOld));
1835       }
1836     }
1837     for (Int_t n = 2; n <= fMaxMoment; n++) {
1838       // 2-particle
1839       cumu2.Add(MakeOutputHist(2, n, "Ref", CumuHistos::kNUA));
1840       if ((fFlags & k3Cor)) {
1841         cumu2.Add(MakeOutputHist(2, n, "DiffA", CumuHistos::kNUA));
1842         cumu2.Add(MakeOutputHist(2, n, "DiffB", CumuHistos::kNUA));
1843       } else {
1844         cumu2.Add(MakeOutputHist(2, n, "Diff", CumuHistos::kNUA));
1845       }
1846     }
1847   }
1848
1849   // Calculating the cumulants
1850   if ((fFlags & k3Cor)) {
1851     Calculate3CorFlow(cumu2, quality, cent, dNdetaRef, dNdetaDiff);
1852   } else {
1853     CalculateReferenceFlow(cumu2, cumu4, quality, cent, dNdetaRef);
1854     CalculateDifferentialFlow(cumu2, cumu4, quality, dNdetaDiff);
1855   }
1856   if ((fFlags & kNUAcorr)) {
1857     SolveCoupledFlowEquations(cumu2, 'r');
1858     if ((fFlags & k3Cor)) {
1859       SolveCoupledFlowEquations(cumu2, 'a');
1860       SolveCoupledFlowEquations(cumu2, 'b');
1861     } else {
1862       SolveCoupledFlowEquations(cumu2, 'd');
1863     }
1864   }
1865
1866   // Add to output for immediate viewing - individual vtx bins are used for final results
1867   AddVertexBins(cumu2, outlist, ((fFlags & kNUAcorr) ? 2 : 0));
1868   if ((fFlags & kStdQC)) AddVertexBins(cumu4, outlist, ((fFlags & kNUAcorr) ? 1 : 0));
1869
1870   // Setup NUA diagnoastics histograms
1871   TList* nualist = (TList*)outlist->FindObject("NUATerms");
1872   if (!nualist) {
1873     nualist = new TList();
1874     nualist->SetName("NUATerms");
1875     outlist->Add(nualist);
1876   }
1877   // Reference
1878   TH2D* nuaRef = (TH2D*)nualist->FindObject(Form("%sReferenceNUA%s", fType.Data(), GetQCType(fFlags)));
1879   TH2D* temp = 0;
1880   if (!nuaRef) {
1881     nuaRef = (TH2D*)fCumuNUARef->Project3D("yz");
1882     nuaRef->Scale(1./fCumuNUARef->GetNbinsX());
1883     nuaRef->SetName(Form("%sReferenceNUA%s", fType.Data(), GetQCType(fFlags)));
1884     nuaRef->SetTitle(Form("%sReferenceNUA%s", fType.Data(), GetQCType(fFlags)));
1885     nualist->Add(nuaRef);
1886   } else {
1887     temp = (TH2D*)fCumuNUARef->Project3D("yz");
1888     temp->Scale(1./fCumuNUARef->GetNbinsX());
1889     nuaRef->Add(temp);
1890     delete temp;
1891   }
1892   // Filling in underflow to make scaling possible in Terminate()
1893   nuaRef->Fill(0., -1., 1.);
1894   // Differential
1895   TH2D* nuaDiff = (TH2D*)nualist->FindObject(Form("%sDifferentialNUA%s", fType.Data(), GetQCType(fFlags)));
1896   if (!nuaDiff) {
1897     nuaDiff = (TH2D*)fCumuNUADiff->Project3D("yz");
1898     nuaDiff->SetName(Form("%sDifferentialNUA%s", fType.Data(), GetQCType(fFlags)));
1899     nuaDiff->SetTitle(Form("%sDifferentialNUA%s", fType.Data(), GetQCType(fFlags)));
1900     nualist->Add(nuaDiff);
1901   } else {
1902     temp = (TH2D*)fCumuNUADiff->Project3D("yz");
1903     nuaDiff->Add(temp);
1904     delete temp;
1905   }
1906   // Filling in underflow to make scaling possible in Terminate()
1907   nuaDiff->Fill(0., -1., 1.);
1908
1909   return;
1910 }
1911 //_____________________________________________________________________
1912 void AliForwardFlowTaskQC::VertexBin::CalculateReferenceFlow(CumuHistos& cumu2h, CumuHistos& cumu4h, TH2I* quality,
1913                                                              TH1D* chist, TH2D* dNdetaRef) const
1914 {
1915   //
1916   //  Calculates the reference flow
1917   //
1918   //  Parameters:
1919   //   cumu2h: CumuHistos object with QC{2} cumulants
1920   //   cumu4h: CumuHistos object with QC{4} cumulants
1921   //   quality: Histogram for success rate of cumulants
1922   //   chist: Centrality histogram
1923   //   dNdetaRef: dN/deta histogram for estimating multiplicity used for ref calculations
1924   //
1925
1926   // Normalizing common NUA hists
1927   for (Int_t cBin = 1; cBin <= fCumuNUARef->GetNbinsY(); cBin++) {
1928     Double_t cent = fCumuNUARef->GetYaxis()->GetBinCenter(cBin);
1929     for (Int_t eBin = 1; eBin <= fCumuNUARef->GetNbinsX(); eBin++) {
1930       Double_t eta = fCumuNUARef->GetXaxis()->GetBinCenter(eBin);
1931       Double_t mult = fCumuNUARef->GetBinContent(eBin, cBin, 0);
1932       if (mult == 0) continue;
1933       for (Int_t qBin = 1; qBin <= fCumuNUARef->GetNbinsZ(); qBin++) {
1934         fCumuNUARef->SetBinContent(eBin, cBin, qBin, fCumuNUARef->GetBinContent(eBin, cBin, qBin)/mult);
1935         fCumuNUARef->SetBinError(eBin, cBin, qBin, fCumuNUARef->GetBinError(eBin, cBin, qBin)/mult);
1936       }
1937       // Fill dN/deta diagnostics
1938       dNdetaRef->Fill(eta, cent, mult);
1939     }
1940   }
1941
1942   // For flow calculations
1943   TH3D* cumuRef = 0;
1944   TH2D* cumu2 = 0;
1945   TH2D* cumu2NUAold = 0;
1946   TH2D* cumu4 = 0;
1947   TH2D* cumu4NUA = 0;
1948   Int_t qualityFactor = ((fFlags & kStdQC) ? 8 : 4); // used for correctly filling in quality hists
1949   // Loop over cumulant histogram for final calculations
1950   for (Int_t n = 2; n <= fMaxMoment; n++) { // Moment loop begins
1951     cumu2 = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNoNUA);
1952     if ((fFlags & kNUAcorr)) {
1953       cumu2NUAold = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNUAOld);
1954     }
1955     if ((fFlags & kStdQC)) {
1956       cumu4 = (TH2D*)cumu4h.Get('r', n, CumuHistos::kNoNUA);
1957       if ((fFlags & kNUAcorr)) cumu4NUA = (TH2D*)cumu4h.Get('r', n, CumuHistos::kNUAOld);
1958     }
1959     cumuRef  = (TH3D*)fCumuHists.Get('r', n);
1960     // Begin loops
1961     for (Int_t cBin = 1; cBin <= cumuRef->GetNbinsY(); cBin++) { // Centrality loop begins
1962       Double_t cent = cumuRef->GetYaxis()->GetBinCenter(cBin);
1963       if (n == 2) chist->Fill(cent, cumuRef->GetBinContent(0, cBin, 0));
1964       if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), n, cent));
1965       for (Int_t etaBin = 1; etaBin <= cumuRef->GetNbinsX(); etaBin++) { // Eta loop begins
1966         Double_t eta = cumuRef->GetXaxis()->GetBinCenter(etaBin);
1967         Double_t refEta = eta;
1968         Int_t refEtaBinA = fCumuNUARef->GetXaxis()->FindBin(refEta);
1969         if ((fFlags & kEtaGap)) refEta = -eta;
1970         Int_t refEtaBinB = fCumuNUARef->GetXaxis()->FindBin(refEta);
1971         // 2-particle reference flow
1972         Double_t w2Two = cumuRef->GetBinContent(refEtaBinA, cBin, kW2Two);
1973         Double_t w2 = cumuRef->GetBinContent(refEtaBinA, cBin, kW2);
1974         if (w2 == 0) continue;
1975         Double_t cosP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberCos(n));
1976         Double_t sinP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberSin(n));
1977         Double_t cosP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(n));
1978         Double_t sinP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(n));
1979         Double_t cos2nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberCos(2*n));
1980         Double_t sin2nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberSin(2*n));
1981         Double_t cos2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(2*n));
1982         Double_t sin2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(2*n));
1983         Double_t two = w2Two / w2;
1984         Double_t qc2 = two;
1985         if (qc2 >= 0) cumu2->Fill(eta, cent, TMath::Sqrt(qc2));
1986
1987         if ((fFlags & kNUAcorr)) {
1988           // Old NUA
1989           // With no eta gap the last two terms are <<cos(phi)>>^2 and <<sin(phi)>>^2,
1990           // with eta gap the different coverage is taken into account.
1991           // The next line covers both cases.
1992           qc2 -= cosP1nPhiA*cosP1nPhiB + sinP1nPhiA*sinP1nPhiB;
1993           // Extra NUA term from 2n cosines and sines
1994           Double_t den = 1-(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB);
1995           if (den != 0) qc2 /= den;
1996           else qc2 = 0;
1997         }
1998         if (qc2 <= 0) {
1999           if (fDebug > 0)
2000             AliInfo(Form("%s: QC_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f - skipping",
2001             fType.Data(), n, qc2, eta, cent));
2002           quality->Fill((n-2)*qualityFactor+2, Int_t(cent));
2003           continue;
2004         }
2005         Double_t vnTwo = TMath::Sqrt(qc2);
2006         if (!TMath::IsNaN(vnTwo)) {
2007           quality->Fill((n-2)*qualityFactor+1, Int_t(cent));
2008           if ((fFlags & kNUAcorr)) cumu2NUAold->Fill(eta, cent, vnTwo);
2009         }
2010
2011         if (!(fFlags & kStdQC)) continue;
2012         // 4-particle reference flow
2013         Double_t w4Four = cumuRef->GetBinContent(refEtaBinA, cBin, kW4Four);
2014         Double_t w4 = cumuRef->GetBinContent(refEtaBinA, cBin, kW4);
2015         Double_t multm1m2 = cumuRef->GetBinContent(refEtaBinA, cBin, k3pWeight);
2016         if (w4 == 0 || multm1m2 == 0) continue;
2017         Double_t cosP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2);
2018         Double_t sinP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2);
2019         Double_t cosP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2phi3m);
2020         Double_t sinP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2phi3m);
2021
2022         cosP1nPhi1P1nPhi2 /= w2;
2023         sinP1nPhi1P1nPhi2 /= w2;
2024         cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2025         sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2026         Double_t four = w4Four / w4;
2027         Double_t qc4 = four-2.*TMath::Power(two,2.);
2028         if (qc4 < 0) cumu4->Fill(eta, cent, TMath::Power(-qc4, 0.25));
2029
2030         if ((fFlags & kNUAcorr)) {
2031            qc4 += - 4.*cosP1nPhiA*cosP1nPhi1M1nPhi2M1nPhi3
2032                   + 4.*sinP1nPhiA*sinP1nPhi1M1nPhi2M1nPhi3-TMath::Power(cosP1nPhi1P1nPhi2,2.)-TMath::Power(sinP1nPhi1P1nPhi2,2.)
2033                   + 4.*cosP1nPhi1P1nPhi2*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
2034                   + 8.*sinP1nPhi1P1nPhi2*sinP1nPhiA*cosP1nPhiA
2035                   + 8.*two*(TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.))
2036                   - 6.*TMath::Power((TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.)),2.);
2037         }
2038         if (qc4 >= 0) {
2039           if (fDebug > 0)
2040             AliInfo(Form("%s: QC_%d{4} = %1.3f for eta = %1.2f and centrality %3.1f - skipping",
2041             fType.Data(), n, qc2, eta, cent));
2042           quality->Fill((n-2)*qualityFactor+6, Int_t(cent));
2043           continue;
2044         }
2045         Double_t vnFour = TMath::Power(-qc4, 0.25);
2046         if (!TMath::IsNaN(vnFour*multm1m2)) {
2047           quality->Fill((n-2)*qualityFactor+5, Int_t(cent));
2048           if ((fFlags & kNUAcorr)) cumu4NUA->Fill(eta, cent, vnFour);
2049         }
2050       } // End of eta
2051     } // End of cent
2052   } // End of moment
2053
2054   return;
2055 }
2056 //_____________________________________________________________________
2057 void AliForwardFlowTaskQC::VertexBin::CalculateDifferentialFlow(CumuHistos& cumu2h, CumuHistos& cumu4h,
2058                                                                 TH2I* quality, TH2D* dNdetaDiff) const
2059 {
2060   //
2061   //  Calculates the differential flow
2062   //
2063   //  Parameters:
2064   //   cumu2h: CumuHistos object with QC{2} cumulants
2065   //   cumu4h: CumuHistos object with QC{4} cumulants
2066   //   quality: Histogram for success rate of cumulants
2067   //   dNdetaDiff: dN/deta histogram for estimating multiplicity used for diff calculations
2068   //
2069
2070   for (Int_t cBin = 1; cBin <= fCumuNUADiff->GetNbinsY(); cBin++) {
2071     Double_t cent = fCumuNUADiff->GetYaxis()->GetBinCenter(cBin);
2072     for (Int_t eBin = 1; eBin <= fCumuNUADiff->GetNbinsX(); eBin++) {
2073       Double_t eta = fCumuNUADiff->GetXaxis()->GetBinCenter(eBin);
2074       Double_t mult = fCumuNUADiff->GetBinContent(eBin, cBin, 0);
2075       if (mult == 0) continue;
2076       for (Int_t qBin = 1; qBin <= fCumuNUADiff->GetNbinsZ(); qBin++) {
2077         fCumuNUADiff->SetBinContent(eBin, cBin, qBin, fCumuNUADiff->GetBinContent(eBin, cBin, qBin)/mult);
2078         fCumuNUADiff->SetBinError(eBin, cBin, qBin, fCumuNUADiff->GetBinError(eBin, cBin, qBin)/mult);
2079       }
2080       dNdetaDiff->Fill(eta, cent, mult);
2081     }
2082   }
2083
2084   // For flow calculations
2085   TH3D* cumuRef = 0;
2086   TH3D* cumuDiff = 0;
2087   TH2D* cumu2 = 0;
2088   TH2D* cumu2NUAold = 0;
2089   TH2D* cumu4 = 0;
2090   TH2D* cumu4NUA = 0;
2091   Int_t qualityFactor = ((fFlags & kStdQC) ? 8 : 4); // used for correctly filling in quality hists
2092   // Loop over cumulant histogram for final calculations
2093   for (Int_t n = 2; n <= fMaxMoment; n++) { // Moment loop begins
2094     cumu2 = (TH2D*)cumu2h.Get('d', n, CumuHistos::kNoNUA);
2095     if ((fFlags & kNUAcorr)) {
2096       cumu2NUAold = (TH2D*)cumu2h.Get('d', n, CumuHistos::kNUAOld);
2097     }
2098     if ((fFlags & kStdQC)) {
2099       cumu4 = (TH2D*)cumu4h.Get('d',n);
2100       if ((fFlags & kNUAcorr)) cumu4NUA = (TH2D*)cumu4h.Get('d', n, CumuHistos::kNUAOld);
2101     }
2102     cumuRef  = (TH3D*)fCumuHists.Get('r',n);
2103     cumuDiff = (TH3D*)fCumuHists.Get('d',n);
2104     for (Int_t cBin = 1; cBin <= cumuDiff->GetNbinsY(); cBin++) { // Centrality loop begins
2105       Double_t cent = cumuDiff->GetYaxis()->GetBinCenter(cBin);
2106       if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), n, cent));
2107       for (Int_t etaBin = 1; etaBin <= cumuDiff->GetNbinsX(); etaBin++) { // Eta loop begins
2108         Double_t eta = cumuDiff->GetXaxis()->GetBinCenter(etaBin);
2109         Double_t refEta = eta;
2110         Int_t refEtaBinA = fCumuNUARef->GetXaxis()->FindBin(refEta);
2111         if ((fFlags & kEtaGap)) refEta = -eta;
2112         Int_t refEtaBinB = fCumuNUARef->GetXaxis()->FindBin(refEta);
2113
2114         // Reference objects
2115         Double_t w2 = cumuRef->GetBinContent(refEtaBinA, cBin, kW2);
2116         if (w2 == 0) continue;
2117         Double_t two = cumuRef->GetBinContent(refEtaBinA, cBin, kW2Two);
2118         two /= w2;
2119         Double_t cosP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberCos(n));
2120         Double_t sinP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberSin(n));
2121         Double_t cosP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(n));
2122         Double_t sinP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(n));
2123         Double_t cos2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(2*n));
2124         Double_t sin2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(2*n));
2125
2126         // 2-particle differential flow
2127         Double_t w2pTwoPrime = cumuDiff->GetBinContent(etaBin, cBin, kW2Two);
2128         Double_t w2p = cumuDiff->GetBinContent(etaBin, cBin, kW2);
2129         if (w2p == 0) continue;
2130         Double_t cosP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(n));
2131         Double_t sinP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(n));
2132         Double_t cos2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(2*n));
2133         Double_t sin2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(2*n));
2134         Double_t twoPrime = w2pTwoPrime / w2p;
2135
2136         Double_t qc2Prime = twoPrime;
2137         cumu2->Fill(eta, cent, qc2Prime);
2138         if ((fFlags & kNUAcorr)) {
2139           // Old nua
2140           qc2Prime -= cosP1nPsi*cosP1nPhiB + sinP1nPsi*sinP1nPhiB;
2141           // Extra NUA term from 2n cosines and sines
2142           qc2Prime /= (1.-(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
2143         }
2144         if (!TMath::IsNaN(qc2Prime)) {
2145           quality->Fill((n-2)*qualityFactor+3, Int_t(cent));
2146           if ((fFlags & kNUAcorr)) cumu2NUAold->Fill(eta, cent, qc2Prime);
2147         }
2148         else
2149           quality->Fill((n-2)*qualityFactor+4, Int_t(cent));
2150         if (fDebug > 1)
2151           AliInfo(Form("%s: QC'_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f",
2152           fType.Data(), n, qc2Prime, eta, cent));
2153
2154         if (!(fFlags & kStdQC)) continue;
2155         // Reference objects
2156         Double_t cosP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2);
2157         Double_t sinP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2);
2158         Double_t cosP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2phi3m);
2159         Double_t sinP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2phi3m);
2160         Double_t multm1m2 = cumuRef->GetBinContent(refEtaBinA, cBin, k3pWeight);
2161         cosP1nPhi1P1nPhi2 /= w2;
2162         sinP1nPhi1P1nPhi2 /= w2;
2163         cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2164         sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2165
2166         // 4-particle differential flow
2167         Double_t w4pFourPrime = cumuDiff->GetBinContent(etaBin, cBin, kW4Four);
2168         Double_t w4p = cumuDiff->GetBinContent(etaBin, cBin, kW4);
2169         Double_t mpqMult = cumuDiff->GetBinContent(etaBin, cBin, k3pWeight);
2170         if (w4p == 0 || mpqMult == 0) continue;
2171         Double_t cosP1nPsi1P1nPhi2 = cumuDiff->GetBinContent(etaBin, cBin, kCosphi1phi2);
2172         Double_t sinP1nPsi1P1nPhi2 = cumuDiff->GetBinContent(etaBin, cBin, kSinphi1phi2);
2173         Double_t cosP1nPsi1M1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kCosphi1phi2phi3m);
2174         Double_t sinP1nPsi1M1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kSinphi1phi2phi3m);
2175         Double_t cosP1nPsi1P1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kCosphi1phi2phi3p);
2176         Double_t sinP1nPsi1P1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kSinphi1phi2phi3p);
2177
2178         cosP1nPsi1P1nPhi2 /= w2p;
2179         sinP1nPsi1P1nPhi2 /= w2p;
2180         cosP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
2181         sinP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
2182         cosP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
2183         sinP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
2184
2185         Double_t fourPrime = w4pFourPrime / w4p;
2186         Double_t qc4Prime = fourPrime-2.*twoPrime*two;
2187         if (cumu4) cumu4->Fill(eta, cent, qc4Prime);
2188
2189         if ((fFlags & kNUAcorr)) {
2190           qc4Prime += - cosP1nPsi*cosP1nPhi1M1nPhi2M1nPhi3
2191                       + sinP1nPsi*sinP1nPhi1M1nPhi2M1nPhi3
2192                       - cosP1nPhiA*cosP1nPsi1M1nPhi2M1nPhi3
2193                       + sinP1nPhiA*sinP1nPsi1M1nPhi2M1nPhi3
2194                       - 2.*cosP1nPhiA*cosP1nPsi1P1nPhi2M1nPhi3
2195                       - 2.*sinP1nPhiA*sinP1nPsi1P1nPhi2M1nPhi3
2196                       - cosP1nPsi1P1nPhi2*cosP1nPhi1P1nPhi2
2197                       - sinP1nPsi1P1nPhi2*sinP1nPhi1P1nPhi2
2198                       + 2.*cosP1nPhi1P1nPhi2*(cosP1nPsi*cosP1nPhiA-sinP1nPsi*sinP1nPhiA)
2199                       + 2.*sinP1nPhi1P1nPhi2*(cosP1nPsi*sinP1nPhiA+sinP1nPsi*cosP1nPhiA)
2200                       + 4.*two*(cosP1nPsi*cosP1nPhiA+sinP1nPsi*sinP1nPhiA)
2201                       + 2.*cosP1nPsi1P1nPhi2*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
2202                       + 4.*sinP1nPsi1P1nPhi2*cosP1nPhiA*sinP1nPhiA
2203                       + 4.*twoPrime*(TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.))
2204                       - 6.*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
2205                       * (cosP1nPsi*cosP1nPhiA-sinP1nPsi*sinP1nPhiA)
2206                       - 12.*cosP1nPhiA*sinP1nPhiA
2207                       * (sinP1nPsi*cosP1nPhiA+cosP1nPsi*sinP1nPhiA);
2208         }
2209 //      Double_t vnFourDiff = - qc4Prime / TMath::Power(-qc4, 0.75);
2210         if (!TMath::IsNaN(qc4Prime*mpqMult)) {
2211           quality->Fill((n-2)*qualityFactor+7, Int_t(cent));
2212           if (cumu4NUA) cumu4NUA->Fill(eta, cent, qc4Prime);
2213         }
2214         else
2215           quality->Fill((n-2)*qualityFactor+8, Int_t(cent));
2216         if (fDebug > 1)
2217           AliInfo(Form("%s: v_%d{4} = %1.3f for eta = %1.2f and centrality %3.1f",
2218           fType.Data(), n, qc4Prime, eta, cent));
2219       } // End of eta loop
2220     } // End of centrality loop
2221   } // End of moment
2222
2223   return;
2224 }
2225 //_____________________________________________________________________
2226 void AliForwardFlowTaskQC::VertexBin::Calculate3CorFlow(CumuHistos& cumu2h, TH2I* quality, TH1D* chist,
2227                                                         TH2D* dNdetaRef, TH2D* dNdetaDiff) const
2228 {
2229   //
2230   //  Calculates the 3 sub flow
2231   //
2232   //  Parameters:
2233   //   cumu2h: CumuHistos object with QC{2} cumulants
2234   //   quality: Histogram for success rate of cumulants
2235   //   chist: Centrality histogram
2236   //   dNdetaDiff: dN/deta histogram for estimating multiplicity used for diff calculations
2237   //
2238
2239   // For flow calculations
2240   TH3D* cumuRef = 0;
2241   TH3D* cumuDiff = 0;
2242   TH2D* cumu2r = 0;
2243   TH2D* cumu2rNUAold = 0;
2244   TH2D* cumu2a = 0;
2245   TH2D* cumu2aNUAold = 0;
2246   TH2D* cumu2b = 0;
2247   TH2D* cumu2bNUAold = 0;
2248   // Loop over cumulant histogram for final calculations
2249   for (Int_t n = 2; n <= fMaxMoment; n++) { // Moment loop begins
2250     cumu2r = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNoNUA);
2251     cumu2a = (TH2D*)cumu2h.Get('a', n, CumuHistos::kNoNUA);
2252     cumu2b = (TH2D*)cumu2h.Get('b', n, CumuHistos::kNoNUA);
2253     if ((fFlags & kNUAcorr)) {
2254       cumu2rNUAold = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNUAOld);
2255       cumu2aNUAold = (TH2D*)cumu2h.Get('a', n, CumuHistos::kNUAOld);
2256       cumu2bNUAold = (TH2D*)cumu2h.Get('b', n, CumuHistos::kNUAOld);
2257     }
2258     cumuRef  = (TH3D*)fCumuHists.Get('r',n);
2259     cumuDiff = (TH3D*)fCumuHists.Get('d',n);
2260     for (Int_t cBin = 1; cBin <= cumuRef->GetNbinsY(); cBin++) { // Centrality loop begins
2261       Double_t cent = cumuRef->GetYaxis()->GetBinCenter(cBin);
2262       if (n == 2) chist->Fill(cent, cumuRef->GetBinContent(0, cBin, 0));
2263       if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), n, cent));
2264       // Here it starts!
2265       Int_t prevLim = 0;
2266       Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2267       Double_t cosP1nPhiA = 0;
2268       Double_t sinP1nPhiA = 0;
2269       Double_t cos2nPhiA = 0;
2270       Double_t sin2nPhiA = 0;
2271       Double_t cosP1nPhiB = 0;
2272       Double_t sinP1nPhiB = 0;
2273       Double_t cos2nPhiB = 0;
2274       Double_t sin2nPhiB = 0;
2275       Double_t multA = 0;
2276       Double_t multB = 0;
2277
2278       for (Int_t etaBin = 1; etaBin <= cumuDiff->GetNbinsX(); etaBin++) { // Eta loop begins
2279         Double_t eta = cumuDiff->GetXaxis()->GetBinCenter(etaBin);
2280         // 2-particle reference flow
2281         Double_t w2Two = cumuRef->GetBinContent(etaBin, cBin, kW2Two);
2282         Double_t w2 = cumuRef->GetBinContent(etaBin, cBin, kW2);
2283         if (w2 == 0) continue;
2284
2285         // Update NUA for new range!
2286         if (fEtaLims[prevLim] < eta) {
2287           GetLimits(prevLim, aLow, aHigh, bLow, bHigh);
2288           prevLim++;
2289           cosP1nPhiA = 0; sinP1nPhiA = 0; cos2nPhiA = 0; sin2nPhiA = 0; multA = 0;
2290           cosP1nPhiB = 0; sinP1nPhiB = 0; cos2nPhiB = 0; sin2nPhiB = 0; multB = 0;
2291           for (Int_t a = aLow; a <= aHigh; a++) {
2292             cosP1nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n));
2293             sinP1nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n));
2294             cos2nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(2*n));
2295             sin2nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(2*n));
2296             multA += fCumuNUARef->GetBinContent(a, cBin, 0);
2297           }
2298           for (Int_t b = bLow; b <= bHigh; b++) {
2299             cosP1nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n));
2300             sinP1nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n));
2301             cos2nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(2*n));
2302             sin2nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(2*n));
2303             multB += fCumuNUARef->GetBinContent(b, cBin, 0);
2304           }
2305           if (multA == 0 || multB == 0) {
2306             AliWarning(Form("Empty NUA values for 3Cor! (%s)", cumuRef->GetName()));
2307             continue;
2308           }
2309           cosP1nPhiA /= multA;
2310           sinP1nPhiA /= multA;
2311           cos2nPhiA /= multA;
2312           sin2nPhiA /= multA;
2313           cosP1nPhiB /= multB;
2314           sinP1nPhiB /= multB;
2315           cos2nPhiB /= multB;
2316           sin2nPhiB /= multB;
2317
2318           dNdetaRef->Fill(eta, cent, multA+multB);
2319         }
2320         Double_t two = w2Two / w2;
2321
2322         Double_t qc2 = two;
2323         if (qc2 >= 0) cumu2r->Fill(eta, cent, TMath::Sqrt(qc2));
2324
2325         if ((fFlags & kNUAcorr)) {
2326           // Old nua
2327           qc2 -= cosP1nPhiA*cosP1nPhiB + sinP1nPhiA*sinP1nPhiB;
2328           // Extra NUA term from 2n cosines and sines
2329           qc2 /= (1-(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB));
2330         }
2331         if (qc2 <= 0) {
2332           if (fDebug > 0)
2333             AliInfo(Form("%s: QC_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f - skipping",
2334             fType.Data(), n, qc2, eta, cent));
2335           quality->Fill((n-2)*4+2, Int_t(cent));
2336           continue;
2337         }
2338         Double_t vnTwo = TMath::Sqrt(qc2);
2339         if (!TMath::IsNaN(vnTwo)) {
2340           quality->Fill((n-2)*4+1, Int_t(cent));
2341           if ((fFlags & kNUAcorr)) cumu2rNUAold->Fill(eta, cent, vnTwo);
2342         }
2343
2344         // 2-particle differential flow
2345         Double_t w2pTwoPrimeA = cumuDiff->GetBinContent(etaBin, cBin, kW2Two);
2346         Double_t w2pA = cumuDiff->GetBinContent(etaBin, cBin, kW2);
2347         Double_t w2pTwoPrimeB = cumuDiff->GetBinContent(etaBin, cBin, kW4Four);
2348         Double_t w2pB = cumuDiff->GetBinContent(etaBin, cBin, kW4);
2349         if (w2pA == 0 || w2pB == 0) continue;
2350         Double_t cosP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(n));
2351         Double_t sinP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(n));
2352         Double_t cos2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(2*n));
2353         Double_t sin2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(2*n));
2354         Double_t mult = fCumuNUADiff->GetBinContent(etaBin, cBin, 0);
2355         if (mult == 0) continue;
2356         cosP1nPsi /= mult;
2357         sinP1nPsi /= mult;
2358         cos2nPsi /= mult;
2359         sin2nPsi /= mult;
2360         Double_t twoPrimeA = w2pTwoPrimeA / w2pA;
2361         Double_t twoPrimeB = w2pTwoPrimeB / w2pB;
2362         dNdetaDiff->Fill(eta, cent, mult);
2363
2364         Double_t qc2PrimeA = twoPrimeA;
2365         Double_t qc2PrimeB = twoPrimeB;
2366         if (qc2PrimeA*qc2PrimeB >= 0) {
2367           cumu2a->Fill(eta, cent, qc2PrimeA);
2368           cumu2b->Fill(eta, cent, qc2PrimeB);
2369         }
2370         if ((fFlags & kNUAcorr)) {
2371           // Old nua
2372           qc2PrimeA -= cosP1nPsi*cosP1nPhiA + sinP1nPsi*sinP1nPhiA;
2373           qc2PrimeB -= cosP1nPsi*cosP1nPhiB + sinP1nPsi*sinP1nPhiB; // Is this OK?
2374           // Extra NUA term from 2n cosines and sines
2375           if (cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA != 1.) qc2PrimeA /= (1.-(cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA));
2376           if (cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB != 1.) qc2PrimeB /= (1.-(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
2377         }
2378         if (!TMath::IsNaN(qc2PrimeA) && !TMath::IsNaN(qc2PrimeB) && qc2 != 0) {
2379         if (qc2PrimeA*qc2PrimeB >= 0) {
2380           quality->Fill((n-2)*4+3, Int_t(cent));
2381           if ((fFlags & kNUAcorr)) cumu2aNUAold->Fill(eta, cent, qc2PrimeA);
2382           if ((fFlags & kNUAcorr)) cumu2bNUAold->Fill(eta, cent, qc2PrimeB);
2383         }
2384       }
2385       else
2386         quality->Fill((n-2)*4+4, Int_t(cent));
2387         if (fDebug > 1)
2388           AliInfo(Form("%s: QC'a_%d{2} = %1.3f, QC'b_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f",
2389           fType.Data(), n, qc2PrimeA, n, qc2PrimeB, eta, cent));
2390       } // End of eta loop
2391     } // End of centrality loop
2392   } // End of moment
2393
2394   return;
2395 }
2396 //_____________________________________________________________________
2397 void AliForwardFlowTaskQC::VertexBin::SolveCoupledFlowEquations(CumuHistos& cumu, const Char_t type) const
2398 {
2399   //
2400   //  Function to solve the coupled flow equations
2401   //  We solve it by using matrix calculations:
2402   //  A*v_n = V => v_n = A^-1*V
2403   //  First we set up a TMatrixD container to make ROOT
2404   //  do the inversions in an efficient way, we multiply the current <<2>> estimates.
2405   //  Then we fill new TH2D's if the corrected <<2>>'s (cumuNUA object).
2406   //
2407   //  Parameters:
2408   //   cumu: CumuHistos object - uncorrected
2409   //   type: Reference ('r') or differential ('d') or ('a' or 'b') for 3 correlator
2410   //
2411
2412   // We start by declaring Matrix and vector objects, as their constructors are quite heavy
2413   TMatrixD mNUA(fMaxMoment-1, fMaxMoment-1);
2414   TVectorD vQC2(fMaxMoment-1);
2415
2416   for (Int_t cBin = 1; cBin <= cumu.Get(type, 2, CumuHistos::kNUAOld)->GetNbinsY(); cBin++) { // cent loop
2417     Double_t cent = cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinCenter(cBin);
2418     for (Int_t eBin = 1; eBin <= cumu.Get(type, 2, CumuHistos::kNUAOld)->GetNbinsX(); eBin++) { // eta loop
2419       Double_t eta = cumu.Get(type, 2, CumuHistos::kNUAOld)->GetXaxis()->GetBinCenter(eBin);
2420       mNUA.Zero(); // reset matrix
2421       vQC2.Zero(); // reset vector
2422       for (Int_t n = 0; n < fMaxMoment-1; n++) { // moment loop
2423         vQC2(n) = static_cast<TH2D*>(cumu.Get(type, n+2, CumuHistos::kNUAOld))->GetBinContent(eBin, cBin);
2424         if (type == 'r' || type == 'R') vQC2(n) *= vQC2(n); // going back to <<2>>
2425         for (Int_t m = 0; m < fMaxMoment-1; m++) { // cross moment
2426           mNUA(n,m) = CalculateNUAMatrixElement(n, m, type, eBin, cBin);
2427         } // End of cross moment loop
2428       } // End of moment loop
2429       // Invert matrix
2430       Double_t det = 0;
2431       mNUA.Invert(&det);
2432       // If determinant is non-zero we go with corrected results
2433       if (det != 0 ) vQC2 = mNUA*vQC2;
2434       else AliWarning(Form("Determinant == 0 - cent: %d-%d, eta: %f, type: '%c', data: %s, vtx: %d-%d%s",
2435                       Int_t(cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinLowEdge(cBin)),
2436                       Int_t(cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinUpEdge(cBin)),
2437                       cumu.Get(type, 2, CumuHistos::kNUAOld)->GetXaxis()->GetBinCenter(eBin),
2438                       type, fType.Data(), fVzMin, fVzMax,
2439                       GetQCType(fFlags, kTRUE)));
2440       // Go back to v_n for ref. keep <<2'>> for diff. flow).
2441       for (Int_t n = 0; n < fMaxMoment-1; n++) {
2442         Double_t vnTwo = 0;
2443         if (type == 'r' || type == 'R')
2444           vnTwo = (vQC2(n) > 0. ? TMath::Sqrt(vQC2(n)) : 0.);
2445         else {
2446           // is really more <<2'>> in this case
2447           vnTwo = vQC2(n);
2448         }
2449         // Fill in corrected v_n
2450         if (vnTwo != 0) static_cast<TH2D*>(cumu.Get(type, n+2, CumuHistos::kNUA))->Fill(eta, cent, vnTwo);
2451       } // End of moment loop
2452     } // End of eta loop
2453   } // End of centrality loop
2454   return;
2455 }
2456 //_____________________________________________________________________
2457 Double_t AliForwardFlowTaskQC::VertexBin::CalculateNUAMatrixElement(Int_t n, Int_t m, Char_t type, Int_t binA, Int_t cBin) const
2458 {
2459   //
2460   //  Calculates the (n,m)-th element in the NUA matrix: 1 if n == m, otherwise:
2461   //               <<cos[(n-m)phi1]>>*<<cos[(n-m)phi2]>> + <<sin[(n-m)phi1]>>*<<sin[(n-m)phi2]>>
2462   //    NUA(n,m) = -----------------------------------------------------------------------------
2463   //                   1 + <<cos[2nphi1]>>*<<cos[2nphi2]>> + <<sin[2nphi1]>>*<<sin[2nphi2]>>
2464   //
2465   //               <<cos[(n+m)phi1]>>*<<cos[(n+m)phi2]>> + <<sin[(n+m)phi1]>>*<<sin[(n+m)phi2]>>
2466   //             + -----------------------------------------------------------------------------
2467   //                   1 + <<cos[2nphi1]>>*<<cos[2nphi2]>> + <<sin[2nphi1]>>*<<sin[2nphi2]>>
2468   //
2469   //  Parameters:
2470   //   n: row
2471   //   m: coumn
2472   //   type: Reference ('r') or differential ('d') or ('a' or 'b')
2473   //   binA: eta bin of phi1
2474   //   cBin: centrality bin
2475   //
2476   //  Return: NUA(n,m)
2477   //
2478   if (n == m) return 1.;
2479   n += 2;
2480   m += 2;
2481
2482   Double_t cosnMmPhi1 = 0, cosnMmPhi2 = 0, sinnMmPhi1 = 0, sinnMmPhi2 = 0;
2483   Double_t cosnPmPhi1 = 0, cosnPmPhi2 = 0, sinnPmPhi1 = 0, sinnPmPhi2 = 0;
2484   Double_t cos2nPhi1 = 0, cos2nPhi2 = 0, sin2nPhi1 = 0, sin2nPhi2 = 0;
2485
2486   // reference flow
2487   if (type == 'r' || type == 'R') {
2488     if ((fFlags & k3Cor)) {
2489       Double_t eta = fCumuNUARef->GetXaxis()->GetBinCenter(binA);
2490       Int_t i = 0;
2491       while (fEtaLims[i] < eta) i++;
2492       Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2493       GetLimits(i-1, aLow, aHigh, bLow, bHigh);
2494       Double_t multA = 0, multB = 0;
2495       for (Int_t a = aLow; a <= aHigh; a++) {
2496         cosnMmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n-m));
2497         sinnMmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n-m));
2498         cosnPmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n+m));
2499         sinnPmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n+m));
2500         cos2nPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(2*n));
2501         sin2nPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(2*n));
2502         multA += fCumuNUARef->GetBinContent(a, cBin, 0);
2503       }
2504        for (Int_t b = bLow; b <= bHigh; b++) {
2505         cosnMmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n-m));
2506         sinnMmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n-m));
2507         cosnPmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n+m));
2508         sinnPmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n+m));
2509         cos2nPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(2*n));
2510         sin2nPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(2*n));
2511         multB += fCumuNUARef->GetBinContent(b, cBin, 0);
2512       }
2513       if (multA == 0 || multB == 0) {
2514         if (fDebug > 0) AliWarning("multA or multB == 0 in matrix elements, aborting NUA");
2515         return 0.;
2516       }
2517       cosnMmPhi1 /= multA;
2518       sinnMmPhi1 /= multA;
2519       cosnPmPhi1 /= multA;
2520       sinnPmPhi1 /= multA;
2521       cos2nPhi1 /= multA;
2522       sin2nPhi1 /= multA;
2523       cosnMmPhi2 /= multB;
2524       sinnMmPhi2 /= multB;
2525       cosnPmPhi2 /= multB;
2526       sinnPmPhi2 /= multB;
2527       cos2nPhi2 /= multB;
2528       sin2nPhi2 /= multB;
2529     } else {
2530       Int_t binB = fCumuNUARef->GetXaxis()->FindBin(-fCumuNUARef->GetXaxis()->GetBinCenter(binA));
2531       cosnMmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2532       sinnMmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2533       cosnPmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2534       sinnPmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2535       cos2nPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2536       sin2nPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2537       cosnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n-m));
2538       sinnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n-m));
2539       cosnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n+m));
2540       sinnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n+m));
2541       cos2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(2*n));
2542       sin2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(2*n));
2543     }
2544   } // differential flow
2545   else if (type == 'd' || type == 'D') {
2546     Int_t binB = fCumuNUARef->GetXaxis()->FindBin(-fCumuNUADiff->GetXaxis()->GetBinCenter(binA));
2547     cosnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2548     sinnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2549     cosnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2550     sinnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2551     cos2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2552     sin2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2553     cosnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n-m));
2554     sinnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n-m));
2555     cosnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n+m));
2556     sinnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n+m));
2557     cos2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(2*n));
2558     sin2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(2*n));
2559   } // 3 correlator part a or b
2560   else if (type == 'a' || type == 'A' || type == 'b' || type == 'B') {
2561     Double_t mult1 = 0, mult2 = 0;
2562     // POIs
2563     cosnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2564     sinnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2565     cosnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2566     sinnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2567     cos2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2568     sin2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2569     mult1 = fCumuNUADiff->GetBinContent(binA, cBin, 0);
2570     // RPs
2571     Double_t eta = fCumuNUADiff->GetXaxis()->GetBinCenter(binA);
2572     Int_t i = 0;
2573     while (fEtaLims[i] < eta) i++;
2574     Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2575     GetLimits(i-1, aLow, aHigh, bLow, bHigh);
2576     Int_t lLow = ((type == 'a' || type == 'A') ? aLow : bLow);
2577     Int_t lHigh = ((type == 'a' || type == 'A') ? aHigh : bHigh);
2578     for (Int_t l = lLow; l <= lHigh; l++) {
2579       cosnMmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(n-m));
2580       sinnMmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(n-m));
2581       cosnPmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(n+m));
2582       sinnPmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(n+m));
2583       cos2nPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(2*n));
2584       sin2nPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(2*n));
2585       mult2 += fCumuNUARef->GetBinContent(l, cBin, 0);
2586     }
2587     if (mult1 == 0 || mult2 == 0) {
2588       if (fDebug > 0) AliWarning("mult1 or mult2 == 0 in matrix elements, aborting NUA");
2589       return 0.;
2590     }
2591     cosnMmPhi1 /= mult1;
2592     sinnMmPhi1 /= mult1;
2593     cosnPmPhi1 /= mult1;
2594     sinnPmPhi1 /= mult1;
2595     cos2nPhi1 /= mult1;
2596     sin2nPhi1 /= mult1;
2597     cosnMmPhi2 /= mult2;
2598     sinnMmPhi2 /= mult2;
2599     cosnPmPhi2 /= mult2;
2600     sinnPmPhi2 /= mult2;
2601     cos2nPhi2 /= mult2;
2602     sin2nPhi2 /= mult2;
2603   }
2604
2605   // Actual calculation
2606   Double_t e = cosnMmPhi1*cosnMmPhi2 + sinnMmPhi1*sinnMmPhi2 + cosnPmPhi1*cosnPmPhi2 + sinnPmPhi1*sinnPmPhi2;
2607   Double_t den = 1 + cos2nPhi1*cos2nPhi2 + sin2nPhi1*sin2nPhi2;
2608   if (den != 0) e /= den;
2609   else return 0.;
2610
2611   return e;
2612 }
2613 //_____________________________________________________________________
2614 void AliForwardFlowTaskQC::VertexBin::AddVertexBins(CumuHistos& cumu, TList* list, UInt_t nNUA) const
2615 {
2616   //
2617   //  Add up vertex bins with flow results
2618   //
2619   //  Parameters:
2620   //   cumu: CumuHistos object with vtxbin results
2621   //   list: Outout list with added results
2622   //   nNUA: # of NUA histograms to loop over
2623   //
2624   TH2D* vtxHist = 0;
2625   TProfile2D* avgProf = 0;
2626   TString name;
2627   Int_t nT = ((fFlags & k3Cor) ? 3 : 2);
2628   Char_t ct = '\0';
2629   for (UInt_t nua = 0; nua <= nNUA; nua++) { // NUA loop
2630     for (Int_t n = 2; n <= fMaxMoment; n++) { // moment loop
2631       for (Int_t t = 0; t < nT; t++) { // type loop (r/d/a/b)
2632         // Find type
2633         switch (t) {
2634           case 0: ct = 'r'; break;
2635           case 1: ct = ((fFlags & k3Cor) ? 'a' : 'd'); break;
2636           case 2: ct = 'b'; break;
2637           default: ct = '\0'; break;
2638         }
2639         vtxHist = static_cast<TH2D*>(cumu.Get(ct,n,nua));
2640         if (!vtxHist) {
2641           AliWarning("VertexBin::AddVertexBins: vtxHist not found!");
2642           continue;
2643         }
2644         name = vtxHist->GetName();
2645         // Strip name of vtx info
2646         Ssiz_t l = name.Last('x')-3;
2647         name.Resize(l);
2648         avgProf = (TProfile2D*)list->FindObject(name.Data());
2649         // if no output profile yet, make one
2650         if (!avgProf) {
2651           avgProf = new TProfile2D(name.Data(), name.Data(),
2652               vtxHist->GetNbinsX(), vtxHist->GetXaxis()->GetXmin(), vtxHist->GetXaxis()->GetXmax(),
2653               vtxHist->GetNbinsY(), vtxHist->GetYaxis()->GetXmin(), vtxHist->GetYaxis()->GetXmax());
2654           if (vtxHist->GetXaxis()->IsVariableBinSize())
2655             avgProf->GetXaxis()->Set(vtxHist->GetNbinsX(), vtxHist->GetXaxis()->GetXbins()->GetArray());
2656           if (vtxHist->GetYaxis()->IsVariableBinSize())
2657             avgProf->GetYaxis()->Set(vtxHist->GetNbinsY(), vtxHist->GetYaxis()->GetXbins()->GetArray());
2658           list->Add(avgProf);
2659         }
2660         // Fill in, cannot be done with Add function.
2661         for (Int_t e = 1; e <= vtxHist->GetNbinsX(); e++) { // eta loop
2662           Double_t eta = vtxHist->GetXaxis()->GetBinCenter(e);
2663           for (Int_t c = 1; c <= vtxHist->GetNbinsY(); c++) { // cent loop
2664             Double_t cent = vtxHist->GetYaxis()->GetBinCenter(c);
2665             Double_t cont = vtxHist->GetBinContent(e, c);
2666             if (cont == 0) continue;
2667             avgProf->Fill(eta, cent, cont);
2668           } // End of cent loop
2669         } //  End of eta loop
2670       } // End of type loop
2671     } // End of moment loop
2672   } // End of nua loop
2673 }
2674 //_____________________________________________________________________
2675 Int_t AliForwardFlowTaskQC::VertexBin::GetBinNumberCos(Int_t n) const
2676 {
2677   //
2678   //  Get the bin number of <<cos(nphi)>>
2679   //
2680   //  Parameters:
2681   //   n: moment
2682   //
2683   //  Return: bin number
2684   //
2685   Int_t bin = 0;
2686   n = TMath::Abs(n);
2687
2688   if (n == 0) bin = fMaxMoment*4-1;
2689   else        bin = n*2-1;
2690
2691   return bin;
2692 }
2693 //_____________________________________________________________________
2694 Int_t AliForwardFlowTaskQC::VertexBin::GetBinNumberSin(Int_t n) const
2695 {
2696   //
2697   //  Get the bin number of <<sin(nphi)>>
2698   //
2699   //  Parameters:
2700   //   n: moment
2701   //
2702   //  Return: bin number
2703   //
2704   Int_t bin = 0;
2705   n = TMath::Abs(n);
2706
2707   if (n == 0) bin = fMaxMoment*4;
2708   else        bin = n*2;
2709
2710   return bin;
2711 }
2712 //_____________________________________________________________________
2713 void AliForwardFlowTaskQC::VertexBin::SetupNUALabels(TAxis* a) const
2714 {
2715   //
2716   //  Setup NUA labels on axis
2717   //
2718   //  Parameters:
2719   //   a: Axis to set up
2720   //
2721   if (a->GetNbins() != GetBinNumberSin()) AliFatal("SetupNUALabels: Wrong number of bins on axis");
2722
2723   Int_t i = 1, j= 1;
2724   while (i <= a->GetNbins()) {
2725     a->SetBinLabel(i++, Form("<<cos(%d#varphi)>>", j));
2726     a->SetBinLabel(i++, Form("<<sin(%d#varphi)>>", j++));
2727   }
2728
2729   return;
2730 }
2731 //_____________________________________________________________________
2732 TH2I* AliForwardFlowTaskQC::VertexBin::MakeQualityHist(const Char_t* name) const
2733 {
2734   //
2735   //  Add a histogram for checking the analysis quality
2736   //
2737   //  Parameters:
2738   //   const char*: name of data type
2739   //
2740   //  Return: histogram for tracking successful calculations
2741   //
2742   Int_t nBins = ((fFlags & kStdQC) ? 8 : 4);
2743   TH2I* quality = new TH2I(name, name, (fMaxMoment-1)*nBins, 1, (fMaxMoment-1)*nBins+1, fCumuNUARef->GetNbinsY(),
2744                            fCumuNUARef->GetYaxis()->GetXmin(), fCumuNUARef->GetYaxis()->GetXmax());
2745   quality->GetYaxis()->Set(fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXbins()->GetArray());
2746   for (Int_t i = 2, j = 1; i <= fMaxMoment; i++) {
2747     quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{2} > 0", i));
2748     quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{2} <= 0", i));
2749     quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{2} > 0", i));
2750     quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{2} <= 0", i));
2751     if ((fFlags & kStdQC)) {
2752       quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{4} < 0", i));
2753       quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{4} >= 0", i));
2754       quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{4} < 0", i));
2755       quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{4} >= 0", i));
2756     }
2757   }
2758
2759   return quality;
2760 }
2761 //_____________________________________________________________________
2762 TH2D* AliForwardFlowTaskQC::VertexBin::MakeOutputHist(Int_t qc, Int_t n, const Char_t* ctype, UInt_t nua) const
2763 {
2764   //
2765   //  Setup a TH2D for the output
2766   //
2767   //  Parameters:
2768   //   qc   # of particle correlations
2769   //   n    flow moment
2770   //   ref  Type: r/d/a/b
2771   //   nua  For nua corrected hists?
2772   //
2773   //  Return: 2D hist for results
2774   //
2775   Bool_t ref = ((ctype[0] == 'R') || (ctype[0] == 'r'));
2776   TAxis* xAxis = (ref ? fCumuNUARef->GetXaxis() : fCumuNUADiff->GetXaxis());
2777   TAxis* yAxis = (ref ? fCumuNUARef->GetYaxis() : fCumuNUADiff->GetYaxis());
2778   TString ntype = "";
2779   switch (nua) {
2780     case CumuHistos::kNoNUA: ntype = "Un"; break;
2781     case CumuHistos::kNUAOld: ntype = "NUAOld"; break;
2782     case CumuHistos::kNUA: ntype = "NUA"; break;
2783     default: break;
2784   }
2785   TH2D* h = new TH2D(Form("%sQC%d_v%d_%s_%sCorr%s_vtx_%d_%d",
2786                         fType.Data(), qc, n, ctype, ntype.Data(),
2787                         GetQCType(fFlags), fVzMin, fVzMax),
2788                       Form("%sQC%d_v%d_%s_%sCorr%s_vtx_%d_%d",
2789                         fType.Data(), qc, n, ctype, ntype.Data(),
2790                         GetQCType(fFlags), fVzMin, fVzMax),
2791                       xAxis->GetNbins(), xAxis->GetXmin(), xAxis->GetXmax(),
2792                       yAxis->GetNbins(), yAxis->GetXmin(), yAxis->GetXmax());
2793   if (xAxis->IsVariableBinSize()) h->GetXaxis()->Set(xAxis->GetNbins(), xAxis->GetXbins()->GetArray());
2794   h->GetYaxis()->Set(yAxis->GetNbins(), yAxis->GetXbins()->GetArray());
2795
2796   return h;
2797 }
2798 //_____________________________________________________________________
2799 void AliForwardFlowTaskQC::PrintFlowSetup() const
2800 {
2801   //
2802   //  Print the setup of the flow task
2803   //
2804   Printf("=======================================================");
2805   Printf("%s::Print", this->IsA()->GetName());
2806   Printf("Forward detector:                 :\t%s", ((fFlowFlags & kFMD) ? "FMD" : "VZERO"));
2807   Printf("Number of bins in vertex axis     :\t%d", fVtxAxis->GetNbins());
2808   Printf("Range of vertex axis              :\t[%3.1f,%3.1f]",
2809                           fVtxAxis->GetXmin(), fVtxAxis->GetXmax());
2810   Printf("Number of bins in centrality axis :\t%d", fCentAxis->GetNbins());
2811   printf("Centrality binning                :\t");
2812   for (Int_t cBin = 1; cBin <= fCentAxis->GetNbins(); cBin++) {
2813     printf("%02d-%02d%% ", Int_t(fCentAxis->GetBinLowEdge(cBin)), Int_t(fCentAxis->GetBinUpEdge(cBin)));
2814     if (cBin == fCentAxis->GetNbins()) printf("\n");
2815     else if (cBin % 4 == 0) printf("\n\t\t\t\t\t");
2816   }
2817   printf("Doing flow analysis for           :\t");
2818   for (Int_t n  = 2; n <= fMaxMoment; n++) printf("v%d ", n);
2819   printf("\n");
2820   TString type = "Standard QC{2} and QC{4} calculations.";
2821   if ((fFlowFlags & kEtaGap)) type = "QC{2} with a rapidity gap.";
2822   if ((fFlowFlags & k3Cor))   type = "QC{2} with 3 correlators.";
2823   if ((fFlowFlags & kTPC) == kTPC)    type.ReplaceAll(".", " with TPC tracks for reference.");
2824   if ((fFlowFlags & kHybrid) == kHybrid) type.ReplaceAll(".", " with hybrid tracks for reference.");
2825   Printf("QC calculation type               :\t%s", type.Data());
2826   Printf("Symmetrize ref. flow wrt. eta = 0 :\t%s", ((fFlowFlags & kSymEta) ? "true" : "false"));
2827   Printf("Apply NUA correction terms        :\t%s", ((fFlowFlags & kNUAcorr) ? "true" : "false"));
2828   Printf("Satellite vertex flag             :\t%s", ((fFlowFlags & kSatVtx) ? "true" : "false"));
2829   Printf("FMD sigma cut:                    :\t%f", fFMDCut);
2830   Printf("SPD sigma cut:                    :\t%f", fSPDCut);
2831   if ((fFlowFlags & kEtaGap) || (fFlowFlags & kTracks))
2832     Printf("Eta gap:                          :\t%f", fEtaGap);
2833   Printf("=======================================================");
2834 }
2835 //_____________________________________________________________________
2836 const Char_t* AliForwardFlowTaskQC::GetQCType(UShort_t flags, Bool_t prependUS)
2837 {
2838   //
2839   //  Get the type of the QC calculations
2840   //  Used for naming of objects in the VertexBin class,
2841   //  important to avoid memory problems when running multiple
2842   //  initializations of the class with different setups
2843   //
2844   //  Parameters:
2845   //   flags: Flow flags for type determination
2846   //   prependUS: Prepend an underscore (_) to the name
2847   //
2848   //  Return: QC calculation type
2849   //
2850   TString type = "";
2851   if      ((flags & kStdQC))  type = "StdQC";
2852   else if ((flags & kEtaGap)) type = "EtaGap";
2853   else if ((flags & k3Cor))   type = "3Cor";
2854   else type = "UNKNOWN";
2855   if (prependUS) type.Prepend("_");
2856   if ((flags & kTPC) == kTPC)    type.Append("TPCTr");
2857   if ((flags & kHybrid) == kHybrid) type.Append("HybTr");
2858
2859   return type.Data();
2860 }
2861 //_____________________________________________________________________
2862 TH1* AliForwardFlowTaskQC::CumuHistos::Get(Char_t t, Int_t n, UInt_t nua) const
2863 {
2864   //
2865   //  Get histogram/profile for type t and moment n
2866   //
2867   //  Parameters:
2868   //   t: type = 'r'/'d'
2869   //   n: moment
2870   //   nua: NUA type
2871   //
2872   n = GetPos(n, nua);
2873   if (n < 0) AliFatal(Form("CumuHistos out of range: (%c,%d)", t, n));
2874
2875   TH1* h = 0;
2876   Int_t pos = -1;
2877   if      (t == 'r' || t == 'R') pos = n;
2878   else if (t == 'd' || t == 'D') pos = n;
2879   else if (t == 'a' || t == 'A') pos = 2*n;
2880   else if (t == 'b' || t == 'B') pos = 2*n+1;
2881   else AliFatal(Form("CumuHistos wrong type: %c", t));
2882
2883   if (t == 'r' || t == 'R') {
2884     if (pos < fRefHists->GetEntries()) {
2885       h = (TH1*)fRefHists->At(pos);
2886     }
2887   } else if (pos < fDiffHists->GetEntries()) {
2888     h = (TH1*)fDiffHists->At(pos);
2889   }
2890   if (!h) AliFatal(Form("No hist found in list %c at pos %d", t, pos));
2891
2892   return h;
2893 }
2894 //_____________________________________________________________________
2895 void AliForwardFlowTaskQC::CumuHistos::ConnectList(TString name, TList* l)
2896 {
2897   //
2898   //  Connect an output list with this object
2899   //
2900   //  Parameters:
2901   //   name: base name
2902   //   l: list to keep outputs in
2903   //
2904   TString ref = name;
2905   ref.ReplaceAll("Cumu","CumuRef");
2906   fRefHists = (TList*)l->FindObject(ref.Data());
2907   if (!fRefHists) {
2908     fRefHists = new TList();
2909     fRefHists->SetName(ref.Data());
2910     l->Add(fRefHists);
2911   } else {
2912     // Check that the list correspond to fMaxMoments
2913     if (fRefHists->GetEntries() != (fMaxMoment-1.)*(fNUA+1))
2914       AliError("CumuHistos::ConnectList Wrong number of hists in ref list,"
2915                "you are doing something wrong!");
2916   }
2917   TString diff = name;
2918   diff.ReplaceAll("Cumu","CumuDiff");
2919   fDiffHists = (TList*)l->FindObject(diff.Data());
2920   if (!fDiffHists) {
2921     fDiffHists = new TList();
2922     fDiffHists->SetName(diff.Data());
2923     l->Add(fDiffHists);
2924   } else {
2925     // Check that the list correspond to fMaxMoment
2926     if ((fDiffHists->GetEntries() != (fMaxMoment-1.)*(fNUA+1)) &&
2927         (fDiffHists->GetEntries() != 2*(fMaxMoment-1.)*(fNUA+1)))
2928       AliError(Form("CumuHistos::ConnectList Wrong number of hists in diff list,"
2929                "you are doing something wrong! (%s)",name.Data()));
2930   }
2931
2932 }
2933 //_____________________________________________________________________
2934 void AliForwardFlowTaskQC::CumuHistos::Add(TH1* h) const
2935 {
2936   //
2937   //  Add a histogram to this objects lists
2938   //
2939   //  Parameters:
2940   //   h: histogram/profile to add
2941   //
2942   TString name = h->GetName();
2943   if (name.Contains("Ref")) {
2944     if (fRefHists) fRefHists->Add(h);
2945     else AliFatal("CumuHistos::Add() - fRefHists does not exist");
2946     // Check that the list correspond to fMaxMoments
2947     if (fRefHists->GetEntries() > (fNUA+1)*(fMaxMoment-1.))
2948       AliError("CumuHistos::Add wrong number of hists in ref list, "
2949                "you are doing something wrong!");
2950   }
2951   else if (name.Contains("Diff")) {
2952     if (fDiffHists) fDiffHists->Add(h);
2953     else AliFatal("CumuHistos::Add() - fDiffHists does not exist");
2954     // Check that the list correspond to fMaxMoment
2955     if (fDiffHists->GetEntries() > 2*(fNUA+1)*(fMaxMoment-1.))
2956       AliError("CumuHistos::Add wrong number of hists in diff list, "
2957          "you are doing something wrong!");
2958   }
2959   return;
2960 }
2961 //_____________________________________________________________________
2962 Int_t AliForwardFlowTaskQC::CumuHistos::GetPos(Int_t n, UInt_t nua) const
2963 {
2964   //
2965   //  Get position in list of moment n objects
2966   //  To take care of different numbering schemes
2967   //
2968   //  Parameters:
2969   //   n: moment
2970   //   nua: # of nua corrections applied
2971   //
2972   //  Return: position #
2973   //
2974   if (n > fMaxMoment) return -1;
2975   else return (n-2)+nua*(fMaxMoment-1);
2976 }
2977 //_____________________________________________________________________
2978 //
2979 //
2980 // EOF