PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.cxx

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