PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.cxx

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