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