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