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