PWGLF/FORWARD/analysis2/AliForwardFlowTaskQC.cxx

   1 //
   2 // Calculate flow in the forward and central regions using the Q cumulants method.
   3 //
   4 // Inputs:
   5 //  - AliAODEvent
   6 //
   7 // Outputs:
   8 //  - AnalysisResults.root
   9 //
  10 #include <TROOT.h>
  11 #include <TSystem.h>
  12 #include <TInterpreter.h>
  13 #include <TChain.h>
  14 #include <TFile.h>
  15 #include <TList.h>
  16 #include <iostream>
  17 #include <TMath.h>
  18 #include <TH3D.h>
  19 #include <TProfile2D.h>
  20 #include <TParameter.h>
  21 #include <TGraph.h>
  22 #include "AliLog.h"
  23 #include "AliForwardFlowTaskQC.h"
  24 #include "AliAnalysisManager.h"
  25 #include "AliAODHandler.h"
  26 #include "AliAODInputHandler.h"
  27 #include "AliAODForwardMult.h"
  28 #include "AliAODCentralMult.h"
  29 #include "AliAODEvent.h"
  30 #include "AliForwardUtil.h"
  31
  32 ClassImp(AliForwardFlowTaskQC)
  33 #if 0
  34 ; // For emacs
  35 #endif
  36
  37 AliForwardFlowTaskQC::AliForwardFlowTaskQC()
  38   : AliAnalysisTaskSE(),
  39     fVtxAxis(),         // Axis to contorl vertex binning
  40     fFMDCut(-1),        // FMD sigma cut
  41     fSPDCut(-1),        // SPD sigma cut
  42     fFlowFlags(kSymEta),// Flow flags
  43     fEtaGap(2.),        // Eta gap value
  44     fBinsFMD(),         // List with FMD flow histos
  45     fBinsSPD(),         // List with SPD flow histos
  46     fSumList(0),        // Event sum list
  47     fOutputList(0),     // Result output list
  48     fAOD(0),            // AOD input event
  49     fV(),               // Flow moments
  50     fVtx(1111),         // Z vertex coordinate
  51     fCent(-1),          // Centrality
  52     fHistCent(),        // Histo for centrality
  53     fHistVertexSel()    // Histo for selected vertices
  54 {
  55   //
  56   // Default constructor
  57   //
  58 }
  59 //_____________________________________________________________________
  60 AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name)
  61   : AliAnalysisTaskSE(name),
  62     fVtxAxis(),         // Axis to contorl vertex binning
  63     fFMDCut(-1),        // FMD sigma cut
  64     fSPDCut(-1),        // SPD sigma cut
  65     fFlowFlags(kSymEta),// Flow flags
  66     fEtaGap(2.),        // Eta gap value
  67     fBinsFMD(),         // List with FMD flow histos
  68     fBinsSPD(),         // List with SPD flow histos
  69     fSumList(0),        // Event sum list
  70     fOutputList(0),     // Result output list
  71     fAOD(0),            // AOD input event
  72     fV(),               // Flow moments
  73     fVtx(1111),         // Z vertex coordinate
  74     fCent(-1),          // Centrality
  75     fHistCent(),        // Histo for centrality
  76     fHistVertexSel()    // Histo for selected vertices
  77 {
  78   //
  79   // Constructor
  80   //
  81   // Parameters:
  82   //  name: Name of task
  83   //
  84   DefineOutput(1, TList::Class());
  85   DefineOutput(2, TList::Class());
  86 }
  87 //_____________________________________________________________________
  88 AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o)
  89   : AliAnalysisTaskSE(o),
  90     fVtxAxis(o.fVtxAxis),              // Axis to contorl vertex binning
  91     fFMDCut(o.fFMDCut),                // FMD sigma cut
  92     fSPDCut(o.fSPDCut),                // SPD sigma cut
  93     fFlowFlags(o.fFlowFlags),          // Flow flags
  94     fEtaGap(o.fEtaGap),                // Eta gap value
  95     fBinsFMD(),                        // List with FMD flow histos
  96     fBinsSPD(),                        // List with SPD flow histos
  97     fSumList(o.fSumList),              // Event sum list
  98     fOutputList(o.fOutputList),        // Result output list
  99     fAOD(o.fAOD),                      // AOD input event
 100     fV(o.fV),                          // Flow moments
 101     fVtx(o.fVtx),                      // Z vertex coordinate
 102     fCent(o.fCent),                    // Centrality
 103     fHistCent(o.fHistCent),            // Histo for centrality
 104     fHistVertexSel(o.fHistVertexSel)   // Histo for selected vertices
 105 {
 106   //
 107   // Copy constructor
 108   //
 109   // Parameters:
 110   //    o Object to copy from
 111   //
 112 }
 113 //_____________________________________________________________________
 114 AliForwardFlowTaskQC&
 115 AliForwardFlowTaskQC::operator=(const AliForwardFlowTaskQC& o)
 116 {
 117   //
 118   // Assignment operator
 119   //
 120   if (&o == this) return *this;
 121   fVtxAxis       = o.fVtxAxis;
 122   fFMDCut        = o.fFMDCut;
 123   fSPDCut        = o.fSPDCut;
 124   fFlowFlags     = o.fFlowFlags;
 125   fEtaGap        = o.fEtaGap;
 126   fSumList       = o.fSumList;
 127   fOutputList    = o.fOutputList;
 128   fAOD           = o.fAOD;
 129   fV             = o.fV;
 130   fVtx           = o.fVtx;
 131   fCent          = o.fCent;
 132   fHistCent      = o.fHistCent;
 133   fHistVertexSel = o.fHistVertexSel;
 134
 135   return *this;
 136 }
 137 //_____________________________________________________________________
 138 void AliForwardFlowTaskQC::UserCreateOutputObjects()
 139 {
 140   //
 141   // Create output objects
 142   //
 143   InitVertexBins();
 144   InitHists();
 145   PrintFlowSetup();
 146
 147   PostData(1, fSumList);
 148   PostData(2, fOutputList);
 149
 150 }
 151 //_____________________________________________________________________
 152 void AliForwardFlowTaskQC::InitVertexBins()
 153 {
 154   //
 155   // Init vertexbin objects for FMD and SPD, and add them to the lists
 156   //
 157   Int_t moment = 0;
 158   for(UShort_t n = 0; n < fV.GetSize(); n++) {
 159     moment = fV.At(n);
 160     for (Int_t v = 1; v <= fVtxAxis->GetNbins(); v++) {
 161       Int_t vL = Int_t(fVtxAxis->GetBinLowEdge(v));
 162       Int_t vH = Int_t(fVtxAxis->GetBinUpEdge(v));
 163       fBinsFMD.Add(new VertexBin(vL, vH, moment, "FMD", fFlowFlags, fFMDCut, fEtaGap));
 164       fBinsSPD.Add(new VertexBin(vL, vH, moment, "SPD", fFlowFlags, fSPDCut, fEtaGap));
 165     }
 166   }
 167 }
 168 //_____________________________________________________________________
 169 void AliForwardFlowTaskQC::InitHists()
 170 {
 171   //
 172   // Init histograms and add vertex bin histograms to the sum list
 173   //
 174
 175   if (!fSumList)
 176     fSumList = new TList();
 177   fSumList->SetName("Sums");
 178   fSumList->SetOwner();
 179
 180   if (!fVtxAxis) fVtxAxis = new TAxis(20, -10, 10);
 181   fVtxAxis->SetName("VtxAxis");
 182   fHistCent         = new TH1D("hCent", "Centralities", 100, 0, 100);
 183   fHistVertexSel    = new TH1D("hVertexSel", "Selected vertices", fVtxAxis->GetNbins(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax());
 184
 185   TList* dList = new TList();
 186   dList->SetName("Diagnostics");
 187 //  dList->Add(fVtxAxis);
 188   dList->Add(fHistCent);
 189   dList->Add(fHistVertexSel);
 190   fSumList->Add(dList);
 191
 192   TIter nextFMD(&fBinsFMD);
 193   VertexBin* bin = 0;
 194   while ((bin = static_cast<VertexBin*>(nextFMD()))) {
 195     bin->AddOutput(fSumList);
 196   }
 197   TIter nextSPD(&fBinsSPD);
 198   while ((bin = static_cast<VertexBin*>(nextSPD()))) {
 199     bin->AddOutput(fSumList);
 200   }
 201 }
 202 //_____________________________________________________________________
 203 void AliForwardFlowTaskQC::UserExec(Option_t */*option*/)
 204 {
 205   //
 206   // Calls the analyze function - called every event
 207   //
 208   // Parameters:
 209   //  option: Not used
 210   //
 211
 212   Analyze();
 213
 214   PostData(1, fSumList);
 215
 216   return;
 217 }
 218 //_____________________________________________________________________
 219 Bool_t AliForwardFlowTaskQC::Analyze()
 220 {
 221   //
 222   // Load FMD and SPD objects from aod tree and call the cumulants
 223   // calculation for the correct vertexbin
 224   //
 225
 226   // Reset data members
 227   fCent = -1;
 228   fVtx = 1111;
 229
 230   // Get input event
 231 //  fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 232   fAOD = dynamic_cast<AliAODEvent*>(AliForwardUtil::GetAODEvent(this));
 233   if (!fAOD) return kFALSE;
 234
 235   const AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("Forward"));
 236   const AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClusters"));
 237   if (!aodfmult) return kFALSE;
 238
 239   // Check event for triggers, get centrality, vtx etc.
 240   if (!CheckEvent(aodfmult)) return kFALSE;
 241   Int_t vtx = fVtxAxis->FindBin(fVtx)-1;
 242
 243   // If everything is OK: get histos and run analysis
 244   const TH2D& fmddNdetadphi = aodfmult->GetHistogram();
 245   if ((fFlowFlags & kEtaGap)) {
 246     FillVtxBinListEtaGap(fBinsFMD, fmddNdetadphi, fmddNdetadphi, vtx);
 247   } else {
 248     FillVtxBinList(fBinsFMD, fmddNdetadphi, vtx);
 249   }
 250
 251   if (aodcmult) {
 252     const TH2D& spddNdetadphi = aodcmult->GetHistogram();
 253     if ((fFlowFlags & kEtaGap)) {
 254       FillVtxBinListEtaGap(fBinsSPD, fmddNdetadphi, spddNdetadphi, vtx);
 255     } else {
 256       FillVtxBinList(fBinsSPD, spddNdetadphi, vtx);
 257     }
 258   }
 259
 260   return kTRUE;
 261 }
 262 //_____________________________________________________________________
 263 Bool_t AliForwardFlowTaskQC::FillVtxBinList(const TList& list, const TH2D& h, Int_t vtx) const
 264 {
 265   //
 266   // Loops over list of VtxBins, fills hists of bins for current vertex
 267   // and runs analysis on those bins
 268   //
 269   // Parameters:
 270   //  list: list of VtxBins
 271   //  h:    dN/detadphi histogram
 272   //  vBin: current vertex bin
 273   //
 274   // return true on success
 275   //
 276   VertexBin* bin = 0;
 277   Int_t i = 0;
 278   Int_t nVtxBins = fVtxAxis->GetNbins();
 279
 280   while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
 281     Bool_t skipFourP = !bin->FillHists(h, fCent, kFillBoth);
 282     bin->CumulantsAccumulate(fCent, skipFourP);
 283     i++;
 284   }
 285
 286   return kTRUE;
 287 }
 288 //_____________________________________________________________________
 289 Bool_t AliForwardFlowTaskQC::FillVtxBinListEtaGap(const TList& list, const TH2D& href,
 290                                                   const TH2D& hdiff, Int_t vtx) const
 291 {
 292   //
 293   // Loops over list of VtxBins, fills hists of bins for current vertex
 294   // and runs analysis on those bins
 295   //
 296   // Parameters:
 297   //  list: list of VtxBins
 298   //  h:    dN/detadphi histogram
 299   //  vBin: current vertex bin
 300   //
 301   // return true on success
 302   //
 303   VertexBin* bin = 0;
 304   Int_t i = 0;
 305   Int_t nVtxBins = fVtxAxis->GetNbins();
 306
 307   while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
 308     bin->FillHists(href, fCent, kFillRef);
 309     bin->FillHists(hdiff, fCent, kFillDiff);
 310     bin->CumulantsAccumulate(fCent);
 311     i++;
 312   }
 313
 314   return kTRUE;
 315 }
 316 //_____________________________________________________________________
 317 void AliForwardFlowTaskQC::Terminate(Option_t */*option*/)
 318 {
 319   //
 320   // Calls the finalize function, done at the end of the analysis
 321   //
 322   // Parameters:
 323   //  option: Not used
 324   //
 325
 326   // Make sure pointers are set to the correct lists
 327   fSumList = dynamic_cast<TList*> (GetOutputData(1));
 328   if(!fSumList) {
 329     AliError("Could not retrieve TList fSumList");
 330     return;
 331   }
 332   if (!fOutputList)
 333     fOutputList = new TList();
 334   fOutputList->SetName("Results");
 335   fOutputList->SetOwner();
 336
 337   if ((fFlowFlags & kEtaGap)) {
 338     TParameter<Double_t>* etaGap = new TParameter<Double_t>("EtaGap", fEtaGap);
 339     fOutputList->Add(etaGap);
 340   }
 341
 342   // We make summary histograms of accepted events.
 343   TList* list = 0;
 344   TH3D* hist = 0;
 345   TH1D* cent = 0;
 346   TObject* o = 0;
 347   for (Int_t i = 1; i < fSumList->GetEntries(); i++) {
 348     list = dynamic_cast<TList*>(fSumList->At(i));
 349     if (!list) continue;
 350     hist = dynamic_cast<TH3D*>(list->At(0));
 351     if (!hist) continue;
 352     const char* histname = hist->GetName();
 353     TString name = "";
 354     for (Int_t j = 0; ; j++) {
 355       if (histname[j] == 'v') break;
 356       name += histname[j];
 357     }
 358     if ((fFlowFlags & kEtaGap)) name += "_etaGap";
 359     cent = (TH1D*)fOutputList->FindObject(name.Data());
 360     if (!cent) {
 361       cent = new TH1D(name.Data(), name.Data(), hist->GetNbinsY(), hist->GetYaxis()->GetXmin(), hist->GetYaxis()->GetXmax());
 362       cent->GetXaxis()->Set(hist->GetNbinsY(), hist->GetYaxis()->GetXbins()->GetArray());
 363       fOutputList->Add(cent);
 364     }
 365     for (Int_t k = 1; k <= hist->GetNbinsY(); k++) {
 366       Double_t centrality = hist->GetYaxis()->GetBinCenter(k);
 367       Double_t events = hist->GetBinContent(0, k, 0);
 368       cent->Fill(centrality, events);
 369     }
 370     o = fOutputList->FindObject(Form("hQCQuality%s", name.Data()));
 371     if (!o) MakeQualityHist(name);
 372   }
 373
 374   // Run finalize on VertexBins
 375   Finalize();
 376
 377   // Collect centralities
 378   MakeCentralityHists(fOutputList);
 379   TList* vtxList = (TList*)fOutputList->FindObject("vtxList");
 380   if (vtxList) MakeCentralityHists(vtxList);
 381   TList* nuaList = (TList*)fOutputList->FindObject("NUATerms");
 382   if (nuaList) MakeCentralityHists(nuaList);
 383
 384   PostData(2, fOutputList);
 385
 386   return;
 387 }
 388 //_____________________________________________________________________
 389 void AliForwardFlowTaskQC::AddFlowMoment(Short_t n)
 390 {
 391     //
 392     // Add a flow moment to be calculated
 393     //
 394     if (n > 10) AliFatal(Form("Too big moment added: %d (bug)", n));
 395     Int_t size = fV.GetSize();
 396     fV.Set(size+1);
 397     fV.AddAt(n, size);
 398 }
 399 //_____________________________________________________________________
 400 void AliForwardFlowTaskQC::Finalize()
 401 {
 402   //
 403   // Finalize command, called by Terminate()
 404   //
 405
 406   // Reinitiate vertex bins if Terminate is called separately!
 407   if (fBinsFMD.GetEntries() == 0) InitVertexBins();
 408
 409   // Iterate over all vertex bins objects and finalize cumulants
 410   // calculations
 411   EndVtxBinList(fBinsFMD);
 412   EndVtxBinList(fBinsSPD);
 413
 414   return;
 415 }
 416 //_____________________________________________________________________
 417 void AliForwardFlowTaskQC::EndVtxBinList(const TList& list) const
 418 {
 419   //
 420   // Loop over VertexBin list and call terminate on each
 421   //
 422   // Parameters:
 423   //  list VertexBin list
 424   //
 425   TIter next(&list);
 426   VertexBin* bin = 0;
 427   while ((bin = static_cast<VertexBin*>(next()))) {
 428     bin->CumulantsTerminate(fSumList, fOutputList);
 429   }
 430   return;
 431 }
 432 // _____________________________________________________________________
 433 void AliForwardFlowTaskQC::MakeCentralityHists(TList* list)
 434 {
 435   //
 436   // Loop over a list containing a TProfile2D with flow results and project
 437   // and project to TH1's in specific centrality bins
 438   //
 439   // Parameters:
 440   //  list Flow results list
 441   //
 442   TProfile2D* hist2D = 0;
 443   TList* centList = 0;
 444   TH1D* hist1D = 0;
 445   TObject* helper = 0;
 446   TIter nextProfile(list);
 447   while ((helper = dynamic_cast<TObject*>(nextProfile()))) {
 448     if (!(hist2D = dynamic_cast<TProfile2D*>(helper))) continue;
 449     for (Int_t cBin = 1; cBin <= hist2D->GetNbinsY(); cBin++) {
 450       Int_t cMin = Int_t(hist2D->GetYaxis()->GetBinLowEdge(cBin));
 451       Int_t cMax = Int_t(hist2D->GetYaxis()->GetBinUpEdge(cBin));
 452       TString name = Form("cent_%d-%d", cMin, cMax);
 453       centList = (TList*)list->FindObject(name.Data());
 454       if (!centList) {
 455         centList = new TList();
 456         centList->SetName(name.Data());
 457         list->Add(centList);
 458       }
 459       hist1D = hist2D->ProjectionX(Form("%s_%s", hist2D->GetName(), name.Data()),
 460                                          cBin, cBin, "E");
 461       hist1D->SetTitle(hist1D->GetName());
 462       centList->Add(hist1D);
 463     }
 464   }
 465 }
 466 // _____________________________________________________________________
 467 Bool_t AliForwardFlowTaskQC::CheckEvent(const AliAODForwardMult* aodfm)
 468 {
 469   //
 470   // Function to check that and AOD event meets the cuts
 471   //
 472   // Parameters:
 473   //  AliAODForwardMult: forward mult object with trigger and vertex info
 474   //
 475   // Returns false if there is no trigger or if the centrality or vertex
 476   // is out of range. Otherwise true.
 477   //
 478
 479   // First check for trigger
 480   if (!CheckTrigger(aodfm)) return kFALSE;
 481
 482   // Then check for centrality
 483   if (!GetCentrality(aodfm)) return kFALSE;
 484
 485   // And finally check for vertex
 486   if (!GetVertex(aodfm)) return kFALSE;
 487
 488   // Ev. accepted - filling diag. hists
 489   fHistCent->Fill(fCent);
 490   fHistVertexSel->Fill(fVtx);
 491
 492   return kTRUE;
 493 }
 494 // _____________________________________________________________________
 495 Bool_t AliForwardFlowTaskQC::CheckTrigger(const AliAODForwardMult* aodfm) const
 496 {
 497   //
 498   // Function to look for a trigger string in the event.
 499   //
 500   // Parameters:
 501   //  AliAODForwardMult: forward mult object with trigger and vertex info
 502   //
 503   // Returns true if offline trigger is present
 504   //
 505   return aodfm->IsTriggerBits(AliAODForwardMult::kOffline);
 506 }
 507 // _____________________________________________________________________
 508 Bool_t AliForwardFlowTaskQC::GetCentrality(const AliAODForwardMult* aodfm)
 509 {
 510   //
 511   // Function to look get centrality of the event.
 512   //
 513   // Parameters:
 514   //  AliAODForwardMult: forward mult object with trigger and vertex info
 515   //
 516   // Returns true if centrality determination is present
 517   //
 518   if (aodfm->HasCentrality()) {
 519     fCent = (Double_t)aodfm->GetCentrality();
 520     if (0. >= fCent || fCent >= 100.) return kFALSE;
 521   }
 522   else fCent = 97.5;
 523
 524   return kTRUE;
 525 }
 526 // _____________________________________________________________________
 527 Bool_t AliForwardFlowTaskQC::GetVertex(const AliAODForwardMult* aodfm)
 528 {
 529   //
 530   // Function to look for vertex determination in the event.
 531   //
 532   // Parameters:
 533   //  AliAODForwardMult: forward mult object with trigger and vertex info
 534   //
 535   // Returns true if vertex is determined
 536   //
 537   fVtx = aodfm->GetIpZ();
 538   if (fVtx < fVtxAxis->GetXmin() || fVtx > fVtxAxis->GetXmax()) return kFALSE;
 539
 540   return kTRUE;
 541 }
 542 //_____________________________________________________________________
 543 void AliForwardFlowTaskQC::MakeQualityHist(const Char_t* name) const {
 544
 545   TH1I* quality = new TH1I(Form("hQCQuality%s", name),
 546                       Form("hQCQuality%s", name),
 547                       fV.GetSize()*8, 1, fV.GetSize()*8+1);
 548   for (Int_t i = 0, j = 1; i < fV.GetSize(); i++) {
 549     quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{2} > 0", fV.At(i)));
 550     quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{2} <= 0", fV.At(i)));
 551     quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{2} > 0", fV.At(i)));
 552     quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{2} <= 0", fV.At(i)));
 553     quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{4} < 0", fV.At(i)));
 554     quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{4} >= 0", fV.At(i)));
 555     quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{4} < 0", fV.At(i)));
 556     quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{4} >= 0", fV.At(i)));
 557   }
 558
 559   fOutputList->Add(quality);
 560 }
 561 //_____________________________________________________________________
 562 AliForwardFlowTaskQC::VertexBin::VertexBin()
 563   : TNamed(),
 564     fMoment(0),      // Flow moment for this vertexbin
 565     fVzMin(0),       // Vertex z-coordinate min
 566     fVzMax(0),       // Vertex z-coordinate max
 567     fType(),         // Data type name e.g., FMD/SPD/FMDTR/SPDTR/MC
 568     fFlags(0),       // Use forward-backward symmetry, if detector allows it
 569     fSigmaCut(-1),   // Sigma cut to remove outlier events
 570     fEtaGap(2.),     // Eta gap value
 571     fCumuRef(),      // Histogram for reference flow
 572     fCumuDiff(),     // Histogram for differential flow
 573     fCumuHist(),     // Sum histogram for cumulants
 574     fdNdedpAcc(),    // Diagnostics histogram to make acc. maps
 575     fOutliers(),     // Histogram for sigma distribution
 576     fDebug()         // Debug level
 577 {
 578   //
 579   // Default constructor
 580   //
 581 }
 582 //_____________________________________________________________________
 583 AliForwardFlowTaskQC::VertexBin::VertexBin(Int_t vLow, Int_t vHigh,
 584                                            UShort_t moment, TString name,
 585                                            UShort_t flags, Double_t cut,
 586                                            Double_t etaGap)
 587   : TNamed("", ""),
 588     fMoment(moment),  // Flow moment for this vertexbin
 589     fVzMin(vLow),     // Vertex z-coordinate min
 590     fVzMax(vHigh),    // Vertex z-coordinate max
 591     fType(name),      // Data type name e.g., FMD/SPD/FMDTR/SPDTR/MC
 592     fFlags(flags),    // Use forward-backward symmetry, if detector allows it
 593     fSigmaCut(cut),   // Sigma cut to remove outlier events
 594     fEtaGap(etaGap),  // Eta gap value
 595     fCumuRef(),       // Histogram for reference flow
 596     fCumuDiff(),      // Histogram for differential flow
 597     fCumuHist(),      // Sum histogram for cumulants
 598     fdNdedpAcc(),     // Diagnostics histogram to make acc. maps
 599     fOutliers(),      // Histogram for sigma distribution
 600     fDebug(0)         // Debug level
 601 {
 602   //
 603   // Constructor
 604   //
 605   // Parameters
 606   //  vLow: min z-coordinate
 607   //  vHigh: max z-coordinate
 608   //  moment: flow moment
 609   //  name: data type name (FMD/SPD/FMDTR/SPDTR/MC)
 610   //  sym: data is symmetric in eta
 611   //
 612   fType.ToUpper();
 613
 614   SetName(Form("%svertexBin%d_%d_%d%s", fType.Data(), moment, vLow, vHigh, ((fFlags & kEtaGap) ? "_etaGap" : "")));
 615   SetTitle(Form("%svertexBin%d_%d_%d%s", fType.Data(), moment, vLow, vHigh, ((fFlags & kEtaGap) ? "_etaGap" : "")));
 616
 617   fDebug = AliAnalysisManager::GetAnalysisManager()->GetDebugLevel();
 618   if (fDebug > 0) Printf("AliForwardFlowTaskQC::VertexBin()\tDebugMode: %d", fDebug);
 619 }
 620 //_____________________________________________________________________
 621 AliForwardFlowTaskQC::VertexBin&
 622 AliForwardFlowTaskQC::VertexBin::operator=(const AliForwardFlowTaskQC::VertexBin& o)
 623 {
 624   //
 625   // Assignment operator
 626   //
 627   // Parameters
 628   //  o: AliForwardFlowTaskQC::VertexBin
 629   //
 630   if (&o == this) return *this;
 631   fType         = o.fType;
 632   fCumuRef      = o.fCumuRef;
 633   fCumuDiff     = o.fCumuDiff;
 634   fCumuHist     = o.fCumuHist;
 635   fdNdedpAcc    = o.fdNdedpAcc;
 636   fOutliers     = o.fOutliers;
 637   fDebug        = o.fDebug;
 638
 639   return *this;
 640 }
 641 //_____________________________________________________________________
 642 void AliForwardFlowTaskQC::VertexBin::AddOutput(TList* outputlist)
 643 {
 644   //
 645   // Add histograms to outputlist
 646   //
 647   // Parameters
 648   //  outputlist: list of histograms
 649   //
 650
 651   // First we try to find an outputlist for this vertexbin
 652   TList* list = (TList*)outputlist->FindObject(Form("%svertex_%d_%d%s", fType.Data(), fVzMin, fVzMax, ((fFlags & kEtaGap) ? "_etaGap" : "")));
 653
 654   // If it doesn't exist we make one
 655   if (!list) {
 656     list = new TList();
 657     list->SetName(Form("%svertex_%d_%d%s", fType.Data(), fVzMin, fVzMax, ((fFlags & kEtaGap) ? "_etaGap" : "")));
 658     outputlist->Add(list);
 659   }
 660
 661   // We initiate the reference histogram
 662   fCumuRef = new TH2D(Form("%s_v%d_%d_%d%s_ref", fType.Data(), fMoment, fVzMin, fVzMax, ((fFlags & kEtaGap) ? "_etaGap" : "")),
 663                         Form("%s_v%d_%d_%d%s_ref", fType.Data(), fMoment, fVzMin, fVzMax, ((fFlags & kEtaGap) ? "_etaGap" : "")),
 664                         2, -6., 6., 8, 0.5, 8.5);
 665   fCumuRef->Sumw2();
 666   //list->Add(fCumuRef);
 667
 668   // We initiate the differential histogram
 669   fCumuDiff = new TH2D(Form("%s_v%d_%d_%d%s_diff", fType.Data(), fMoment, fVzMin, fVzMax, ((fFlags & kEtaGap) ? "_etaGap" : "")),
 670                        Form("%s_v%d_%d_%d%s_diff", fType.Data(), fMoment, fVzMin, fVzMax, ((fFlags & kEtaGap) ? "_etaGap" : "")),
 671                        48, -6., 6., 8, 0.5, 8.5);
 672   fCumuDiff->Sumw2();
 673   //list->Add(fCumuDiff);
 674
 675   // Initiate the cumulant sum histogram
 676   fCumuHist = new TH3D(Form("%sv%d_vertex_%d_%d%s_cumu", fType.Data(), fMoment, fVzMin, fVzMax, ((fFlags & kEtaGap) ? "_etaGap" : "")),
 677                        Form("%sv%d_vertex_%d_%d%s_cumu", fType.Data(), fMoment, fVzMin, fVzMax, ((fFlags & kEtaGap) ? "_etaGap" : "")),
 678                        48, -6., 6., 20, 0., 100., 29, 0.5, 29.5);
 679   fCumuHist->Sumw2();
 680   SetupCentAxis(fCumuHist->GetYaxis());
 681
 682   list->Add(fCumuHist);
 683
 684   // We check for diagnostics histograms (only done per type and moment, not vertexbin)
 685   // If they are not found we create them.
 686   TList* dList = (TList*)outputlist->FindObject("Diagnostics");
 687   if (!dList) AliFatal("No diagnostics list found, what kind of game are you running here?!?!");
 688
 689   // Acceptance hists are shared over all moments
 690   fdNdedpAcc = (TH2F*)dList->FindObject(Form("h%sdNdedpAcc_%d_%d%s", fType.Data(), fVzMin, fVzMax, ((fFlags & kEtaGap) ? "_etaGap" : "")));
 691   if (!fdNdedpAcc) {
 692     fdNdedpAcc = new TH2F(Form("h%sdNdedpAcc_%d_%d%s", fType.Data(), fVzMin, fVzMax, ((fFlags & kEtaGap) ? "_etaGap" : "")),
 693                           Form("%s acceptance map for %d cm < v_{z} < %d cm", fType.Data(), fVzMin, fVzMax),
 694                           48, -6, 6, 20, 0, TMath::TwoPi());
 695     fdNdedpAcc->Sumw2();
 696     dList->Add(fdNdedpAcc);
 697   }
 698
 699   if (!(fFlags & kEtaGap)) {
 700     fOutliers = new TH2F(Form("hOutliers_%s_v%d_%d_%d", fType.Data(), fMoment, fVzMin, fVzMax),
 701                         Form("Maximum #sigma from mean N_{ch} pr. bin - %s v_{%d}, %d < v_{z} < %d",
 702                         fType.Data(), fMoment, fVzMin, fVzMax),
 703                         20, 0., 100., 500, 0., (fType.Contains("MC") ? 15. : 5.));
 704     dList->Add(fOutliers);
 705   }
 706
 707 }
 708 //_____________________________________________________________________
 709 Bool_t AliForwardFlowTaskQC::VertexBin::FillHists(const TH2D& dNdetadphi, Double_t cent, EFillFlow mode)
 710 {
 711   //
 712   // Fill reference and differential eta-histograms
 713   //
 714   // Parameters:
 715   //  dNdetadphi: 2D histogram with input data
 716   //
 717
 718   if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
 719
 720   Bool_t useEvent = kTRUE;
 721
 722   // Fist we reset histograms
 723   if ((mode & kFillRef)) fCumuRef->Reset();
 724   fCumuDiff->Reset();
 725
 726   // Numbers to cut away bad events and acceptance.
 727   Double_t runAvg = 0;
 728   Double_t max = 0;
 729   Int_t nInAvg = 0;
 730   Double_t avgSqr = 0;
 731   Int_t nBins = (dNdetadphi.GetNbinsX() * 6) / (fCumuDiff->GetNbinsX() * 5);
 732   Int_t nInBin = 0;
 733   Int_t nCurBin = 0, nPrevBin = 0;
 734   Int_t nCurRefBin = 0, nPrevRefBin = 0;
 735   Int_t nBadBins = 0;
 736   Bool_t firstBin = kFALSE;
 737   Double_t limit = 9999.;
 738   Bool_t mc = fType.Contains("MC");
 739
 740   // Then we loop over the input and calculate sum cos(k*n*phi)
 741   // and fill it in the reference and differential histograms
 742   Double_t eta, phi, weight;
 743   Double_t dQnRe = 0, dQ2nRe = 0, dQnIm = 0, dQ2nIm = 0;
 744
 745   for (Int_t etaBin = 1; etaBin <= dNdetadphi.GetNbinsX(); etaBin++) {
 746     eta = dNdetadphi.GetXaxis()->GetBinCenter(etaBin);
 747     nCurBin = fCumuDiff->GetXaxis()->FindBin(eta);
 748     nCurRefBin = fCumuRef->GetXaxis()->FindBin(eta);
 749     // If we have moved to a new bin in the flow hist, and less than half the eta
 750     // region has been covered by it we cut it away.
 751     if (nPrevBin == 0) nPrevBin = nCurBin;
 752     if (nPrevRefBin == 0) nPrevRefBin = nCurRefBin;
 753     if (nCurBin != nPrevBin) {
 754       if (nInBin <= nBins/2) {
 755         for (Int_t qBin = 1; qBin <= fCumuDiff->GetNbinsY(); qBin++) {
 756           Double_t removeContent = fCumuDiff->GetBinContent(nPrevBin, qBin);
 757           Double_t removeEta = fCumuDiff->GetXaxis()->GetBinCenter(nPrevBin);
 758           if (nCurRefBin != nPrevRefBin) {
 759             if (!(fFlags & kEtaGap)) {
 760               fCumuRef->Fill(removeEta, qBin, -removeContent);
 761               if ((fFlags & kSymEta)) {
 762                 fCumuRef->Fill(-1.*removeEta, qBin, -removeContent);
 763               }
 764             }
 765             if ((fFlags & kEtaGap)) {
 766               switch(qBin) {
 767                 case kHmultA: fCumuRef->Fill(removeEta, qBin, removeContent);  break;
 768                 case kHQnReA: fCumuRef->Fill(removeEta, qBin, removeContent);  break;
 769                 case kHQnImA: fCumuRef->Fill(removeEta, qBin, removeContent);  break;
 770                 case kHmultB: fCumuRef->Fill(-removeEta, qBin, removeContent); break;
 771                 case kHQnReB: fCumuRef->Fill(-removeEta, qBin, removeContent); break;
 772                 case kHQnImB: fCumuRef->Fill(-removeEta, qBin, removeContent); break;
 773                 default: break;
 774               }
 775             }
 776           }
 777           fCumuDiff->SetBinContent(nPrevBin, qBin, 0);
 778           fCumuDiff->SetBinError(nPrevBin, qBin, 0);
 779         }
 780       }
 781       nInBin = 0;
 782       nPrevBin = nCurBin;
 783       if (nCurRefBin != nPrevRefBin) nPrevRefBin = nCurRefBin;
 784     }
 785
 786     Bool_t data = kFALSE;
 787     for (Int_t phiBin = 0; phiBin <= dNdetadphi.GetNbinsY(); phiBin++) {
 788       if (phiBin == 0) {
 789         if (dNdetadphi.GetBinContent(etaBin, phiBin) == 0) break;
 790         else data = kTRUE;
 791         if ((fFlags & kEtaGap) && (mode & kFillRef) && TMath::Abs(eta) < fEtaGap) break;
 792         if ((fFlags & kEtaGap) && (mode & kFillRef) && TMath::Abs(eta) > TMath::Abs(limit)) break;
 793         if (data && !firstBin) {
 794           limit = dNdetadphi.GetXaxis()->GetBinLowEdge(etaBin);
 795           firstBin = kTRUE;
 796         }
 797         continue;
 798       }
 799       phi = dNdetadphi.GetYaxis()->GetBinCenter(phiBin);
 800       if ((fFlags & kEtaGap) && !mc && (phiBin == 7 || phiBin == 14 || phiBin == 15 || phiBin == 20)) continue;
 801       weight = dNdetadphi.GetBinContent(etaBin, phiBin);
 802
 803       // We calculate the average Nch per. bin
 804       avgSqr += weight*weight;
 805       runAvg += weight;
 806       nInAvg++;
 807       if (weight == 0) continue;
 808       if (weight > max) max = weight;
 809
 810       dQnRe = weight*TMath::Cos(fMoment*phi);
 811       dQnIm = weight*TMath::Sin(fMoment*phi);
 812       dQ2nRe = weight*TMath::Cos(2.*fMoment*phi);
 813       dQ2nIm = weight*TMath::Sin(2.*fMoment*phi);
 814
 815       // if we do not have an eta gap we fill the ref sums in eta
 816       if (!(fFlags & kEtaGap)) {
 817         if ((mode & kFillRef)){
 818           fCumuRef->Fill(eta, kHmultA, weight);
 819           fCumuRef->Fill(eta, kHQnReA, dQnRe);
 820           fCumuRef->Fill(eta, kHQnImA, dQnIm);
 821           fCumuRef->Fill(eta, kHmultB, weight);
 822           fCumuRef->Fill(eta, kHQnReB, dQnRe);
 823           fCumuRef->Fill(eta, kHQnImB, dQnIm);
 824
 825           fCumuRef->Fill(eta, kHQ2nRe, dQ2nRe);
 826           fCumuRef->Fill(eta, kHQ2nIm, dQ2nIm);
 827         }
 828       }
 829       // if we have an eta gap or we symmetrize around eta = 0
 830       // we fill in -eta
 831       if ((fFlags & kEtaGap)){
 832         if ((mode & kFillRef)){
 833           fCumuRef->Fill(eta, kHmultA, weight);
 834           fCumuRef->Fill(eta, kHQnReA, dQnRe);
 835           fCumuRef->Fill(eta, kHQnImA, dQnIm);
 836
 837           fCumuRef->Fill(-eta, kHmultB, weight);
 838           fCumuRef->Fill(-eta, kHQnReB, dQnRe);
 839           fCumuRef->Fill(-eta, kHQnImB, dQnIm);
 840         }
 841       }
 842
 843       if ((fFlags & kSymEta)){
 844         if ((mode & kFillRef)){
 845           fCumuRef->Fill(-eta, kHmultA, weight);
 846           fCumuRef->Fill(-eta, kHQnReA, dQnRe);
 847           fCumuRef->Fill(-eta, kHQnImA, dQnIm);
 848           fCumuRef->Fill(-eta, kHmultB, weight);
 849           fCumuRef->Fill(-eta, kHQnReB, dQnRe);
 850           fCumuRef->Fill(-eta, kHQnImB, dQnIm);
 851
 852           fCumuRef->Fill(-eta, kHQ2nRe, dQ2nRe);
 853           fCumuRef->Fill(-eta, kHQ2nIm, dQ2nIm);
 854         }
 855       }
 856
 857       // If we fill diff flow, we always fill it in eta
 858       // This is filled always, to be able to remove edge effects
 859       // It is reset both for kFillRef and kFillDiff to make the eta
 860       // gap analysis work.
 861       fCumuDiff->Fill(eta, kHmultA, weight);
 862       fCumuDiff->Fill(eta, kHQnReA, dQnRe);
 863       fCumuDiff->Fill(eta, kHQnImA, dQnIm);
 864       fCumuDiff->Fill(eta, kHQ2nRe, dQ2nRe);
 865       fCumuDiff->Fill(eta, kHQ2nIm, dQ2nIm);
 866
 867       // Fill acc. map
 868       if ((mode & kFillDiff)) fdNdedpAcc->Fill(eta, phi, weight);
 869     }
 870     if (data) {
 871       nInBin++;
 872     }
 873     // Outlier cut calculations
 874     if (nInAvg > 3 && !(fFlags & kEtaGap)) {
 875       runAvg /= nInAvg;
 876       avgSqr /= nInAvg;
 877       Double_t stdev = TMath::Sqrt(nInAvg/(nInAvg-1))*TMath::Sqrt(avgSqr - runAvg*runAvg);
 878       Double_t nSigma = (stdev == 0 ? 0 : (max-runAvg)/stdev);
 879       if (fSigmaCut > 0. && nSigma >= fSigmaCut && cent <= 60) nBadBins++;
 880       else nBadBins = 0;
 881       fOutliers->Fill(cent, nSigma);
 882       // We still finish the loop, for fOutliers to make sense,
 883       // but we do no keep the event for analysis
 884       if (nBadBins > 3) useEvent = kFALSE;
 885     }
 886     runAvg = 0;
 887     avgSqr = 0;
 888     nInAvg = 0;
 889     max = 0;
 890   }
 891
 892   return useEvent;
 893 }
 894 //_____________________________________________________________________
 895 void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate(Double_t cent, Bool_t skipFourP)
 896 {
 897   //
 898   // Calculate the Q cumulant of order fMoment
 899   //
 900   // Parameters:
 901   //  cent: Centrality of event
 902   //
 903
 904   if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
 905
 906   // We create the objects needed for the analysis
 907   Double_t dQnReA = 0, dQnImA = 0, multA;
 908   Double_t dQnReB = 0, dQnImB = 0, multB;
 909   Double_t dQ2nRe = 0, dQ2nIm = 0;
 910   Double_t pnRe = 0, p2nRe = 0, qnRe = 0, q2nRe = 0, pnIm = 0, p2nIm = 0, qnIm = 0, q2nIm = 0;
 911   Double_t two = 0, four = 0, twoPrime = 0, fourPrime = 0;
 912   Double_t cosPhi1Phi2 = 0, cosPhi1Phi2Phi3m = 0;
 913   Double_t sinPhi1Phi2 = 0, sinPhi1Phi2Phi3m = 0;
 914   Double_t cosPsi1Phi2 = 0, cosPsi1Phi2Phi3m = 0, cosPsi1Phi2Phi3p = 0;
 915   Double_t sinPsi1Phi2 = 0, sinPsi1Phi2Phi3m = 0, sinPsi1Phi2Phi3p = 0;
 916   Double_t eta = 0;
 917   Double_t mp = 0, mq = 0;
 918   Double_t w2 = 0, w4 = 0, w2p = 0, w4p = 0;
 919   Int_t refEtaBin = 0;
 920
 921   // We loop over the data 1 time!
 922   for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
 923     eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
 924     refEtaBin = fCumuRef->GetXaxis()->FindBin(eta);
 925     // The values for each individual etaBin bins are reset
 926     multA = 0;
 927     multB = 0;
 928     dQnReA = 0;
 929     dQnImA = 0;
 930     dQnReB = 0;
 931     dQnImB = 0;
 932     dQ2nRe = 0;
 933     dQ2nIm = 0;
 934
 935     mp = 0;
 936     pnRe = 0;
 937     pnIm = 0;
 938     p2nRe = 0;
 939     p2nIm = 0;
 940
 941     mq = 0;
 942     qnRe = 0;
 943     qnIm = 0;
 944     q2nRe = 0;
 945     q2nIm = 0;
 946
 947     // Reference flow
 948     multA  = fCumuRef->GetBinContent(refEtaBin, kHmultA);
 949     multB  = fCumuRef->GetBinContent(refEtaBin, kHmultB);
 950     dQnReA = fCumuRef->GetBinContent(refEtaBin, kHQnReA);
 951     dQnImA = fCumuRef->GetBinContent(refEtaBin, kHQnImA);
 952     dQnReB = fCumuRef->GetBinContent(refEtaBin, kHQnReB);
 953     dQnImB = fCumuRef->GetBinContent(refEtaBin, kHQnImB);
 954     dQ2nRe = fCumuRef->GetBinContent(refEtaBin, kHQ2nRe);
 955     dQ2nIm = fCumuRef->GetBinContent(refEtaBin, kHQ2nIm);
 956
 957     // For each etaBin bin the necessary values for differential flow
 958     // is calculated. .
 959     mp = fCumuDiff->GetBinContent(etaBin, kHmultA);
 960     pnRe = fCumuDiff->GetBinContent(etaBin, kHQnReA);
 961     pnIm = fCumuDiff->GetBinContent(etaBin, kHQnImA);
 962     p2nRe = fCumuDiff->GetBinContent(etaBin, kHQ2nRe);
 963     p2nIm = fCumuDiff->GetBinContent(etaBin, kHQ2nIm);
 964
 965     if (multA <= 3) continue;
 966     if (mp == 0) continue;
 967
 968     // The reference flow is calculated
 969     // 2-particle
 970     if (!(fFlags & kEtaGap)) {
 971       w2 = multA * (multA - 1.);
 972       two = dQnReA*dQnReA + dQnImA*dQnImA - multA;
 973     } else {
 974       w2 = multA * multB;
 975       two = dQnReA*dQnReB + dQnImA*dQnImB;
 976     }
 977
 978     fCumuHist->Fill(eta, cent, kW2Two, two);
 979     fCumuHist->Fill(eta, cent, kW2, w2);
 980
 981     fCumuHist->Fill(eta, cent, kQnReA, dQnReA);
 982     fCumuHist->Fill(eta, cent, kQnImA, dQnImA);
 983     fCumuHist->Fill(eta, cent, kMA, multA);
 984
 985     fCumuHist->Fill(eta, cent, kQnReB, dQnReB);
 986     fCumuHist->Fill(eta, cent, kQnImB, dQnImB);
 987     fCumuHist->Fill(eta, cent, kMB, multB);
 988
 989     // Differential flow calculations for each eta bin is done:
 990     // 2-particle differential flow
 991     if (!(fFlags & kEtaGap)) {
 992       mq = mp;
 993       qnRe = pnRe;
 994       qnIm = pnIm;
 995       q2nRe = p2nRe;
 996       q2nIm = p2nIm;
 997     }
 998
 999     w2p = mp * multB - mq;
1000     twoPrime = pnRe*dQnReB + pnIm*dQnImB - mq;
1001
1002     fCumuHist->Fill(eta, cent, kw2two, twoPrime);
1003     fCumuHist->Fill(eta, cent, kw2, w2p);
1004
1005     fCumuHist->Fill(eta, cent, kpnRe, pnRe);
1006     fCumuHist->Fill(eta, cent, kpnIm, pnIm);
1007     fCumuHist->Fill(eta, cent, kmp, mp);
1008
1009     if ((fFlags & kEtaGap) || skipFourP) continue;
1010     // The reference flow is calculated
1011     // 4-particle
1012     w4 = multA * (multA - 1.) * (multA - 2.) * (multA - 3.);
1013
1014     four = 2.*multA*(multA-3.) + TMath::Power((TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.)),2.)
1015              -4.*(multA-2.)*(TMath::Power(dQnReA,2.) + TMath::Power(dQnImA,2.))
1016              -2.*(TMath::Power(dQnReA,2.)*dQ2nRe+2.*dQnReA*dQnImA*dQ2nIm-TMath::Power(dQnImA,2.)*dQ2nRe)
1017              +(TMath::Power(dQ2nRe,2.)+TMath::Power(dQ2nIm,2.));
1018
1019     fCumuHist->Fill(eta, cent, kW4Four, four);
1020     fCumuHist->Fill(eta, cent, kW4, w4);
1021
1022     cosPhi1Phi2 = dQnReA*dQnReA - dQnImA*dQnImA - dQ2nRe;
1023     sinPhi1Phi2 = 2.*dQnReA*dQnImA - dQ2nIm;
1024
1025     cosPhi1Phi2Phi3m = dQnReA*(TMath::Power(dQnReA,2)+TMath::Power(dQnImA,2))-dQnReA*dQ2nRe-dQnImA*dQ2nIm-2.*(multA-1)*dQnReA;
1026
1027     sinPhi1Phi2Phi3m = -dQnImA*(TMath::Power(dQnReA,2)+TMath::Power(dQnImA,2))+dQnReA*dQ2nIm-dQnImA*dQ2nRe+2.*(multA-1)*dQnImA;
1028
1029     fCumuHist->Fill(eta, cent, kCosphi1phi2, cosPhi1Phi2);
1030     fCumuHist->Fill(eta, cent, kSinphi1phi2, sinPhi1Phi2);
1031     fCumuHist->Fill(eta, cent, kCosphi1phi2phi3m, cosPhi1Phi2Phi3m);
1032     fCumuHist->Fill(eta, cent, kSinphi1phi2phi3m, sinPhi1Phi2Phi3m);
1033     fCumuHist->Fill(eta, cent, kMm1m2, multA*(multA-1.)*(multA-2.));
1034
1035     // Differential flow calculations for each eta bin bin is done:
1036     // 4-particle differential flow
1037     w4p = (mp * multA - 3.*mq)*(multA - 1.)*(multA - 2.);
1038
1039     fourPrime = (TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.))*(pnRe*dQnReA+pnIm*dQnImA)
1040                       - q2nRe*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))
1041                       - 2.*q2nIm*dQnReA*dQnImA
1042                       - pnRe*(dQnReA*dQ2nRe+dQnImA*dQ2nIm)
1043                       + pnIm*(dQnImA*dQ2nRe-dQnReA*dQ2nIm)
1044                       - 2.*multA*(pnRe*dQnReA+pnIm*dQnImA)
1045                       - 2.*(TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.))*mq
1046                       + 6.*(qnRe*dQnReA+qnIm*dQnImA)
1047                       + 1.*(q2nRe*dQ2nRe+q2nIm*dQ2nIm)
1048                       + 2.*(pnRe*dQnReA+pnIm*dQnImA)
1049                       + 2.*mq*multA
1050                       - 6.*mq;
1051
1052     fCumuHist->Fill(eta, cent, kw4four, fourPrime);
1053     fCumuHist->Fill(eta, cent, kw4, w4p);
1054
1055     cosPsi1Phi2 = pnRe*dQnReA - pnIm*dQnImA - q2nRe;
1056     sinPsi1Phi2 = pnRe*dQnImA + pnIm*dQnReA - q2nIm;
1057
1058     cosPsi1Phi2Phi3p = pnRe*(TMath::Power(dQnImA,2.)+TMath::Power(dQnReA,2.)-multA)
1059                           - 1.*(q2nRe*dQnReA+q2nIm*dQnImA)
1060                           - mq*dQnReA+2.*qnRe;
1061
1062     sinPsi1Phi2Phi3p = pnIm*(TMath::Power(dQnImA,2.)+TMath::Power(dQnReA,2.)-multA)
1063                           - 1.*(q2nIm*dQnReA-q2nRe*dQnImA)
1064                           - mq*dQnImA+2.*qnIm;
1065
1066     cosPsi1Phi2Phi3m = pnRe*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))+2.*pnIm*dQnReA*dQnImA
1067                           - 1.*(pnRe*dQ2nRe+pnIm*dQ2nIm)
1068                           - 2.*mq*dQnReA+2.*qnRe;
1069
1070     sinPsi1Phi2Phi3m = pnIm*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))-2.*pnRe*dQnReA*dQnImA
1071                           - 1.*(pnIm*dQ2nRe-pnRe*dQ2nIm)
1072                           + 2.*mq*dQnImA-2.*qnIm;
1073
1074     fCumuHist->Fill(eta, cent, kCospsi1phi2, cosPsi1Phi2);
1075     fCumuHist->Fill(eta, cent, kSinpsi1phi2, sinPsi1Phi2);
1076     fCumuHist->Fill(eta, cent, kCospsi1phi2phi3m, cosPsi1Phi2Phi3m);
1077     fCumuHist->Fill(eta, cent, kSinpsi1phi2phi3m, sinPsi1Phi2Phi3m);
1078     fCumuHist->Fill(eta, cent, kmpmq, (mp*multA-2.*mq)*(multA-1.));
1079     fCumuHist->Fill(eta, cent, kCospsi1phi2phi3p, cosPsi1Phi2Phi3p);
1080     fCumuHist->Fill(eta, cent, kSinpsi1phi2phi3p, sinPsi1Phi2Phi3p);
1081   }
1082   // Event count
1083   fCumuHist->Fill(-7., cent, -0.5, 1.);
1084
1085 }
1086 //_____________________________________________________________________
1087 void AliForwardFlowTaskQC::VertexBin::CumulantsTerminate(TList* inlist, TList* outlist)
1088 {
1089   //
1090   //  Finalizes the Q cumulant calculations
1091   //
1092   //  Parameters:
1093   //   inlist: input sumlist
1094   //   outlist: output result list
1095   //
1096
1097   // Re-find cumulants hist if Terminate is called separately
1098   if (!fCumuHist) {
1099     TList* list = (TList*)inlist->FindObject(Form("%svertex_%d_%d%s", fType.Data(), fVzMin, fVzMax, ((fFlags & kEtaGap) ? "_etaGap" : "")));
1100     fCumuHist = (TH3D*)list->FindObject(Form("%sv%d_vertex_%d_%d%s_cumu", fType.Data(), fMoment,
1101                                         fVzMin, fVzMax, ((fFlags & kEtaGap) ? "_etaGap" : "")));
1102   }
1103
1104   // Create outputs
1105   TList* vtxList = (TList*)outlist->FindObject("vtxList");
1106   if (!vtxList) {
1107     vtxList = new TList();
1108     vtxList->SetName("vtxList");
1109     outlist->Add(vtxList);
1110   }
1111   TList* nualist = (TList*)outlist->FindObject("NUATerms");
1112   if (!nualist) {
1113     nualist = new TList();
1114     nualist->SetName("NUATerms");
1115     outlist->Add(nualist);
1116   }
1117
1118   TH1I* quality = (TH1I*)outlist->FindObject(Form("hQCQuality%s%s", fType.Data(), ((fFlags & kEtaGap) ? "_etaGap" : "")));
1119
1120   // Differential flow
1121   TProfile2D* cumu2Sum = (TProfile2D*)outlist->FindObject(Form("%sQC2_v%d_unCorr%s", fType.Data(), fMoment, ((fFlags & kEtaGap) ? "_etaGap" : "")));
1122   TProfile2D* cumu4Sum = (TProfile2D*)outlist->FindObject(Form("%sQC4_v%d_unCorr%s", fType.Data(), fMoment, ((fFlags & kEtaGap) ? "_etaGap" : "")));
1123   if (!cumu2Sum) {
1124     cumu2Sum = new TProfile2D(Form("%sQC2_v%d_unCorr%s", fType.Data(), fMoment, ((fFlags & kEtaGap) ? "_etaGap" : "")),
1125                            Form("%sQC2_v%d_unCorr%s", fType.Data(), fMoment, ((fFlags & kEtaGap) ? "_etaGap" : "")),
1126               fCumuHist->GetNbinsX(), fCumuHist->GetXaxis()->GetXmin(), fCumuHist->GetXaxis()->GetXmax(),
1127               fCumuHist->GetNbinsY(), fCumuHist->GetYaxis()->GetXmin(), fCumuHist->GetYaxis()->GetXmax());
1128     SetupCentAxis(cumu2Sum->GetYaxis());
1129     outlist->Add(cumu2Sum);
1130   }
1131   if (!cumu4Sum && !(fFlags & kEtaGap)) {
1132     cumu4Sum = new TProfile2D(Form("%sQC4_v%d_unCorr%s", fType.Data(), fMoment, ((fFlags & kEtaGap) ? "_etaGap" : "")),
1133                            Form("%sQC4_v%d_unCorr%s", fType.Data(), fMoment, ((fFlags & kEtaGap) ? "_etaGap" : "")),
1134               fCumuHist->GetNbinsX(), fCumuHist->GetXaxis()->GetXmin(), fCumuHist->GetXaxis()->GetXmax(),
1135               fCumuHist->GetNbinsY(), fCumuHist->GetYaxis()->GetXmin(), fCumuHist->GetYaxis()->GetXmax());
1136     SetupCentAxis(cumu4Sum->GetYaxis());
1137     outlist->Add(cumu4Sum);
1138   }
1139
1140   TProfile2D* cumu2 = new TProfile2D(Form("%sQC2_v%d_unCorr%s_vtx%3.1f", fType.Data(), fMoment,
1141                          ((fFlags & kEtaGap) ? "_etaGap" : ""), (fVzMin+fVzMax)/2.),
1142                          Form("%sQC2_v%d_unCorr%s_vtx%3.1f", fType.Data(), fMoment,
1143                          ((fFlags & kEtaGap) ? "_etaGap" : ""), (fVzMin+fVzMax)/2.),
1144             fCumuHist->GetNbinsX(), fCumuHist->GetXaxis()->GetXmin(), fCumuHist->GetXaxis()->GetXmax(),
1145             fCumuHist->GetNbinsY(), fCumuHist->GetYaxis()->GetXmin(), fCumuHist->GetYaxis()->GetXmax());
1146   SetupCentAxis(cumu2->GetYaxis());
1147   vtxList->Add(cumu2);
1148   TProfile2D* cumu4 = 0;
1149   if (!(fFlags & kEtaGap)) {
1150     cumu4 = new TProfile2D(Form("%sQC4_v%d_unCorr%s_vtx%3.1f", fType.Data(), fMoment,
1151                             ((fFlags & kEtaGap) ? "_etaGap" : ""), (fVzMin+fVzMax)/2.),
1152                             Form("%sQC4_v%d_unCorr%s_vtx%3.1f", fType.Data(), fMoment,
1153                             ((fFlags & kEtaGap) ? "_etaGap" : ""), (fVzMin+fVzMax)/2.),
1154               fCumuHist->GetNbinsX(), fCumuHist->GetXaxis()->GetXmin(), fCumuHist->GetXaxis()->GetXmax(),
1155               fCumuHist->GetNbinsY(), fCumuHist->GetYaxis()->GetXmin(), fCumuHist->GetYaxis()->GetXmax());
1156     SetupCentAxis(cumu4->GetYaxis());
1157     vtxList->Add(cumu4);
1158   }
1159
1160   // Reference flow
1161   TProfile2D* cumu2Ref = (TProfile2D*)outlist->FindObject(Form("%sQC2_v%d_ref_unCorr%s", fType.Data(), fMoment, ((fFlags & kEtaGap) ? "_etaGap" : "")));
1162   TProfile2D* cumu4Ref = (TProfile2D*)outlist->FindObject(Form("%sQC4_v%d_ref_unCorr%s", fType.Data(), fMoment, ((fFlags & kEtaGap) ? "_etaGap" : "")));
1163   if (!cumu2Ref) {
1164     cumu2Ref = new TProfile2D(Form("%sQC2_v%d_ref_unCorr%s", fType.Data(), fMoment, ((fFlags & kEtaGap) ? "_etaGap" : "")),
1165                            Form("%sQC2_v%d_ref_unCorr%s", fType.Data(), fMoment, ((fFlags & kEtaGap) ? "_etaGap" : "")),
1166               fCumuHist->GetNbinsX(), fCumuHist->GetXaxis()->GetXmin(), fCumuHist->GetXaxis()->GetXmax(),
1167               fCumuHist->GetNbinsY(), fCumuHist->GetYaxis()->GetXmin(), fCumuHist->GetYaxis()->GetXmax());
1168     SetupCentAxis(cumu2Ref->GetYaxis());
1169     outlist->Add(cumu2Ref);
1170   }
1171   if (!cumu4Ref && !(fFlags & kEtaGap)) {
1172     cumu4Ref = new TProfile2D(Form("%sQC4_v%d_ref_unCorr%s", fType.Data(), fMoment, ((fFlags & kEtaGap) ? "_etaGap" : "")),
1173                            Form("%sQC4_v%d_ref_unCorr%s", fType.Data(), fMoment, ((fFlags & kEtaGap) ? "_etaGap" : "")),
1174               fCumuHist->GetNbinsX(), fCumuHist->GetXaxis()->GetXmin(), fCumuHist->GetXaxis()->GetXmax(),
1175               fCumuHist->GetNbinsY(), fCumuHist->GetYaxis()->GetXmin(), fCumuHist->GetYaxis()->GetXmax());
1176     SetupCentAxis(cumu4Ref->GetYaxis());
1177     outlist->Add(cumu4Ref);
1178   }
1179
1180   // NUA terms
1181   TProfile2D* cosPhi = (TProfile2D*)nualist->FindObject(Form("%sCosPhi_v%d_unCorr%s", fType.Data(), fMoment,
1182                             ((fFlags & kEtaGap) ? "_etaGap" : "")));
1183   if (!cosPhi) {
1184     cosPhi = new TProfile2D(Form("%sCosPhi_v%d_unCorr%s", fType.Data(), fMoment,
1185                             ((fFlags & kEtaGap) ? "_etaGap" : "")),
1186                                   Form("%sCosPhi_v%d_unCorr%s", fType.Data(), fMoment,
1187                             ((fFlags & kEtaGap) ? "_etaGap" : "")),
1188                             48, -6, 6,
1189                             fCumuHist->GetNbinsY(), fCumuHist->GetYaxis()->GetXmin(), fCumuHist->GetYaxis()->GetXmax());
1190     SetupCentAxis(cosPhi->GetYaxis());
1191     nualist->Add(cosPhi);
1192   }
1193
1194   TProfile2D* cosPsi = (TProfile2D*)nualist->FindObject(Form("%sCosPsi_v%d_unCorr%s", fType.Data(), fMoment,
1195                             ((fFlags & kEtaGap) ? "_etaGap" : "")));
1196   if (!cosPsi) {
1197     cosPsi = new TProfile2D(Form("%sCosPsi_v%d_unCorr%s", fType.Data(), fMoment,
1198                             ((fFlags & kEtaGap) ? "_etaGap" : "")),
1199                                   Form("%sCosPsi_v%d_unCorr%s", fType.Data(), fMoment,
1200                             ((fFlags & kEtaGap) ? "_etaGap" : "")),
1201                             48, -6, 6,
1202                             fCumuHist->GetNbinsY(), fCumuHist->GetYaxis()->GetXmin(), fCumuHist->GetYaxis()->GetXmax());
1203     SetupCentAxis(cosPsi->GetYaxis());
1204     nualist->Add(cosPsi);
1205   }
1206
1207   TProfile2D* sinPhi = (TProfile2D*)nualist->FindObject(Form("%sSinPhi_v%d_unCorr%s", fType.Data(), fMoment,
1208                             ((fFlags & kEtaGap) ? "_etaGap" : "")));
1209   if (!sinPhi) {
1210     sinPhi = new TProfile2D(Form("%sSinPhi_v%d_unCorr%s", fType.Data(), fMoment,
1211                             ((fFlags & kEtaGap) ? "_etaGap" : "")),
1212                                   Form("%sSinPhi_v%d_unCorr%s", fType.Data(), fMoment,
1213                             ((fFlags & kEtaGap) ? "_etaGap" : "")),
1214                             48, -6, 6,
1215                             fCumuHist->GetNbinsY(), fCumuHist->GetYaxis()->GetXmin(), fCumuHist->GetYaxis()->GetXmax());
1216     SetupCentAxis(sinPhi->GetYaxis());
1217     nualist->Add(sinPhi);
1218   }
1219
1220   TProfile2D* sinPsi = (TProfile2D*)nualist->FindObject(Form("%sSinPsi_v%d_unCorr%s", fType.Data(), fMoment,
1221                             ((fFlags & kEtaGap) ? "_etaGap" : "")));
1222   if (!sinPsi) {
1223     sinPsi = new TProfile2D(Form("%sSinPsi_v%d_unCorr%s", fType.Data(), fMoment,
1224                             ((fFlags & kEtaGap) ? "_etaGap" : "")),
1225                                   Form("%sSinPsi_v%d_unCorr%s", fType.Data(), fMoment,
1226                             ((fFlags & kEtaGap) ? "_etaGap" : "")),
1227                             48, -6, 6,
1228                             fCumuHist->GetNbinsY(), fCumuHist->GetYaxis()->GetXmin(), fCumuHist->GetYaxis()->GetXmax());
1229     SetupCentAxis(sinPsi->GetYaxis());
1230     nualist->Add(sinPsi);
1231   }
1232
1233   // For flow calculations
1234   Double_t two = 0, qc2 = 0, vnTwo = 0, four = 0, qc4 = 0, vnFour = 0;
1235   Double_t twoPrime = 0, qc2Prime = 0, vnTwoDiff = 0, fourPrime = 0, qc4Prime = 0, vnFourDiff = 0;
1236   Double_t w2 = 0, w4 = 0, w2p = 0, w4p = 0;
1237   Double_t w2Two = 0, w2pTwoPrime = 0, w4Four = 0, w4pFourPrime = 0;
1238   Double_t cosP1nPhiA = 0, sinP1nPhiA = 0, multA = 0;
1239   Double_t cosP1nPhiB = 0, sinP1nPhiB = 0, multB = 0;
1240   Double_t cosP1nPhi1P1nPhi2 = 0, sinP1nPhi1P1nPhi2 = 0;
1241   Double_t cosP1nPhi1M1nPhi2M1nPhi3 = 0, sinP1nPhi1M1nPhi2M1nPhi3 = 0, multm1m2 = 0;
1242   Double_t cosP1nPsi = 0, sinP1nPsi = 0, mp = 0, cosP1nPsi1P1nPhi2 = 0, sinP1nPsi1P1nPhi2 = 0;
1243   Double_t cosP1nPsi1M1nPhi2M1nPhi3 = 0, sinP1nPsi1M1nPhi2M1nPhi3 = 0, mpqMult = 0;
1244   Double_t cosP1nPsi1P1nPhi2M1nPhi3 = 0, sinP1nPsi1P1nPhi2M1nPhi3 = 0;
1245
1246   // Loop over cumulant histogram for final calculations
1247   // Centrality loop
1248   for (Int_t cBin = 1; cBin <= fCumuHist->GetNbinsY(); cBin++) {
1249     // for weighted avg.
1250     Double_t cent = fCumuHist->GetYaxis()->GetBinCenter(cBin);
1251     if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), fMoment, cent));
1252     // Eta loop
1253     for (Int_t etaBin = 1; etaBin <= fCumuHist->GetNbinsX(); etaBin++) {
1254       Double_t eta = fCumuHist->GetXaxis()->GetBinCenter(etaBin);
1255       // 2-particle reference flow
1256       w2Two = fCumuHist->GetBinContent(etaBin, cBin, kW2Two);
1257       w2 = fCumuHist->GetBinContent(etaBin, cBin, kW2);
1258       multA = fCumuHist->GetBinContent(etaBin, cBin, kMA);
1259       multB = fCumuHist->GetBinContent(etaBin, cBin, kMB);
1260       if (w2 == 0 || multA == 0 || multB == 0) continue;
1261       cosP1nPhiA = fCumuHist->GetBinContent(etaBin, cBin, kQnReA);
1262       sinP1nPhiA = fCumuHist->GetBinContent(etaBin, cBin, kQnImA);
1263       cosP1nPhiB = fCumuHist->GetBinContent(etaBin, cBin, kQnReB);
1264       sinP1nPhiB = fCumuHist->GetBinContent(etaBin, cBin, kQnImB);
1265
1266       cosP1nPhiA /= multA;
1267       sinP1nPhiA /= multA;
1268       cosP1nPhiB /= multB;
1269       sinP1nPhiB /= multB;
1270       two = w2Two / w2;
1271       // With no eta gap the last two terms are <<cos(phi)>>^2 and <<sin(phi)>>^2,
1272       // with eta gap the different coverage is taken into account.
1273       // The next line covers both cases.
1274       qc2 = two - cosP1nPhiA*cosP1nPhiB - sinP1nPhiA*sinP1nPhiB;
1275       if (qc2 <= 0) {
1276         if (fDebug > 0) AliInfo(Form("%s: QC_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f - skipping", fType.Data(), fMoment, qc2, eta, cent));
1277         quality->Fill((fMoment-2)*8+2);
1278         continue;
1279       }
1280        vnTwo = TMath::Sqrt(qc2);
1281       if (!TMath::IsNaN(vnTwo*multA)) {
1282         quality->Fill((fMoment-2)*8+1);
1283         cumu2Ref->Fill(eta, cent, vnTwo);
1284       }
1285
1286       // 2-particle differential flow
1287       w2pTwoPrime = fCumuHist->GetBinContent(etaBin, cBin, kw2two);
1288       w2p = fCumuHist->GetBinContent(etaBin, cBin, kw2);
1289       mp = fCumuHist->GetBinContent(etaBin, cBin, kmp);
1290       if (w2p == 0 || mp == 0) continue;
1291       cosP1nPsi = fCumuHist->GetBinContent(etaBin, cBin, kpnRe);
1292       sinP1nPsi = fCumuHist->GetBinContent(etaBin, cBin, kpnIm);
1293
1294       cosP1nPsi /= mp;
1295       sinP1nPsi /= mp;
1296       twoPrime = w2pTwoPrime / w2p;
1297       qc2Prime = twoPrime - sinP1nPsi*sinP1nPhiB - cosP1nPsi*cosP1nPhiB;
1298
1299       cosPhi->Fill(eta, cent, cosP1nPhiB);
1300       cosPsi->Fill(eta, cent, cosP1nPsi);
1301       sinPhi->Fill(eta, cent, sinP1nPhiB);
1302       sinPsi->Fill(eta, cent, sinP1nPsi);
1303
1304       vnTwoDiff = qc2Prime / TMath::Sqrt(qc2);
1305       if (!TMath::IsNaN(vnTwoDiff*mp)) {
1306         cumu2->Fill(eta, cent, vnTwoDiff);
1307         quality->Fill((fMoment-2)*8+3);
1308       }
1309       else
1310         quality->Fill((fMoment-2)*8+4);
1311
1312       if (fDebug > 1) AliInfo(Form("%s: v_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f", fType.Data(), fMoment, vnTwoDiff, eta, cent));
1313       if ((fFlags & kEtaGap)) continue;
1314
1315       // 4-particle reference flow
1316       w4Four = fCumuHist->GetBinContent(etaBin, cBin, kW4Four);
1317       w4 = fCumuHist->GetBinContent(etaBin, cBin, kW4);
1318       multm1m2 = fCumuHist->GetBinContent(etaBin, cBin, kMm1m2);
1319       if (w4 == 0 || multm1m2 == 0) continue;
1320       cosP1nPhi1P1nPhi2 = fCumuHist->GetBinContent(etaBin, cBin, kCosphi1phi2);
1321       sinP1nPhi1P1nPhi2 = fCumuHist->GetBinContent(etaBin, cBin, kSinphi1phi2);
1322       cosP1nPhi1M1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, cBin, kCosphi1phi2phi3m);
1323       sinP1nPhi1M1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, cBin, kSinphi1phi2phi3m);
1324
1325       cosP1nPhi1P1nPhi2 /= w2;
1326       sinP1nPhi1P1nPhi2 /= w2;
1327       cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
1328       sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
1329       four = w4Four / w4;
1330       qc4 = four-2.*TMath::Power(two,2.)
1331          - 4.*cosP1nPhiA*cosP1nPhi1M1nPhi2M1nPhi3
1332          + 4.*sinP1nPhiA*sinP1nPhi1M1nPhi2M1nPhi3-TMath::Power(cosP1nPhi1P1nPhi2,2.)-TMath::Power(sinP1nPhi1P1nPhi2,2.)
1333          + 4.*cosP1nPhi1P1nPhi2*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
1334          + 8.*sinP1nPhi1P1nPhi2*sinP1nPhiA*cosP1nPhiA
1335          + 8.*two*(TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.))
1336          - 6.*TMath::Power((TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.)),2.);
1337
1338       if (qc4 >= 0) {
1339         if (fDebug > 0) AliInfo(Form("%s: QC_%d{4} = %1.3f for eta = %1.2f and centrality %3.1f - skipping", fType.Data(), fMoment, qc2, eta, cent));
1340         quality->Fill((fMoment-2)*8+6);
1341         continue;
1342       }
1343       vnFour = TMath::Power(-qc4, 0.25);
1344       if (!TMath::IsNaN(vnFour*multm1m2)) {
1345         quality->Fill((fMoment-2)*8+5);
1346         cumu4Ref->Fill(eta, cent, vnFour);
1347       }
1348
1349       // 4-particle differential flow
1350       w4pFourPrime = fCumuHist->GetBinContent(etaBin, cBin, kw4four);
1351       w4p = fCumuHist->GetBinContent(etaBin, cBin, kw4);
1352       mpqMult = fCumuHist->GetBinContent(etaBin, cBin, kmpmq);
1353       if (w4p == 0 || mpqMult == 0) continue;
1354       cosP1nPsi1P1nPhi2 = fCumuHist->GetBinContent(etaBin, cBin, kCospsi1phi2);
1355       sinP1nPsi1P1nPhi2 = fCumuHist->GetBinContent(etaBin, cBin, kSinpsi1phi2);
1356       cosP1nPsi1M1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, cBin, kCospsi1phi2phi3m);
1357       sinP1nPsi1M1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, cBin, kSinpsi1phi2phi3m);
1358       cosP1nPsi1P1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, cBin, kCospsi1phi2phi3p);
1359       sinP1nPsi1P1nPhi2M1nPhi3 = fCumuHist->GetBinContent(etaBin, cBin, kSinpsi1phi2phi3p);
1360
1361       cosP1nPsi1P1nPhi2 /= w2p;
1362       sinP1nPsi1P1nPhi2 /= w2p;
1363       cosP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
1364       sinP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
1365       cosP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
1366       sinP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
1367
1368       fourPrime = w4pFourPrime / w4p;
1369
1370       qc4Prime = fourPrime-2.*twoPrime*two
1371                 - cosP1nPsi*cosP1nPhi1M1nPhi2M1nPhi3
1372                 + sinP1nPsi*sinP1nPhi1M1nPhi2M1nPhi3
1373                 - cosP1nPhiA*cosP1nPsi1M1nPhi2M1nPhi3
1374                 + sinP1nPhiA*sinP1nPsi1M1nPhi2M1nPhi3
1375                 - 2.*cosP1nPhiA*cosP1nPsi1P1nPhi2M1nPhi3
1376                 - 2.*sinP1nPhiA*sinP1nPsi1P1nPhi2M1nPhi3
1377                 - cosP1nPsi1P1nPhi2*cosP1nPhi1P1nPhi2
1378                 - sinP1nPsi1P1nPhi2*sinP1nPhi1P1nPhi2
1379                 + 2.*cosP1nPhi1P1nPhi2*(cosP1nPsi*cosP1nPhiA-sinP1nPsi*sinP1nPhiA)
1380                 + 2.*sinP1nPhi1P1nPhi2*(cosP1nPsi*sinP1nPhiA+sinP1nPsi*cosP1nPhiA)
1381                 + 4.*two*(cosP1nPsi*cosP1nPhiA+sinP1nPsi*sinP1nPhiA)
1382                 + 2.*cosP1nPsi1P1nPhi2*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
1383                 + 4.*sinP1nPsi1P1nPhi2*cosP1nPhiA*sinP1nPhiA
1384                 + 4.*twoPrime*(TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.))
1385                 - 6.*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
1386                 * (cosP1nPsi*cosP1nPhiA-sinP1nPsi*sinP1nPhiA)
1387                 - 12.*cosP1nPhiA*sinP1nPhiA
1388                 * (sinP1nPsi*cosP1nPhiA+cosP1nPsi*sinP1nPhiA);
1389
1390       vnFourDiff = - qc4Prime / TMath::Power(-qc4, 0.75);
1391       if (!TMath::IsNaN(vnFourDiff*mpqMult) && vnFourDiff > 0) {
1392         cumu4->Fill(eta, cent, vnFourDiff);
1393         quality->Fill((fMoment-2)*8+7);
1394       }
1395       else
1396         quality->Fill((fMoment-2)*8+8);
1397
1398       if (fDebug > 1) AliInfo(Form("%s: v_%d{4} = %1.3f for eta = %1.2f and centrality %3.1f", fType.Data(), fMoment, vnFourDiff, eta, cent));
1399     } // End of eta loop
1400   } // End of centrality loop
1401
1402   cumu2Sum->Add(cumu2);
1403   if (!(fFlags & kEtaGap)) cumu4Sum->Add(cumu4);
1404
1405   return;
1406 }
1407 //_____________________________________________________________________
1408 void AliForwardFlowTaskQC::VertexBin::SetupCentAxis(TAxis* axis)
1409 {
1410   //
1411   // Setup centrality axis for histogram
1412   //
1413   // Parameters:
1414   //  axis: centrality axis
1415   //
1416   if (!axis) {
1417     AliError("Null pointer passed for axis");
1418     return;
1419   }
1420
1421   if ((fFlags & kSatVtx)) {
1422     Double_t cent[3] = {0, 40, 100};
1423     axis->Set(2, cent);
1424   } else {
1425     Double_t cent[13] = {0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 100 };
1426 //    Double_t cent[13] = {0, 2.5, 15, 25, 50, 60, 70, 80, 90, 100, 110, 115, 120 };
1427     axis->Set(12, cent);
1428   }
1429
1430   return;
1431 }
1432 //_____________________________________________________________________
1433 void AliForwardFlowTaskQC::PrintFlowSetup() const
1434 {
1435   //
1436   // Print the setup of the flow task
1437   //
1438   Printf("AliForwardFlowTaskQC::Print");
1439   Printf("Number of bins in vertex axis   :\t%d", fVtxAxis->GetNbins());
1440   Printf("Range of vertex axis            :\t[%3.1f,%3.1f]",
1441                           fVtxAxis->GetXmin(), fVtxAxis->GetXmax());
1442   printf("Doing flow analysis for         :\t");
1443   for (Int_t n  = 0; n < fV.GetSize(); n++) printf("v%d ", fV.At(n));
1444   printf("\n");
1445   Printf("Satellite vertex flag           :\t%s", ((fFlowFlags & kSatVtx) ? "true" : "false"));
1446   Printf("Symmetrize ref. flow wrt eta = 0:\t%s", ((fFlowFlags & kSymEta) ? "true" : "false"));
1447   Printf("Use an eta-gap for ref. flow    :\t%s", ((fFlowFlags & kEtaGap) ? "true" : "false"));
1448   Printf("FMD sigma cut:                  :\t%f", fFMDCut);
1449   Printf("SPD sigma cut:                  :\t%f", fSPDCut);
1450   if ((fFlowFlags & kEtaGap))
1451     Printf("Eta gap:                        :\t%f", fEtaGap);
1452
1453 }
1454 //_____________________________________________________________________
1455 //
1456 //
1457 // EOF