PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.cxx

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