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, bit = 0;
1326       if ((fFlags & kTPC) == kTPC)    pTMin = 0.2, pTMax = 5., etaMin = -0.8, etaMax = 0.8, minNCl = 70, bit = 128;
1327       if ((fFlags & kHybrid) == kHybrid) pTMin = 0.2, pTMax = 5., etaMin = -0.8, etaMax = 0.8, minNCl = 70, bit = 272;
1328
1329       AliAODTrack* aodTr = (AliAODTrack*)tr;
1330       if (aodTr->GetID() > -1) continue;
1331       if (!aodTr->TestFilterBit(bit) || aodTr->Pt() > pTMax || aodTr->Pt() < pTMin ||
1332         aodTr->Eta() > etaMax || aodTr->Eta() < etaMin || aodTr->GetTPCNcls() < minNCl) continue;
1333     }
1334
1335 //    if (tr->GetTPCsignal() < tpcdEdx) continue;
1336     // Track accepted
1337     Double_t eta = tr->Eta();
1338     if (((fFlags & kSPD) || (fFlags & kEtaGap)) && TMath::Abs(eta) < fEtaGap) continue;
1339     Double_t phi = tr->Phi();
1340 //    Double_t pT = tr->Pt();
1341 //    AliAODMCHeader* head = static_cast<AliAODMCHeader*>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
1342 //    Double_t rp = head->GetReactionPlaneAngle();
1343 //    Double_t b = head->GetImpactParameter();
1344     Double_t weight = 1.;//AliForwardFlowWeights::Instance().CalcWeight(eta, pT, phi, 0, rp, b);
1345
1346     if ((mode & kFillRef)) {
1347       fCumuRef->Fill(eta, 0., weight);// mult goes in underflowbin - no visual, but not needed?
1348       fdNdedpRefAcc->Fill(eta, phi, weight);
1349     }
1350     if ((mode & kFillDiff)) {
1351       fCumuDiff->Fill(eta, 0., weight);
1352       fdNdedpDiffAcc->Fill(eta, phi, weight);
1353     }
1354     for (Int_t n = 1; n <= 2*fMaxMoment; n++) {
1355       Double_t cosBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberCos(n));
1356       Double_t sinBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberSin(n));
1357       Double_t cosnPhi = weight*TMath::Cos(n*phi);
1358       Double_t sinnPhi = weight*TMath::Sin(n*phi);
1359       // fill ref
1360       if ((mode & kFillRef)) {
1361         fCumuRef->Fill(eta, cosBin, cosnPhi);
1362         fCumuRef->Fill(eta, sinBin, sinnPhi);
1363       }
1364       // fill diff
1365       if ((mode & kFillDiff)) {
1366         fCumuDiff->Fill(eta, cosBin, cosnPhi);
1367         fCumuDiff->Fill(eta, sinBin, sinnPhi);
1368       }
1369     } // End of NUA loop
1370   } // End of track loop
1371   return kTRUE;
1372 }
1373 //_____________________________________________________________________
1374 void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate(Double_t cent)
1375 {
1376   //
1377   //  Calculate the Q cumulant up to order fMaxMoment
1378   //
1379   //  Parameters:
1380   //   cent: centrality of event
1381   //
1382   if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1383
1384   // Fill out NUA hists
1385   for (Int_t etaBin = 1; etaBin <= fCumuRef->GetNbinsX(); etaBin++) {
1386     Double_t eta = fCumuRef->GetXaxis()->GetBinCenter(etaBin);
1387     if (fCumuRef->GetBinContent(etaBin, 0) == 0) continue;
1388     if ((fFlags & kTracks) && (fFlags && kSPD) && !(fFlags & kEtaGap)) eta = -eta;
1389     for (Int_t qBin = 0; qBin <= fCumuRef->GetNbinsY(); qBin++) {
1390       fCumuNUARef->Fill(eta, cent, Double_t(qBin), fCumuRef->GetBinContent(etaBin, qBin));
1391     }
1392   }
1393   for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1394     Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1395     if (fCumuDiff->GetBinContent(etaBin, 0) == 0) continue;
1396     for (Int_t qBin = 0; qBin <= fCumuDiff->GetNbinsY(); qBin++) {
1397       fCumuNUADiff->Fill(eta, cent, Double_t(qBin), fCumuDiff->GetBinContent(etaBin, qBin));
1398     }
1399   }
1400
1401   // We create the objects needed for the analysis
1402   TH3D* cumuRef = 0;
1403   TH3D* cumuDiff = 0;
1404   // For each n we loop over the hists
1405   for (Int_t n = 2; n <= fMaxMoment; n++) {
1406     cumuRef  = (TH3D*)fCumuHists.Get('r',n);
1407     cumuDiff = (TH3D*)fCumuHists.Get('d',n);
1408     Int_t prevRefEtaBin = 0;
1409
1410     // Per mom. quantities
1411     Double_t dQnReA = 0, dQnImA = 0, multA = 0;
1412     Double_t dQnReB = 0, dQnImB = 0, multB = 0;
1413     Double_t dQ2nReA = 0, dQ2nImA = 0;
1414     Double_t two = 0, w2 = 0, four = 0, w4 = 0;
1415     Double_t cosPhi1Phi2 = 0, cosPhi1Phi2Phi3m = 0;
1416     Double_t sinPhi1Phi2 = 0, sinPhi1Phi2Phi3m = 0;
1417     for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1418       Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1419       Double_t refEta = eta;
1420       if ((fFlags & kTracks) && (fFlags && kSPD) && !(fFlags & kEtaGap)) refEta = -eta;
1421       Int_t refEtaBinA = fCumuRef->GetXaxis()->FindBin(refEta);
1422       if ((fFlags & kEtaGap)) refEta = -eta;
1423       Int_t refEtaBinB = fCumuRef->GetXaxis()->FindBin(refEta);
1424       if (refEtaBinA != prevRefEtaBin) {
1425         prevRefEtaBin = refEtaBinA;
1426         // Reference flow
1427         multA  = fCumuRef->GetBinContent(refEtaBinA, 0);
1428         dQnReA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberCos(n));
1429         dQnImA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberSin(n));
1430         dQ2nReA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberCos(n*2));
1431         dQ2nImA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberSin(n*2));
1432
1433         multB  = fCumuRef->GetBinContent(refEtaBinB, 0);
1434         dQnReB = fCumuRef->GetBinContent(refEtaBinB, GetBinNumberCos(n));
1435         dQnImB = fCumuRef->GetBinContent(refEtaBinB, GetBinNumberSin(n));
1436
1437         if (multA <= 3 || multB <= 3) return;
1438         // The reference flow is calculated
1439         // 2-particle
1440         if ((fFlags & kStdQC)) {
1441           w2 = multA * (multA - 1.);
1442           two = dQnReA*dQnReA + dQnImA*dQnImA - multA;
1443         } else {
1444           w2 = multA * multB;
1445           two = dQnReA*dQnReB + dQnImA*dQnImB;
1446         }
1447         cumuRef->Fill(eta, cent, kW2Two, two);
1448         cumuRef->Fill(eta, cent, kW2, w2);
1449
1450         // The reference flow is calculated
1451         // 4-particle
1452         if ((fFlags & kStdQC)) {
1453           w4 = multA * (multA - 1.) * (multA - 2.) * (multA - 3.);
1454
1455           four = 2.*multA*(multA-3.) + TMath::Power((TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.)),2.)
1456                    -4.*(multA-2.)*(TMath::Power(dQnReA,2.) + TMath::Power(dQnImA,2.))
1457                    -2.*(TMath::Power(dQnReA,2.)*dQ2nReA+2.*dQnReA*dQnImA*dQ2nImA-TMath::Power(dQnImA,2.)*dQ2nReA)
1458                    +(TMath::Power(dQ2nReA,2.)+TMath::Power(dQ2nImA,2.));
1459
1460           cumuRef->Fill(eta, cent, kW4Four, four);
1461           cumuRef->Fill(eta, cent, kW4, w4);
1462
1463           // NUA
1464           cosPhi1Phi2 = dQnReA*dQnReA - dQnImA*dQnImA - dQ2nReA;
1465           sinPhi1Phi2 = 2.*dQnReA*dQnImA - dQ2nImA;
1466
1467           cosPhi1Phi2Phi3m = dQnReA*(TMath::Power(dQnReA,2)+TMath::Power(dQnImA,2))
1468                               -dQnReA*dQ2nReA-dQnImA*dQ2nImA-2.*(multA-1)*dQnReA;
1469
1470           sinPhi1Phi2Phi3m = -dQnImA*(TMath::Power(dQnReA,2)+TMath::Power(dQnImA,2))
1471                               +dQnReA*dQ2nImA-dQnImA*dQ2nReA+2.*(multA-1)*dQnImA;
1472
1473           cumuRef->Fill(eta, cent, kCosphi1phi2, cosPhi1Phi2);
1474           cumuRef->Fill(eta, cent, kSinphi1phi2, sinPhi1Phi2);
1475           cumuRef->Fill(eta, cent, kCosphi1phi2phi3m, cosPhi1Phi2Phi3m);
1476           cumuRef->Fill(eta, cent, kSinphi1phi2phi3m, sinPhi1Phi2Phi3m);
1477           cumuRef->Fill(eta, cent, k3pWeight, multA*(multA-1.)*(multA-2.));
1478         } // End of QC{4}
1479       } // End of reference flow
1480       // For each etaBin bin the necessary values for differential flow is calculated
1481       Double_t mp = fCumuDiff->GetBinContent(etaBin, 0);
1482       Double_t pnRe = fCumuDiff->GetBinContent(etaBin, GetBinNumberCos(n));
1483       Double_t pnIm = fCumuDiff->GetBinContent(etaBin, GetBinNumberSin(n));
1484       Double_t p2nRe = fCumuDiff->GetBinContent(etaBin, GetBinNumberCos(n*2));
1485       Double_t p2nIm = fCumuDiff->GetBinContent(etaBin, GetBinNumberSin(n*2));
1486       if (mp == 0) continue;
1487       Double_t mq = 0;
1488       Double_t qnRe = 0;
1489       Double_t qnIm = 0;
1490       Double_t q2nRe = 0;
1491       Double_t q2nIm = 0;
1492
1493       // Differential flow calculations for each eta bin is done:
1494       // 2-particle differential flow
1495       if ((fFlags & kStdQC) && (!(fFlags & kTracks) || ((fFlags & kTracks) && (fFlags & kMC) && !(fFlags & kSPD) && TMath::Abs(eta) < 0.75))) {
1496         mq = mp;
1497         qnRe = pnRe;
1498         qnIm = pnIm;
1499         q2nRe = p2nRe;
1500         q2nIm = p2nIm;
1501       }
1502
1503       Double_t w2p = mp * multB - mq;
1504       Double_t twoPrime = pnRe*dQnReB + pnIm*dQnImB - mq;
1505
1506       cumuDiff->Fill(eta, cent, kW2Two, twoPrime);
1507       cumuDiff->Fill(eta, cent, kW2, w2p);
1508
1509       if ((fFlags & kEtaGap)) continue;
1510       // Differential flow calculations for each eta bin bin is done:
1511       // 4-particle differential flow
1512       Double_t w4p = (mp * multA - 3.*mq)*(multA - 1.)*(multA - 2.);
1513
1514       Double_t fourPrime = (TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.))*(pnRe*dQnReA+pnIm*dQnImA)
1515                           - q2nRe*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))
1516                           - 2.*q2nIm*dQnReA*dQnImA
1517                           - pnRe*(dQnReA*dQ2nReA+dQnImA*dQ2nImA)
1518                           + pnIm*(dQnImA*dQ2nReA-dQnReA*dQ2nImA)
1519                           - 2.*multA*(pnRe*dQnReA+pnIm*dQnImA)
1520                           - 2.*(TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.))*mq
1521                           + 6.*(qnRe*dQnReA+qnIm*dQnImA)
1522                           + 1.*(q2nRe*dQ2nReA+q2nIm*dQ2nImA)
1523                           + 2.*(pnRe*dQnReA+pnIm*dQnImA)
1524                           + 2.*mq*multA
1525                           - 6.*mq;
1526
1527       cumuDiff->Fill(eta, cent, kW4Four, fourPrime);
1528       cumuDiff->Fill(eta, cent, kW4, w4p);
1529
1530       // NUA
1531       Double_t cosPsi1Phi2 = pnRe*dQnReA - pnIm*dQnImA - q2nRe;
1532       Double_t sinPsi1Phi2 = pnRe*dQnImA + pnIm*dQnReA - q2nIm;
1533
1534       Double_t cosPsi1Phi2Phi3p = pnRe*(TMath::Power(dQnImA,2.)+TMath::Power(dQnReA,2.)-multA)
1535                             - 1.*(q2nRe*dQnReA+q2nIm*dQnImA)
1536                             - mq*dQnReA+2.*qnRe;
1537
1538       Double_t sinPsi1Phi2Phi3p = pnIm*(TMath::Power(dQnImA,2.)+TMath::Power(dQnReA,2.)-multA)
1539                             - 1.*(q2nIm*dQnReA-q2nRe*dQnImA)
1540                             - mq*dQnImA+2.*qnIm;
1541
1542       Double_t cosPsi1Phi2Phi3m = pnRe*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))+2.*pnIm*dQnReA*dQnImA
1543                             - 1.*(pnRe*dQ2nReA+pnIm*dQ2nImA)
1544                             - 2.*mq*dQnReA+2.*qnRe;
1545
1546       Double_t sinPsi1Phi2Phi3m = pnIm*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))-2.*pnRe*dQnReA*dQnImA
1547                             - 1.*(pnIm*dQ2nReA-pnRe*dQ2nImA)
1548                             + 2.*mq*dQnImA-2.*qnIm;
1549
1550       cumuDiff->Fill(eta, cent, kCosphi1phi2, cosPsi1Phi2);
1551       cumuDiff->Fill(eta, cent, kSinphi1phi2, sinPsi1Phi2);
1552       cumuDiff->Fill(eta, cent, kCosphi1phi2phi3m, cosPsi1Phi2Phi3m);
1553       cumuDiff->Fill(eta, cent, kSinphi1phi2phi3m, sinPsi1Phi2Phi3m);
1554       cumuDiff->Fill(eta, cent, k3pWeight, (mp*multA-2.*mq)*(multA-1.));
1555       cumuDiff->Fill(eta, cent, kCosphi1phi2phi3p, cosPsi1Phi2Phi3p);
1556       cumuDiff->Fill(eta, cent, kSinphi1phi2phi3p, sinPsi1Phi2Phi3p);
1557     } // End of eta loop
1558     // Event count
1559     cumuRef->Fill(-7., cent, -0.5, 1.);
1560   } // End of moment loop
1561   return;
1562 }
1563 //_____________________________________________________________________
1564 void AliForwardFlowTaskQC::VertexBin::GetLimits(Int_t bin, Int_t& aLow, Int_t& aHigh,
1565                                                 Int_t& bLow, Int_t& bHigh) const
1566 {
1567   //
1568   //  Get the limits for the 3 correlator method
1569   //
1570   //  Parameters:
1571   //   bin  : reference bin #
1572   //   aLow : Lowest bin to be included in v_A calculations
1573   //   aHigh: Highest bin to be included in v_A calculations
1574   //   bLow : Lowest bin to be included in v_B calculations
1575   //   bHigh: Highest bin to be included in v_B calculations
1576   //
1577   if ((fFlags & kFMD)) {
1578     switch(bin) {
1579       case 0:
1580         aLow = 14; aHigh = 15;
1581         bLow = 20; bHigh = 22;
1582         break;
1583       case 1:
1584         aLow = 16; aHigh = 16;
1585         bLow = 21; bHigh = 22;
1586         break;
1587       case 2:
1588         aLow =  6; aHigh =  7;
1589         bLow = 21; bHigh = 22;
1590         break;
1591       case 3:
1592         aLow =  6; aHigh =  7;
1593         bLow = 12; bHigh = 12;
1594         break;
1595       case 4:
1596         aLow =  6; aHigh =  8;
1597         bLow = 13; bHigh = 14;
1598         break;
1599       default:
1600         AliFatal(Form("No limits for this eta region! (%d)", bin));
1601     }
1602   }
1603   else if ((fFlags & kVZERO)) {
1604     switch(bin) {
1605       case 0:
1606         aLow =  6; aHigh = 13;
1607         bLow = 17; bHigh = 18;
1608         break;
1609       case 1:
1610         aLow =  6; aHigh =  9;
1611         bLow = 17; bHigh = 18;
1612         break;
1613       case 2:
1614         aLow =  2; aHigh =  3;
1615         bLow = 17; bHigh = 18;
1616         break;
1617       case 3:
1618         aLow =  2; aHigh =  3;
1619         bLow =  6; bHigh =  9;
1620         break;
1621       case 4:
1622         aLow =  2; aHigh =  3;
1623         bLow =  6; bHigh = 13;
1624         break;
1625       default:
1626         AliFatal(Form("No limits for this eta region! (%d)", bin));
1627     }
1628   }
1629   // Try to catch cases where fEtaLimits and these values do not correspond to each other
1630   if (aHigh > fCumuNUARef->GetNbinsX() || bHigh > fCumuNUARef->GetNbinsX())
1631     AliFatal(Form("Limits outside vtx range! (%d) - aHigh = %d, bHigh = %d, Nbins = %d",
1632       bin, aHigh, bHigh, fCumuNUARef->GetNbinsX()));
1633 }
1634 //_____________________________________________________________________
1635 void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate3Cor(Double_t cent)
1636 {
1637   //
1638   //  Calculate the Q cumulant up to order fMaxMoment
1639   //
1640   //  Parameters:
1641   //   cent: centrality of event
1642   //
1643   if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1644
1645   // Fill out NUA hists
1646   for (Int_t etaBin = 1; etaBin <= fCumuRef->GetNbinsX(); etaBin++) {
1647     Double_t eta = fCumuRef->GetXaxis()->GetBinCenter(etaBin);
1648     if (fCumuRef->GetBinContent(etaBin, 0) == 0) continue;
1649     for (Int_t qBin = 0; qBin <= fCumuRef->GetNbinsY(); qBin++) {
1650       fCumuNUARef->Fill(eta, cent, Double_t(qBin), fCumuRef->GetBinContent(etaBin, qBin));
1651     }
1652   }
1653   for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1654     Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1655     if (fCumuDiff->GetBinContent(etaBin, 0) == 0) continue;
1656     for (Int_t qBin = 0; qBin <= fCumuDiff->GetNbinsY(); qBin++) {
1657       fCumuNUADiff->Fill(eta, cent, Double_t(qBin), fCumuDiff->GetBinContent(etaBin, qBin));
1658     }
1659   }
1660
1661   // We create the objects needed for the analysis
1662   TH3D* cumuRef = 0;
1663   TH3D* cumuDiff = 0;
1664   // For each n we loop over the hists
1665   for (Int_t n = 2; n <= fMaxMoment; n++) {
1666     cumuRef  = (TH3D*)fCumuHists.Get('r',n);
1667     cumuDiff = (TH3D*)fCumuHists.Get('d',n);
1668
1669     // Per mom. quantities
1670     Int_t prevLim = 0;
1671     Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
1672     Double_t dQnReA = 0, dQnImA = 0, multA = 0;
1673     Double_t dQnReB = 0, dQnImB = 0, multB = 0;
1674     Double_t two = 0, w2 = 0;
1675     for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1676       Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1677       if (fEtaLims[prevLim] < eta) {
1678         GetLimits(prevLim, aLow, aHigh, bLow, bHigh);
1679         prevLim++;
1680         multA = 0; dQnReA = 0; dQnImA = 0;
1681         multB = 0; dQnReB = 0; dQnImB = 0;
1682         // Reference flow
1683         for (Int_t a = aLow; a <= aHigh; a++) {
1684           multA  += fCumuRef->GetBinContent(a, 0);
1685           dQnReA += fCumuRef->GetBinContent(a, GetBinNumberCos(n));
1686           dQnImA += fCumuRef->GetBinContent(a, GetBinNumberSin(n));
1687         }
1688         for (Int_t b = bLow; b <= bHigh; b++) {
1689           multB  += fCumuRef->GetBinContent(b, 0);
1690           dQnReB += fCumuRef->GetBinContent(b, GetBinNumberCos(n));
1691           dQnImB += fCumuRef->GetBinContent(b, GetBinNumberSin(n));
1692         }
1693         // The reference flow is calculated
1694         // 2-particle
1695         w2 = multA * multB;
1696         two = dQnReA*dQnReB + dQnImA*dQnImB;
1697       } // End of reference flow
1698       cumuRef->Fill(eta, cent, kW2Two, two);
1699       cumuRef->Fill(eta, cent, kW2, w2);
1700
1701       // For each etaBin bin the necessary values for differential flow is calculated
1702       Double_t mp = fCumuDiff->GetBinContent(etaBin, 0);
1703       Double_t pnRe = fCumuDiff->GetBinContent(etaBin, GetBinNumberCos(n));
1704       Double_t pnIm = fCumuDiff->GetBinContent(etaBin, GetBinNumberSin(n));
1705       if (mp == 0) continue;
1706
1707       // Differential flow calculations for each eta bin is done:
1708       // 2-particle differential flow
1709       Double_t w2pA = mp * multA;
1710       Double_t twoPrimeA = pnRe*dQnReA + pnIm*dQnImA;
1711       cumuDiff->Fill(eta, cent, kW2Two, twoPrimeA);
1712       cumuDiff->Fill(eta, cent, kW2, w2pA);
1713
1714       Double_t w2pB = mp * multB;
1715       Double_t twoPrimeB = pnRe*dQnReB + pnIm*dQnImB;
1716       cumuDiff->Fill(eta, cent, kW4Four, twoPrimeB);
1717       cumuDiff->Fill(eta, cent, kW4, w2pB);
1718      } // End of eta loop
1719     // Event count
1720     cumuRef->Fill(-7., cent, -0.5, 1.);
1721   } // End of moment loop
1722   return;
1723
1724 }
1725 //_____________________________________________________________________
1726 void AliForwardFlowTaskQC::VertexBin::CumulantsTerminate(TList* inlist, TList* outlist)
1727 {
1728   //
1729   //  Finalizes the Q cumulant calculations
1730   //
1731   //  Parameters:
1732   //   inlist: input sumlist
1733   //   outlist: output result list
1734   //
1735
1736   // Re-find cumulants hist if Terminate is called separately
1737   if (!fCumuHists.IsConnected()) {
1738     TList* list = (TList*)inlist->FindObject(Form("%svertex_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1739     fCumuHists.ConnectList(Form("%sCumu_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)), list);
1740
1741     if (!fCumuNUARef)
1742       fCumuNUARef = (TH3D*)list->FindObject(Form("%s_vertex_%d_%d%s_cumuNUARef", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1743     if (!fCumuNUADiff)
1744       fCumuNUADiff = (TH3D*)list->FindObject(Form("%s_vertex_%d_%d%s_cumuNUADiff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1745   }
1746   // Clone to avoid normalization problems when redoing terminate locally
1747   fCumuNUARef = (TH3D*)fCumuNUARef->Clone(Form("%s_vertex_%d_%d%s_cumuNUARefNorm", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1748   fCumuNUADiff = (TH3D*)fCumuNUADiff->Clone(Form("%s_vertex_%d_%d%s_cumuNUADiffNorm", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1749
1750   // Diagnostics histograms
1751   TH2I* quality = (TH2I*)outlist->FindObject(Form("hQCQuality%s%s", fType.Data(), GetQCType(fFlags)));
1752   if (!quality) {
1753     quality = MakeQualityHist(Form("hQCQuality%s%s", fType.Data(), GetQCType(fFlags)));
1754     outlist->Add(quality);
1755   }
1756   TH1D* cent = (TH1D*)outlist->FindObject(Form("%s%s_cent", fType.Data(), GetQCType(fFlags)));
1757   if (!cent) {
1758     cent = new TH1D(Form("%s%s_cent", fType.Data(), GetQCType(fFlags)),
1759                     Form("%s%s_cent", fType.Data(), GetQCType(fFlags)),
1760                 fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXmin(), fCumuNUARef->GetYaxis()->GetXmax());
1761     cent->GetXaxis()->Set(fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXbins()->GetArray());
1762     outlist->Add(cent);
1763   }
1764   TH2D* dNdetaRef = (TH2D*)outlist->FindObject(Form("%s%s_dNdetaRef", fType.Data(), GetQCType(fFlags)));
1765   if (!dNdetaRef) {
1766     dNdetaRef = new TH2D(Form("%s%s_dNdetaRef", fType.Data(), GetQCType(fFlags)),
1767                                Form("%s%s_dNdetaRef", fType.Data(), GetQCType(fFlags)),
1768                 fCumuNUARef->GetNbinsX(), fCumuNUARef->GetXaxis()->GetXmin(), fCumuNUARef->GetXaxis()->GetXmax(),
1769                 fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXmin(), fCumuNUARef->GetYaxis()->GetXmax());
1770     dNdetaRef->GetYaxis()->Set(fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXbins()->GetArray());
1771     dNdetaRef->Sumw2();
1772     outlist->Add(dNdetaRef);
1773   }
1774   TH2D* dNdetaDiff = (TH2D*)outlist->FindObject(Form("%s%s_dNdetaDiff", fType.Data(), GetQCType(fFlags)));
1775   if (!dNdetaDiff) {
1776     dNdetaDiff = new TH2D(Form("%s%s_dNdetaDiff", fType.Data(), GetQCType(fFlags)),
1777                                 Form("%s%s_dNdetaDiff", fType.Data(), GetQCType(fFlags)),
1778                 fCumuNUADiff->GetNbinsX(), fCumuNUADiff->GetXaxis()->GetXmin(), fCumuNUADiff->GetXaxis()->GetXmax(),
1779                 fCumuNUADiff->GetNbinsY(), fCumuNUADiff->GetYaxis()->GetXmin(), fCumuNUADiff->GetYaxis()->GetXmax());
1780     dNdetaDiff->GetYaxis()->Set(fCumuNUADiff->GetNbinsY(), fCumuNUADiff->GetYaxis()->GetXbins()->GetArray());
1781     dNdetaDiff->Sumw2();
1782     outlist->Add(dNdetaDiff);
1783   }
1784
1785   // Setting up outputs
1786   // Create output lists and diagnostics
1787   TList* vtxList = (TList*)outlist->FindObject("vtxList");
1788   if (!vtxList) {
1789     vtxList = new TList();
1790     vtxList->SetName("vtxList");
1791     outlist->Add(vtxList);
1792   }
1793   vtxList->Add(fCumuNUARef);
1794   vtxList->Add(fCumuNUADiff);
1795
1796   // Setup output profiles
1797   CumuHistos cumu2(fMaxMoment, ((fFlags & kNUAcorr) ? 2 : 0));
1798   CumuHistos cumu4(fMaxMoment, ((fFlags & kNUAcorr) ? 1 : 0));
1799
1800   cumu2.ConnectList(Form("%sQC2_Cumu%s_vtx_%d_%d", fType.Data(), GetQCType(fFlags), fVzMin, fVzMax), vtxList);
1801   if ((fFlags & kStdQC))
1802     cumu4.ConnectList(Form("%sQC4_Cumu%s_vtx_%d_%d", fType.Data(), GetQCType(fFlags), fVzMin, fVzMax), vtxList);
1803
1804   for (Int_t n = 2; n <= fMaxMoment; n++) {
1805     // 2-particle
1806     cumu2.Add(MakeOutputHist(2, n, "Ref", CumuHistos::kNoNUA));
1807     if ((fFlags & k3Cor)){
1808       cumu2.Add(MakeOutputHist(2, n, "DiffA", CumuHistos::kNoNUA));
1809       cumu2.Add(MakeOutputHist(2, n, "DiffB", CumuHistos::kNoNUA));
1810     } else {
1811       cumu2.Add(MakeOutputHist(2, n, "Diff", CumuHistos::kNoNUA));
1812     }
1813     // 4-particle
1814     if ((fFlags & kStdQC)) {
1815       cumu4.Add(MakeOutputHist(4, n, "Ref", CumuHistos::kNoNUA));
1816       cumu4.Add(MakeOutputHist(4, n, "Diff", CumuHistos::kNoNUA));
1817     }
1818   } // End of v_n result loop
1819   // NUA corrected
1820   if ((fFlags & kNUAcorr)) {
1821     for (Int_t n = 2; n <= fMaxMoment; n++) {
1822       // 2-particle
1823       cumu2.Add(MakeOutputHist(2, n, "Ref", CumuHistos::kNUAOld));
1824       if ((fFlags & k3Cor)) {
1825         cumu2.Add(MakeOutputHist(2, n, "DiffA", CumuHistos::kNUAOld));
1826         cumu2.Add(MakeOutputHist(2, n, "DiffB", CumuHistos::kNUAOld));
1827       } else {
1828         cumu2.Add(MakeOutputHist(2, n, "Diff", CumuHistos::kNUAOld));
1829       }
1830       // 4-particle
1831       if ((fFlags & kStdQC)) {
1832         cumu4.Add(MakeOutputHist(4, n, "Ref", CumuHistos::kNUAOld));
1833         cumu4.Add(MakeOutputHist(4, n, "Diff", CumuHistos::kNUAOld));
1834       }
1835     }
1836     for (Int_t n = 2; n <= fMaxMoment; n++) {
1837       // 2-particle
1838       cumu2.Add(MakeOutputHist(2, n, "Ref", CumuHistos::kNUA));
1839       if ((fFlags & k3Cor)) {
1840         cumu2.Add(MakeOutputHist(2, n, "DiffA", CumuHistos::kNUA));
1841         cumu2.Add(MakeOutputHist(2, n, "DiffB", CumuHistos::kNUA));
1842       } else {
1843         cumu2.Add(MakeOutputHist(2, n, "Diff", CumuHistos::kNUA));
1844       }
1845     }
1846   }
1847
1848   // Calculating the cumulants
1849   if ((fFlags & k3Cor)) {
1850     Calculate3CorFlow(cumu2, quality, cent, dNdetaRef, dNdetaDiff);
1851   } else {
1852     CalculateReferenceFlow(cumu2, cumu4, quality, cent, dNdetaRef);
1853     CalculateDifferentialFlow(cumu2, cumu4, quality, dNdetaDiff);
1854   }
1855   if ((fFlags & kNUAcorr)) {
1856     SolveCoupledFlowEquations(cumu2, 'r');
1857     if ((fFlags & k3Cor)) {
1858       SolveCoupledFlowEquations(cumu2, 'a');
1859       SolveCoupledFlowEquations(cumu2, 'b');
1860     } else {
1861       SolveCoupledFlowEquations(cumu2, 'd');
1862     }
1863   }
1864
1865   // Add to output for immediate viewing - individual vtx bins are used for final results
1866   AddVertexBins(cumu2, outlist, ((fFlags & kNUAcorr) ? 2 : 0));
1867   if ((fFlags & kStdQC)) AddVertexBins(cumu4, outlist, ((fFlags & kNUAcorr) ? 1 : 0));
1868
1869   // Setup NUA diagnoastics histograms
1870   TList* nualist = (TList*)outlist->FindObject("NUATerms");
1871   if (!nualist) {
1872     nualist = new TList();
1873     nualist->SetName("NUATerms");
1874     outlist->Add(nualist);
1875   }
1876   // Reference
1877   TH2D* nuaRef = (TH2D*)nualist->FindObject(Form("%sReferenceNUA%s", fType.Data(), GetQCType(fFlags)));
1878   TH2D* temp = 0;
1879   if (!nuaRef) {
1880     nuaRef = (TH2D*)fCumuNUARef->Project3D("yz");
1881     nuaRef->Scale(1./fCumuNUARef->GetNbinsX());
1882     nuaRef->SetName(Form("%sReferenceNUA%s", fType.Data(), GetQCType(fFlags)));
1883     nuaRef->SetTitle(Form("%sReferenceNUA%s", fType.Data(), GetQCType(fFlags)));
1884     nualist->Add(nuaRef);
1885   } else {
1886     temp = (TH2D*)fCumuNUARef->Project3D("yz");
1887     temp->Scale(1./fCumuNUARef->GetNbinsX());
1888     nuaRef->Add(temp);
1889     delete temp;
1890   }
1891   // Filling in underflow to make scaling possible in Terminate()
1892   nuaRef->Fill(0., -1., 1.);
1893   // Differential
1894   TH2D* nuaDiff = (TH2D*)nualist->FindObject(Form("%sDifferentialNUA%s", fType.Data(), GetQCType(fFlags)));
1895   if (!nuaDiff) {
1896     nuaDiff = (TH2D*)fCumuNUADiff->Project3D("yz");
1897     nuaDiff->SetName(Form("%sDifferentialNUA%s", fType.Data(), GetQCType(fFlags)));
1898     nuaDiff->SetTitle(Form("%sDifferentialNUA%s", fType.Data(), GetQCType(fFlags)));
1899     nualist->Add(nuaDiff);
1900   } else {
1901     temp = (TH2D*)fCumuNUADiff->Project3D("yz");
1902     nuaDiff->Add(temp);
1903     delete temp;
1904   }
1905   // Filling in underflow to make scaling possible in Terminate()
1906   nuaDiff->Fill(0., -1., 1.);
1907
1908   return;
1909 }
1910 //_____________________________________________________________________
1911 void AliForwardFlowTaskQC::VertexBin::CalculateReferenceFlow(CumuHistos& cumu2h, CumuHistos& cumu4h, TH2I* quality,
1912                                                              TH1D* chist, TH2D* dNdetaRef) const
1913 {
1914   //
1915   //  Calculates the reference flow
1916   //
1917   //  Parameters:
1918   //   cumu2h: CumuHistos object with QC{2} cumulants
1919   //   cumu4h: CumuHistos object with QC{4} cumulants
1920   //   quality: Histogram for success rate of cumulants
1921   //   chist: Centrality histogram
1922   //   dNdetaRef: dN/deta histogram for estimating multiplicity used for ref calculations
1923   //
1924
1925   // Normalizing common NUA hists
1926   for (Int_t cBin = 1; cBin <= fCumuNUARef->GetNbinsY(); cBin++) {
1927     Double_t cent = fCumuNUARef->GetYaxis()->GetBinCenter(cBin);
1928     for (Int_t eBin = 1; eBin <= fCumuNUARef->GetNbinsX(); eBin++) {
1929       Double_t eta = fCumuNUARef->GetXaxis()->GetBinCenter(eBin);
1930       Double_t mult = fCumuNUARef->GetBinContent(eBin, cBin, 0);
1931       if (mult == 0) continue;
1932       for (Int_t qBin = 1; qBin <= fCumuNUARef->GetNbinsZ(); qBin++) {
1933         fCumuNUARef->SetBinContent(eBin, cBin, qBin, fCumuNUARef->GetBinContent(eBin, cBin, qBin)/mult);
1934         fCumuNUARef->SetBinError(eBin, cBin, qBin, fCumuNUARef->GetBinError(eBin, cBin, qBin)/mult);
1935       }
1936       // Fill dN/deta diagnostics
1937       dNdetaRef->Fill(eta, cent, mult);
1938     }
1939   }
1940
1941   // For flow calculations
1942   TH3D* cumuRef = 0;
1943   TH2D* cumu2 = 0;
1944   TH2D* cumu2NUAold = 0;
1945   TH2D* cumu4 = 0;
1946   TH2D* cumu4NUA = 0;
1947   Int_t qualityFactor = ((fFlags & kStdQC) ? 8 : 4); // used for correctly filling in quality hists
1948   // Loop over cumulant histogram for final calculations
1949   for (Int_t n = 2; n <= fMaxMoment; n++) { // Moment loop begins
1950     cumu2 = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNoNUA);
1951     if ((fFlags & kNUAcorr)) {
1952       cumu2NUAold = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNUAOld);
1953     }
1954     if ((fFlags & kStdQC)) {
1955       cumu4 = (TH2D*)cumu4h.Get('r', n, CumuHistos::kNoNUA);
1956       if ((fFlags & kNUAcorr)) cumu4NUA = (TH2D*)cumu4h.Get('r', n, CumuHistos::kNUAOld);
1957     }
1958     cumuRef  = (TH3D*)fCumuHists.Get('r', n);
1959     // Begin loops
1960     for (Int_t cBin = 1; cBin <= cumuRef->GetNbinsY(); cBin++) { // Centrality loop begins
1961       Double_t cent = cumuRef->GetYaxis()->GetBinCenter(cBin);
1962       if (n == 2) chist->Fill(cent, cumuRef->GetBinContent(0, cBin, 0));
1963       if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), n, cent));
1964       for (Int_t etaBin = 1; etaBin <= cumuRef->GetNbinsX(); etaBin++) { // Eta loop begins
1965         Double_t eta = cumuRef->GetXaxis()->GetBinCenter(etaBin);
1966         Double_t refEta = eta;
1967         Int_t refEtaBinA = fCumuNUARef->GetXaxis()->FindBin(refEta);
1968         if ((fFlags & kEtaGap)) refEta = -eta;
1969         Int_t refEtaBinB = fCumuNUARef->GetXaxis()->FindBin(refEta);
1970         // 2-particle reference flow
1971         Double_t w2Two = cumuRef->GetBinContent(refEtaBinA, cBin, kW2Two);
1972         Double_t w2 = cumuRef->GetBinContent(refEtaBinA, cBin, kW2);
1973         if (w2 == 0) continue;
1974         Double_t cosP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberCos(n));
1975         Double_t sinP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberSin(n));
1976         Double_t cosP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(n));
1977         Double_t sinP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(n));
1978         Double_t cos2nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberCos(2*n));
1979         Double_t sin2nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberSin(2*n));
1980         Double_t cos2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(2*n));
1981         Double_t sin2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(2*n));
1982         Double_t two = w2Two / w2;
1983         Double_t qc2 = two;
1984         if (qc2 >= 0) cumu2->Fill(eta, cent, TMath::Sqrt(qc2));
1985
1986         if ((fFlags & kNUAcorr)) {
1987           // Old NUA
1988           // With no eta gap the last two terms are <<cos(phi)>>^2 and <<sin(phi)>>^2,
1989           // with eta gap the different coverage is taken into account.
1990           // The next line covers both cases.
1991           qc2 -= cosP1nPhiA*cosP1nPhiB + sinP1nPhiA*sinP1nPhiB;
1992           // Extra NUA term from 2n cosines and sines
1993           Double_t den = 1-(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB);
1994           if (den != 0) qc2 /= den;
1995           else qc2 = 0;
1996         }
1997         if (qc2 <= 0) {
1998           if (fDebug > 0)
1999             AliInfo(Form("%s: QC_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f - skipping",
2000             fType.Data(), n, qc2, eta, cent));
2001           quality->Fill((n-2)*qualityFactor+2, Int_t(cent));
2002           continue;
2003         }
2004         Double_t vnTwo = TMath::Sqrt(qc2);
2005         if (!TMath::IsNaN(vnTwo)) {
2006           quality->Fill((n-2)*qualityFactor+1, Int_t(cent));
2007           if ((fFlags & kNUAcorr)) cumu2NUAold->Fill(eta, cent, vnTwo);
2008         }
2009
2010         if (!(fFlags & kStdQC)) continue;
2011         // 4-particle reference flow
2012         Double_t w4Four = cumuRef->GetBinContent(refEtaBinA, cBin, kW4Four);
2013         Double_t w4 = cumuRef->GetBinContent(refEtaBinA, cBin, kW4);
2014         Double_t multm1m2 = cumuRef->GetBinContent(refEtaBinA, cBin, k3pWeight);
2015         if (w4 == 0 || multm1m2 == 0) continue;
2016         Double_t cosP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2);
2017         Double_t sinP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2);
2018         Double_t cosP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2phi3m);
2019         Double_t sinP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2phi3m);
2020
2021         cosP1nPhi1P1nPhi2 /= w2;
2022         sinP1nPhi1P1nPhi2 /= w2;
2023         cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2024         sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2025         Double_t four = w4Four / w4;
2026         Double_t qc4 = four-2.*TMath::Power(two,2.);
2027         if (qc4 < 0) cumu4->Fill(eta, cent, TMath::Power(-qc4, 0.25));
2028
2029         if ((fFlags & kNUAcorr)) {
2030            qc4 += - 4.*cosP1nPhiA*cosP1nPhi1M1nPhi2M1nPhi3
2031                   + 4.*sinP1nPhiA*sinP1nPhi1M1nPhi2M1nPhi3-TMath::Power(cosP1nPhi1P1nPhi2,2.)-TMath::Power(sinP1nPhi1P1nPhi2,2.)
2032                   + 4.*cosP1nPhi1P1nPhi2*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
2033                   + 8.*sinP1nPhi1P1nPhi2*sinP1nPhiA*cosP1nPhiA
2034                   + 8.*two*(TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.))
2035                   - 6.*TMath::Power((TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.)),2.);
2036         }
2037         if (qc4 >= 0) {
2038           if (fDebug > 0)
2039             AliInfo(Form("%s: QC_%d{4} = %1.3f for eta = %1.2f and centrality %3.1f - skipping",
2040             fType.Data(), n, qc2, eta, cent));
2041           quality->Fill((n-2)*qualityFactor+6, Int_t(cent));
2042           continue;
2043         }
2044         Double_t vnFour = TMath::Power(-qc4, 0.25);
2045         if (!TMath::IsNaN(vnFour*multm1m2)) {
2046           quality->Fill((n-2)*qualityFactor+5, Int_t(cent));
2047           if ((fFlags & kNUAcorr)) cumu4NUA->Fill(eta, cent, vnFour);
2048         }
2049       } // End of eta
2050     } // End of cent
2051   } // End of moment
2052
2053   return;
2054 }
2055 //_____________________________________________________________________
2056 void AliForwardFlowTaskQC::VertexBin::CalculateDifferentialFlow(CumuHistos& cumu2h, CumuHistos& cumu4h,
2057                                                                 TH2I* quality, TH2D* dNdetaDiff) const
2058 {
2059   //
2060   //  Calculates the differential flow
2061   //
2062   //  Parameters:
2063   //   cumu2h: CumuHistos object with QC{2} cumulants
2064   //   cumu4h: CumuHistos object with QC{4} cumulants
2065   //   quality: Histogram for success rate of cumulants
2066   //   dNdetaDiff: dN/deta histogram for estimating multiplicity used for diff calculations
2067   //
2068
2069   for (Int_t cBin = 1; cBin <= fCumuNUADiff->GetNbinsY(); cBin++) {
2070     Double_t cent = fCumuNUADiff->GetYaxis()->GetBinCenter(cBin);
2071     for (Int_t eBin = 1; eBin <= fCumuNUADiff->GetNbinsX(); eBin++) {
2072       Double_t eta = fCumuNUADiff->GetXaxis()->GetBinCenter(eBin);
2073       Double_t mult = fCumuNUADiff->GetBinContent(eBin, cBin, 0);
2074       if (mult == 0) continue;
2075       for (Int_t qBin = 1; qBin <= fCumuNUADiff->GetNbinsZ(); qBin++) {
2076         fCumuNUADiff->SetBinContent(eBin, cBin, qBin, fCumuNUADiff->GetBinContent(eBin, cBin, qBin)/mult);
2077         fCumuNUADiff->SetBinError(eBin, cBin, qBin, fCumuNUADiff->GetBinError(eBin, cBin, qBin)/mult);
2078       }
2079       dNdetaDiff->Fill(eta, cent, mult);
2080     }
2081   }
2082
2083   // For flow calculations
2084   TH3D* cumuRef = 0;
2085   TH3D* cumuDiff = 0;
2086   TH2D* cumu2 = 0;
2087   TH2D* cumu2NUAold = 0;
2088   TH2D* cumu4 = 0;
2089   TH2D* cumu4NUA = 0;
2090   Int_t qualityFactor = ((fFlags & kStdQC) ? 8 : 4); // used for correctly filling in quality hists
2091   // Loop over cumulant histogram for final calculations
2092   for (Int_t n = 2; n <= fMaxMoment; n++) { // Moment loop begins
2093     cumu2 = (TH2D*)cumu2h.Get('d', n, CumuHistos::kNoNUA);
2094     if ((fFlags & kNUAcorr)) {
2095       cumu2NUAold = (TH2D*)cumu2h.Get('d', n, CumuHistos::kNUAOld);
2096     }
2097     if ((fFlags & kStdQC)) {
2098       cumu4 = (TH2D*)cumu4h.Get('d',n);
2099       if ((fFlags & kNUAcorr)) cumu4NUA = (TH2D*)cumu4h.Get('d', n, CumuHistos::kNUAOld);
2100     }
2101     cumuRef  = (TH3D*)fCumuHists.Get('r',n);
2102     cumuDiff = (TH3D*)fCumuHists.Get('d',n);
2103     for (Int_t cBin = 1; cBin <= cumuDiff->GetNbinsY(); cBin++) { // Centrality loop begins
2104       Double_t cent = cumuDiff->GetYaxis()->GetBinCenter(cBin);
2105       if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), n, cent));
2106       for (Int_t etaBin = 1; etaBin <= cumuDiff->GetNbinsX(); etaBin++) { // Eta loop begins
2107         Double_t eta = cumuDiff->GetXaxis()->GetBinCenter(etaBin);
2108         Double_t refEta = eta;
2109         Int_t refEtaBinA = fCumuNUARef->GetXaxis()->FindBin(refEta);
2110         if ((fFlags & kEtaGap)) refEta = -eta;
2111         Int_t refEtaBinB = fCumuNUARef->GetXaxis()->FindBin(refEta);
2112
2113         // Reference objects
2114         Double_t w2 = cumuRef->GetBinContent(refEtaBinA, cBin, kW2);
2115         if (w2 == 0) continue;
2116         Double_t two = cumuRef->GetBinContent(refEtaBinA, cBin, kW2Two);
2117         two /= w2;
2118         Double_t cosP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberCos(n));
2119         Double_t sinP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberSin(n));
2120         Double_t cosP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(n));
2121         Double_t sinP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(n));
2122         Double_t cos2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(2*n));
2123         Double_t sin2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(2*n));
2124
2125         // 2-particle differential flow
2126         Double_t w2pTwoPrime = cumuDiff->GetBinContent(etaBin, cBin, kW2Two);
2127         Double_t w2p = cumuDiff->GetBinContent(etaBin, cBin, kW2);
2128         if (w2p == 0) continue;
2129         Double_t cosP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(n));
2130         Double_t sinP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(n));
2131         Double_t cos2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(2*n));
2132         Double_t sin2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(2*n));
2133         Double_t twoPrime = w2pTwoPrime / w2p;
2134
2135         Double_t qc2Prime = twoPrime;
2136         cumu2->Fill(eta, cent, qc2Prime);
2137         if ((fFlags & kNUAcorr)) {
2138           // Old nua
2139           qc2Prime -= cosP1nPsi*cosP1nPhiB + sinP1nPsi*sinP1nPhiB;
2140           // Extra NUA term from 2n cosines and sines
2141           qc2Prime /= (1.-(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
2142         }
2143         if (!TMath::IsNaN(qc2Prime)) {
2144           quality->Fill((n-2)*qualityFactor+3, Int_t(cent));
2145           if ((fFlags & kNUAcorr)) cumu2NUAold->Fill(eta, cent, qc2Prime);
2146         }
2147         else
2148           quality->Fill((n-2)*qualityFactor+4, Int_t(cent));
2149         if (fDebug > 1)
2150           AliInfo(Form("%s: QC'_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f",
2151           fType.Data(), n, qc2Prime, eta, cent));
2152
2153         if (!(fFlags & kStdQC)) continue;
2154         // Reference objects
2155         Double_t cosP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2);
2156         Double_t sinP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2);
2157         Double_t cosP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2phi3m);
2158         Double_t sinP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2phi3m);
2159         Double_t multm1m2 = cumuRef->GetBinContent(refEtaBinA, cBin, k3pWeight);
2160         cosP1nPhi1P1nPhi2 /= w2;
2161         sinP1nPhi1P1nPhi2 /= w2;
2162         cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2163         sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2164
2165         // 4-particle differential flow
2166         Double_t w4pFourPrime = cumuDiff->GetBinContent(etaBin, cBin, kW4Four);
2167         Double_t w4p = cumuDiff->GetBinContent(etaBin, cBin, kW4);
2168         Double_t mpqMult = cumuDiff->GetBinContent(etaBin, cBin, k3pWeight);
2169         if (w4p == 0 || mpqMult == 0) continue;
2170         Double_t cosP1nPsi1P1nPhi2 = cumuDiff->GetBinContent(etaBin, cBin, kCosphi1phi2);
2171         Double_t sinP1nPsi1P1nPhi2 = cumuDiff->GetBinContent(etaBin, cBin, kSinphi1phi2);
2172         Double_t cosP1nPsi1M1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kCosphi1phi2phi3m);
2173         Double_t sinP1nPsi1M1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kSinphi1phi2phi3m);
2174         Double_t cosP1nPsi1P1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kCosphi1phi2phi3p);
2175         Double_t sinP1nPsi1P1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kSinphi1phi2phi3p);
2176
2177         cosP1nPsi1P1nPhi2 /= w2p;
2178         sinP1nPsi1P1nPhi2 /= w2p;
2179         cosP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
2180         sinP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
2181         cosP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
2182         sinP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
2183
2184         Double_t fourPrime = w4pFourPrime / w4p;
2185         Double_t qc4Prime = fourPrime-2.*twoPrime*two;
2186         if (cumu4) cumu4->Fill(eta, cent, qc4Prime);
2187
2188         if ((fFlags & kNUAcorr)) {
2189           qc4Prime += - cosP1nPsi*cosP1nPhi1M1nPhi2M1nPhi3
2190                       + sinP1nPsi*sinP1nPhi1M1nPhi2M1nPhi3
2191                       - cosP1nPhiA*cosP1nPsi1M1nPhi2M1nPhi3
2192                       + sinP1nPhiA*sinP1nPsi1M1nPhi2M1nPhi3
2193                       - 2.*cosP1nPhiA*cosP1nPsi1P1nPhi2M1nPhi3
2194                       - 2.*sinP1nPhiA*sinP1nPsi1P1nPhi2M1nPhi3
2195                       - cosP1nPsi1P1nPhi2*cosP1nPhi1P1nPhi2
2196                       - sinP1nPsi1P1nPhi2*sinP1nPhi1P1nPhi2
2197                       + 2.*cosP1nPhi1P1nPhi2*(cosP1nPsi*cosP1nPhiA-sinP1nPsi*sinP1nPhiA)
2198                       + 2.*sinP1nPhi1P1nPhi2*(cosP1nPsi*sinP1nPhiA+sinP1nPsi*cosP1nPhiA)
2199                       + 4.*two*(cosP1nPsi*cosP1nPhiA+sinP1nPsi*sinP1nPhiA)
2200                       + 2.*cosP1nPsi1P1nPhi2*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
2201                       + 4.*sinP1nPsi1P1nPhi2*cosP1nPhiA*sinP1nPhiA
2202                       + 4.*twoPrime*(TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.))
2203                       - 6.*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
2204                       * (cosP1nPsi*cosP1nPhiA-sinP1nPsi*sinP1nPhiA)
2205                       - 12.*cosP1nPhiA*sinP1nPhiA
2206                       * (sinP1nPsi*cosP1nPhiA+cosP1nPsi*sinP1nPhiA);
2207         }
2208 //      Double_t vnFourDiff = - qc4Prime / TMath::Power(-qc4, 0.75);
2209         if (!TMath::IsNaN(qc4Prime*mpqMult)) {
2210           quality->Fill((n-2)*qualityFactor+7, Int_t(cent));
2211           if (cumu4NUA) cumu4NUA->Fill(eta, cent, qc4Prime);
2212         }
2213         else
2214           quality->Fill((n-2)*qualityFactor+8, Int_t(cent));
2215         if (fDebug > 1)
2216           AliInfo(Form("%s: v_%d{4} = %1.3f for eta = %1.2f and centrality %3.1f",
2217           fType.Data(), n, qc4Prime, eta, cent));
2218       } // End of eta loop
2219     } // End of centrality loop
2220   } // End of moment
2221
2222   return;
2223 }
2224 //_____________________________________________________________________
2225 void AliForwardFlowTaskQC::VertexBin::Calculate3CorFlow(CumuHistos& cumu2h, TH2I* quality, TH1D* chist,
2226                                                         TH2D* dNdetaRef, TH2D* dNdetaDiff) const
2227 {
2228   //
2229   //  Calculates the 3 sub flow
2230   //
2231   //  Parameters:
2232   //   cumu2h: CumuHistos object with QC{2} cumulants
2233   //   quality: Histogram for success rate of cumulants
2234   //   chist: Centrality histogram
2235   //   dNdetaDiff: dN/deta histogram for estimating multiplicity used for diff calculations
2236   //
2237
2238   // For flow calculations
2239   TH3D* cumuRef = 0;
2240   TH3D* cumuDiff = 0;
2241   TH2D* cumu2r = 0;
2242   TH2D* cumu2rNUAold = 0;
2243   TH2D* cumu2a = 0;
2244   TH2D* cumu2aNUAold = 0;
2245   TH2D* cumu2b = 0;
2246   TH2D* cumu2bNUAold = 0;
2247   // Loop over cumulant histogram for final calculations
2248   for (Int_t n = 2; n <= fMaxMoment; n++) { // Moment loop begins
2249     cumu2r = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNoNUA);
2250     cumu2a = (TH2D*)cumu2h.Get('a', n, CumuHistos::kNoNUA);
2251     cumu2b = (TH2D*)cumu2h.Get('b', n, CumuHistos::kNoNUA);
2252     if ((fFlags & kNUAcorr)) {
2253       cumu2rNUAold = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNUAOld);
2254       cumu2aNUAold = (TH2D*)cumu2h.Get('a', n, CumuHistos::kNUAOld);
2255       cumu2bNUAold = (TH2D*)cumu2h.Get('b', n, CumuHistos::kNUAOld);
2256     }
2257     cumuRef  = (TH3D*)fCumuHists.Get('r',n);
2258     cumuDiff = (TH3D*)fCumuHists.Get('d',n);
2259     for (Int_t cBin = 1; cBin <= cumuRef->GetNbinsY(); cBin++) { // Centrality loop begins
2260       Double_t cent = cumuRef->GetYaxis()->GetBinCenter(cBin);
2261       if (n == 2) chist->Fill(cent, cumuRef->GetBinContent(0, cBin, 0));
2262       if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), n, cent));
2263       // Here it starts!
2264       Int_t prevLim = 0;
2265       Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2266       Double_t cosP1nPhiA = 0;
2267       Double_t sinP1nPhiA = 0;
2268       Double_t cos2nPhiA = 0;
2269       Double_t sin2nPhiA = 0;
2270       Double_t cosP1nPhiB = 0;
2271       Double_t sinP1nPhiB = 0;
2272       Double_t cos2nPhiB = 0;
2273       Double_t sin2nPhiB = 0;
2274       Double_t multA = 0;
2275       Double_t multB = 0;
2276
2277       for (Int_t etaBin = 1; etaBin <= cumuDiff->GetNbinsX(); etaBin++) { // Eta loop begins
2278         Double_t eta = cumuDiff->GetXaxis()->GetBinCenter(etaBin);
2279         // 2-particle reference flow
2280         Double_t w2Two = cumuRef->GetBinContent(etaBin, cBin, kW2Two);
2281         Double_t w2 = cumuRef->GetBinContent(etaBin, cBin, kW2);
2282         if (w2 == 0) continue;
2283
2284         // Update NUA for new range!
2285         if (fEtaLims[prevLim] < eta) {
2286           GetLimits(prevLim, aLow, aHigh, bLow, bHigh);
2287           prevLim++;
2288           cosP1nPhiA = 0; sinP1nPhiA = 0; cos2nPhiA = 0; sin2nPhiA = 0; multA = 0;
2289           cosP1nPhiB = 0; sinP1nPhiB = 0; cos2nPhiB = 0; sin2nPhiB = 0; multB = 0;
2290           for (Int_t a = aLow; a <= aHigh; a++) {
2291             cosP1nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n));
2292             sinP1nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n));
2293             cos2nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(2*n));
2294             sin2nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(2*n));
2295             multA += fCumuNUARef->GetBinContent(a, cBin, 0);
2296           }
2297           for (Int_t b = bLow; b <= bHigh; b++) {
2298             cosP1nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n));
2299             sinP1nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n));
2300             cos2nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(2*n));
2301             sin2nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(2*n));
2302             multB += fCumuNUARef->GetBinContent(b, cBin, 0);
2303           }
2304           if (multA == 0 || multB == 0) {
2305             AliWarning(Form("Empty NUA values for 3Cor! (%s)", cumuRef->GetName()));
2306             continue;
2307           }
2308           cosP1nPhiA /= multA;
2309           sinP1nPhiA /= multA;
2310           cos2nPhiA /= multA;
2311           sin2nPhiA /= multA;
2312           cosP1nPhiB /= multB;
2313           sinP1nPhiB /= multB;
2314           cos2nPhiB /= multB;
2315           sin2nPhiB /= multB;
2316
2317           dNdetaRef->Fill(eta, cent, multA+multB);
2318         }
2319         Double_t two = w2Two / w2;
2320
2321         Double_t qc2 = two;
2322         if (qc2 >= 0) cumu2r->Fill(eta, cent, TMath::Sqrt(qc2));
2323
2324         if ((fFlags & kNUAcorr)) {
2325           // Old nua
2326           qc2 -= cosP1nPhiA*cosP1nPhiB + sinP1nPhiA*sinP1nPhiB;
2327           // Extra NUA term from 2n cosines and sines
2328           qc2 /= (1-(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB));
2329         }
2330         if (qc2 <= 0) {
2331           if (fDebug > 0)
2332             AliInfo(Form("%s: QC_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f - skipping",
2333             fType.Data(), n, qc2, eta, cent));
2334           quality->Fill((n-2)*4+2, Int_t(cent));
2335           continue;
2336         }
2337         Double_t vnTwo = TMath::Sqrt(qc2);
2338         if (!TMath::IsNaN(vnTwo)) {
2339           quality->Fill((n-2)*4+1, Int_t(cent));
2340           if ((fFlags & kNUAcorr)) cumu2rNUAold->Fill(eta, cent, vnTwo);
2341         }
2342
2343         // 2-particle differential flow
2344         Double_t w2pTwoPrimeA = cumuDiff->GetBinContent(etaBin, cBin, kW2Two);
2345         Double_t w2pA = cumuDiff->GetBinContent(etaBin, cBin, kW2);
2346         Double_t w2pTwoPrimeB = cumuDiff->GetBinContent(etaBin, cBin, kW4Four);
2347         Double_t w2pB = cumuDiff->GetBinContent(etaBin, cBin, kW4);
2348         if (w2pA == 0 || w2pB == 0) continue;
2349         Double_t cosP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(n));
2350         Double_t sinP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(n));
2351         Double_t cos2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(2*n));
2352         Double_t sin2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(2*n));
2353         Double_t mult = fCumuNUADiff->GetBinContent(etaBin, cBin, 0);
2354         if (mult == 0) continue;
2355         cosP1nPsi /= mult;
2356         sinP1nPsi /= mult;
2357         cos2nPsi /= mult;
2358         sin2nPsi /= mult;
2359         Double_t twoPrimeA = w2pTwoPrimeA / w2pA;
2360         Double_t twoPrimeB = w2pTwoPrimeB / w2pB;
2361         dNdetaDiff->Fill(eta, cent, mult);
2362
2363         Double_t qc2PrimeA = twoPrimeA;
2364         Double_t qc2PrimeB = twoPrimeB;
2365         if (qc2PrimeA*qc2PrimeB >= 0) {
2366           cumu2a->Fill(eta, cent, qc2PrimeA);
2367           cumu2b->Fill(eta, cent, qc2PrimeB);
2368         }
2369         if ((fFlags & kNUAcorr)) {
2370           // Old nua
2371           qc2PrimeA -= cosP1nPsi*cosP1nPhiA + sinP1nPsi*sinP1nPhiA;
2372           qc2PrimeB -= cosP1nPsi*cosP1nPhiB + sinP1nPsi*sinP1nPhiB; // Is this OK?
2373           // Extra NUA term from 2n cosines and sines
2374           if (cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA != 1.) qc2PrimeA /= (1.-(cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA));
2375           if (cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB != 1.) qc2PrimeB /= (1.-(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
2376         }
2377         if (!TMath::IsNaN(qc2PrimeA) && !TMath::IsNaN(qc2PrimeB) && qc2 != 0) {
2378         if (qc2PrimeA*qc2PrimeB >= 0) {
2379           quality->Fill((n-2)*4+3, Int_t(cent));
2380           if ((fFlags & kNUAcorr)) cumu2aNUAold->Fill(eta, cent, qc2PrimeA);
2381           if ((fFlags & kNUAcorr)) cumu2bNUAold->Fill(eta, cent, qc2PrimeB);
2382         }
2383       }
2384       else
2385         quality->Fill((n-2)*4+4, Int_t(cent));
2386         if (fDebug > 1)
2387           AliInfo(Form("%s: QC'a_%d{2} = %1.3f, QC'b_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f",
2388           fType.Data(), n, qc2PrimeA, n, qc2PrimeB, eta, cent));
2389       } // End of eta loop
2390     } // End of centrality loop
2391   } // End of moment
2392
2393   return;
2394 }
2395 //_____________________________________________________________________
2396 void AliForwardFlowTaskQC::VertexBin::SolveCoupledFlowEquations(CumuHistos& cumu, const Char_t type) const
2397 {
2398   //
2399   //  Function to solve the coupled flow equations
2400   //  We solve it by using matrix calculations:
2401   //  A*v_n = V => v_n = A^-1*V
2402   //  First we set up a TMatrixD container to make ROOT
2403   //  do the inversions in an efficient way, we multiply the current <<2>> estimates.
2404   //  Then we fill new TH2D's if the corrected <<2>>'s (cumuNUA object).
2405   //
2406   //  Parameters:
2407   //   cumu: CumuHistos object - uncorrected
2408   //   type: Reference ('r') or differential ('d') or ('a' or 'b') for 3 correlator
2409   //
2410
2411   // We start by declaring Matrix and vector objects, as their constructors are quite heavy
2412   TMatrixD mNUA(fMaxMoment-1, fMaxMoment-1);
2413   TVectorD vQC2(fMaxMoment-1);
2414
2415   for (Int_t cBin = 1; cBin <= cumu.Get(type, 2, CumuHistos::kNUAOld)->GetNbinsY(); cBin++) { // cent loop
2416     Double_t cent = cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinCenter(cBin);
2417     for (Int_t eBin = 1; eBin <= cumu.Get(type, 2, CumuHistos::kNUAOld)->GetNbinsX(); eBin++) { // eta loop
2418       Double_t eta = cumu.Get(type, 2, CumuHistos::kNUAOld)->GetXaxis()->GetBinCenter(eBin);
2419       mNUA.Zero(); // reset matrix
2420       vQC2.Zero(); // reset vector
2421       for (Int_t n = 0; n < fMaxMoment-1; n++) { // moment loop
2422         vQC2(n) = static_cast<TH2D*>(cumu.Get(type, n+2, CumuHistos::kNUAOld))->GetBinContent(eBin, cBin);
2423         if (type == 'r' || type == 'R') vQC2(n) *= vQC2(n); // going back to <<2>>
2424         for (Int_t m = 0; m < fMaxMoment-1; m++) { // cross moment
2425           mNUA(n,m) = CalculateNUAMatrixElement(n, m, type, eBin, cBin);
2426         } // End of cross moment loop
2427       } // End of moment loop
2428       // Invert matrix
2429       Double_t det = 0;
2430       mNUA.Invert(&det);
2431       // If determinant is non-zero we go with corrected results
2432       if (det != 0 ) vQC2 = mNUA*vQC2;
2433       else AliWarning(Form("Determinant == 0 - cent: %d-%d, eta: %f, type: '%c', data: %s, vtx: %d-%d%s",
2434                       Int_t(cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinLowEdge(cBin)),
2435                       Int_t(cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinUpEdge(cBin)),
2436                       cumu.Get(type, 2, CumuHistos::kNUAOld)->GetXaxis()->GetBinCenter(eBin),
2437                       type, fType.Data(), fVzMin, fVzMax,
2438                       GetQCType(fFlags, kTRUE)));
2439       // Go back to v_n for ref. keep <<2'>> for diff. flow).
2440       for (Int_t n = 0; n < fMaxMoment-1; n++) {
2441         Double_t vnTwo = 0;
2442         if (type == 'r' || type == 'R')
2443           vnTwo = (vQC2(n) > 0. ? TMath::Sqrt(vQC2(n)) : 0.);
2444         else {
2445           // is really more <<2'>> in this case
2446           vnTwo = vQC2(n);
2447         }
2448         // Fill in corrected v_n
2449         if (vnTwo != 0) static_cast<TH2D*>(cumu.Get(type, n+2, CumuHistos::kNUA))->Fill(eta, cent, vnTwo);
2450       } // End of moment loop
2451     } // End of eta loop
2452   } // End of centrality loop
2453   return;
2454 }
2455 //_____________________________________________________________________
2456 Double_t AliForwardFlowTaskQC::VertexBin::CalculateNUAMatrixElement(Int_t n, Int_t m, Char_t type, Int_t binA, Int_t cBin) const
2457 {
2458   //
2459   //  Calculates the (n,m)-th element in the NUA matrix: 1 if n == m, otherwise:
2460   //               <<cos[(n-m)phi1]>>*<<cos[(n-m)phi2]>> + <<sin[(n-m)phi1]>>*<<sin[(n-m)phi2]>>
2461   //    NUA(n,m) = -----------------------------------------------------------------------------
2462   //                   1 + <<cos[2nphi1]>>*<<cos[2nphi2]>> + <<sin[2nphi1]>>*<<sin[2nphi2]>>
2463   //
2464   //               <<cos[(n+m)phi1]>>*<<cos[(n+m)phi2]>> + <<sin[(n+m)phi1]>>*<<sin[(n+m)phi2]>>
2465   //             + -----------------------------------------------------------------------------
2466   //                   1 + <<cos[2nphi1]>>*<<cos[2nphi2]>> + <<sin[2nphi1]>>*<<sin[2nphi2]>>
2467   //
2468   //  Parameters:
2469   //   n: row
2470   //   m: coumn
2471   //   type: Reference ('r') or differential ('d') or ('a' or 'b')
2472   //   binA: eta bin of phi1
2473   //   cBin: centrality bin
2474   //
2475   //  Return: NUA(n,m)
2476   //
2477   if (n == m) return 1.;
2478   n += 2;
2479   m += 2;
2480
2481   Double_t cosnMmPhi1 = 0, cosnMmPhi2 = 0, sinnMmPhi1 = 0, sinnMmPhi2 = 0;
2482   Double_t cosnPmPhi1 = 0, cosnPmPhi2 = 0, sinnPmPhi1 = 0, sinnPmPhi2 = 0;
2483   Double_t cos2nPhi1 = 0, cos2nPhi2 = 0, sin2nPhi1 = 0, sin2nPhi2 = 0;
2484
2485   // reference flow
2486   if (type == 'r' || type == 'R') {
2487     if ((fFlags & k3Cor)) {
2488       Double_t eta = fCumuNUARef->GetXaxis()->GetBinCenter(binA);
2489       Int_t i = 0;
2490       while (fEtaLims[i] < eta) i++;
2491       Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2492       GetLimits(i-1, aLow, aHigh, bLow, bHigh);
2493       Double_t multA = 0, multB = 0;
2494       for (Int_t a = aLow; a <= aHigh; a++) {
2495         cosnMmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n-m));
2496         sinnMmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n-m));
2497         cosnPmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n+m));
2498         sinnPmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n+m));
2499         cos2nPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(2*n));
2500         sin2nPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(2*n));
2501         multA += fCumuNUARef->GetBinContent(a, cBin, 0);
2502       }
2503        for (Int_t b = bLow; b <= bHigh; b++) {
2504         cosnMmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n-m));
2505         sinnMmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n-m));
2506         cosnPmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n+m));
2507         sinnPmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n+m));
2508         cos2nPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(2*n));
2509         sin2nPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(2*n));
2510         multB += fCumuNUARef->GetBinContent(b, cBin, 0);
2511       }
2512       if (multA == 0 || multB == 0) {
2513         if (fDebug > 0) AliWarning("multA or multB == 0 in matrix elements, aborting NUA");
2514         return 0.;
2515       }
2516       cosnMmPhi1 /= multA;
2517       sinnMmPhi1 /= multA;
2518       cosnPmPhi1 /= multA;
2519       sinnPmPhi1 /= multA;
2520       cos2nPhi1 /= multA;
2521       sin2nPhi1 /= multA;
2522       cosnMmPhi2 /= multB;
2523       sinnMmPhi2 /= multB;
2524       cosnPmPhi2 /= multB;
2525       sinnPmPhi2 /= multB;
2526       cos2nPhi2 /= multB;
2527       sin2nPhi2 /= multB;
2528     } else {
2529       Int_t binB = fCumuNUARef->GetXaxis()->FindBin(-fCumuNUARef->GetXaxis()->GetBinCenter(binA));
2530       cosnMmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2531       sinnMmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2532       cosnPmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2533       sinnPmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2534       cos2nPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2535       sin2nPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2536       cosnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n-m));
2537       sinnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n-m));
2538       cosnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n+m));
2539       sinnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n+m));
2540       cos2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(2*n));
2541       sin2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(2*n));
2542     }
2543   } // differential flow
2544   else if (type == 'd' || type == 'D') {
2545     Int_t binB = fCumuNUARef->GetXaxis()->FindBin(-fCumuNUADiff->GetXaxis()->GetBinCenter(binA));
2546     cosnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2547     sinnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2548     cosnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2549     sinnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2550     cos2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2551     sin2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2552     cosnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n-m));
2553     sinnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n-m));
2554     cosnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n+m));
2555     sinnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n+m));
2556     cos2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(2*n));
2557     sin2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(2*n));
2558   } // 3 correlator part a or b
2559   else if (type == 'a' || type == 'A' || type == 'b' || type == 'B') {
2560     Double_t mult1 = 0, mult2 = 0;
2561     // POIs
2562     cosnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2563     sinnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2564     cosnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2565     sinnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2566     cos2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2567     sin2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2568     mult1 = fCumuNUADiff->GetBinContent(binA, cBin, 0);
2569     // RPs
2570     Double_t eta = fCumuNUADiff->GetXaxis()->GetBinCenter(binA);
2571     Int_t i = 0;
2572     while (fEtaLims[i] < eta) i++;
2573     Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2574     GetLimits(i-1, aLow, aHigh, bLow, bHigh);
2575     Int_t lLow = ((type == 'a' || type == 'A') ? aLow : bLow);
2576     Int_t lHigh = ((type == 'a' || type == 'A') ? aHigh : bHigh);
2577     for (Int_t l = lLow; l <= lHigh; l++) {
2578       cosnMmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(n-m));
2579       sinnMmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(n-m));
2580       cosnPmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(n+m));
2581       sinnPmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(n+m));
2582       cos2nPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(2*n));
2583       sin2nPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(2*n));
2584       mult2 += fCumuNUARef->GetBinContent(l, cBin, 0);
2585     }
2586     if (mult1 == 0 || mult2 == 0) {
2587       if (fDebug > 0) AliWarning("mult1 or mult2 == 0 in matrix elements, aborting NUA");
2588       return 0.;
2589     }
2590     cosnMmPhi1 /= mult1;
2591     sinnMmPhi1 /= mult1;
2592     cosnPmPhi1 /= mult1;
2593     sinnPmPhi1 /= mult1;
2594     cos2nPhi1 /= mult1;
2595     sin2nPhi1 /= mult1;
2596     cosnMmPhi2 /= mult2;
2597     sinnMmPhi2 /= mult2;
2598     cosnPmPhi2 /= mult2;
2599     sinnPmPhi2 /= mult2;
2600     cos2nPhi2 /= mult2;
2601     sin2nPhi2 /= mult2;
2602   }
2603
2604   // Actual calculation
2605   Double_t e = cosnMmPhi1*cosnMmPhi2 + sinnMmPhi1*sinnMmPhi2 + cosnPmPhi1*cosnPmPhi2 + sinnPmPhi1*sinnPmPhi2;
2606   Double_t den = 1 + cos2nPhi1*cos2nPhi2 + sin2nPhi1*sin2nPhi2;
2607   if (den != 0) e /= den;
2608   else return 0.;
2609
2610   return e;
2611 }
2612 //_____________________________________________________________________
2613 void AliForwardFlowTaskQC::VertexBin::AddVertexBins(CumuHistos& cumu, TList* list, UInt_t nNUA) const
2614 {
2615   //
2616   //  Add up vertex bins with flow results
2617   //
2618   //  Parameters:
2619   //   cumu: CumuHistos object with vtxbin results
2620   //   list: Outout list with added results
2621   //   nNUA: # of NUA histograms to loop over
2622   //
2623   TH2D* vtxHist = 0;
2624   TProfile2D* avgProf = 0;
2625   TString name;
2626   Int_t nT = ((fFlags & k3Cor) ? 3 : 2);
2627   Char_t ct = '\0';
2628   for (UInt_t nua = 0; nua <= nNUA; nua++) { // NUA loop
2629     for (Int_t n = 2; n <= fMaxMoment; n++) { // moment loop
2630       for (Int_t t = 0; t < nT; t++) { // type loop (r/d/a/b)
2631         // Find type
2632         switch (t) {
2633           case 0: ct = 'r'; break;
2634           case 1: ct = ((fFlags & k3Cor) ? 'a' : 'd'); break;
2635           case 2: ct = 'b'; break;
2636           default: ct = '\0'; break;
2637         }
2638         vtxHist = static_cast<TH2D*>(cumu.Get(ct,n,nua));
2639         if (!vtxHist) {
2640           AliWarning("VertexBin::AddVertexBins: vtxHist not found!");
2641           continue;
2642         }
2643         name = vtxHist->GetName();
2644         // Strip name of vtx info
2645         Ssiz_t l = name.Last('x')-3;
2646         name.Resize(l);
2647         avgProf = (TProfile2D*)list->FindObject(name.Data());
2648         // if no output profile yet, make one
2649         if (!avgProf) {
2650           avgProf = new TProfile2D(name.Data(), name.Data(),
2651               vtxHist->GetNbinsX(), vtxHist->GetXaxis()->GetXmin(), vtxHist->GetXaxis()->GetXmax(),
2652               vtxHist->GetNbinsY(), vtxHist->GetYaxis()->GetXmin(), vtxHist->GetYaxis()->GetXmax());
2653           if (vtxHist->GetXaxis()->IsVariableBinSize())
2654             avgProf->GetXaxis()->Set(vtxHist->GetNbinsX(), vtxHist->GetXaxis()->GetXbins()->GetArray());
2655           if (vtxHist->GetYaxis()->IsVariableBinSize())
2656             avgProf->GetYaxis()->Set(vtxHist->GetNbinsY(), vtxHist->GetYaxis()->GetXbins()->GetArray());
2657           list->Add(avgProf);
2658         }
2659         // Fill in, cannot be done with Add function.
2660         for (Int_t e = 1; e <= vtxHist->GetNbinsX(); e++) { // eta loop
2661           Double_t eta = vtxHist->GetXaxis()->GetBinCenter(e);
2662           for (Int_t c = 1; c <= vtxHist->GetNbinsY(); c++) { // cent loop
2663             Double_t cent = vtxHist->GetYaxis()->GetBinCenter(c);
2664             Double_t cont = vtxHist->GetBinContent(e, c);
2665             if (cont == 0) continue;
2666             avgProf->Fill(eta, cent, cont);
2667           } // End of cent loop
2668         } //  End of eta loop
2669       } // End of type loop
2670     } // End of moment loop
2671   } // End of nua loop
2672 }
2673 //_____________________________________________________________________
2674 Int_t AliForwardFlowTaskQC::VertexBin::GetBinNumberCos(Int_t n) const
2675 {
2676   //
2677   //  Get the bin number of <<cos(nphi)>>
2678   //
2679   //  Parameters:
2680   //   n: moment
2681   //
2682   //  Return: bin number
2683   //
2684   Int_t bin = 0;
2685   n = TMath::Abs(n);
2686
2687   if (n == 0) bin = fMaxMoment*4-1;
2688   else        bin = n*2-1;
2689
2690   return bin;
2691 }
2692 //_____________________________________________________________________
2693 Int_t AliForwardFlowTaskQC::VertexBin::GetBinNumberSin(Int_t n) const
2694 {
2695   //
2696   //  Get the bin number of <<sin(nphi)>>
2697   //
2698   //  Parameters:
2699   //   n: moment
2700   //
2701   //  Return: bin number
2702   //
2703   Int_t bin = 0;
2704   n = TMath::Abs(n);
2705
2706   if (n == 0) bin = fMaxMoment*4;
2707   else        bin = n*2;
2708
2709   return bin;
2710 }
2711 //_____________________________________________________________________
2712 void AliForwardFlowTaskQC::VertexBin::SetupNUALabels(TAxis* a) const
2713 {
2714   //
2715   //  Setup NUA labels on axis
2716   //
2717   //  Parameters:
2718   //   a: Axis to set up
2719   //
2720   if (a->GetNbins() != GetBinNumberSin()) AliFatal("SetupNUALabels: Wrong number of bins on axis");
2721
2722   Int_t i = 1, j= 1;
2723   while (i <= a->GetNbins()) {
2724     a->SetBinLabel(i++, Form("<<cos(%d#varphi)>>", j));
2725     a->SetBinLabel(i++, Form("<<sin(%d#varphi)>>", j++));
2726   }
2727
2728   return;
2729 }
2730 //_____________________________________________________________________
2731 TH2I* AliForwardFlowTaskQC::VertexBin::MakeQualityHist(const Char_t* name) const
2732 {
2733   //
2734   //  Add a histogram for checking the analysis quality
2735   //
2736   //  Parameters:
2737   //   const char*: name of data type
2738   //
2739   //  Return: histogram for tracking successful calculations
2740   //
2741   Int_t nBins = ((fFlags & kStdQC) ? 8 : 4);
2742   TH2I* quality = new TH2I(name, name, (fMaxMoment-1)*nBins, 1, (fMaxMoment-1)*nBins+1, fCumuNUARef->GetNbinsY(),
2743                            fCumuNUARef->GetYaxis()->GetXmin(), fCumuNUARef->GetYaxis()->GetXmax());
2744   quality->GetYaxis()->Set(fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXbins()->GetArray());
2745   for (Int_t i = 2, j = 1; i <= fMaxMoment; i++) {
2746     quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{2} > 0", 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     if ((fFlags & kStdQC)) {
2751       quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{4} < 0", i));
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     }
2756   }
2757
2758   return quality;
2759 }
2760 //_____________________________________________________________________
2761 TH2D* AliForwardFlowTaskQC::VertexBin::MakeOutputHist(Int_t qc, Int_t n, const Char_t* ctype, UInt_t nua) const
2762 {
2763   //
2764   //  Setup a TH2D for the output
2765   //
2766   //  Parameters:
2767   //   qc   # of particle correlations
2768   //   n    flow moment
2769   //   ref  Type: r/d/a/b
2770   //   nua  For nua corrected hists?
2771   //
2772   //  Return: 2D hist for results
2773   //
2774   Bool_t ref = ((ctype[0] == 'R') || (ctype[0] == 'r'));
2775   TAxis* xAxis = (ref ? fCumuNUARef->GetXaxis() : fCumuNUADiff->GetXaxis());
2776   TAxis* yAxis = (ref ? fCumuNUARef->GetYaxis() : fCumuNUADiff->GetYaxis());
2777   TString ntype = "";
2778   switch (nua) {
2779     case CumuHistos::kNoNUA: ntype = "Un"; break;
2780     case CumuHistos::kNUAOld: ntype = "NUAOld"; break;
2781     case CumuHistos::kNUA: ntype = "NUA"; break;
2782     default: break;
2783   }
2784   TH2D* h = new TH2D(Form("%sQC%d_v%d_%s_%sCorr%s_vtx_%d_%d",
2785                         fType.Data(), qc, n, ctype, ntype.Data(),
2786                         GetQCType(fFlags), fVzMin, fVzMax),
2787                       Form("%sQC%d_v%d_%s_%sCorr%s_vtx_%d_%d",
2788                         fType.Data(), qc, n, ctype, ntype.Data(),
2789                         GetQCType(fFlags), fVzMin, fVzMax),
2790                       xAxis->GetNbins(), xAxis->GetXmin(), xAxis->GetXmax(),
2791                       yAxis->GetNbins(), yAxis->GetXmin(), yAxis->GetXmax());
2792   if (xAxis->IsVariableBinSize()) h->GetXaxis()->Set(xAxis->GetNbins(), xAxis->GetXbins()->GetArray());
2793   h->GetYaxis()->Set(yAxis->GetNbins(), yAxis->GetXbins()->GetArray());
2794
2795   return h;
2796 }
2797 //_____________________________________________________________________
2798 void AliForwardFlowTaskQC::PrintFlowSetup() const
2799 {
2800   //
2801   //  Print the setup of the flow task
2802   //
2803   Printf("=======================================================");
2804   Printf("%s::Print", this->IsA()->GetName());
2805   Printf("Forward detector:                 :\t%s", ((fFlowFlags & kFMD) ? "FMD" : "VZERO"));
2806   Printf("Number of bins in vertex axis     :\t%d", fVtxAxis->GetNbins());
2807   Printf("Range of vertex axis              :\t[%3.1f,%3.1f]",
2808                           fVtxAxis->GetXmin(), fVtxAxis->GetXmax());
2809   Printf("Number of bins in centrality axis :\t%d", fCentAxis->GetNbins());
2810   printf("Centrality binning                :\t");
2811   for (Int_t cBin = 1; cBin <= fCentAxis->GetNbins(); cBin++) {
2812     printf("%02d-%02d%% ", Int_t(fCentAxis->GetBinLowEdge(cBin)), Int_t(fCentAxis->GetBinUpEdge(cBin)));
2813     if (cBin == fCentAxis->GetNbins()) printf("\n");
2814     else if (cBin % 4 == 0) printf("\n\t\t\t\t\t");
2815   }
2816   printf("Doing flow analysis for           :\t");
2817   for (Int_t n  = 2; n <= fMaxMoment; n++) printf("v%d ", n);
2818   printf("\n");
2819   TString type = "Standard QC{2} and QC{4} calculations.";
2820   if ((fFlowFlags & kEtaGap)) type = "QC{2} with a rapidity gap.";
2821   if ((fFlowFlags & k3Cor))   type = "QC{2} with 3 correlators.";
2822   if ((fFlowFlags & kTPC) == kTPC)    type.ReplaceAll(".", " with TPC tracks for reference.");
2823   if ((fFlowFlags & kHybrid) == kHybrid) type.ReplaceAll(".", " with hybrid tracks for reference.");
2824   Printf("QC calculation type               :\t%s", type.Data());
2825   Printf("Symmetrize ref. flow wrt. eta = 0 :\t%s", ((fFlowFlags & kSymEta) ? "true" : "false"));
2826   Printf("Apply NUA correction terms        :\t%s", ((fFlowFlags & kNUAcorr) ? "true" : "false"));
2827   Printf("Satellite vertex flag             :\t%s", ((fFlowFlags & kSatVtx) ? "true" : "false"));
2828   Printf("FMD sigma cut:                    :\t%f", fFMDCut);
2829   Printf("SPD sigma cut:                    :\t%f", fSPDCut);
2830   if ((fFlowFlags & kEtaGap) || (fFlowFlags & kTracks))
2831     Printf("Eta gap:                          :\t%f", fEtaGap);
2832   Printf("=======================================================");
2833 }
2834 //_____________________________________________________________________
2835 const Char_t* AliForwardFlowTaskQC::GetQCType(UShort_t flags, Bool_t prependUS)
2836 {
2837   //
2838   //  Get the type of the QC calculations
2839   //  Used for naming of objects in the VertexBin class,
2840   //  important to avoid memory problems when running multiple
2841   //  initializations of the class with different setups
2842   //
2843   //  Parameters:
2844   //   flags: Flow flags for type determination
2845   //   prependUS: Prepend an underscore (_) to the name
2846   //
2847   //  Return: QC calculation type
2848   //
2849   TString type = "";
2850   if      ((flags & kStdQC))  type = "StdQC";
2851   else if ((flags & kEtaGap)) type = "EtaGap";
2852   else if ((flags & k3Cor))   type = "3Cor";
2853   else type = "UNKNOWN";
2854   if (prependUS) type.Prepend("_");
2855   if ((flags & kTPC) == kTPC)    type.Append("TPCTr");
2856   if ((flags & kHybrid) == kHybrid) type.Append("HybTr");
2857
2858   return type.Data();
2859 }
2860 //_____________________________________________________________________
2861 TH1* AliForwardFlowTaskQC::CumuHistos::Get(Char_t t, Int_t n, UInt_t nua) const
2862 {
2863   //
2864   //  Get histogram/profile for type t and moment n
2865   //
2866   //  Parameters:
2867   //   t: type = 'r'/'d'
2868   //   n: moment
2869   //   nua: NUA type
2870   //
2871   n = GetPos(n, nua);
2872   if (n < 0) AliFatal(Form("CumuHistos out of range: (%c,%d)", t, n));
2873
2874   TH1* h = 0;
2875   Int_t pos = -1;
2876   if      (t == 'r' || t == 'R') pos = n;
2877   else if (t == 'd' || t == 'D') pos = n;
2878   else if (t == 'a' || t == 'A') pos = 2*n;
2879   else if (t == 'b' || t == 'B') pos = 2*n+1;
2880   else AliFatal(Form("CumuHistos wrong type: %c", t));
2881
2882   if (t == 'r' || t == 'R') {
2883     if (pos < fRefHists->GetEntries()) {
2884       h = (TH1*)fRefHists->At(pos);
2885     }
2886   } else if (pos < fDiffHists->GetEntries()) {
2887     h = (TH1*)fDiffHists->At(pos);
2888   }
2889   if (!h) AliFatal(Form("No hist found in list %c at pos %d", t, pos));
2890
2891   return h;
2892 }
2893 //_____________________________________________________________________
2894 void AliForwardFlowTaskQC::CumuHistos::ConnectList(TString name, TList* l)
2895 {
2896   //
2897   //  Connect an output list with this object
2898   //
2899   //  Parameters:
2900   //   name: base name
2901   //   l: list to keep outputs in
2902   //
2903   TString ref = name;
2904   ref.ReplaceAll("Cumu","CumuRef");
2905   fRefHists = (TList*)l->FindObject(ref.Data());
2906   if (!fRefHists) {
2907     fRefHists = new TList();
2908     fRefHists->SetName(ref.Data());
2909     l->Add(fRefHists);
2910   } else {
2911     // Check that the list correspond to fMaxMoments
2912     if (fRefHists->GetEntries() != (fMaxMoment-1.)*(fNUA+1))
2913       AliError("CumuHistos::ConnectList Wrong number of hists in ref list,"
2914                "you are doing something wrong!");
2915   }
2916   TString diff = name;
2917   diff.ReplaceAll("Cumu","CumuDiff");
2918   fDiffHists = (TList*)l->FindObject(diff.Data());
2919   if (!fDiffHists) {
2920     fDiffHists = new TList();
2921     fDiffHists->SetName(diff.Data());
2922     l->Add(fDiffHists);
2923   } else {
2924     // Check that the list correspond to fMaxMoment
2925     if ((fDiffHists->GetEntries() != (fMaxMoment-1.)*(fNUA+1)) &&
2926         (fDiffHists->GetEntries() != 2*(fMaxMoment-1.)*(fNUA+1)))
2927       AliError(Form("CumuHistos::ConnectList Wrong number of hists in diff list,"
2928                "you are doing something wrong! (%s)",name.Data()));
2929   }
2930
2931 }
2932 //_____________________________________________________________________
2933 void AliForwardFlowTaskQC::CumuHistos::Add(TH1* h) const
2934 {
2935   //
2936   //  Add a histogram to this objects lists
2937   //
2938   //  Parameters:
2939   //   h: histogram/profile to add
2940   //
2941   TString name = h->GetName();
2942   if (name.Contains("Ref")) {
2943     if (fRefHists) fRefHists->Add(h);
2944     else AliFatal("CumuHistos::Add() - fRefHists does not exist");
2945     // Check that the list correspond to fMaxMoments
2946     if (fRefHists->GetEntries() > (fNUA+1)*(fMaxMoment-1.))
2947       AliError("CumuHistos::Add wrong number of hists in ref list, "
2948                "you are doing something wrong!");
2949   }
2950   else if (name.Contains("Diff")) {
2951     if (fDiffHists) fDiffHists->Add(h);
2952     else AliFatal("CumuHistos::Add() - fDiffHists does not exist");
2953     // Check that the list correspond to fMaxMoment
2954     if (fDiffHists->GetEntries() > 2*(fNUA+1)*(fMaxMoment-1.))
2955       AliError("CumuHistos::Add wrong number of hists in diff list, "
2956          "you are doing something wrong!");
2957   }
2958   return;
2959 }
2960 //_____________________________________________________________________
2961 Int_t AliForwardFlowTaskQC::CumuHistos::GetPos(Int_t n, UInt_t nua) const
2962 {
2963   //
2964   //  Get position in list of moment n objects
2965   //  To take care of different numbering schemes
2966   //
2967   //  Parameters:
2968   //   n: moment
2969   //   nua: # of nua corrections applied
2970   //
2971   //  Return: position #
2972   //
2973   if (n > fMaxMoment) return -1;
2974   else return (n-2)+nua*(fMaxMoment-1);
2975 }
2976 //_____________________________________________________________________
2977 //
2978 //
2979 // EOF