PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.cxx

   1
   2 //
   3 // Calculate flow in the forward and central regions using the Q cumulants method.
   4 //
   5 // Inputs:
   6 //  - AliAODEvent
   7 //
   8 // Outputs:
   9 //  - AnalysisResults.root or forward_flow.root
  10 //
  11 #include <TROOT.h>
  12 #include <TSystem.h>
  13 #include <TInterpreter.h>
  14 #include <TChain.h>
  15 #include <TFile.h>
  16 #include <TList.h>
  17 #include <TMath.h>
  18 #include <TH3D.h>
  19 #include <TProfile2D.h>
  20 #include <TParameter.h>
  21 #include <TMatrixD.h>
  22 #include <TVectorD.h>
  23 #include <TGraph.h>
  24 #include "AliLog.h"
  25 #include "AliForwardFlowTaskQC.h"
  26 #include "AliAnalysisManager.h"
  27 #include "AliAODHandler.h"
  28 #include "AliAODInputHandler.h"
  29 #include "AliAODForwardMult.h"
  30 #include "AliAODCentralMult.h"
  31 #include "AliAODEvent.h"
  32 #include "AliForwardUtil.h"
  33 #include "AliVVZERO.h"
  34 #include "AliAODVertex.h"
  35 #include "AliCentrality.h"
  36 #include "AliESDEvent.h"
  37 #include "AliVTrack.h"
  38 #include "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(0x0),       // 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|kSPD, fSPDCut, fEtaGap));
 212     }
 213     else if ((fFlowFlags & kVZERO)) {
 214       fBinsForward.Add(new VertexBin(vL, vH, fMaxMoment, "VZERO", fFlowFlags, 0, fEtaGap));
 215       if ((fFlowFlags & kEtaGap)) fBinsCentral.Add(new VertexBin(vL, vH, fMaxMoment, "SPD-VZERO", fFlowFlags|kNUAcorr|kSPD, 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   AliVVZERO* vzero = GetVZERO();
 326   if ((fFlowFlags & kVZERO)) {
 327     if (vzero) {
 328       fHistdNdedpV0.Reset();
 329       FillVZEROHist(vzero);
 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) && !vzero) {
 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 AliVVZERO* AliForwardFlowTaskQC::GetVZERO() const
 864 {
 865   AliVVZERO* vzero = 0;
 866   // Get input type
 867   UShort_t input = AliForwardUtil::CheckForAOD();
 868   switch (input) {
 869     // If AOD input, simply get the track array from the event
 870     case 1: vzero = (AliVVZERO*)fAOD->GetVZEROData();
 871             break;
 872     case 2: {
 873     // If ESD input get event, apply track cuts
 874               AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 875               if (!esd) return 0;
 876               vzero = (AliVVZERO*)esd->GetVZEROData();
 877               break;
 878             }
 879     default: AliFatal("Neither ESD or AOD input. This should never happen");
 880             break;
 881   }
 882   return vzero;
 883 }
 884 // _____________________________________________________________________
 885 void AliForwardFlowTaskQC::FillVZEROHist(AliVVZERO* vzero)
 886 {
 887   //
 888   //  Loops over VZERO data object and fill up d^2N/detadphi histogram for flow analysis
 889   //
 890   //  Parameters:
 891   //   vzero: VZERO AOD data object
 892   //
 893   Int_t ring = 0;
 894   Int_t bin = 0;
 895   Double_t eta = 0;
 896   // Sort of right for 2010 data, do not use for 2011!
 897   Double_t eq[64] = { 1.43536, 1.45727, 1.44993, 1.30051, 1.17425, 1.2335, 1.22247, 1.14362,
 898                       1.14647, 1.25208, 1.17681, 1.21642, 1.16604, 1.05532, 1.03212, 1.1032,
 899                       1.22941, 1.36986, 1.14652, 1.20056, 0.927086, 1.10809, 1.03343, 1.29472,
 900                       1.21204, 1.29217, 1.2003, 2.10382, 1.28513, 1.40558, 1.25784, 1.21848,
 901                       0.475162, 0.50421, 0.503617, 0.512471, 0.515276, 0.39831, 0.415199, 0.444664,
 902                       0.521922, 0.785915, 0.703658, 0.832479, 0.77461, 0.73129, 0.778697, 0.710265,
 903                       0.89686, 0.967688, 0.974225, 0.873445, 0.811096, 0.828493, 0.889609, 0.586056,
 904                       1.15877, 0.954656, 0.914557, 0.979028, 1.04907, 0.748518, 0.928043, 0.98175 };
 905   for (Int_t i = 0; i < 64; i++) {
 906     if (i % 8 == 0) {
 907       ring++;
 908       bin = (ring < 5 ? ring+1 : 15-ring);
 909       eta = fHistdNdedpV0.GetXaxis()->GetBinCenter(bin);
 910       fHistdNdedpV0.SetBinContent(bin, 0, 1);
 911     }
 912     Float_t amp = vzero->GetMultiplicity(i);
 913     amp /= eq[i];
 914     Double_t phi = TMath::Pi()/8.+TMath::TwoPi()*i/8.;
 915     while (phi > TMath::TwoPi()) phi -= TMath::TwoPi();
 916     fHistdNdedpV0.Fill(eta, phi, amp);
 917   }
 918 }
 919 //_____________________________________________________________________
 920 AliForwardFlowTaskQC::VertexBin::VertexBin()
 921   : TNamed(),
 922     fMaxMoment(0),   // Max flow moment for this vertexbin
 923     fVzMin(0),       // Vertex z-coordinate min [cm]
 924     fVzMax(0),       // Vertex z-coordinate max [cm]
 925     fType(),         // Data type name e.g., FMD/SPD/FMDTR/SPDTR/MC
 926     fFlags(0),       // Flow flags
 927     fSigmaCut(-1),   // Sigma cut to remove outlier events
 928     fEtaGap(-1),     // Eta gap value
 929     fEtaLims(),      // Limits for binning in 3Cor method
 930     fCumuRef(),      // Histogram for reference flow
 931     fCumuDiff(),     // Histogram for differential flow
 932     fCumuHists(0,0), // CumuHists object for keeping track of results
 933     fCumuNUARef(),   // Histogram for ref NUA terms
 934     fCumuNUADiff(),  // Histogram for diff NUA terms
 935     fdNdedpRefAcc(), // Diagnostics histogram for acc. maps
 936     fdNdedpDiffAcc(),// Diagnostics histogram for acc. maps
 937     fOutliers(),     // Histogram for sigma distribution
 938     fDebug()         // Debug level
 939 {
 940   //
 941   //  Default constructor
 942   //
 943 }
 944 //_____________________________________________________________________
 945 AliForwardFlowTaskQC::VertexBin::VertexBin(Int_t vLow, Int_t vHigh,
 946                                            UShort_t moment, TString name,
 947                                            UShort_t flags, Double_t cut,
 948                                            Double_t etaGap)
 949   : TNamed("", ""),
 950     fMaxMoment(moment),  // Max flow moment for this vertexbin
 951     fVzMin(vLow),        // Vertex z-coordinate min [cm]
 952     fVzMax(vHigh),       // Vertex z-coordinate max [cm]
 953     fType(name),         // Data type name e.g., FMD/SPD/FMDTR/SPDTR/MC
 954     fFlags(flags),       // Flow flags
 955     fSigmaCut(cut),      // Sigma cut to remove outlier events
 956     fEtaGap(etaGap),     // Eta gap value
 957     fEtaLims(),          // Limits for binning in 3Cor method
 958     fCumuRef(),          // Histogram for reference flow
 959     fCumuDiff(),         // Histogram for differential flow
 960     fCumuHists(moment,0),// CumuHists object for keeping track of results
 961     fCumuNUARef(),       // Histogram for ref NUA terms
 962     fCumuNUADiff(),      // Histogram for diff NUA terms
 963     fdNdedpRefAcc(),     // Diagnostics histogram for acc. maps
 964     fdNdedpDiffAcc(),    // Diagnostics histogram for acc. maps
 965     fOutliers(),         // Histogram for sigma distribution
 966     fDebug(0)            // Debug level
 967 {
 968   //
 969   //  Constructor
 970   //
 971   //  Parameters
 972   //   vLow: min z-coordinate
 973   //   vHigh: max z-coordinate
 974   //   moment: max flow moment
 975   //   name: data type name (FMD/SPD/FMDTR/SPDTR/MC)
 976   //   flags: flow flags
 977   //   cut: sigma cut
 978   //   etaGap: eta-gap value
 979   //
 980   fType.ToUpper();
 981
 982   SetName(Form("%svertexBin%d_%d_%d%s", fType.Data(), moment, vLow, vHigh, GetQCType(fFlags)));
 983   SetTitle(Form("%svertexBin%d_%d_%d%s", fType.Data(), moment, vLow, vHigh, GetQCType(fFlags)));
 984
 985   fDebug = AliAnalysisManager::GetAnalysisManager()->GetDebugLevel();
 986   if (fDebug > 0) Printf("AliForwardFlowTaskQC::VertexBin()\tDebugMode: %d", fDebug);
 987
 988   // Set limits for 3 correlator method
 989   if ((fFlags & kFMD)) {
 990     fEtaLims[0] = -6.;
 991     fEtaLims[1] = -1.5;
 992     fEtaLims[2] = -0.5;
 993     fEtaLims[3] = 2.;
 994     fEtaLims[4] = 3.;
 995     fEtaLims[5] = 6.;
 996   } else if ((fFlags & kVZERO)) {
 997     fEtaLims[0] = -6;
 998     fEtaLims[1] = -2.7;
 999     fEtaLims[2] = -2.0;
1000     fEtaLims[3] = 2.0;
1001     fEtaLims[4] = 3.9;
1002     fEtaLims[5] = 6;
1003   }
1004 }
1005 //_____________________________________________________________________
1006 AliForwardFlowTaskQC::VertexBin&
1007 AliForwardFlowTaskQC::VertexBin::operator=(const AliForwardFlowTaskQC::VertexBin& o)
1008 {
1009   //
1010   //  Assignment operator
1011   //
1012   //  Parameters
1013   //   o: AliForwardFlowTaskQC::VertexBin
1014   //
1015   if (&o == this) return *this;
1016   fMaxMoment     = o.fMaxMoment;
1017   fVzMin         = o.fVzMin;
1018   fVzMax         = o.fVzMax;
1019   fType          = o.fType;
1020   fFlags         = o.fFlags;
1021   fSigmaCut      = o.fSigmaCut;
1022   fEtaGap        = o.fEtaGap;
1023   fCumuRef       = o.fCumuRef;
1024   fCumuDiff      = o.fCumuDiff;
1025   fCumuHists     = o.fCumuHists;
1026   fCumuNUARef    = o.fCumuNUARef;
1027   fCumuNUADiff   = o.fCumuNUADiff;
1028   fdNdedpRefAcc  = o.fdNdedpRefAcc;
1029   fdNdedpDiffAcc = o.fdNdedpDiffAcc;
1030   fOutliers      = o.fOutliers;
1031   fDebug         = o.fDebug;
1032   for (UInt_t i = 0; i < sizeof(fEtaLims)/sizeof(Double_t); i++) fEtaLims[i] = o.fEtaLims[i];
1033
1034   return *this;
1035 }
1036 //_____________________________________________________________________
1037 void AliForwardFlowTaskQC::VertexBin::AddOutput(TList* outputlist, TAxis* centAxis)
1038 {
1039   //
1040   //  Add histograms to outputlist
1041   //
1042   //  Parameters
1043   //   outputlist: list of histograms
1044   //   centAxis: centrality axis
1045   //
1046
1047   // First we try to find an outputlist for this vertexbin
1048   TList* list = (TList*)outputlist->FindObject(Form("%svertex_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1049   // If it doesn't exist we make one
1050   if (!list) {
1051     list = new TList();
1052     list->SetName(Form("%svertex_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1053     outputlist->Add(list);
1054   }
1055
1056   // Get bin numbers and binning defined
1057   Int_t nHBins = GetBinNumberSin();
1058   Int_t nEtaBins = 48;
1059   if ((fFlags & k3Cor)) {
1060     if ((fFlags & kFMD)) nEtaBins = 24;
1061     else if ((fFlags & kVZERO)) nEtaBins = 19;
1062   }
1063   else if ((fFlags & kVZERO) && (fFlags & kEtaGap)) nEtaBins = 19;
1064   else if ((fFlags & kVZERO)) nEtaBins = 11;
1065
1066   Double_t vzeroBins[12] = { -6, -3.7, -3.2, -2.7, -2.2, -1.7,
1067                              2.8, 3.4,  3.9,  4.5,  5.1, 6 };
1068   Double_t vzeroBins2[20] = { -6, -3.7, -3.2, -2.7, -2.2, // VZERO
1069                              -2.0, -1.5, -1.0, -0.5 , 0., 0.5, 1.0, 1.5, 2.0, // SPD
1070                             2.8, 3.4,  3.9,  4.5,  5.1, 6 }; // VZERO
1071
1072   Int_t nRefBins = nEtaBins; // needs to be something as default
1073   if ((fFlags & kStdQC)) {
1074     if ((fFlags & kSymEta) && !((fFlags & kTPC) && (fFlags & kSPD))) nRefBins = 1;
1075     else nRefBins = 2;
1076   } else if ((fFlags & kEtaGap )) {
1077     nRefBins = 2;
1078   } else if ((fFlags & k3Cor)) {
1079     nRefBins = 24;
1080   }
1081
1082   // We initiate the reference histogram
1083   fCumuRef = new TH2D(Form("%s_%d_%d%s_ref", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1084                       Form("%s_%d_%d%s_ref", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1085                         nRefBins, -6., 6.,
1086                         nHBins, 0.5, nHBins+0.5);
1087   if ((fFlags & kVZERO) && (fFlags & k3Cor)) fCumuRef->GetXaxis()->Set(nEtaBins, vzeroBins2);
1088   SetupNUALabels(fCumuRef->GetYaxis());
1089   list->Add(fCumuRef);
1090   // We initiate the differential histogram
1091   fCumuDiff = new TH2D(Form("%s_%d_%d%s_diff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1092                        Form("%s_%d_%d%s_diff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1093                        nEtaBins, -6., 6., nHBins, 0.5, nHBins+0.5);
1094   if ((fFlags & kVZERO)) {
1095     if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1096       fCumuDiff->GetXaxis()->Set(nEtaBins, vzeroBins2);
1097     else fCumuDiff->GetXaxis()->Set(nEtaBins, vzeroBins);
1098   }
1099   SetupNUALabels(fCumuDiff->GetYaxis());
1100   list->Add(fCumuDiff);
1101
1102   // Cumulants sum hists
1103   Int_t cBins = centAxis->GetNbins();
1104   fCumuHists.ConnectList(Form("%sCumu_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)), list);
1105   TH3D* cumuHist = 0;
1106   Int_t nC2Bins = ((fFlags & kEtaGap) || (fFlags & k3Cor) ? kW2 : k3pWeight);
1107   Int_t nC4Bins = ((fFlags & kEtaGap) ? kW2 : ((fFlags & k3Cor) ? kW4 : kSinphi1phi2phi3p));
1108   for (Int_t n = 2; n <= fMaxMoment; n++) {
1109     // Initiate the ref cumulant sum histogram
1110     cumuHist = new TH3D(Form("%sv%d_vertex_%d_%d%s_cumuRef", fType.Data(), n, fVzMin, fVzMax, GetQCType(fFlags)),
1111                         Form("%sv%d_vertex_%d_%d%s_cumuRef", fType.Data(), n, fVzMin, fVzMax, GetQCType(fFlags)),
1112                         nRefBins, -6., 6.,
1113                         cBins, 0., 100., nC2Bins, 0.5, nC2Bins+0.5);
1114     if ((fFlags & kVZERO) && (fFlags & k3Cor)) cumuHist->GetXaxis()->Set(nEtaBins, vzeroBins2);
1115     cumuHist->GetYaxis()->Set(cBins, centAxis->GetXbins()->GetArray());
1116     fCumuHists.Add(cumuHist);
1117     // Initiate the diff cumulant sum histogram
1118     cumuHist = new TH3D(Form("%sv%d_vertex_%d_%d%s_cumuDiff", fType.Data(), n, fVzMin, fVzMax, GetQCType(fFlags)),
1119                         Form("%sv%d_vertex_%d_%d%s_cumuDiff", fType.Data(), n, fVzMin, fVzMax, GetQCType(fFlags)),
1120                         nEtaBins, -6., 6., cBins, 0., 100., nC4Bins, 0.5, nC4Bins+0.5);
1121     if ((fFlags & kVZERO)) {
1122       if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1123         cumuHist->GetXaxis()->Set(nEtaBins, vzeroBins2);
1124       else cumuHist->GetXaxis()->Set(nEtaBins, vzeroBins);
1125     }
1126     cumuHist->GetYaxis()->Set(cBins, centAxis->GetXbins()->GetArray());
1127     fCumuHists.Add(cumuHist);
1128   }
1129
1130   // Common NUA histograms
1131   fCumuNUARef = new TH3D(Form("%s_vertex_%d_%d%s_cumuNUARef", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1132                          Form("%s_vertex_%d_%d%s_cumuNUARef", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1133                          nRefBins, -6., 6.,
1134                           cBins, 0., 100., nHBins, 0.5, nHBins+0.5);
1135   if ((fFlags & kVZERO) && (fFlags & k3Cor)) fCumuNUARef->GetXaxis()->Set(nEtaBins, vzeroBins2);
1136   fCumuNUARef->GetYaxis()->Set(cBins, centAxis->GetXbins()->GetArray());
1137   SetupNUALabels(fCumuNUARef->GetZaxis());
1138   fCumuNUARef->Sumw2();
1139   list->Add(fCumuNUARef);
1140
1141   fCumuNUADiff = new TH3D(Form("%s_vertex_%d_%d%s_cumuNUADiff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1142                           Form("%s_vertex_%d_%d%s_cumuNUADiff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1143                           nEtaBins, -6., 6., cBins, 0., 100., nHBins, 0.5, nHBins+0.5);
1144   if ((fFlags & kVZERO)) {
1145     if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1146       fCumuNUADiff->GetXaxis()->Set(nEtaBins, vzeroBins2);
1147     else fCumuNUADiff->GetXaxis()->Set(nEtaBins, vzeroBins);
1148   }
1149   fCumuNUADiff->GetYaxis()->Set(cBins, centAxis->GetXbins()->GetArray());
1150   SetupNUALabels(fCumuNUADiff->GetZaxis());
1151   fCumuNUADiff->Sumw2();
1152   list->Add(fCumuNUADiff);
1153
1154   // We create diagnostic histograms.
1155   TList* dList = (TList*)outputlist->FindObject("Diagnostics");
1156   if (!dList) AliFatal("No diagnostics list found");
1157
1158   // Acceptance hist
1159   Int_t nPhiBins = ((fFlags & kFMD) ? 20 : 8);
1160   fdNdedpRefAcc = new TH2F(Form("h%sdNdedpRefAcc_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1161     Form("%s reference flow acceptance map for %d cm < v_{z} < %d cm", fType.Data(), fVzMin, fVzMax),
1162     nEtaBins, -6., 6., nPhiBins, 0, TMath::TwoPi());
1163   if ((fFlags & kVZERO)) {
1164     if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1165       fdNdedpRefAcc->GetXaxis()->Set(nEtaBins, vzeroBins2);
1166     else fdNdedpRefAcc->GetXaxis()->Set(nEtaBins, vzeroBins);
1167   }
1168   dList->Add(fdNdedpRefAcc);
1169
1170   fdNdedpDiffAcc = new TH2F(Form("h%sdNdedpDiffAcc_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1171     Form("%s differential flow acceptance map for %d cm < v_{z} < %d cm", fType.Data(), fVzMin, fVzMax),
1172     nEtaBins, -6., 6., nPhiBins, 0, TMath::TwoPi());
1173   if ((fFlags & kVZERO)) {
1174     if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1175       fdNdedpDiffAcc->GetXaxis()->Set(nEtaBins, vzeroBins2);
1176     else fdNdedpDiffAcc->GetXaxis()->Set(nEtaBins, vzeroBins);
1177   }
1178   dList->Add(fdNdedpDiffAcc);
1179
1180   fOutliers = new TH2F(Form("hOutliers_%s_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1181                        Form("Maximum #sigma from mean N_{ch} pr. bin - %s, %d < v_{z} < %d",
1182                        fType.Data(), fVzMin, fVzMax),
1183                        20, 0., 100., 500, 0., ((fFlags & kMC) ? 15. : 5.));
1184   dList->Add(fOutliers);
1185
1186   return;
1187 }
1188 //_____________________________________________________________________
1189 Bool_t AliForwardFlowTaskQC::VertexBin::FillHists(TH2D& dNdetadphi, Double_t cent, UShort_t mode)
1190 {
1191   //
1192   //  Fill reference and differential eta-histograms
1193   //
1194   //  Parameters:
1195   //   dNdetadphi: 2D histogram with input data
1196   //   cent: centrality
1197   //   mode: filling mode: kFillRef/kFillDiff/kFillBoth
1198   //
1199   if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1200   Bool_t useEvent = kTRUE;
1201
1202   // Fist we reset histograms
1203   if ((mode & kReset)) {
1204     if ((mode & kFillRef))  fCumuRef->Reset();
1205     if ((mode & kFillDiff)) fCumuDiff->Reset();
1206   }
1207
1208   // Then we loop over the input and calculate sum cos(k*n*phi)
1209   // and fill it in the reference and differential histograms
1210   Int_t nBadBins = 0;
1211   Double_t limit = 9999.;
1212   for (Int_t etaBin = 1; etaBin <= dNdetadphi.GetNbinsX(); etaBin++) {
1213     Double_t eta = dNdetadphi.GetXaxis()->GetBinCenter(etaBin);
1214     // Numbers to cut away bad events
1215     Double_t runAvg = 0;
1216     Double_t max = 0;
1217     Int_t nInAvg = 0;
1218     Double_t avgSqr = 0;
1219     for (Int_t phiBin = 0; phiBin <= dNdetadphi.GetNbinsY(); phiBin++) {
1220       // Check for acceptance
1221       if (phiBin == 0) {
1222         if (dNdetadphi.GetBinContent(etaBin, 0) == 0) break;
1223         // Central limit for eta gap break for reference flow
1224         if ((fFlags & kEtaGap) && (mode & kFillRef) &&
1225              TMath::Abs(eta) < fEtaGap) break;
1226         // Backward and forward eta gap break for reference flow
1227         if ((fFlags & kEtaGap) && (mode & kFillRef) && TMath::Abs(eta) > TMath::Abs(limit)) break;
1228         if ((fFlags & kStdQC) && (fFlags & kMC)) {
1229           if (!(fFlags & kSPD) && TMath::Abs(eta) < 1.75) break;
1230           if ((fFlags & kSPD) && TMath::Abs(eta) > 2.00) break;
1231         }
1232         if (limit > 1e3) limit = dNdetadphi.GetXaxis()->GetBinLowEdge(etaBin);
1233         continue;
1234       } // End of phiBin == 0
1235       Double_t phi = dNdetadphi.GetYaxis()->GetBinCenter(phiBin);
1236       Double_t weight = dNdetadphi.GetBinContent(etaBin, phiBin);
1237
1238       // We calculate the average Nch per. bin
1239       avgSqr += weight*weight;
1240       runAvg += weight;
1241       nInAvg++;
1242       if (weight == 0) continue;
1243       if (weight > max) max = weight;
1244
1245       // Fill into Cos() and Sin() hists
1246       if ((mode & kFillRef)) {
1247         fCumuRef->Fill(eta, 0., weight);// mult goes in underflowbin - no visual, but not needed?
1248         fdNdedpRefAcc->Fill(eta, phi, weight);
1249       }
1250       if ((mode & kFillDiff)) {
1251         fCumuDiff->Fill(eta, 0., weight);
1252         fdNdedpDiffAcc->Fill(eta, phi, weight);
1253       }
1254       for (Int_t n = 1; n <= 2*fMaxMoment; n++) {
1255         Double_t cosBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberCos(n));
1256         Double_t sinBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberSin(n));
1257         Double_t cosnPhi = weight*TMath::Cos(n*phi);
1258         Double_t sinnPhi = weight*TMath::Sin(n*phi);
1259         // fill ref
1260         if ((mode & kFillRef)) {
1261           fCumuRef->Fill(eta, cosBin, cosnPhi);
1262           fCumuRef->Fill(eta, sinBin, sinnPhi);
1263         }
1264         // fill diff
1265         if ((mode & kFillDiff)) {
1266           fCumuDiff->Fill(eta, cosBin, cosnPhi);
1267           fCumuDiff->Fill(eta, sinBin, sinnPhi);
1268         }
1269       } // End of NUA loop
1270     } // End of phi loop
1271     // Outlier cut calculations
1272     if (nInAvg > 0) {
1273       runAvg /= nInAvg;
1274       avgSqr /= nInAvg;
1275       Double_t stdev = (nInAvg > 1 ? TMath::Sqrt(nInAvg/(nInAvg-1))*TMath::Sqrt(avgSqr - runAvg*runAvg) : 0);
1276       Double_t nSigma = (stdev == 0 ? 0 : (max-runAvg)/stdev);
1277       if (fSigmaCut > 0. && nSigma >= fSigmaCut && cent < 60) nBadBins++;
1278       else nBadBins = 0;
1279       fOutliers->Fill(cent, nSigma);
1280       // We still finish the loop, for fOutliers to make sense,
1281       // but we do no keep the event for analysis
1282       if (nBadBins > 3) useEvent = kFALSE;
1283     }
1284   } // End of eta bin
1285
1286   return useEvent;
1287 }
1288 //_____________________________________________________________________
1289 Bool_t AliForwardFlowTaskQC::VertexBin::FillTracks(TObjArray* trList, UShort_t mode)
1290 {
1291   //
1292   //  Fill reference and differential eta-histograms
1293   //
1294   //  Parameters:
1295   //   trList: list of tracks
1296   //   mode: filling mode: kFillRef/kFillDiff/kFillBoth
1297   //
1298   if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1299
1300   // Fist we reset histograms
1301   if ((mode & kReset)) {
1302     if ((mode & kFillRef))  fCumuRef->Reset();
1303     if ((mode & kFillDiff)) fCumuDiff->Reset();
1304   }
1305
1306   UShort_t input = AliForwardUtil::CheckForAOD();
1307   // Then we loop over the input and calculate sum cos(k*n*phi)
1308   // and fill it in the reference and differential histograms
1309   Int_t nTr = trList->GetEntries();
1310   if (nTr == 0) return kFALSE;
1311   AliVTrack* tr = 0;
1312   AliAODTrack* aodTr = 0;
1313   // Cuts for AOD tracks (have already been applied to ESD tracks)
1314   const Double_t pTMin = 0.5, pTMax = 20., etaMin = -0.8, etaMax = 0.8, minNCl = 50;
1315   for (Int_t i = 0; i < nTr; i++) { // track loop
1316     tr = (AliVTrack*)trList->At(i);
1317     if (!tr) continue;
1318     if (input == 1) { // If AOD input
1319       // A dynamic cast would be more safe here, but this is faster...
1320       aodTr = (AliAODTrack*)tr;
1321       if (aodTr->GetID() > -1) continue;
1322       if (!aodTr->TestFilterBit(128) || !aodTr->Pt() > pTMax || aodTr->Pt() < pTMin ||
1323         aodTr->Eta() > etaMax || aodTr->Eta() < etaMin || aodTr->GetTPCNcls() < minNCl) continue;
1324     }
1325     // Track accepted
1326     Double_t eta = tr->Eta();
1327 //    if ((fFlags & kSPD) && TMath::Abs(eta) < 0.4) continue;
1328     Double_t phi = tr->Phi();
1329     if ((mode & kFillRef)) {
1330       fCumuRef->Fill(eta, 0.);// mult goes in underflowbin - no visual, but not needed?
1331       fdNdedpRefAcc->Fill(eta, phi);
1332     }
1333     if ((mode & kFillDiff)) {
1334       fCumuDiff->Fill(eta, 0);
1335       fdNdedpDiffAcc->Fill(eta, phi);
1336     }
1337     for (Int_t n = 1; n <= 2*fMaxMoment; n++) {
1338       Double_t cosBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberCos(n));
1339       Double_t sinBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberSin(n));
1340       Double_t cosnPhi = TMath::Cos(n*phi);
1341       Double_t sinnPhi = TMath::Sin(n*phi);
1342       // fill ref
1343       if ((mode & kFillRef)) {
1344         fCumuRef->Fill(eta, cosBin, cosnPhi);
1345         fCumuRef->Fill(eta, sinBin, sinnPhi);
1346       }
1347       // fill diff
1348       if ((mode & kFillDiff)) {
1349         fCumuDiff->Fill(eta, cosBin, cosnPhi);
1350         fCumuDiff->Fill(eta, sinBin, sinnPhi);
1351       }
1352     } // End of NUA loop
1353   } // End of track loop
1354   return kTRUE;
1355 }
1356 //_____________________________________________________________________
1357 void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate(Double_t cent)
1358 {
1359   //
1360   //  Calculate the Q cumulant up to order fMaxMoment
1361   //
1362   //  Parameters:
1363   //   cent: centrality of event
1364   //
1365   if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1366
1367   // Fill out NUA hists
1368   for (Int_t etaBin = 1; etaBin <= fCumuRef->GetNbinsX(); etaBin++) {
1369     Double_t eta = fCumuRef->GetXaxis()->GetBinCenter(etaBin);
1370     if (fCumuRef->GetBinContent(etaBin, 0) == 0) continue;
1371     if ((fFlags & kTPC) && (fFlags && kSPD)) eta = -eta;
1372     for (Int_t qBin = 0; qBin <= fCumuRef->GetNbinsY(); qBin++) {
1373       fCumuNUARef->Fill(eta, cent, Double_t(qBin), fCumuRef->GetBinContent(etaBin, qBin));
1374     }
1375   }
1376   for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1377     Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1378     if (fCumuDiff->GetBinContent(etaBin, 0) == 0) continue;
1379     for (Int_t qBin = 0; qBin <= fCumuDiff->GetNbinsY(); qBin++) {
1380       fCumuNUADiff->Fill(eta, cent, Double_t(qBin), fCumuDiff->GetBinContent(etaBin, qBin));
1381     }
1382   }
1383
1384   // We create the objects needed for the analysis
1385   TH3D* cumuRef = 0;
1386   TH3D* cumuDiff = 0;
1387   // For each n we loop over the hists
1388   for (Int_t n = 2; n <= fMaxMoment; n++) {
1389     cumuRef  = (TH3D*)fCumuHists.Get('r',n);
1390     cumuDiff = (TH3D*)fCumuHists.Get('d',n);
1391     Int_t prevRefEtaBin = 0;
1392
1393     // Per mom. quantities
1394     Double_t dQnReA = 0, dQnImA = 0, multA = 0;
1395     Double_t dQnReB = 0, dQnImB = 0, multB = 0;
1396     Double_t dQ2nReA = 0, dQ2nImA = 0;
1397     Double_t two = 0, w2 = 0, four = 0, w4 = 0;
1398     Double_t cosPhi1Phi2 = 0, cosPhi1Phi2Phi3m = 0;
1399     Double_t sinPhi1Phi2 = 0, sinPhi1Phi2Phi3m = 0;
1400     for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1401       Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1402       Double_t refEta = eta;
1403       if ((fFlags & kTPC) && (fFlags && kSPD)) refEta = -eta;
1404       Int_t refEtaBinA = fCumuRef->GetXaxis()->FindBin(refEta);
1405       if ((fFlags & kEtaGap)) refEta = -eta;
1406       Int_t refEtaBinB = fCumuRef->GetXaxis()->FindBin(refEta);
1407       if (refEtaBinA != prevRefEtaBin) {
1408         prevRefEtaBin = refEtaBinA;
1409         // Reference flow
1410         multA  = fCumuRef->GetBinContent(refEtaBinA, 0);
1411         dQnReA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberCos(n));
1412         dQnImA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberSin(n));
1413         dQ2nReA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberCos(n*2));
1414         dQ2nImA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberSin(n*2));
1415
1416         multB  = fCumuRef->GetBinContent(refEtaBinB, 0);
1417         dQnReB = fCumuRef->GetBinContent(refEtaBinB, GetBinNumberCos(n));
1418         dQnImB = fCumuRef->GetBinContent(refEtaBinB, GetBinNumberSin(n));
1419
1420         if (multA <= 3 || multB <= 3) return;
1421         // The reference flow is calculated
1422         // 2-particle
1423         if ((fFlags & kStdQC)) {
1424           w2 = multA * (multA - 1.);
1425           two = dQnReA*dQnReA + dQnImA*dQnImA - multA;
1426         } else {
1427           w2 = multA * multB;
1428           two = dQnReA*dQnReB + dQnImA*dQnImB;
1429         }
1430         cumuRef->Fill(eta, cent, kW2Two, two);
1431         cumuRef->Fill(eta, cent, kW2, w2);
1432
1433         // The reference flow is calculated
1434         // 4-particle
1435         if ((fFlags & kStdQC)) {
1436           w4 = multA * (multA - 1.) * (multA - 2.) * (multA - 3.);
1437
1438           four = 2.*multA*(multA-3.) + TMath::Power((TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.)),2.)
1439                    -4.*(multA-2.)*(TMath::Power(dQnReA,2.) + TMath::Power(dQnImA,2.))
1440                    -2.*(TMath::Power(dQnReA,2.)*dQ2nReA+2.*dQnReA*dQnImA*dQ2nImA-TMath::Power(dQnImA,2.)*dQ2nReA)
1441                    +(TMath::Power(dQ2nReA,2.)+TMath::Power(dQ2nImA,2.));
1442
1443           cumuRef->Fill(eta, cent, kW4Four, four);
1444           cumuRef->Fill(eta, cent, kW4, w4);
1445
1446           // NUA
1447           cosPhi1Phi2 = dQnReA*dQnReA - dQnImA*dQnImA - dQ2nReA;
1448           sinPhi1Phi2 = 2.*dQnReA*dQnImA - dQ2nImA;
1449
1450           cosPhi1Phi2Phi3m = dQnReA*(TMath::Power(dQnReA,2)+TMath::Power(dQnImA,2))
1451                               -dQnReA*dQ2nReA-dQnImA*dQ2nImA-2.*(multA-1)*dQnReA;
1452
1453           sinPhi1Phi2Phi3m = -dQnImA*(TMath::Power(dQnReA,2)+TMath::Power(dQnImA,2))
1454                               +dQnReA*dQ2nImA-dQnImA*dQ2nReA+2.*(multA-1)*dQnImA;
1455
1456           cumuRef->Fill(eta, cent, kCosphi1phi2, cosPhi1Phi2);
1457           cumuRef->Fill(eta, cent, kSinphi1phi2, sinPhi1Phi2);
1458           cumuRef->Fill(eta, cent, kCosphi1phi2phi3m, cosPhi1Phi2Phi3m);
1459           cumuRef->Fill(eta, cent, kSinphi1phi2phi3m, sinPhi1Phi2Phi3m);
1460           cumuRef->Fill(eta, cent, k3pWeight, multA*(multA-1.)*(multA-2.));
1461         } // End of QC{4}
1462       } // End of reference flow
1463       // For each etaBin bin the necessary values for differential flow is calculated
1464       Double_t mp = fCumuDiff->GetBinContent(etaBin, 0);
1465       Double_t pnRe = fCumuDiff->GetBinContent(etaBin, GetBinNumberCos(n));
1466       Double_t pnIm = fCumuDiff->GetBinContent(etaBin, GetBinNumberSin(n));
1467       Double_t p2nRe = fCumuDiff->GetBinContent(etaBin, GetBinNumberCos(n*2));
1468       Double_t p2nIm = fCumuDiff->GetBinContent(etaBin, GetBinNumberSin(n*2));
1469       if (mp == 0) continue;
1470       Double_t mq = 0;
1471       Double_t qnRe = 0;
1472       Double_t qnIm = 0;
1473       Double_t q2nRe = 0;
1474       Double_t q2nIm = 0;
1475
1476       // Differential flow calculations for each eta bin is done:
1477       // 2-particle differential flow
1478       if ((fFlags & kStdQC) && !(fFlags & kTPC)) {
1479         mq = mp;
1480         qnRe = pnRe;
1481         qnIm = pnIm;
1482         q2nRe = p2nRe;
1483         q2nIm = p2nIm;
1484       }
1485
1486       Double_t w2p = mp * multB - mq;
1487       Double_t twoPrime = pnRe*dQnReB + pnIm*dQnImB - mq;
1488
1489       cumuDiff->Fill(eta, cent, kW2Two, twoPrime);
1490       cumuDiff->Fill(eta, cent, kW2, w2p);
1491
1492       if ((fFlags & kEtaGap)) continue;
1493       // Differential flow calculations for each eta bin bin is done:
1494       // 4-particle differential flow
1495       Double_t w4p = (mp * multA - 3.*mq)*(multA - 1.)*(multA - 2.);
1496
1497       Double_t fourPrime = (TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.))*(pnRe*dQnReA+pnIm*dQnImA)
1498                           - q2nRe*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))
1499                           - 2.*q2nIm*dQnReA*dQnImA
1500                           - pnRe*(dQnReA*dQ2nReA+dQnImA*dQ2nImA)
1501                           + pnIm*(dQnImA*dQ2nReA-dQnReA*dQ2nImA)
1502                           - 2.*multA*(pnRe*dQnReA+pnIm*dQnImA)
1503                           - 2.*(TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.))*mq
1504                           + 6.*(qnRe*dQnReA+qnIm*dQnImA)
1505                           + 1.*(q2nRe*dQ2nReA+q2nIm*dQ2nImA)
1506                           + 2.*(pnRe*dQnReA+pnIm*dQnImA)
1507                           + 2.*mq*multA
1508                           - 6.*mq;
1509
1510       cumuDiff->Fill(eta, cent, kW4Four, fourPrime);
1511       cumuDiff->Fill(eta, cent, kW4, w4p);
1512
1513       // NUA
1514       Double_t cosPsi1Phi2 = pnRe*dQnReA - pnIm*dQnImA - q2nRe;
1515       Double_t sinPsi1Phi2 = pnRe*dQnImA + pnIm*dQnReA - q2nIm;
1516
1517       Double_t cosPsi1Phi2Phi3p = pnRe*(TMath::Power(dQnImA,2.)+TMath::Power(dQnReA,2.)-multA)
1518                             - 1.*(q2nRe*dQnReA+q2nIm*dQnImA)
1519                             - mq*dQnReA+2.*qnRe;
1520
1521       Double_t sinPsi1Phi2Phi3p = pnIm*(TMath::Power(dQnImA,2.)+TMath::Power(dQnReA,2.)-multA)
1522                             - 1.*(q2nIm*dQnReA-q2nRe*dQnImA)
1523                             - mq*dQnImA+2.*qnIm;
1524
1525       Double_t cosPsi1Phi2Phi3m = pnRe*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))+2.*pnIm*dQnReA*dQnImA
1526                             - 1.*(pnRe*dQ2nReA+pnIm*dQ2nImA)
1527                             - 2.*mq*dQnReA+2.*qnRe;
1528
1529       Double_t sinPsi1Phi2Phi3m = pnIm*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))-2.*pnRe*dQnReA*dQnImA
1530                             - 1.*(pnIm*dQ2nReA-pnRe*dQ2nImA)
1531                             + 2.*mq*dQnImA-2.*qnIm;
1532
1533       cumuDiff->Fill(eta, cent, kCosphi1phi2, cosPsi1Phi2);
1534       cumuDiff->Fill(eta, cent, kSinphi1phi2, sinPsi1Phi2);
1535       cumuDiff->Fill(eta, cent, kCosphi1phi2phi3m, cosPsi1Phi2Phi3m);
1536       cumuDiff->Fill(eta, cent, kSinphi1phi2phi3m, sinPsi1Phi2Phi3m);
1537       cumuDiff->Fill(eta, cent, k3pWeight, (mp*multA-2.*mq)*(multA-1.));
1538       cumuDiff->Fill(eta, cent, kCosphi1phi2phi3p, cosPsi1Phi2Phi3p);
1539       cumuDiff->Fill(eta, cent, kSinphi1phi2phi3p, sinPsi1Phi2Phi3p);
1540     } // End of eta loop
1541     // Event count
1542     cumuRef->Fill(-7., cent, -0.5, 1.);
1543   } // End of moment loop
1544   return;
1545 }
1546 //_____________________________________________________________________
1547 void AliForwardFlowTaskQC::VertexBin::GetLimits(Int_t bin, Int_t& aLow, Int_t& aHigh,
1548                                                 Int_t& bLow, Int_t& bHigh) const
1549 {
1550   //
1551   //  Get the limits for the 3 correlator method
1552   //
1553   //  Parameters:
1554   //   bin  : reference bin #
1555   //   aLow : Lowest bin to be included in v_A calculations
1556   //   aHigh: Highest bin to be included in v_A calculations
1557   //   bLow : Lowest bin to be included in v_B calculations
1558   //   bHigh: Highest bin to be included in v_B calculations
1559   //
1560   if ((fFlags & kFMD)) {
1561     switch(bin) {
1562       case 0:
1563         aLow = 14; aHigh = 15;
1564         bLow = 20; bHigh = 22;
1565         break;
1566       case 1:
1567         aLow = 16; aHigh = 16;
1568         bLow = 21; bHigh = 22;
1569         break;
1570       case 2:
1571         aLow =  6; aHigh =  7;
1572         bLow = 21; bHigh = 22;
1573         break;
1574       case 3:
1575         aLow =  6; aHigh =  7;
1576         bLow = 12; bHigh = 12;
1577         break;
1578       case 4:
1579         aLow =  6; aHigh =  8;
1580         bLow = 13; bHigh = 14;
1581         break;
1582       default:
1583         AliFatal(Form("No limits for this eta region! (%d)", bin));
1584     }
1585   }
1586   else if ((fFlags & kVZERO)) {
1587     switch(bin) {
1588       case 0:
1589         aLow =  6; aHigh = 13;
1590         bLow = 17; bHigh = 18;
1591         break;
1592       case 1:
1593         aLow =  6; aHigh =  9;
1594         bLow = 17; bHigh = 18;
1595         break;
1596       case 2:
1597         aLow =  2; aHigh =  3;
1598         bLow = 17; bHigh = 18;
1599         break;
1600       case 3:
1601         aLow =  2; aHigh =  3;
1602         bLow =  6; bHigh =  9;
1603         break;
1604       case 4:
1605         aLow =  2; aHigh =  3;
1606         bLow =  6; bHigh = 13;
1607         break;
1608       default:
1609         AliFatal(Form("No limits for this eta region! (%d)", bin));
1610     }
1611   }
1612   // Try to catch cases where fEtaLimits and these values do not correspond to each other
1613   if (aHigh > fCumuNUARef->GetNbinsX() || bHigh > fCumuNUARef->GetNbinsX())
1614     AliFatal(Form("Limits outside vtx range! (%d) - aHigh = %d, bHigh = %d, Nbins = %d",
1615       bin, aHigh, bHigh, fCumuNUARef->GetNbinsX()));
1616 }
1617 //_____________________________________________________________________
1618 void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate3Cor(Double_t cent)
1619 {
1620   //
1621   //  Calculate the Q cumulant up to order fMaxMoment
1622   //
1623   //  Parameters:
1624   //   cent: centrality of event
1625   //
1626   if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1627
1628   // Fill out NUA hists
1629   for (Int_t etaBin = 1; etaBin <= fCumuRef->GetNbinsX(); etaBin++) {
1630     Double_t eta = fCumuRef->GetXaxis()->GetBinCenter(etaBin);
1631     if (fCumuRef->GetBinContent(etaBin, 0) == 0) continue;
1632     for (Int_t qBin = 0; qBin <= fCumuRef->GetNbinsY(); qBin++) {
1633       fCumuNUARef->Fill(eta, cent, Double_t(qBin), fCumuRef->GetBinContent(etaBin, qBin));
1634     }
1635   }
1636   for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1637     Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1638     if (fCumuDiff->GetBinContent(etaBin, 0) == 0) continue;
1639     for (Int_t qBin = 0; qBin <= fCumuDiff->GetNbinsY(); qBin++) {
1640       fCumuNUADiff->Fill(eta, cent, Double_t(qBin), fCumuDiff->GetBinContent(etaBin, qBin));
1641     }
1642   }
1643
1644   // We create the objects needed for the analysis
1645   TH3D* cumuRef = 0;
1646   TH3D* cumuDiff = 0;
1647   // For each n we loop over the hists
1648   for (Int_t n = 2; n <= fMaxMoment; n++) {
1649     cumuRef  = (TH3D*)fCumuHists.Get('r',n);
1650     cumuDiff = (TH3D*)fCumuHists.Get('d',n);
1651
1652     // Per mom. quantities
1653     Int_t prevLim = 0;
1654     Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
1655     Double_t dQnReA = 0, dQnImA = 0, multA = 0;
1656     Double_t dQnReB = 0, dQnImB = 0, multB = 0;
1657     Double_t two = 0, w2 = 0;
1658     for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1659       Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1660       if (fEtaLims[prevLim] < eta) {
1661         GetLimits(prevLim, aLow, aHigh, bLow, bHigh);
1662         prevLim++;
1663         multA = 0; dQnReA = 0; dQnImA = 0;
1664         multB = 0; dQnReB = 0; dQnImB = 0;
1665         // Reference flow
1666         for (Int_t a = aLow; a <= aHigh; a++) {
1667           multA  += fCumuRef->GetBinContent(a, 0);
1668           dQnReA += fCumuRef->GetBinContent(a, GetBinNumberCos(n));
1669           dQnImA += fCumuRef->GetBinContent(a, GetBinNumberSin(n));
1670         }
1671         for (Int_t b = bLow; b <= bHigh; b++) {
1672           multB  += fCumuRef->GetBinContent(b, 0);
1673           dQnReB += fCumuRef->GetBinContent(b, GetBinNumberCos(n));
1674           dQnImB += fCumuRef->GetBinContent(b, GetBinNumberSin(n));
1675         }
1676         // The reference flow is calculated
1677         // 2-particle
1678         w2 = multA * multB;
1679         two = dQnReA*dQnReB + dQnImA*dQnImB;
1680       } // End of reference flow
1681       cumuRef->Fill(eta, cent, kW2Two, two);
1682       cumuRef->Fill(eta, cent, kW2, w2);
1683
1684       // For each etaBin bin the necessary values for differential flow is calculated
1685       Double_t mp = fCumuDiff->GetBinContent(etaBin, 0);
1686       Double_t pnRe = fCumuDiff->GetBinContent(etaBin, GetBinNumberCos(n));
1687       Double_t pnIm = fCumuDiff->GetBinContent(etaBin, GetBinNumberSin(n));
1688       if (mp == 0) continue;
1689
1690       // Differential flow calculations for each eta bin is done:
1691       // 2-particle differential flow
1692       Double_t w2pA = mp * multA;
1693       Double_t twoPrimeA = pnRe*dQnReA + pnIm*dQnImA;
1694       cumuDiff->Fill(eta, cent, kW2Two, twoPrimeA);
1695       cumuDiff->Fill(eta, cent, kW2, w2pA);
1696
1697       Double_t w2pB = mp * multB;
1698       Double_t twoPrimeB = pnRe*dQnReB + pnIm*dQnImB;
1699       cumuDiff->Fill(eta, cent, kW4Four, twoPrimeB);
1700       cumuDiff->Fill(eta, cent, kW4, w2pB);
1701      } // End of eta loop
1702     // Event count
1703     cumuRef->Fill(-7., cent, -0.5, 1.);
1704   } // End of moment loop
1705   return;
1706
1707 }
1708 //_____________________________________________________________________
1709 void AliForwardFlowTaskQC::VertexBin::CumulantsTerminate(TList* inlist, TList* outlist)
1710 {
1711   //
1712   //  Finalizes the Q cumulant calculations
1713   //
1714   //  Parameters:
1715   //   inlist: input sumlist
1716   //   outlist: output result list
1717   //
1718
1719   // Re-find cumulants hist if Terminate is called separately
1720   if (!fCumuHists.IsConnected()) {
1721     TList* list = (TList*)inlist->FindObject(Form("%svertex_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1722     fCumuHists.ConnectList(Form("%sCumu_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)), list);
1723
1724     if (!fCumuNUARef)
1725       fCumuNUARef = (TH3D*)list->FindObject(Form("%s_vertex_%d_%d%s_cumuNUARef", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1726     if (!fCumuNUADiff)
1727       fCumuNUADiff = (TH3D*)list->FindObject(Form("%s_vertex_%d_%d%s_cumuNUADiff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1728   }
1729   // Clone to avoid normalization problems when redoing terminate locally
1730   fCumuNUARef = (TH3D*)fCumuNUARef->Clone(Form("%s_vertex_%d_%d%s_cumuNUARefNorm", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1731   fCumuNUADiff = (TH3D*)fCumuNUADiff->Clone(Form("%s_vertex_%d_%d%s_cumuNUADiffNorm", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1732
1733   // Diagnostics histograms
1734   TH2I* quality = (TH2I*)outlist->FindObject(Form("hQCQuality%s%s", fType.Data(), GetQCType(fFlags)));
1735   if (!quality) {
1736     quality = MakeQualityHist(Form("hQCQuality%s%s", fType.Data(), GetQCType(fFlags)));
1737     outlist->Add(quality);
1738   }
1739   TH1D* cent = (TH1D*)outlist->FindObject(Form("%s%s_cent", fType.Data(), GetQCType(fFlags)));
1740   if (!cent) {
1741     cent = new TH1D(Form("%s%s_cent", fType.Data(), GetQCType(fFlags)),
1742                     Form("%s%s_cent", fType.Data(), GetQCType(fFlags)),
1743                 fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXmin(), fCumuNUARef->GetYaxis()->GetXmax());
1744     cent->GetXaxis()->Set(fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXbins()->GetArray());
1745     outlist->Add(cent);
1746   }
1747   TH2D* dNdetaRef = (TH2D*)outlist->FindObject(Form("%s%s_dNdetaRef", fType.Data(), GetQCType(fFlags)));
1748   if (!dNdetaRef) {
1749     dNdetaRef = new TH2D(Form("%s%s_dNdetaRef", fType.Data(), GetQCType(fFlags)),
1750                                Form("%s%s_dNdetaRef", fType.Data(), GetQCType(fFlags)),
1751                 fCumuNUARef->GetNbinsX(), fCumuNUARef->GetXaxis()->GetXmin(), fCumuNUARef->GetXaxis()->GetXmax(),
1752                 fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXmin(), fCumuNUARef->GetYaxis()->GetXmax());
1753     dNdetaRef->GetYaxis()->Set(fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXbins()->GetArray());
1754     dNdetaRef->Sumw2();
1755     outlist->Add(dNdetaRef);
1756   }
1757   TH2D* dNdetaDiff = (TH2D*)outlist->FindObject(Form("%s%s_dNdetaDiff", fType.Data(), GetQCType(fFlags)));
1758   if (!dNdetaDiff) {
1759     dNdetaDiff = new TH2D(Form("%s%s_dNdetaDiff", fType.Data(), GetQCType(fFlags)),
1760                                 Form("%s%s_dNdetaDiff", fType.Data(), GetQCType(fFlags)),
1761                 fCumuNUADiff->GetNbinsX(), fCumuNUADiff->GetXaxis()->GetXmin(), fCumuNUADiff->GetXaxis()->GetXmax(),
1762                 fCumuNUADiff->GetNbinsY(), fCumuNUADiff->GetYaxis()->GetXmin(), fCumuNUADiff->GetYaxis()->GetXmax());
1763     dNdetaDiff->GetYaxis()->Set(fCumuNUADiff->GetNbinsY(), fCumuNUADiff->GetYaxis()->GetXbins()->GetArray());
1764     dNdetaDiff->Sumw2();
1765     outlist->Add(dNdetaDiff);
1766   }
1767
1768   // Setting up outputs
1769   // Create output lists and diagnostics
1770   TList* vtxList = (TList*)outlist->FindObject("vtxList");
1771   if (!vtxList) {
1772     vtxList = new TList();
1773     vtxList->SetName("vtxList");
1774     outlist->Add(vtxList);
1775   }
1776   vtxList->Add(fCumuNUARef);
1777   vtxList->Add(fCumuNUADiff);
1778
1779   // Setup output profiles
1780   CumuHistos cumu2(fMaxMoment, ((fFlags & kNUAcorr) ? 2 : 0));
1781   CumuHistos cumu4(fMaxMoment, ((fFlags & kNUAcorr) ? 1 : 0));
1782
1783   cumu2.ConnectList(Form("%sQC2_Cumu%s_vtx_%d_%d", fType.Data(), GetQCType(fFlags), fVzMin, fVzMax), vtxList);
1784   if ((fFlags & kStdQC))
1785     cumu4.ConnectList(Form("%sQC4_Cumu%s_vtx_%d_%d", fType.Data(), GetQCType(fFlags), fVzMin, fVzMax), vtxList);
1786
1787   for (Int_t n = 2; n <= fMaxMoment; n++) {
1788     // 2-particle
1789     cumu2.Add(MakeOutputHist(2, n, "Ref", CumuHistos::kNoNUA));
1790     if ((fFlags & k3Cor)){
1791       cumu2.Add(MakeOutputHist(2, n, "DiffA", CumuHistos::kNoNUA));
1792       cumu2.Add(MakeOutputHist(2, n, "DiffB", CumuHistos::kNoNUA));
1793     } else {
1794       cumu2.Add(MakeOutputHist(2, n, "Diff", CumuHistos::kNoNUA));
1795     }
1796     // 4-particle
1797     if ((fFlags & kStdQC)) {
1798       cumu4.Add(MakeOutputHist(4, n, "Ref", CumuHistos::kNoNUA));
1799       cumu4.Add(MakeOutputHist(4, n, "Diff", CumuHistos::kNoNUA));
1800     }
1801   } // End of v_n result loop
1802   // NUA corrected
1803   if ((fFlags & kNUAcorr)) {
1804     for (Int_t n = 2; n <= fMaxMoment; n++) {
1805       // 2-particle
1806       cumu2.Add(MakeOutputHist(2, n, "Ref", CumuHistos::kNUAOld));
1807       if ((fFlags & k3Cor)) {
1808         cumu2.Add(MakeOutputHist(2, n, "DiffA", CumuHistos::kNUAOld));
1809         cumu2.Add(MakeOutputHist(2, n, "DiffB", CumuHistos::kNUAOld));
1810       } else {
1811         cumu2.Add(MakeOutputHist(2, n, "Diff", CumuHistos::kNUAOld));
1812       }
1813       // 4-particle
1814       if ((fFlags & kStdQC)) {
1815         cumu4.Add(MakeOutputHist(4, n, "Ref", CumuHistos::kNUAOld));
1816         cumu4.Add(MakeOutputHist(4, n, "Diff", CumuHistos::kNUAOld));
1817       }
1818     }
1819     for (Int_t n = 2; n <= fMaxMoment; n++) {
1820       // 2-particle
1821       cumu2.Add(MakeOutputHist(2, n, "Ref", CumuHistos::kNUA));
1822       if ((fFlags & k3Cor)) {
1823         cumu2.Add(MakeOutputHist(2, n, "DiffA", CumuHistos::kNUA));
1824         cumu2.Add(MakeOutputHist(2, n, "DiffB", CumuHistos::kNUA));
1825       } else {
1826         cumu2.Add(MakeOutputHist(2, n, "Diff", CumuHistos::kNUA));
1827       }
1828     }
1829   }
1830
1831   // Calculating the cumulants
1832   if ((fFlags & k3Cor)) {
1833     Calculate3CorFlow(cumu2, quality, cent, dNdetaRef, dNdetaDiff);
1834   } else {
1835     CalculateReferenceFlow(cumu2, cumu4, quality, cent, dNdetaRef);
1836     CalculateDifferentialFlow(cumu2, cumu4, quality, dNdetaDiff);
1837   }
1838   if ((fFlags & kNUAcorr)) {
1839     SolveCoupledFlowEquations(cumu2, 'r');
1840     if ((fFlags & k3Cor)) {
1841       SolveCoupledFlowEquations(cumu2, 'a');
1842       SolveCoupledFlowEquations(cumu2, 'b');
1843     } else {
1844       SolveCoupledFlowEquations(cumu2, 'd');
1845     }
1846   }
1847
1848   // Add to output for immediate viewing - individual vtx bins are used for final results
1849   AddVertexBins(cumu2, outlist, ((fFlags & kNUAcorr) ? 2 : 0));
1850   if ((fFlags & kStdQC)) AddVertexBins(cumu4, outlist, ((fFlags & kNUAcorr) ? 1 : 0));
1851
1852   // Setup NUA diagnoastics histograms
1853   TList* nualist = (TList*)outlist->FindObject("NUATerms");
1854   if (!nualist) {
1855     nualist = new TList();
1856     nualist->SetName("NUATerms");
1857     outlist->Add(nualist);
1858   }
1859   // Reference
1860   TH2D* nuaRef = (TH2D*)nualist->FindObject(Form("%sReferenceNUA%s", fType.Data(), GetQCType(fFlags)));
1861   TH2D* temp = 0;
1862   if (!nuaRef) {
1863     nuaRef = (TH2D*)fCumuNUARef->Project3D("yz");
1864     nuaRef->Scale(1./fCumuNUARef->GetNbinsX());
1865     nuaRef->SetName(Form("%sReferenceNUA%s", fType.Data(), GetQCType(fFlags)));
1866     nuaRef->SetTitle(Form("%sReferenceNUA%s", fType.Data(), GetQCType(fFlags)));
1867     nualist->Add(nuaRef);
1868   } else {
1869     temp = (TH2D*)fCumuNUARef->Project3D("yz");
1870     temp->Scale(1./fCumuNUARef->GetNbinsX());
1871     nuaRef->Add(temp);
1872     delete temp;
1873   }
1874   // Filling in underflow to make scaling possible in Terminate()
1875   nuaRef->Fill(0., -1., 1.);
1876   // Differential
1877   TH2D* nuaDiff = (TH2D*)nualist->FindObject(Form("%sDifferentialNUA%s", fType.Data(), GetQCType(fFlags)));
1878   if (!nuaDiff) {
1879     nuaDiff = (TH2D*)fCumuNUADiff->Project3D("yz");
1880     nuaDiff->SetName(Form("%sDifferentialNUA%s", fType.Data(), GetQCType(fFlags)));
1881     nuaDiff->SetTitle(Form("%sDifferentialNUA%s", fType.Data(), GetQCType(fFlags)));
1882     nualist->Add(nuaDiff);
1883   } else {
1884     temp = (TH2D*)fCumuNUADiff->Project3D("yz");
1885     nuaDiff->Add(temp);
1886     delete temp;
1887   }
1888   // Filling in underflow to make scaling possible in Terminate()
1889   nuaDiff->Fill(0., -1., 1.);
1890
1891   return;
1892 }
1893 //_____________________________________________________________________
1894 void AliForwardFlowTaskQC::VertexBin::CalculateReferenceFlow(CumuHistos& cumu2h, CumuHistos& cumu4h, TH2I* quality,
1895                                                              TH1D* chist, TH2D* dNdetaRef) const
1896 {
1897   //
1898   //  Calculates the reference flow
1899   //
1900   //  Parameters:
1901   //   cumu2h: CumuHistos object with QC{2} cumulants
1902   //   cumu4h: CumuHistos object with QC{4} cumulants
1903   //   quality: Histogram for success rate of cumulants
1904   //   chist: Centrality histogram
1905   //   dNdetaRef: dN/deta histogram for estimating multiplicity used for ref calculations
1906   //
1907
1908   // Normalizing common NUA hists
1909   for (Int_t cBin = 1; cBin <= fCumuNUARef->GetNbinsY(); cBin++) {
1910     Double_t cent = fCumuNUARef->GetYaxis()->GetBinCenter(cBin);
1911     for (Int_t eBin = 1; eBin <= fCumuNUARef->GetNbinsX(); eBin++) {
1912       Double_t eta = fCumuNUARef->GetXaxis()->GetBinCenter(eBin);
1913       Double_t mult = fCumuNUARef->GetBinContent(eBin, cBin, 0);
1914       if (mult == 0) continue;
1915       for (Int_t qBin = 1; qBin <= fCumuNUARef->GetNbinsZ(); qBin++) {
1916         fCumuNUARef->SetBinContent(eBin, cBin, qBin, fCumuNUARef->GetBinContent(eBin, cBin, qBin)/mult);
1917         fCumuNUARef->SetBinError(eBin, cBin, qBin, fCumuNUARef->GetBinError(eBin, cBin, qBin)/mult);
1918       }
1919       // Fill dN/deta diagnostics
1920       dNdetaRef->Fill(eta, cent, mult);
1921     }
1922   }
1923
1924   // For flow calculations
1925   TH3D* cumuRef = 0;
1926   TH2D* cumu2 = 0;
1927   TH2D* cumu2NUAold = 0;
1928   TH2D* cumu4 = 0;
1929   TH2D* cumu4NUA = 0;
1930   Int_t qualityFactor = ((fFlags & kStdQC) ? 8 : 4); // used for correctly filling in quality hists
1931   // Loop over cumulant histogram for final calculations
1932   for (Int_t n = 2; n <= fMaxMoment; n++) { // Moment loop begins
1933     cumu2 = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNoNUA);
1934     if ((fFlags & kNUAcorr)) {
1935       cumu2NUAold = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNUAOld);
1936     }
1937     if ((fFlags & kStdQC)) {
1938       cumu4 = (TH2D*)cumu4h.Get('r', n, CumuHistos::kNoNUA);
1939       if ((fFlags & kNUAcorr)) cumu4NUA = (TH2D*)cumu4h.Get('r', n, CumuHistos::kNUAOld);
1940     }
1941     cumuRef  = (TH3D*)fCumuHists.Get('r', n);
1942     // Begin loops
1943     for (Int_t cBin = 1; cBin <= cumuRef->GetNbinsY(); cBin++) { // Centrality loop begins
1944       Double_t cent = cumuRef->GetYaxis()->GetBinCenter(cBin);
1945       if (n == 2) chist->Fill(cent, cumuRef->GetBinContent(0, cBin, 0));
1946       if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), n, cent));
1947       for (Int_t etaBin = 1; etaBin <= cumuRef->GetNbinsX(); etaBin++) { // Eta loop begins
1948         Double_t eta = cumuRef->GetXaxis()->GetBinCenter(etaBin);
1949         Double_t refEta = eta;
1950         Int_t refEtaBinA = fCumuNUARef->GetXaxis()->FindBin(refEta);
1951         if ((fFlags & kEtaGap)) refEta = -eta;
1952         Int_t refEtaBinB = fCumuNUARef->GetXaxis()->FindBin(refEta);
1953         // 2-particle reference flow
1954         Double_t w2Two = cumuRef->GetBinContent(refEtaBinA, cBin, kW2Two);
1955         Double_t w2 = cumuRef->GetBinContent(refEtaBinA, cBin, kW2);
1956         if (w2 == 0) continue;
1957         Double_t cosP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberCos(n));
1958         Double_t sinP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberSin(n));
1959         Double_t cosP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(n));
1960         Double_t sinP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(n));
1961         Double_t cos2nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberCos(2*n));
1962         Double_t sin2nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberSin(2*n));
1963         Double_t cos2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(2*n));
1964         Double_t sin2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(2*n));
1965         Double_t two = w2Two / w2;
1966         Double_t qc2 = two;
1967         if (qc2 >= 0) cumu2->Fill(eta, cent, TMath::Sqrt(qc2));
1968
1969         if ((fFlags & kNUAcorr)) {
1970           // Old NUA
1971           // With no eta gap the last two terms are <<cos(phi)>>^2 and <<sin(phi)>>^2,
1972           // with eta gap the different coverage is taken into account.
1973           // The next line covers both cases.
1974           qc2 -= cosP1nPhiA*cosP1nPhiB + sinP1nPhiA*sinP1nPhiB;
1975           // Extra NUA term from 2n cosines and sines
1976           Double_t den = 1-(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB);
1977           if (den != 0) qc2 /= den;
1978           else qc2 = 0;
1979         }
1980         if (qc2 <= 0) {
1981           if (fDebug > 0)
1982             AliInfo(Form("%s: QC_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f - skipping",
1983             fType.Data(), n, qc2, eta, cent));
1984           quality->Fill((n-2)*qualityFactor+2, Int_t(cent));
1985           continue;
1986         }
1987         Double_t vnTwo = TMath::Sqrt(qc2);
1988         if (!TMath::IsNaN(vnTwo)) {
1989           quality->Fill((n-2)*qualityFactor+1, Int_t(cent));
1990           if ((fFlags & kNUAcorr)) cumu2NUAold->Fill(eta, cent, vnTwo);
1991         }
1992
1993         if (!(fFlags & kStdQC)) continue;
1994         // 4-particle reference flow
1995         Double_t w4Four = cumuRef->GetBinContent(refEtaBinA, cBin, kW4Four);
1996         Double_t w4 = cumuRef->GetBinContent(refEtaBinA, cBin, kW4);
1997         Double_t multm1m2 = cumuRef->GetBinContent(refEtaBinA, cBin, k3pWeight);
1998         if (w4 == 0 || multm1m2 == 0) continue;
1999         Double_t cosP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2);
2000         Double_t sinP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2);
2001         Double_t cosP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2phi3m);
2002         Double_t sinP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2phi3m);
2003
2004         cosP1nPhi1P1nPhi2 /= w2;
2005         sinP1nPhi1P1nPhi2 /= w2;
2006         cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2007         sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2008         Double_t four = w4Four / w4;
2009         Double_t qc4 = four-2.*TMath::Power(two,2.);
2010         if (qc4 < 0) cumu4->Fill(eta, cent, TMath::Power(-qc4, 0.25));
2011
2012         if ((fFlags & kNUAcorr)) {
2013            qc4 += - 4.*cosP1nPhiA*cosP1nPhi1M1nPhi2M1nPhi3
2014                   + 4.*sinP1nPhiA*sinP1nPhi1M1nPhi2M1nPhi3-TMath::Power(cosP1nPhi1P1nPhi2,2.)-TMath::Power(sinP1nPhi1P1nPhi2,2.)
2015                   + 4.*cosP1nPhi1P1nPhi2*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
2016                   + 8.*sinP1nPhi1P1nPhi2*sinP1nPhiA*cosP1nPhiA
2017                   + 8.*two*(TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.))
2018                   - 6.*TMath::Power((TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.)),2.);
2019         }
2020         if (qc4 >= 0) {
2021           if (fDebug > 0)
2022             AliInfo(Form("%s: QC_%d{4} = %1.3f for eta = %1.2f and centrality %3.1f - skipping",
2023             fType.Data(), n, qc2, eta, cent));
2024           quality->Fill((n-2)*qualityFactor+6, Int_t(cent));
2025           continue;
2026         }
2027         Double_t vnFour = TMath::Power(-qc4, 0.25);
2028         if (!TMath::IsNaN(vnFour*multm1m2)) {
2029           quality->Fill((n-2)*qualityFactor+5, Int_t(cent));
2030           if ((fFlags & kNUAcorr)) cumu4NUA->Fill(eta, cent, vnFour);
2031         }
2032       } // End of eta
2033     } // End of cent
2034   } // End of moment
2035
2036   return;
2037 }
2038 //_____________________________________________________________________
2039 void AliForwardFlowTaskQC::VertexBin::CalculateDifferentialFlow(CumuHistos& cumu2h, CumuHistos& cumu4h,
2040                                                                 TH2I* quality, TH2D* dNdetaDiff) const
2041 {
2042   //
2043   //  Calculates the differential flow
2044   //
2045   //  Parameters:
2046   //   cumu2h: CumuHistos object with QC{2} cumulants
2047   //   cumu4h: CumuHistos object with QC{4} cumulants
2048   //   quality: Histogram for success rate of cumulants
2049   //   dNdetaDiff: dN/deta histogram for estimating multiplicity used for diff calculations
2050   //
2051
2052   for (Int_t cBin = 1; cBin <= fCumuNUADiff->GetNbinsY(); cBin++) {
2053     Double_t cent = fCumuNUADiff->GetYaxis()->GetBinCenter(cBin);
2054     for (Int_t eBin = 1; eBin <= fCumuNUADiff->GetNbinsX(); eBin++) {
2055       Double_t eta = fCumuNUADiff->GetXaxis()->GetBinCenter(eBin);
2056       Double_t mult = fCumuNUADiff->GetBinContent(eBin, cBin, 0);
2057       if (mult == 0) continue;
2058       for (Int_t qBin = 1; qBin <= fCumuNUADiff->GetNbinsZ(); qBin++) {
2059         fCumuNUADiff->SetBinContent(eBin, cBin, qBin, fCumuNUADiff->GetBinContent(eBin, cBin, qBin)/mult);
2060         fCumuNUADiff->SetBinError(eBin, cBin, qBin, fCumuNUADiff->GetBinError(eBin, cBin, qBin)/mult);
2061       }
2062       dNdetaDiff->Fill(eta, cent, mult);
2063     }
2064   }
2065
2066   // For flow calculations
2067   TH3D* cumuRef = 0;
2068   TH3D* cumuDiff = 0;
2069   TH2D* cumu2 = 0;
2070   TH2D* cumu2NUAold = 0;
2071   TH2D* cumu4 = 0;
2072   TH2D* cumu4NUA = 0;
2073   Int_t qualityFactor = ((fFlags & kStdQC) ? 8 : 4); // used for correctly filling in quality hists
2074   // Loop over cumulant histogram for final calculations
2075   for (Int_t n = 2; n <= fMaxMoment; n++) { // Moment loop begins
2076     cumu2 = (TH2D*)cumu2h.Get('d', n, CumuHistos::kNoNUA);
2077     if ((fFlags & kNUAcorr)) {
2078       cumu2NUAold = (TH2D*)cumu2h.Get('d', n, CumuHistos::kNUAOld);
2079     }
2080     if ((fFlags & kStdQC)) {
2081       cumu4 = (TH2D*)cumu4h.Get('d',n);
2082       if ((fFlags & kNUAcorr)) cumu4NUA = (TH2D*)cumu4h.Get('d', n, CumuHistos::kNUAOld);
2083     }
2084     cumuRef  = (TH3D*)fCumuHists.Get('r',n);
2085     cumuDiff = (TH3D*)fCumuHists.Get('d',n);
2086     for (Int_t cBin = 1; cBin <= cumuDiff->GetNbinsY(); cBin++) { // Centrality loop begins
2087       Double_t cent = cumuDiff->GetYaxis()->GetBinCenter(cBin);
2088       if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), n, cent));
2089       for (Int_t etaBin = 1; etaBin <= cumuDiff->GetNbinsX(); etaBin++) { // Eta loop begins
2090         Double_t eta = cumuDiff->GetXaxis()->GetBinCenter(etaBin);
2091         Double_t refEta = eta;
2092         Int_t refEtaBinA = fCumuNUARef->GetXaxis()->FindBin(refEta);
2093         if ((fFlags & kEtaGap)) refEta = -eta;
2094         Int_t refEtaBinB = fCumuNUARef->GetXaxis()->FindBin(refEta);
2095
2096         // Reference objects
2097         Double_t w2 = cumuRef->GetBinContent(refEtaBinA, cBin, kW2);
2098         if (w2 == 0) continue;
2099         Double_t two = cumuRef->GetBinContent(refEtaBinA, cBin, kW2Two);
2100         two /= w2;
2101         Double_t cosP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberCos(n));
2102         Double_t sinP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberSin(n));
2103         Double_t cosP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(n));
2104         Double_t sinP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(n));
2105         Double_t cos2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(2*n));
2106         Double_t sin2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(2*n));
2107
2108         // 2-particle differential flow
2109         Double_t w2pTwoPrime = cumuDiff->GetBinContent(etaBin, cBin, kW2Two);
2110         Double_t w2p = cumuDiff->GetBinContent(etaBin, cBin, kW2);
2111         if (w2p == 0) continue;
2112         Double_t cosP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(n));
2113         Double_t sinP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(n));
2114         Double_t cos2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(2*n));
2115         Double_t sin2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(2*n));
2116         Double_t twoPrime = w2pTwoPrime / w2p;
2117
2118         Double_t qc2Prime = twoPrime;
2119         cumu2->Fill(eta, cent, qc2Prime);
2120         if ((fFlags & kNUAcorr)) {
2121           // Old nua
2122           qc2Prime -= cosP1nPsi*cosP1nPhiB + sinP1nPsi*sinP1nPhiB;
2123           // Extra NUA term from 2n cosines and sines
2124           qc2Prime /= (1.-(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
2125         }
2126         if (!TMath::IsNaN(qc2Prime)) {
2127           quality->Fill((n-2)*qualityFactor+3, Int_t(cent));
2128           if ((fFlags & kNUAcorr)) cumu2NUAold->Fill(eta, cent, qc2Prime);
2129         }
2130         else
2131           quality->Fill((n-2)*qualityFactor+4, Int_t(cent));
2132         if (fDebug > 1)
2133           AliInfo(Form("%s: QC'_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f",
2134           fType.Data(), n, qc2Prime, eta, cent));
2135
2136         if (!(fFlags & kStdQC)) continue;
2137         // Reference objects
2138         Double_t cosP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2);
2139         Double_t sinP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2);
2140         Double_t cosP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2phi3m);
2141         Double_t sinP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2phi3m);
2142         Double_t multm1m2 = cumuRef->GetBinContent(refEtaBinA, cBin, k3pWeight);
2143         cosP1nPhi1P1nPhi2 /= w2;
2144         sinP1nPhi1P1nPhi2 /= w2;
2145         cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2146         sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2147
2148         // 4-particle differential flow
2149         Double_t w4pFourPrime = cumuDiff->GetBinContent(etaBin, cBin, kW4Four);
2150         Double_t w4p = cumuDiff->GetBinContent(etaBin, cBin, kW4);
2151         Double_t mpqMult = cumuDiff->GetBinContent(etaBin, cBin, k3pWeight);
2152         if (w4p == 0 || mpqMult == 0) continue;
2153         Double_t cosP1nPsi1P1nPhi2 = cumuDiff->GetBinContent(etaBin, cBin, kCosphi1phi2);
2154         Double_t sinP1nPsi1P1nPhi2 = cumuDiff->GetBinContent(etaBin, cBin, kSinphi1phi2);
2155         Double_t cosP1nPsi1M1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kCosphi1phi2phi3m);
2156         Double_t sinP1nPsi1M1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kSinphi1phi2phi3m);
2157         Double_t cosP1nPsi1P1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kCosphi1phi2phi3p);
2158         Double_t sinP1nPsi1P1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kSinphi1phi2phi3p);
2159
2160         cosP1nPsi1P1nPhi2 /= w2p;
2161         sinP1nPsi1P1nPhi2 /= w2p;
2162         cosP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
2163         sinP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
2164         cosP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
2165         sinP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
2166
2167         Double_t fourPrime = w4pFourPrime / w4p;
2168         Double_t qc4Prime = fourPrime-2.*twoPrime*two;
2169         if (cumu4) cumu4->Fill(eta, cent, qc4Prime);
2170
2171         if ((fFlags & kNUAcorr)) {
2172           qc4Prime += - cosP1nPsi*cosP1nPhi1M1nPhi2M1nPhi3
2173                       + sinP1nPsi*sinP1nPhi1M1nPhi2M1nPhi3
2174                       - cosP1nPhiA*cosP1nPsi1M1nPhi2M1nPhi3
2175                       + sinP1nPhiA*sinP1nPsi1M1nPhi2M1nPhi3
2176                       - 2.*cosP1nPhiA*cosP1nPsi1P1nPhi2M1nPhi3
2177                       - 2.*sinP1nPhiA*sinP1nPsi1P1nPhi2M1nPhi3
2178                       - cosP1nPsi1P1nPhi2*cosP1nPhi1P1nPhi2
2179                       - sinP1nPsi1P1nPhi2*sinP1nPhi1P1nPhi2
2180                       + 2.*cosP1nPhi1P1nPhi2*(cosP1nPsi*cosP1nPhiA-sinP1nPsi*sinP1nPhiA)
2181                       + 2.*sinP1nPhi1P1nPhi2*(cosP1nPsi*sinP1nPhiA+sinP1nPsi*cosP1nPhiA)
2182                       + 4.*two*(cosP1nPsi*cosP1nPhiA+sinP1nPsi*sinP1nPhiA)
2183                       + 2.*cosP1nPsi1P1nPhi2*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
2184                       + 4.*sinP1nPsi1P1nPhi2*cosP1nPhiA*sinP1nPhiA
2185                       + 4.*twoPrime*(TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.))
2186                       - 6.*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
2187                       * (cosP1nPsi*cosP1nPhiA-sinP1nPsi*sinP1nPhiA)
2188                       - 12.*cosP1nPhiA*sinP1nPhiA
2189                       * (sinP1nPsi*cosP1nPhiA+cosP1nPsi*sinP1nPhiA);
2190         }
2191 //      Double_t vnFourDiff = - qc4Prime / TMath::Power(-qc4, 0.75);
2192         if (!TMath::IsNaN(qc4Prime*mpqMult)) {
2193           quality->Fill((n-2)*qualityFactor+7, Int_t(cent));
2194           if (cumu4NUA) cumu4NUA->Fill(eta, cent, qc4Prime);
2195         }
2196         else
2197           quality->Fill((n-2)*qualityFactor+8, Int_t(cent));
2198         if (fDebug > 1)
2199           AliInfo(Form("%s: v_%d{4} = %1.3f for eta = %1.2f and centrality %3.1f",
2200           fType.Data(), n, qc4Prime, eta, cent));
2201       } // End of eta loop
2202     } // End of centrality loop
2203   } // End of moment
2204
2205   return;
2206 }
2207 //_____________________________________________________________________
2208 void AliForwardFlowTaskQC::VertexBin::Calculate3CorFlow(CumuHistos& cumu2h, TH2I* quality, TH1D* chist,
2209                                                         TH2D* dNdetaRef, TH2D* dNdetaDiff) const
2210 {
2211   //
2212   //  Calculates the 3 sub flow
2213   //
2214   //  Parameters:
2215   //   cumu2h: CumuHistos object with QC{2} cumulants
2216   //   quality: Histogram for success rate of cumulants
2217   //   chist: Centrality histogram
2218   //   dNdetaDiff: dN/deta histogram for estimating multiplicity used for diff calculations
2219   //
2220
2221   // For flow calculations
2222   TH3D* cumuRef = 0;
2223   TH3D* cumuDiff = 0;
2224   TH2D* cumu2r = 0;
2225   TH2D* cumu2rNUAold = 0;
2226   TH2D* cumu2a = 0;
2227   TH2D* cumu2aNUAold = 0;
2228   TH2D* cumu2b = 0;
2229   TH2D* cumu2bNUAold = 0;
2230   // Loop over cumulant histogram for final calculations
2231   for (Int_t n = 2; n <= fMaxMoment; n++) { // Moment loop begins
2232     cumu2r = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNoNUA);
2233     cumu2a = (TH2D*)cumu2h.Get('a', n, CumuHistos::kNoNUA);
2234     cumu2b = (TH2D*)cumu2h.Get('b', n, CumuHistos::kNoNUA);
2235     if ((fFlags & kNUAcorr)) {
2236       cumu2rNUAold = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNUAOld);
2237       cumu2aNUAold = (TH2D*)cumu2h.Get('a', n, CumuHistos::kNUAOld);
2238       cumu2bNUAold = (TH2D*)cumu2h.Get('b', n, CumuHistos::kNUAOld);
2239     }
2240     cumuRef  = (TH3D*)fCumuHists.Get('r',n);
2241     cumuDiff = (TH3D*)fCumuHists.Get('d',n);
2242     for (Int_t cBin = 1; cBin <= cumuRef->GetNbinsY(); cBin++) { // Centrality loop begins
2243       Double_t cent = cumuRef->GetYaxis()->GetBinCenter(cBin);
2244       if (n == 2) chist->Fill(cent, cumuRef->GetBinContent(0, cBin, 0));
2245       if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), n, cent));
2246       // Here it starts!
2247       Int_t prevLim = 0;
2248       Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2249       Double_t cosP1nPhiA = 0;
2250       Double_t sinP1nPhiA = 0;
2251       Double_t cos2nPhiA = 0;
2252       Double_t sin2nPhiA = 0;
2253       Double_t cosP1nPhiB = 0;
2254       Double_t sinP1nPhiB = 0;
2255       Double_t cos2nPhiB = 0;
2256       Double_t sin2nPhiB = 0;
2257       Double_t multA = 0;
2258       Double_t multB = 0;
2259
2260       for (Int_t etaBin = 1; etaBin <= cumuDiff->GetNbinsX(); etaBin++) { // Eta loop begins
2261         Double_t eta = cumuDiff->GetXaxis()->GetBinCenter(etaBin);
2262         // 2-particle reference flow
2263         Double_t w2Two = cumuRef->GetBinContent(etaBin, cBin, kW2Two);
2264         Double_t w2 = cumuRef->GetBinContent(etaBin, cBin, kW2);
2265         if (w2 == 0) continue;
2266
2267         // Update NUA for new range!
2268         if (fEtaLims[prevLim] < eta) {
2269           GetLimits(prevLim, aLow, aHigh, bLow, bHigh);
2270           prevLim++;
2271           cosP1nPhiA = 0; sinP1nPhiA = 0; cos2nPhiA = 0; sin2nPhiA = 0; multA = 0;
2272           cosP1nPhiB = 0; sinP1nPhiB = 0; cos2nPhiB = 0; sin2nPhiB = 0; multB = 0;
2273           for (Int_t a = aLow; a <= aHigh; a++) {
2274             cosP1nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n));
2275             sinP1nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n));
2276             cos2nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(2*n));
2277             sin2nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(2*n));
2278             multA += fCumuNUARef->GetBinContent(a, cBin, 0);
2279           }
2280           for (Int_t b = bLow; b <= bHigh; b++) {
2281             cosP1nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n));
2282             sinP1nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n));
2283             cos2nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(2*n));
2284             sin2nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(2*n));
2285             multB += fCumuNUARef->GetBinContent(b, cBin, 0);
2286           }
2287           if (multA == 0 || multB == 0) AliFatal("Empty NUA values for 3Cor!");
2288           cosP1nPhiA /= multA;
2289           sinP1nPhiA /= multA;
2290           cos2nPhiA /= multA;
2291           sin2nPhiA /= multA;
2292           cosP1nPhiB /= multB;
2293           sinP1nPhiB /= multB;
2294           cos2nPhiB /= multB;
2295           sin2nPhiB /= multB;
2296
2297           dNdetaRef->Fill(eta, cent, multA+multB);
2298         }
2299         Double_t two = w2Two / w2;
2300
2301         Double_t qc2 = two;
2302         if (qc2 >= 0) cumu2r->Fill(eta, cent, TMath::Sqrt(qc2));
2303
2304         if ((fFlags & kNUAcorr)) {
2305           // Old nua
2306           qc2 -= cosP1nPhiA*cosP1nPhiB + sinP1nPhiA*sinP1nPhiB;
2307           // Extra NUA term from 2n cosines and sines
2308           qc2 /= (1-(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB));
2309         }
2310         if (qc2 <= 0) {
2311           if (fDebug > 0)
2312             AliInfo(Form("%s: QC_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f - skipping",
2313             fType.Data(), n, qc2, eta, cent));
2314           quality->Fill((n-2)*4+2, Int_t(cent));
2315           continue;
2316         }
2317         Double_t vnTwo = TMath::Sqrt(qc2);
2318         if (!TMath::IsNaN(vnTwo)) {
2319           quality->Fill((n-2)*4+1, Int_t(cent));
2320           if ((fFlags & kNUAcorr)) cumu2rNUAold->Fill(eta, cent, vnTwo);
2321         }
2322
2323         // 2-particle differential flow
2324         Double_t w2pTwoPrimeA = cumuDiff->GetBinContent(etaBin, cBin, kW2Two);
2325         Double_t w2pA = cumuDiff->GetBinContent(etaBin, cBin, kW2);
2326         Double_t w2pTwoPrimeB = cumuDiff->GetBinContent(etaBin, cBin, kW4Four);
2327         Double_t w2pB = cumuDiff->GetBinContent(etaBin, cBin, kW4);
2328         if (w2pA == 0 || w2pB == 0) continue;
2329         Double_t cosP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(n));
2330         Double_t sinP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(n));
2331         Double_t cos2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(2*n));
2332         Double_t sin2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(2*n));
2333         Double_t mult = fCumuNUADiff->GetBinContent(etaBin, cBin, 0);
2334         if (mult == 0) continue;
2335         cosP1nPsi /= mult;
2336         sinP1nPsi /= mult;
2337         cos2nPsi /= mult;
2338         sin2nPsi /= mult;
2339         Double_t twoPrimeA = w2pTwoPrimeA / w2pA;
2340         Double_t twoPrimeB = w2pTwoPrimeB / w2pB;
2341         dNdetaDiff->Fill(eta, cent, mult);
2342
2343         Double_t qc2PrimeA = twoPrimeA;
2344         Double_t qc2PrimeB = twoPrimeB;
2345         if (qc2PrimeA*qc2PrimeB >= 0) {
2346           cumu2a->Fill(eta, cent, qc2PrimeA);
2347           cumu2b->Fill(eta, cent, qc2PrimeB);
2348         }
2349         if ((fFlags & kNUAcorr)) {
2350           // Old nua
2351           qc2PrimeA -= cosP1nPsi*cosP1nPhiA + sinP1nPsi*sinP1nPhiA;
2352           qc2PrimeB -= cosP1nPsi*cosP1nPhiB + sinP1nPsi*sinP1nPhiB; // Is this OK?
2353           // Extra NUA term from 2n cosines and sines
2354           qc2PrimeA /= (1.-(cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA));
2355           qc2PrimeB /= (1.-(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
2356         }
2357         if (!TMath::IsNaN(qc2PrimeA) && !TMath::IsNaN(qc2PrimeB) && qc2 != 0) {
2358         if (qc2PrimeA*qc2PrimeB >= 0) {
2359           quality->Fill((n-2)*4+3, Int_t(cent));
2360           if ((fFlags & kNUAcorr)) cumu2aNUAold->Fill(eta, cent, qc2PrimeA);
2361           if ((fFlags & kNUAcorr)) cumu2bNUAold->Fill(eta, cent, qc2PrimeB);
2362         }
2363       }
2364       else
2365         quality->Fill((n-2)*4+4, Int_t(cent));
2366         if (fDebug > 1)
2367           AliInfo(Form("%s: QC'a_%d{2} = %1.3f, QC'b_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f",
2368           fType.Data(), n, qc2PrimeA, n, qc2PrimeB, eta, cent));
2369       } // End of eta loop
2370     } // End of centrality loop
2371   } // End of moment
2372
2373   return;
2374 }
2375 //_____________________________________________________________________
2376 void AliForwardFlowTaskQC::VertexBin::SolveCoupledFlowEquations(CumuHistos& cumu, const Char_t type) const
2377 {
2378   //
2379   //  Function to solve the coupled flow equations
2380   //  We solve it by using matrix calculations:
2381   //  A*v_n = V => v_n = A^-1*V
2382   //  First we set up a TMatrixD container to make ROOT
2383   //  do the inversions in an efficient way, we multiply the current <<2>> estimates.
2384   //  Then we fill new TH2D's if the corrected <<2>>'s (cumuNUA object).
2385   //
2386   //  Parameters:
2387   //   cumu: CumuHistos object - uncorrected
2388   //   type: Reference ('r') or differential ('d') or ('a' or 'b') for 3 correlator
2389   //
2390
2391   // We start by declaring Matrix and vector objects, as their constructors are quite heavy
2392   TMatrixD mNUA(fMaxMoment-1, fMaxMoment-1);
2393   TVectorD vQC2(fMaxMoment-1);
2394
2395   for (Int_t cBin = 1; cBin <= cumu.Get(type, 2, CumuHistos::kNUAOld)->GetNbinsY(); cBin++) { // cent loop
2396     Double_t cent = cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinCenter(cBin);
2397     for (Int_t eBin = 1; eBin <= cumu.Get(type, 2, CumuHistos::kNUAOld)->GetNbinsX(); eBin++) { // eta loop
2398       Double_t eta = cumu.Get(type, 2, CumuHistos::kNUAOld)->GetXaxis()->GetBinCenter(eBin);
2399       mNUA.Zero(); // reset matrix
2400       vQC2.Zero(); // reset vector
2401       for (Int_t n = 0; n < fMaxMoment-1; n++) { // moment loop
2402         vQC2(n) = static_cast<TH2D*>(cumu.Get(type, n+2, CumuHistos::kNUAOld))->GetBinContent(eBin, cBin);
2403         if (type == 'r' || type == 'R') vQC2(n) *= vQC2(n); // going back to <<2>>
2404         for (Int_t m = 0; m < fMaxMoment-1; m++) { // cross moment
2405           mNUA(n,m) = CalculateNUAMatrixElement(n, m, type, eBin, cBin);
2406         } // End of cross moment loop
2407       } // End of moment loop
2408       // Invert matrix
2409       Double_t det = 0;
2410       mNUA.Invert(&det);
2411       // If determinant is non-zero we go with corrected results
2412       if (det != 0 ) vQC2 = mNUA*vQC2;
2413       else AliWarning(Form("Determinant == 0 - cent: %d-%d, eta: %f, type: '%c', data: %s, vtx: %d-%d%s",
2414                       Int_t(cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinLowEdge(cBin)),
2415                       Int_t(cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinUpEdge(cBin)),
2416                       cumu.Get(type, 2, CumuHistos::kNUAOld)->GetXaxis()->GetBinCenter(eBin),
2417                       type, fType.Data(), fVzMin, fVzMax,
2418                       GetQCType(fFlags, kTRUE)));
2419       // Go back to v_n for ref. keep <<2'>> for diff. flow).
2420       for (Int_t n = 0; n < fMaxMoment-1; n++) {
2421         Double_t vnTwo = 0;
2422         if (type == 'r' || type == 'R')
2423           vnTwo = (vQC2(n) > 0. ? TMath::Sqrt(vQC2(n)) : 0.);
2424         else {
2425           // is really more <<2'>> in this case
2426           vnTwo = vQC2(n);
2427         }
2428         // Fill in corrected v_n
2429         if (vnTwo != 0) static_cast<TH2D*>(cumu.Get(type, n+2, CumuHistos::kNUA))->Fill(eta, cent, vnTwo);
2430       } // End of moment loop
2431     } // End of eta loop
2432   } // End of centrality loop
2433   return;
2434 }
2435 //_____________________________________________________________________
2436 Double_t AliForwardFlowTaskQC::VertexBin::CalculateNUAMatrixElement(Int_t n, Int_t m, Char_t type, Int_t binA, Int_t cBin) const
2437 {
2438   //
2439   //  Calculates the (n,m)-th element in the NUA matrix: 1 if n == m, otherwise:
2440   //               <<cos[(n-m)phi1]>>*<<cos[(n-m)phi2]>> + <<sin[(n-m)phi1]>>*<<sin[(n-m)phi2]>>
2441   //    NUA(n,m) = -----------------------------------------------------------------------------
2442   //                   1 + <<cos[2nphi1]>>*<<cos[2nphi2]>> + <<sin[2nphi1]>>*<<sin[2nphi2]>>
2443   //
2444   //               <<cos[(n+m)phi1]>>*<<cos[(n+m)phi2]>> + <<sin[(n+m)phi1]>>*<<sin[(n+m)phi2]>>
2445   //             + -----------------------------------------------------------------------------
2446   //                   1 + <<cos[2nphi1]>>*<<cos[2nphi2]>> + <<sin[2nphi1]>>*<<sin[2nphi2]>>
2447   //
2448   //  Parameters:
2449   //   n: row
2450   //   m: coumn
2451   //   type: Reference ('r') or differential ('d') or ('a' or 'b')
2452   //   binA: eta bin of phi1
2453   //   cBin: centrality bin
2454   //
2455   //  Return: NUA(n,m)
2456   //
2457   if (n == m) return 1.;
2458   n += 2;
2459   m += 2;
2460
2461   Double_t cosnMmPhi1 = 0, cosnMmPhi2 = 0, sinnMmPhi1 = 0, sinnMmPhi2 = 0;
2462   Double_t cosnPmPhi1 = 0, cosnPmPhi2 = 0, sinnPmPhi1 = 0, sinnPmPhi2 = 0;
2463   Double_t cos2nPhi1 = 0, cos2nPhi2 = 0, sin2nPhi1 = 0, sin2nPhi2 = 0;
2464
2465   // reference flow
2466   if (type == 'r' || type == 'R') {
2467     if ((fFlags & k3Cor)) {
2468       Double_t eta = fCumuNUARef->GetXaxis()->GetBinCenter(binA);
2469       Int_t i = 0;
2470       while (fEtaLims[i] < eta) i++;
2471       Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2472       GetLimits(i-1, aLow, aHigh, bLow, bHigh);
2473       Double_t multA = 0, multB = 0;
2474       for (Int_t a = aLow; a <= aHigh; a++) {
2475         cosnMmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n-m));
2476         sinnMmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n-m));
2477         cosnPmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n+m));
2478         sinnPmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n+m));
2479         cos2nPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(2*n));
2480         sin2nPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(2*n));
2481         multA += fCumuNUARef->GetBinContent(a, cBin, 0);
2482       }
2483        for (Int_t b = bLow; b <= bHigh; b++) {
2484         cosnMmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n-m));
2485         sinnMmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n-m));
2486         cosnPmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n+m));
2487         sinnPmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n+m));
2488         cos2nPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(2*n));
2489         sin2nPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(2*n));
2490         multB += fCumuNUARef->GetBinContent(b, cBin, 0);
2491       }
2492       if (multA == 0 || multB == 0) {
2493         if (fDebug > 0) AliWarning("multA or multB == 0 in matrix elements, aborting NUA");
2494         return 0.;
2495       }
2496       cosnMmPhi1 /= multA;
2497       sinnMmPhi1 /= multA;
2498       cosnPmPhi1 /= multA;
2499       sinnPmPhi1 /= multA;
2500       cos2nPhi1 /= multA;
2501       sin2nPhi1 /= multA;
2502       cosnMmPhi2 /= multB;
2503       sinnMmPhi2 /= multB;
2504       cosnPmPhi2 /= multB;
2505       sinnPmPhi2 /= multB;
2506       cos2nPhi2 /= multB;
2507       sin2nPhi2 /= multB;
2508     } else {
2509       Int_t binB = fCumuNUARef->GetXaxis()->FindBin(-fCumuNUARef->GetXaxis()->GetBinCenter(binA));
2510       cosnMmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2511       sinnMmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2512       cosnPmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2513       sinnPmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2514       cos2nPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2515       sin2nPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2516       cosnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n-m));
2517       sinnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n-m));
2518       cosnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n+m));
2519       sinnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n+m));
2520       cos2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(2*n));
2521       sin2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(2*n));
2522     }
2523   } // differential flow
2524   else if (type == 'd' || type == 'D') {
2525     Int_t binB = fCumuNUARef->GetXaxis()->FindBin(-fCumuNUADiff->GetXaxis()->GetBinCenter(binA));
2526     cosnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2527     sinnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2528     cosnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2529     sinnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2530     cos2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2531     sin2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2532     cosnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n-m));
2533     sinnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n-m));
2534     cosnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n+m));
2535     sinnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n+m));
2536     cos2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(2*n));
2537     sin2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(2*n));
2538   } // 3 correlator part a or b
2539   else if (type == 'a' || type == 'A' || type == 'b' || type == 'B') {
2540     Double_t mult1 = 0, mult2 = 0;
2541     // POIs
2542     cosnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2543     sinnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2544     cosnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2545     sinnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2546     cos2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2547     sin2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2548     mult1 = fCumuNUADiff->GetBinContent(binA, cBin, 0);
2549     // RPs
2550     Double_t eta = fCumuNUADiff->GetXaxis()->GetBinCenter(binA);
2551     Int_t i = 0;
2552     while (fEtaLims[i] < eta) i++;
2553     Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2554     GetLimits(i-1, aLow, aHigh, bLow, bHigh);
2555     Int_t lLow = ((type == 'a' || type == 'A') ? aLow : bLow);
2556     Int_t lHigh = ((type == 'a' || type == 'A') ? aHigh : bHigh);
2557     for (Int_t l = lLow; l <= lHigh; l++) {
2558       cosnMmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(n-m));
2559       sinnMmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(n-m));
2560       cosnPmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(n+m));
2561       sinnPmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(n+m));
2562       cos2nPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(2*n));
2563       sin2nPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(2*n));
2564       mult2 += fCumuNUARef->GetBinContent(l, cBin, 0);
2565     }
2566     if (mult1 == 0 || mult2 == 0) {
2567       if (fDebug > 0) AliWarning("mult1 or mult2 == 0 in matrix elements, aborting NUA");
2568       return 0.;
2569     }
2570     cosnMmPhi1 /= mult1;
2571     sinnMmPhi1 /= mult1;
2572     cosnPmPhi1 /= mult1;
2573     sinnPmPhi1 /= mult1;
2574     cos2nPhi1 /= mult1;
2575     sin2nPhi1 /= mult1;
2576     cosnMmPhi2 /= mult2;
2577     sinnMmPhi2 /= mult2;
2578     cosnPmPhi2 /= mult2;
2579     sinnPmPhi2 /= mult2;
2580     cos2nPhi2 /= mult2;
2581     sin2nPhi2 /= mult2;
2582   }
2583
2584   // Actual calculation
2585   Double_t e = cosnMmPhi1*cosnMmPhi2 + sinnMmPhi1*sinnMmPhi2 + cosnPmPhi1*cosnPmPhi2 + sinnPmPhi1*sinnPmPhi2;
2586   Double_t den = 1 + cos2nPhi1*cos2nPhi2 + sin2nPhi1*sin2nPhi2;
2587   if (den != 0) e /= den;
2588   else return 0.;
2589
2590   return e;
2591 }
2592 //_____________________________________________________________________
2593 void AliForwardFlowTaskQC::VertexBin::AddVertexBins(CumuHistos& cumu, TList* list, UInt_t nNUA) const
2594 {
2595   //
2596   //  Add up vertex bins with flow results
2597   //
2598   //  Parameters:
2599   //   cumu: CumuHistos object with vtxbin results
2600   //   list: Outout list with added results
2601   //   nNUA: # of NUA histograms to loop over
2602   //
2603   TH2D* vtxHist = 0;
2604   TProfile2D* avgProf = 0;
2605   TString name;
2606   Int_t nT = ((fFlags & k3Cor) ? 3 : 2);
2607   Char_t ct = '\0';
2608   for (UInt_t nua = 0; nua <= nNUA; nua++) { // NUA loop
2609     for (Int_t n = 2; n <= fMaxMoment; n++) { // moment loop
2610       for (Int_t t = 0; t < nT; t++) { // type loop (r/d/a/b)
2611         // Find type
2612         switch (t) {
2613           case 0: ct = 'r'; break;
2614           case 1: ct = ((fFlags & k3Cor) ? 'a' : 'd'); break;
2615           case 2: ct = 'b'; break;
2616           default: ct = '\0'; break;
2617         }
2618         vtxHist = static_cast<TH2D*>(cumu.Get(ct,n,nua));
2619         if (!vtxHist) {
2620           AliWarning("VertexBin::AddVertexBins: vtxHist not found!");
2621           continue;
2622         }
2623         name = vtxHist->GetName();
2624         // Strip name of vtx info
2625         Ssiz_t l = name.Last('x')-3;
2626         name.Resize(l);
2627         avgProf = (TProfile2D*)list->FindObject(name.Data());
2628         // if no output profile yet, make one
2629         if (!avgProf) {
2630           avgProf = new TProfile2D(name.Data(), name.Data(),
2631               vtxHist->GetNbinsX(), vtxHist->GetXaxis()->GetXmin(), vtxHist->GetXaxis()->GetXmax(),
2632               vtxHist->GetNbinsY(), vtxHist->GetYaxis()->GetXmin(), vtxHist->GetYaxis()->GetXmax());
2633           if (vtxHist->GetXaxis()->IsVariableBinSize())
2634             avgProf->GetXaxis()->Set(vtxHist->GetNbinsX(), vtxHist->GetXaxis()->GetXbins()->GetArray());
2635           if (vtxHist->GetYaxis()->IsVariableBinSize())
2636             avgProf->GetYaxis()->Set(vtxHist->GetNbinsY(), vtxHist->GetYaxis()->GetXbins()->GetArray());
2637           list->Add(avgProf);
2638         }
2639         // Fill in, cannot be done with Add function.
2640         for (Int_t e = 1; e <= vtxHist->GetNbinsX(); e++) { // eta loop
2641           Double_t eta = vtxHist->GetXaxis()->GetBinCenter(e);
2642           for (Int_t c = 1; c <= vtxHist->GetNbinsY(); c++) { // cent loop
2643             Double_t cent = vtxHist->GetYaxis()->GetBinCenter(c);
2644             Double_t cont = vtxHist->GetBinContent(e, c);
2645             if (cont == 0) continue;
2646             avgProf->Fill(eta, cent, cont);
2647           } // End of cent loop
2648         } //  End of eta loop
2649       } // End of type loop
2650     } // End of moment loop
2651   } // End of nua loop
2652 }
2653 //_____________________________________________________________________
2654 Int_t AliForwardFlowTaskQC::VertexBin::GetBinNumberCos(Int_t n) const
2655 {
2656   //
2657   //  Get the bin number of <<cos(nphi)>>
2658   //
2659   //  Parameters:
2660   //   n: moment
2661   //
2662   //  Return: bin number
2663   //
2664   Int_t bin = 0;
2665   n = TMath::Abs(n);
2666
2667   if (n == 0) bin = fMaxMoment*4-1;
2668   else        bin = n*2-1;
2669
2670   return bin;
2671 }
2672 //_____________________________________________________________________
2673 Int_t AliForwardFlowTaskQC::VertexBin::GetBinNumberSin(Int_t n) const
2674 {
2675   //
2676   //  Get the bin number of <<sin(nphi)>>
2677   //
2678   //  Parameters:
2679   //   n: moment
2680   //
2681   //  Return: bin number
2682   //
2683   Int_t bin = 0;
2684   n = TMath::Abs(n);
2685
2686   if (n == 0) bin = fMaxMoment*4;
2687   else        bin = n*2;
2688
2689   return bin;
2690 }
2691 //_____________________________________________________________________
2692 void AliForwardFlowTaskQC::VertexBin::SetupNUALabels(TAxis* a) const
2693 {
2694   //
2695   //  Setup NUA labels on axis
2696   //
2697   //  Parameters:
2698   //   a: Axis to set up
2699   //
2700   if (a->GetNbins() != GetBinNumberSin()) AliFatal("SetupNUALabels: Wrong number of bins on axis");
2701
2702   Int_t i = 1, j= 1;
2703   while (i <= a->GetNbins()) {
2704     a->SetBinLabel(i++, Form("<<cos(%d#varphi)>>", j));
2705     a->SetBinLabel(i++, Form("<<sin(%d#varphi)>>", j++));
2706   }
2707
2708   return;
2709 }
2710 //_____________________________________________________________________
2711 TH2I* AliForwardFlowTaskQC::VertexBin::MakeQualityHist(const Char_t* name) const
2712 {
2713   //
2714   //  Add a histogram for checking the analysis quality
2715   //
2716   //  Parameters:
2717   //   const char*: name of data type
2718   //
2719   //  Return: histogram for tracking successful calculations
2720   //
2721   Int_t nBins = ((fFlags & kStdQC) ? 8 : 4);
2722   TH2I* quality = new TH2I(name, name, (fMaxMoment-1)*nBins, 1, (fMaxMoment-1)*nBins+1, fCumuNUARef->GetNbinsY(),
2723                            fCumuNUARef->GetYaxis()->GetXmin(), fCumuNUARef->GetYaxis()->GetXmax());
2724   quality->GetYaxis()->Set(fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXbins()->GetArray());
2725   for (Int_t i = 2, j = 1; i <= fMaxMoment; i++) {
2726     quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{2} > 0", i));
2727     quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{2} <= 0", i));
2728     quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{2} > 0", i));
2729     quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{2} <= 0", i));
2730     if ((fFlags & kStdQC)) {
2731       quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{4} < 0", i));
2732       quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{4} >= 0", i));
2733       quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{4} < 0", i));
2734       quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{4} >= 0", i));
2735     }
2736   }
2737
2738   return quality;
2739 }
2740 //_____________________________________________________________________
2741 TH2D* AliForwardFlowTaskQC::VertexBin::MakeOutputHist(Int_t qc, Int_t n, const Char_t* ctype, UInt_t nua) const
2742 {
2743   //
2744   //  Setup a TH2D for the output
2745   //
2746   //  Parameters:
2747   //   qc   # of particle correlations
2748   //   n    flow moment
2749   //   ref  Type: r/d/a/b
2750   //   nua  For nua corrected hists?
2751   //
2752   //  Return: 2D hist for results
2753   //
2754   Bool_t ref = ((ctype[0] == 'R') || (ctype[0] == 'r'));
2755   TAxis* xAxis = (ref ? fCumuNUARef->GetXaxis() : fCumuNUADiff->GetXaxis());
2756   TAxis* yAxis = (ref ? fCumuNUARef->GetYaxis() : fCumuNUADiff->GetYaxis());
2757   TString ntype = "";
2758   switch (nua) {
2759     case CumuHistos::kNoNUA: ntype = "Un"; break;
2760     case CumuHistos::kNUAOld: ntype = "NUAOld"; break;
2761     case CumuHistos::kNUA: ntype = "NUA"; break;
2762     default: break;
2763   }
2764   TH2D* h = new TH2D(Form("%sQC%d_v%d_%s_%sCorr%s_vtx_%d_%d",
2765                         fType.Data(), qc, n, ctype, ntype.Data(),
2766                         GetQCType(fFlags), fVzMin, fVzMax),
2767                       Form("%sQC%d_v%d_%s_%sCorr%s_vtx_%d_%d",
2768                         fType.Data(), qc, n, ctype, ntype.Data(),
2769                         GetQCType(fFlags), fVzMin, fVzMax),
2770                       xAxis->GetNbins(), xAxis->GetXmin(), xAxis->GetXmax(),
2771                       yAxis->GetNbins(), yAxis->GetXmin(), yAxis->GetXmax());
2772   if (xAxis->IsVariableBinSize()) h->GetXaxis()->Set(xAxis->GetNbins(), xAxis->GetXbins()->GetArray());
2773   h->GetYaxis()->Set(yAxis->GetNbins(), yAxis->GetXbins()->GetArray());
2774
2775   return h;
2776 }
2777 //_____________________________________________________________________
2778 void AliForwardFlowTaskQC::PrintFlowSetup() const
2779 {
2780   //
2781   //  Print the setup of the flow task
2782   //
2783   Printf("=======================================================");
2784   Printf("%s::Print", this->IsA()->GetName());
2785   Printf("Forward detector:                 :\t%s", ((fFlowFlags & kFMD) ? "FMD" : "VZERO"));
2786   Printf("Number of bins in vertex axis     :\t%d", fVtxAxis->GetNbins());
2787   Printf("Range of vertex axis              :\t[%3.1f,%3.1f]",
2788                           fVtxAxis->GetXmin(), fVtxAxis->GetXmax());
2789   Printf("Number of bins in centrality axis :\t%d", fCentAxis->GetNbins());
2790   printf("Centrality binning                :\t");
2791   for (Int_t cBin = 1; cBin <= fCentAxis->GetNbins(); cBin++) {
2792     printf("%02d-%02d%% ", Int_t(fCentAxis->GetBinLowEdge(cBin)), Int_t(fCentAxis->GetBinUpEdge(cBin)));
2793     if (cBin == fCentAxis->GetNbins()) printf("\n");
2794     else if (cBin % 4 == 0) printf("\n\t\t\t\t\t");
2795   }
2796   printf("Doing flow analysis for           :\t");
2797   for (Int_t n  = 2; n <= fMaxMoment; n++) printf("v%d ", n);
2798   printf("\n");
2799   TString type = "Standard QC{2} and QC{4} calculations.";
2800   if ((fFlowFlags & kTPC)) type.ReplaceAll(".", " with TPC tracks for reference.");
2801   if ((fFlowFlags & kEtaGap)) type = "QC{2} with a rapidity gap.";
2802   if ((fFlowFlags & k3Cor))   type = "QC{2} with 3 correlators.";
2803   Printf("QC calculation type               :\t%s", type.Data());
2804   Printf("Symmetrize ref. flow wrt. eta = 0 :\t%s", ((fFlowFlags & kSymEta) ? "true" : "false"));
2805   Printf("Apply NUA correction terms        :\t%s", ((fFlowFlags & kNUAcorr) ? "true" : "false"));
2806   Printf("Satellite vertex flag             :\t%s", ((fFlowFlags & kSatVtx) ? "true" : "false"));
2807   Printf("FMD sigma cut:                    :\t%f", fFMDCut);
2808   Printf("SPD sigma cut:                    :\t%f", fSPDCut);
2809   if ((fFlowFlags & kEtaGap))
2810     Printf("Eta gap:                          :\t%f", fEtaGap);
2811   Printf("=======================================================");
2812 }
2813 //_____________________________________________________________________
2814 const Char_t* AliForwardFlowTaskQC::GetQCType(UShort_t flags, Bool_t prependUS)
2815 {
2816   //
2817   //  Get the type of the QC calculations
2818   //  Used for naming of objects in the VertexBin class,
2819   //  important to avoid memory problems when running multiple
2820   //  initializations of the class with different setups
2821   //
2822   //  Parameters:
2823   //   flags: Flow flags for type determination
2824   //   prependUS: Prepend an underscore (_) to the name
2825   //
2826   //  Return: QC calculation type
2827   //
2828   TString type = "";
2829   if      ((flags & kStdQC))  type = "StdQC";
2830   else if ((flags & kEtaGap)) type = "EtaGap";
2831   else if ((flags & k3Cor))   type = "3Cor";
2832   else type = "UNKNOWN";
2833   if (prependUS) type.Prepend("_");
2834   if ((flags & kTPC)) type.Append("TPCTr");
2835   return type.Data();
2836 }
2837 //_____________________________________________________________________
2838 TH1* AliForwardFlowTaskQC::CumuHistos::Get(Char_t t, Int_t n, UInt_t nua) const
2839 {
2840   //
2841   //  Get histogram/profile for type t and moment n
2842   //
2843   //  Parameters:
2844   //   t: type = 'r'/'d'
2845   //   n: moment
2846   //   nua: NUA type
2847   //
2848   n = GetPos(n, nua);
2849   if (n < 0) AliFatal(Form("CumuHistos out of range: (%c,%d)", t, n));
2850
2851   TH1* h = 0;
2852   Int_t pos = -1;
2853   if      (t == 'r' || t == 'R') pos = n;
2854   else if (t == 'd' || t == 'D') pos = n;
2855   else if (t == 'a' || t == 'A') pos = 2*n;
2856   else if (t == 'b' || t == 'B') pos = 2*n+1;
2857   else AliFatal(Form("CumuHistos wrong type: %c", t));
2858
2859   if (t == 'r' || t == 'R') {
2860     if (pos < fRefHists->GetEntries()) {
2861       h = (TH1*)fRefHists->At(pos);
2862     }
2863   } else if (pos < fDiffHists->GetEntries()) {
2864     h = (TH1*)fDiffHists->At(pos);
2865   }
2866   if (!h) AliFatal(Form("No hist found in list %c at pos %d", t, pos));
2867
2868   return h;
2869 }
2870 //_____________________________________________________________________
2871 void AliForwardFlowTaskQC::CumuHistos::ConnectList(TString name, TList* l)
2872 {
2873   //
2874   //  Connect an output list with this object
2875   //
2876   //  Parameters:
2877   //   name: base name
2878   //   l: list to keep outputs in
2879   //
2880   TString ref = name;
2881   ref.ReplaceAll("Cumu","CumuRef");
2882   fRefHists = (TList*)l->FindObject(ref.Data());
2883   if (!fRefHists) {
2884     fRefHists = new TList();
2885     fRefHists->SetName(ref.Data());
2886     l->Add(fRefHists);
2887   } else {
2888     // Check that the list correspond to fMaxMoments
2889     if (fRefHists->GetEntries() != (fMaxMoment-1.)*(fNUA+1))
2890       AliError("CumuHistos::ConnectList Wrong number of hists in ref list,"
2891                "you are doing something wrong!");
2892   }
2893   TString diff = name;
2894   diff.ReplaceAll("Cumu","CumuDiff");
2895   fDiffHists = (TList*)l->FindObject(diff.Data());
2896   if (!fDiffHists) {
2897     fDiffHists = new TList();
2898     fDiffHists->SetName(diff.Data());
2899     l->Add(fDiffHists);
2900   } else {
2901     // Check that the list correspond to fMaxMoment
2902     if ((fDiffHists->GetEntries() != (fMaxMoment-1.)*(fNUA+1)) &&
2903         (fDiffHists->GetEntries() != 2*(fMaxMoment-1.)*(fNUA+1)))
2904       AliError(Form("CumuHistos::ConnectList Wrong number of hists in diff list,"
2905                "you are doing something wrong! (%s)",name.Data()));
2906   }
2907
2908 }
2909 //_____________________________________________________________________
2910 void AliForwardFlowTaskQC::CumuHistos::Add(TH1* h) const
2911 {
2912   //
2913   //  Add a histogram to this objects lists
2914   //
2915   //  Parameters:
2916   //   h: histogram/profile to add
2917   //
2918   TString name = h->GetName();
2919   if (name.Contains("Ref")) {
2920     if (fRefHists) fRefHists->Add(h);
2921     else AliFatal("CumuHistos::Add() - fRefHists does not exist");
2922     // Check that the list correspond to fMaxMoments
2923     if (fRefHists->GetEntries() > (fNUA+1)*(fMaxMoment-1.))
2924       AliError("CumuHistos::Add wrong number of hists in ref list, "
2925                "you are doing something wrong!");
2926   }
2927   else if (name.Contains("Diff")) {
2928     if (fDiffHists) fDiffHists->Add(h);
2929     else AliFatal("CumuHistos::Add() - fDiffHists does not exist");
2930     // Check that the list correspond to fMaxMoment
2931     if (fDiffHists->GetEntries() > 2*(fNUA+1)*(fMaxMoment-1.))
2932       AliError("CumuHistos::Add wrong number of hists in diff list, "
2933          "you are doing something wrong!");
2934   }
2935   return;
2936 }
2937 //_____________________________________________________________________
2938 Int_t AliForwardFlowTaskQC::CumuHistos::GetPos(Int_t n, UInt_t nua) const
2939 {
2940   //
2941   //  Get position in list of moment n objects
2942   //  To take care of different numbering schemes
2943   //
2944   //  Parameters:
2945   //   n: moment
2946   //   nua: # of nua corrections applied
2947   //
2948   //  Return: position #
2949   //
2950   if (n > fMaxMoment) return -1;
2951   else return (n-2)+nua*(fMaxMoment-1);
2952 }
2953 //_____________________________________________________________________
2954 //
2955 //
2956 // EOF