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