PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.cxx

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