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