AliPhysics  aaf9c62 (aaf9c62)
AliForwardFlowTaskQC.cxx
Go to the documentation of this file.
1 
2 //
3 // Calculate flow in the forward and central regions using the Q cumulants method.
4 //
5 // Inputs:
6 // - AliAODEvent
7 //
8 // Outputs:
9 // - AnalysisResults.root or forward_flow.root
10 //
11 #include <TROOT.h>
12 #include <TSystem.h>
13 #include <TInterpreter.h>
14 #include <TChain.h>
15 #include <TFile.h>
16 #include <TList.h>
17 #include <TMath.h>
18 #include <TH3D.h>
19 #include <TProfile2D.h>
20 #include <TParameter.h>
21 #include <TMatrixD.h>
22 #include <TVectorD.h>
23 #include <TGraph.h>
24 #include "AliLog.h"
25 #include "AliForwardFlowTaskQC.h"
26 #include "AliAnalysisManager.h"
27 #include "AliAODHandler.h"
28 #include "AliAODInputHandler.h"
29 #include "AliAODForwardMult.h"
30 #include "AliAODCentralMult.h"
31 #include "AliAODEvent.h"
32 #include "AliForwardUtil.h"
33 #include "AliVVZERO.h"
34 #include "AliAODVertex.h"
35 #include "AliCentrality.h"
36 #include "AliESDEvent.h"
37 #include "AliVTrack.h"
38 #include "AliESDtrack.h"
39 #include "AliAODTrack.h"
40 #include "AliAnalysisFilter.h"
41 
42 ClassImp(AliForwardFlowTaskQC)
43 #if 0
44 ; // For emacs
45 #endif
46 
49  fVtxAxis(), // Axis to control vertex binning
50  fCentAxis(), // Axis to control centrality/multiplicity binning
51  fFMDCut(-1), // FMD sigma cut
52  fSPDCut(-1), // SPD sigma cut
53  fFlowFlags(0), // Flow flags
54  fEtaGap(0.), // Eta gap value
55  fBinsForward(), // List with forward flow hists
56  fBinsCentral(), // List with central flow hists
57  fSumList(0), // Event sum list
58  fOutputList(0), // Result output list
59  fAOD(0), // AOD input event
60  fTrackCuts(0), // ESD track cuts
61  fMaxMoment(0), // Max flow moment
62  fVtx(1111), // Z vertex coordinate
63  fCent(-1), // Centrality
64  fHistdNdedpV0(), // Hist for v0
65  fHistdNdedp3Cor(), // Hist for combining detectors
66  fHistFMDSPDCorr(), // FMD SPD correlation
67  fHistCent(), // Hist for centrality
68  fHistVertexSel(), // Hist for selected vertices
69  fHistEventSel() // Hist for event selection
70 {
71  //
72  // Default constructor
73  //
74 }
75 //_____________________________________________________________________
77  : AliAnalysisTaskSE(name),
78  fVtxAxis(), // Axis to control vertex binning
79  fCentAxis(), // Axis to control centrality/multiplicity binning
80  fFMDCut(-1), // FMD sigma cut
81  fSPDCut(-1), // SPD sigma cut
82  fFlowFlags(0x0), // Flow flags
83  fEtaGap(0.), // Eta gap value
84  fBinsForward(), // List with forward flow hists
85  fBinsCentral(), // List with central flow hists
86  fSumList(0), // Event sum list
87  fOutputList(0), // Result output list
88  fAOD(0), // AOD input event
89  fTrackCuts(0), // ESD track cuts
90  fMaxMoment(4), // Max flow moment
91  fVtx(1111), // Z vertex coordinate
92  fCent(-1), // Centrality
93  fHistdNdedpV0(), // Histo for v0
94  fHistdNdedp3Cor(), // Histo for combining detectors
95  fHistFMDSPDCorr(), // FMD SPD correlation
96  fHistCent(), // Hist for centrality
97  fHistVertexSel(), // Hist for selected vertices
98  fHistEventSel() // Hist for event selection
99 {
100  //
101  // Constructor
102  //
103  // Parameters:
104  // name: Name of task
105  //
106  DefineOutput(1, TList::Class());
107  DefineOutput(2, TList::Class());
108 }
109 //_____________________________________________________________________
111  : AliAnalysisTaskSE(o),
112  fVtxAxis(o.fVtxAxis), // Axis to control vertex binning
113  fCentAxis(o.fCentAxis), // Array to control centrality/multiplicity binning
114  fFMDCut(o.fFMDCut), // FMD sigma cut
115  fSPDCut(o.fSPDCut), // SPD sigma cut
116  fFlowFlags(o.fFlowFlags), // Flow flags
117  fEtaGap(o.fEtaGap), // Eta gap value
118  fBinsForward(), // List with forward flow hists
119  fBinsCentral(), // List with central flow hists
120  fSumList(o.fSumList), // Event sum list
121  fOutputList(o.fOutputList), // Result output list
122  fAOD(o.fAOD), // AOD input event
123  fTrackCuts(o.fTrackCuts), // ESD track cuts
124  fMaxMoment(o.fMaxMoment), // Flow moments
125  fVtx(o.fVtx), // Z vertex coordinate
126  fCent(o.fCent), // Centrality
127  fHistdNdedpV0(o.fHistdNdedpV0), // Histo for v0
128  fHistdNdedp3Cor(o.fHistdNdedp3Cor),// Histo for combining detectors
129  fHistFMDSPDCorr(o.fHistFMDSPDCorr),// FMD SPD correlation
130  fHistCent(o.fHistCent), // Hist for centrality
131  fHistVertexSel(o.fHistVertexSel), // Hist for selected vertices
132  fHistEventSel(o.fHistEventSel) // Hist for event selection
133 {
134  //
135  // Copy constructor
136  //
137  // Parameters:
138  // o: Object to copy from
139  //
140 }
141 //_____________________________________________________________________
144 {
145  //
146  // Assignment operator
147  //
148  if (&o == this) return *this;
149  fVtxAxis = o.fVtxAxis;
150  fCentAxis = o.fCentAxis;
151  fFMDCut = o.fFMDCut;
152  fSPDCut = o.fSPDCut;
154  fEtaGap = o.fEtaGap;
155  fSumList = o.fSumList;
157  fAOD = o.fAOD;
160  fVtx = o.fVtx;
161  fCent = o.fCent;
165  fHistCent = o.fHistCent;
168 
169  return *this;
170 }
171 //____________________________________________________________________
172 namespace {
173  Double_t GetEdge(const TString& str, Int_t start, Int_t end)
174  {
175  TString sub(str(start, end));
176  return sub.Atof();
177  }
178  Bool_t ExtractBins(const TString& spec, TArrayD& edges)
179  {
180  TArrayD tmp(200);
181  Int_t start = 0;
182  Int_t cnt = 0;
183  for (Int_t i=1; i<spec.Length(); i++) {
184  if (spec[i] == '-' || spec[i] == ':') {
185  Double_t c = GetEdge(spec, start, i);
186  if (cnt > 0 && c < tmp[cnt-1]) {
187  Warning("ExtractBins",
188  "Invalid edge @ %d: %f (< %f)", cnt, c, tmp[cnt-1]);
189  tmp.Set(0);
190  return false;
191  }
192  tmp[cnt] = c;
193  i++;
194  start = i;
195  cnt++;
196  }
197  }
198  if (start+1 != spec.Length()) {
199  Double_t c = GetEdge(spec, start, spec.Length());
200  tmp[cnt] = c;
201  cnt++;
202  }
203  edges.Set(cnt, tmp.GetArray());
204  return true;
205  }
206 }
207 
208 //________________________________________________________________________
209 void
211 {
212  if (!bins || bins[0] == '\0') return;
213 
214  TString spec(bins);
215  if (spec.EqualTo("none", TString::kIgnoreCase))
216  return;
217 
218  TArrayD edges;
219  if (spec.EqualTo("default", TString::kIgnoreCase) ||
220  spec.EqualTo("pbpb", TString::kIgnoreCase)) {
221  // 1 2 3 4 5 6 7 8 9 10 11
222  Double_t tmp[] = { 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 100 };
223  edges.Set(11, tmp);
224  }
225  else if (spec.EqualTo("ppb", TString::kIgnoreCase) ||
226  spec.EqualTo("pbp", TString::kIgnoreCase)) {
227  // 1 2 3 4 5 6 7 8
228  Double_t tmp[] = { 0, 5, 10, 20, 40, 60, 80, 100 };
229  edges.Set(8, tmp);
230  }
231  else {
232  ExtractBins(spec, edges);
233  }
234  TAxis* centAxis = new TAxis(edges.GetSize()-1, edges.GetArray());
235  SetCentralityAxis(centAxis);
236 }
237 //_____________________________________________________________________
239 {
240  //
241  // Set flow flags, making sure the detector setup is right
242  //
243  // Parameters:
244  // flags: Flow flags
245  //
246  if ((flags & kFMD) && (flags & kVZERO))
247  AliFatal("Cannot do analysis on more than one forward detector!");
248  else if (!(flags & kFMD) && !(flags & kVZERO))
249  AliFatal("You need to add a forward detector!");
250  else fFlowFlags = flags;
251 }
252 //_____________________________________________________________________
254 {
255  //
256  // Create output objects
257  //
258  InitVertexBins();
259  InitHists();
260  if ((fFlowFlags & kTracks) && !fTrackCuts) AliFatal("No track cuts set!");
261  PrintFlowSetup();
262 
263  PostData(1, fSumList);
264 }
265 //_____________________________________________________________________
267 {
268  //
269  // Init vertexbin objects for forward and central detectors, and add them to the lists
270  //
271  for (Int_t v = 1; v <= fVtxAxis->GetNbins(); v++) {
272  Int_t vL = Int_t(fVtxAxis->GetBinLowEdge(v));
273  Int_t vH = Int_t(fVtxAxis->GetBinUpEdge(v));
274  if ((fFlowFlags & kFMD)) {
275  fBinsForward.Add(new VertexBin(vL, vH, fMaxMoment, "FMD", fFlowFlags, fFMDCut, fEtaGap));
276  if (!(fFlowFlags & k3Cor))
277  fBinsCentral.Add(new VertexBin(vL, vH, fMaxMoment, "SPD-FMD", fFlowFlags|kNUAcorr|kSPD, fSPDCut, fEtaGap));
278  }
279  else if ((fFlowFlags & kVZERO)) {
280  fBinsForward.Add(new VertexBin(vL, vH, fMaxMoment, "VZERO", fFlowFlags, 0, fEtaGap));
281  if ((fFlowFlags & kEtaGap) && !(fFlowFlags & kTracks))
282  fBinsCentral.Add(new VertexBin(vL, vH, fMaxMoment, "SPD-VZERO", fFlowFlags|kNUAcorr|kSPD, fSPDCut, fEtaGap));
283  }
284  }
285 }
286 //_____________________________________________________________________
288 {
289  //
290  // Init histograms and add vertex bin histograms to the sum list
291  //
292  if (!fSumList)
293  fSumList = new TList();
294  fSumList->SetName("Sums");
295  fSumList->SetOwner();
296 
297  if (!fVtxAxis) fVtxAxis = new TAxis(20, -10, 10);
298  fVtxAxis->SetName("VtxAxis");
299  if (!fCentAxis) fCentAxis = new TAxis(20, 0, 100);
300  fVtxAxis->SetName("CentAxis");
301 
302  fHistCent = new TH1D("hCent", "Centralities", 100, 0, 100);
303  fHistVertexSel = new TH1D("hVertexSel", "Selected vertices", fVtxAxis->GetNbins(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax());
304  fHistEventSel = new TH1I("hEventSel", "Event Selection", kOK, 0.5, kOK+0.5);
305  fHistEventSel->GetXaxis()->SetBinLabel(kNoEvent, "No AOD event");
306  fHistEventSel->GetXaxis()->SetBinLabel(kNoForward, "No forward det");
307  fHistEventSel->GetXaxis()->SetBinLabel(kNoCentral, "No central det");
308  fHistEventSel->GetXaxis()->SetBinLabel(kNoTrigger, "Not triggered");
309  fHistEventSel->GetXaxis()->SetBinLabel(kNoCent, "No centrality");
310  fHistEventSel->GetXaxis()->SetBinLabel(kInvCent, "Centrality outside range");
311  fHistEventSel->GetXaxis()->SetBinLabel(kNoVtx, "No vertex");
312  fHistEventSel->GetXaxis()->SetBinLabel(kInvVtx, "Vtx outside range");
313  fHistEventSel->GetXaxis()->SetBinLabel(kOK, "OK!");
314 
315  fHistFMDSPDCorr = new TH2D("hFMDSPDCorr", "hFMDSPCCorr", 200, 0., 20000., 200, 0, 7500);
316 
317  TList* dList = new TList();
318  dList->SetName("Diagnostics");
319  dList->Add(fHistCent);
320  dList->Add(fHistVertexSel);
321  dList->Add(fHistEventSel);
322  dList->Add(fHistFMDSPDCorr);
323  fSumList->Add(dList);
324 
325  fHistdNdedp3Cor = TH2D(Form("hdNdedpCombined_%s", GetQCType(fFlowFlags)), Form("hdNdedpCombined_%s", GetQCType(fFlowFlags)),
326  200, -4., 6., 20, 0., TMath::TwoPi());
327  if ((fFlowFlags & kVZERO)) {
328  Double_t bins[12] = { -6, -3.7, -3.2, -2.7, -2.2, -1.7,
329  2.8, 3.4, 3.9, 4.5, 5.1, 6 };
330  fHistdNdedpV0 = TH2D(Form("hdNdedpv0%s", GetQCType(fFlowFlags)), Form("hdNdedpv0%s", GetQCType(fFlowFlags)),
331  11, -6, 6, 8, 0, TMath::TwoPi());
332  fHistdNdedpV0.GetXaxis()->Set(11, bins);
333  if ((fFlowFlags & k3Cor)) {
334  Double_t bins2[20] = { -6, -3.7, -3.2, -2.7, -2.2, // VZERO
335  -2.0, -1.5, -1.0, -0.5 , 0., 0.5, 1.0, 1.5, 2.0, // SPD
336  2.8, 3.4, 3.9, 4.5, 5.1, 6 }; // VZERO
337  fHistdNdedp3Cor.GetXaxis()->Set(19, bins2);
338  fHistdNdedp3Cor.GetYaxis()->Set(8, 0., TMath::TwoPi());
339  }
340  }
341 
342  TIter nextForward(&fBinsForward);
343  VertexBin* bin = 0;
344  while ((bin = static_cast<VertexBin*>(nextForward()))) {
346  }
347  TIter nextCentral(&fBinsCentral);
348  while ((bin = static_cast<VertexBin*>(nextCentral()))) {
350  }
351 }
352 //_____________________________________________________________________
354 {
355  //
356  // Calls the analyze function - called every event
357  //
358  // Parameters:
359  // option: Not used
360  //
361 
362  // Reset data members
363  fCent = -1;
364  fVtx = 1111;
365 
366  Analyze();
367 
368  PostData(1, fSumList);
369 
370  return;
371 }
372 //_____________________________________________________________________
374 {
375  //
376  // Load forward and central detector objects from aod tree and call the
377  // cumulants calculation for the correct vertex bin
378  //
379  // Return: true on success
380  //
381 
382  // Get input event
383  fAOD = dynamic_cast<AliAODEvent*>(AliForwardUtil::GetAODEvent(this));
384  if (!fAOD) {
385  fHistEventSel->Fill(kNoEvent);
386  return kFALSE;
387  }
388 
389  // Get detector objects
390  AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("Forward"));
391  AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClusters"));
392  AliVVZERO* vzero = GetVZERO();
393  if ((fFlowFlags & kVZERO)) {
394  if (vzero) {
395  fHistdNdedpV0.Reset();
396  FillVZEROHist(vzero);
397  }
398  }
399 
400  // We make sure that the necessary forward object is there
401  if ((fFlowFlags & kFMD) && !aodfmult) {
402  fHistEventSel->Fill(kNoForward);
403  return kFALSE;
404  }
405  else if ((fFlowFlags & kVZERO) && !vzero) {
406  fHistEventSel->Fill(kNoForward);
407  return kFALSE;
408  }
409  if (!aodcmult) fHistEventSel->Fill(kNoCentral);
410 
411  // Check event for triggers, get centrality, vtx etc.
412  if (!CheckEvent(aodfmult)) return kFALSE;
413  Int_t vtx = fVtxAxis->FindBin(fVtx)-1;
414 
415  // Then we assign a reference to the forward histogram of interest
416  TH2D& forwarddNdedp = ((fFlowFlags & kFMD) ? aodfmult->GetHistogram() : fHistdNdedpV0);
417  TH2D& spddNdedp = aodcmult->GetHistogram();
418  if ((fFlowFlags & kStdQC)) {
419  FillVtxBinList(fBinsForward, forwarddNdedp, vtx);
420  } else if ((fFlowFlags & kEtaGap)) {
421  FillVtxBinListEtaGap(fBinsForward, forwarddNdedp, forwarddNdedp, vtx);
422  }
423  // At the moment only clusters are supported for the central region (some day add tracks?)
424  // So no extra checks necessary
425  if (aodcmult) {
426  if ((fFlowFlags & kStdQC)) {
427  FillVtxBinList(fBinsCentral, spddNdedp, vtx);
428  } else if ((fFlowFlags & kEtaGap)) {
429  FillVtxBinListEtaGap(fBinsCentral, forwarddNdedp, spddNdedp, vtx);
430  } else if ((fFlowFlags & k3Cor)) {
431  FillVtxBinList3Cor(fBinsForward, spddNdedp, forwarddNdedp, vtx);
432  }
433  // Diagnostics
434  if (aodfmult) {
435  Double_t totForward = forwarddNdedp.Integral(1, forwarddNdedp.GetNbinsX(), 1, forwarddNdedp.GetNbinsY());
436  Double_t totSPD = spddNdedp.Integral(1, spddNdedp.GetNbinsX(), 1, spddNdedp.GetNbinsY());
437  fHistFMDSPDCorr->Fill(totForward, totSPD);
438  }
439  }
440 
441  return kTRUE;
442 }
443 //_____________________________________________________________________
444 void AliForwardFlowTaskQC::FillVtxBinList(const TList& list, TH2D& h, Int_t vtx, UShort_t flags) const
445 {
446  //
447  // Loops over list of VtxBins, fills hists of bins for current vertex
448  // and runs analysis on those bins
449  //
450  // Parameters:
451  // list: list of VtxBins
452  // h: dN/detadphi histogram
453  // vtx: current vertex bin
454  // flags: extra flags to handle calculations.
455  //
456  // Note: The while loop is used in this function and the next 2 for historical reasons,
457  // as originally each moment had it's own VertexBin object.
458  VertexBin* bin = 0;
459  Int_t i = 0;
460  Int_t nVtxBins = fVtxAxis->GetNbins();
461 
462  while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
463  i++;
464  // If no tracks do things normally
465  if (!(fFlowFlags & kTracks) || (flags & kMC)) {
466  if (!bin->FillHists(h, fCent, kFillBoth|flags|kReset)) continue;
467  }
468  // if tracks things are more complicated
469  else if ((fFlowFlags & kTracks)) {
470  if (!FillTracks(bin, kFillRef|kReset|flags)) continue;
471  if (!bin->FillHists(h, fCent, kFillDiff|kReset|flags)) continue;
472  }
474  }
475 
476  return;
477 }
478 //_____________________________________________________________________
480  TH2D& hdiff, Int_t vtx, UShort_t flags) const
481 {
482  //
483  // Loops over list of VtxBins, fills hists of bins for current vertex
484  // and runs analysis on those bins
485  //
486  // Parameters:
487  // list: list of VtxBins
488  // href: dN/detadphi histogram for ref. flow
489  // hdiff: dN/detadphi histogram for diff. flow
490  // vBin: current vertex bin
491  // flags: extra flags to handle calculations.
492  //
493  VertexBin* bin = 0;
494  Int_t i = 0;
495  Int_t nVtxBins = fVtxAxis->GetNbins();
496 
497  while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
498  i++;
499  if (!(fFlowFlags & kTracks) || (flags & kMC)) {
500  if(!bin->FillHists(href, fCent, kFillRef|flags|kReset)) continue;
501  }
502  else if ((fFlowFlags & kTracks)) {
503  if (!FillTracks(bin, kFillRef|kReset|flags)) continue;
504  }
505  if (!bin->FillHists(hdiff, fCent, kFillDiff|kReset|flags)) continue;
507  }
508 
509  return;
510 }
511 //_____________________________________________________________________
513  TH2D& hfwd, Int_t vtx, UShort_t flags)
514 {
515  //
516  // Loops over list of VtxBins, fills hists of bins for current vertex
517  // and runs analysis on those bins
518  //
519  // Parameters:
520  // list: list of VtxBins
521  // hcent: dN/detadphi histogram for central coverage
522  // hfwd: dN/detadphi histogram for forward coverage
523  // vBin: current vertex bin
524  // flags: extra flags to handle calculations.
525  //
526  VertexBin* bin = 0;
527  Int_t i = 0;
528  Int_t nVtxBins = fVtxAxis->GetNbins();
529 
530  TH2D& h = CombineHists(hcent, hfwd);
531 
532  while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
533  i++;
534  if (!bin->FillHists(h, fCent, kFillBoth|flags|kReset)) continue;
536  }
537 
538  return;
539 }
540 //_____________________________________________________________________
542 {
543  //
544  // Combines a forward and central d^2N/detadphi histogram.
545  // At some point it might need a flag to choose which histogram gets
546  // priority when there is an overlap, at the moment the average is chosen
547  //
548  // Parameters:
549  // hcent: Central barrel detector
550  // hfwd: Forward detector
551  //
552  // Return: reference to combined hist
553  //
554 
555  // If hists are the same (MC input) don't do anything
556  if (&hcent == &hfwd) return hcent;
557 
558  fHistdNdedp3Cor.Reset();
559  // FMD, SPD input
560  if ((fFlowFlags & kFMD)) {
561  for (Int_t e = 1; e <= fHistdNdedp3Cor.GetNbinsX(); e++) {
562  Double_t eta = fHistdNdedp3Cor.GetXaxis()->GetBinCenter(e);
563  Bool_t fwdCov = (hfwd.GetBinContent(e, 0) != 0);
564  Bool_t centCov = (hcent.GetBinContent(e, 0) != 0);
565  if (!fwdCov && !centCov) continue;
566  else fHistdNdedp3Cor.SetBinContent(e, 0, 1);
567  for (Int_t p = 1; p <= fHistdNdedp3Cor.GetNbinsY(); p++) {
568  Double_t phi = fHistdNdedp3Cor.GetYaxis()->GetBinCenter(p);
569  Int_t n = 0;
570  Double_t cont = 0.;
571  if (fwdCov) {
572  cont += hfwd.GetBinContent(e, p);
573  n++;
574  }
575  if (centCov) {
576  cont += hcent.GetBinContent(e, p);
577  n++;
578  }
579  if (cont == 0 || n == 0) continue;
580  cont /= n;
581  fHistdNdedp3Cor.Fill(eta, phi, cont);
582  }
583  }
584  // VZERO, SPD input, here we do not average but cut to avoid
585  // (too much) overlap.
586  } else if ((fFlowFlags & kVZERO)) {
587  // VZERO loop
588  for (Int_t eV = 1; eV <= hfwd.GetNbinsX(); eV++) {
589  Double_t eta = hfwd.GetXaxis()->GetBinLowEdge(eV)+0.1;
590  if (hfwd.GetBinContent(eV, 0) == 0) continue;
591  else {
592  Int_t he = fHistdNdedp3Cor.GetXaxis()->FindBin(eta);
593  fHistdNdedp3Cor.SetBinContent(he, 0, 1);
594  }
595  for (Int_t p = 1; p <= hfwd.GetNbinsY(); p++) {
596  Double_t phi = hfwd.GetYaxis()->GetBinCenter(p);
597  Double_t cont = hfwd.GetBinContent(eV, p);
598  fHistdNdedp3Cor.Fill(eta, phi, cont);
599  }
600  }
601  // SPD loop
602  Int_t eSs = hcent.GetXaxis()->FindBin(-1.99);
603  Int_t eSe = hcent.GetXaxis()->FindBin(1.99);
604  for (Int_t eS = eSs; eS <= eSe; eS++) {
605  Double_t eta = hcent.GetXaxis()->GetBinCenter(eS);
606  if (hcent.GetBinContent(eS, 0) == 0) continue;
607  else {
608  Int_t he = fHistdNdedp3Cor.GetXaxis()->FindBin(eta);
609  fHistdNdedp3Cor.SetBinContent(he, 0, 1);
610  }
611  for (Int_t p = 1; p <= hcent.GetNbinsY(); p++) {
612  Double_t phi = hcent.GetYaxis()->GetBinCenter(p);
613  Double_t cont = hcent.GetBinContent(eS, p);
614  fHistdNdedp3Cor.Fill(eta, phi, cont);
615  }
616  }
617  }
618  return fHistdNdedp3Cor;
619 }
620 //_____________________________________________________________________
622 {
623  //
624  // Get TPC tracks to use for reference flow.
625  //
626  // Return: TObjArray with tracks
627  //
628  TObjArray* trList = 0;
629  AliESDEvent* esdEv = 0;
630  if (AliForwardUtil::CheckForAOD() == 1) // AOD tracks
631  trList = static_cast<TObjArray*>(fAOD->GetTracks());
632  else
633  esdEv = dynamic_cast<AliESDEvent*>(InputEvent());
634 
635  Bool_t useEvent = bin->FillTracks(trList, esdEv, fTrackCuts, mode);
636  return useEvent;
637 }
638 //_____________________________________________________________________
640 {
641  //
642  // Calls the finalize function, done at the end of the analysis
643  //
644  // Parameters:
645  // option: Not used
646  //
647 
648  // Make sure pointers are set to the correct lists
649  fSumList = dynamic_cast<TList*> (GetOutputData(1));
650  if(!fSumList) {
651  AliError("Could not retrieve TList fSumList");
652  return;
653  }
654  if (!fOutputList)
655  fOutputList = new TList();
656  fOutputList->SetName("Results");
657  fOutputList->SetOwner();
658 
659  if ((fFlowFlags & kEtaGap) || (fFlowFlags & kTracks)) {
660  TParameter<Double_t>* etaGap = new TParameter<Double_t>("EtaGap", fEtaGap);
661  fOutputList->Add(etaGap);
662  }
663  // We only add axes in terminate, as TAxis object do not merge well,
664  // and so we get a mess when running on the grid if we put them in the sum list...
665  fVtxAxis->SetName("VtxAxis");
666  fOutputList->Add(fVtxAxis);
667  fCentAxis->SetName("CentAxis");
668  fOutputList->Add(fCentAxis);
669 
670  // Run finalize on VertexBins
671  Finalize();
672 
673  // Loop over output to get dN/deta hists - used for diagnostics
674  TIter next(fOutputList);
675  TObject* o = 0;
676  TString name;
677  TH2D* dNdeta = 0;
678  TH1D* cent = 0;
679  while ((o = next())) {
680  name = o->GetName();
681  if (name.Contains("dNdeta")) {
682  dNdeta = dynamic_cast<TH2D*>(o);
683  name.ReplaceAll("dNdeta", "cent");
684  name.ReplaceAll("Ref", "");
685  name.ReplaceAll("Diff", "");
686  cent = dynamic_cast<TH1D*>(fOutputList->FindObject(name.Data()));
687  if (!dNdeta || !cent) continue;
688  for (Int_t cBin = 1; cBin <= dNdeta->GetNbinsY(); cBin++) {
689  Double_t nEvents = cent->GetBinContent(cBin);
690  if (nEvents == 0) continue;
691  for (Int_t eBin = 1; eBin <= dNdeta->GetNbinsX(); eBin++) {
692  dNdeta->SetBinContent(eBin, cBin, dNdeta->GetBinContent(eBin, cBin)/nEvents);
693  dNdeta->SetBinError(eBin, cBin, dNdeta->GetBinError(eBin, cBin)/nEvents);
694  }
695  }
696  }
697  }
698 
699  // Loop over output and make 1D projections for fast look at results
701  TList* vtxList = (TList*)fOutputList->FindObject("vtxList");
702  if (vtxList) MakeCentralityHists(vtxList);
703  TList* nuaList = (TList*)fOutputList->FindObject("NUATerms");
704  TIter nextNua(nuaList);
705  o = 0;
706  TH2D* h = 0;
707  while ((o = nextNua())) {
708  if (!(h = dynamic_cast<TH2D*>(o))) continue;
709  Double_t evts = h->GetBinContent(0, 0);
710  if (evts != 0) h->Scale(1./evts);
711  }
712  if (nuaList) MakeCentralityHists(nuaList);
713 
714  PostData(2, fOutputList);
715 
716  return;
717 }
718 //_____________________________________________________________________
720 {
721  //
722  // Finalize command, called by Terminate()
723  //
724 
725  // Reinitiate vertex bins if Terminate is called separately!
726  if (fBinsForward.GetEntries() == 0) InitVertexBins();
727 
728  // Iterate over all vertex bins objects and finalize cumulants
729  // calculations
732 
733  return;
734 }
735 //_____________________________________________________________________
737 {
738  //
739  // Loop over VertexBin list and call terminate on each
740  //
741  // Parameters:
742  // list: VertexBin list
743  //
744  TIter next(&list);
745  VertexBin* bin = 0;
746  while ((bin = static_cast<VertexBin*>(next()))) {
748  }
749  return;
750 }
751 // _____________________________________________________________________
753 {
754  //
755  // Loop over a list containing a TH2D with flow results
756  // and project to TH1's in specific centrality bins
757  //
758  // Parameters:
759  // list: Flow results list
760  //
761  TH2D* hist2D = 0;
762  TList* centList = 0;
763  TH1D* hist1D = 0;
764  TObject* helper = 0;
765  TIter nextHist(list);
766  while ((helper = dynamic_cast<TObject*>(nextHist()))) {
767  if (!(hist2D = dynamic_cast<TH2D*>(helper))) continue;
768  for (Int_t cBin = 1; cBin <= hist2D->GetNbinsY(); cBin++) {
769  Int_t cMin = Int_t(hist2D->GetYaxis()->GetBinLowEdge(cBin));
770  Int_t cMax = Int_t(hist2D->GetYaxis()->GetBinUpEdge(cBin));
771  TString name = Form("cent_%d-%d", cMin, cMax);
772  centList = (TList*)list->FindObject(name.Data());
773  if (!centList) {
774  centList = new TList();
775  centList->SetName(name.Data());
776  list->Add(centList);
777  }
778  hist1D = hist2D->ProjectionX(Form("%s_%s", hist2D->GetName(), name.Data()),
779  cBin, cBin, "E");
780  hist1D->SetTitle(hist1D->GetName());
781  centList->Add(hist1D);
782  }
783  }
784 }
785 // _____________________________________________________________________
787 {
788  //
789  // Function to check that an AOD event meets the cuts
790  //
791  // Parameters:
792  // AliAODForwardMult: forward mult object with trigger and vertex info
793  //
794  // Return: false if there is no trigger or if the centrality or vertex
795  // is out of range. Otherwise true.
796  //
797 
798  // First check for trigger
799  if (!CheckTrigger(aodfm)) {
800  fHistEventSel->Fill(kNoTrigger);
801  return kFALSE;
802  }
803 
804  // Then check for centrality
805  if (!GetCentrality(aodfm)) {
806  return kFALSE;
807  }
808 
809  // And finally check for vertex
810  if (!GetVertex(aodfm)) {
811  return kFALSE;
812  }
813 
814  // Ev. accepted - filling diag. hists
815  fHistCent->Fill(fCent);
816  fHistVertexSel->Fill(fVtx);
817  fHistEventSel->Fill(kOK);
818 
819  return kTRUE;
820 }
821 // _____________________________________________________________________
823 {
824  //
825  // Function to look for a trigger string in the event.
826  // First check for info in forward mult object, if not there, use the AOD header
827  //
828  // Parameters:
829  // AliAODForwardMult: forward mult object with trigger and vertex info
830  //
831  // Return: true if offline trigger is present
832  //
833  if (aodfm) return aodfm->IsTriggerBits(AliAODForwardMult::kOffline);
834  // this may need to be changed for 2011 data to handle kCentral and so on...
835  else return (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
836  ->IsEventSelected() & AliVEvent::kMB);
837 }
838 // _____________________________________________________________________
840 {
841  //
842  // Function to look get centrality of the event.
843  // First check for info in forward mult object, if not there, use the AOD header
844  //
845  // Parameters:
846  // AliAODForwardMult: forward mult object with trigger and vertex info
847  //
848  // Return: true if centrality determination is present
849  //
850  if (aodfm) {
851  if (aodfm->HasCentrality()) {
852  fCent = (Double_t)aodfm->GetCentrality();
853  if (fCentAxis->GetXmin() > fCent || fCent >= fCentAxis->GetXmax()) {
854  fHistEventSel->Fill(kInvCent);
855  return kFALSE;
856  }
857  }
858  else {
859  fCent = 97.5;
860  fHistEventSel->Fill(kNoCent);
861  }
862  return kTRUE;
863  } else {
864  AliCentrality* aodCent = fAOD->GetCentrality();
865  if (aodCent) {
866  fCent = (Double_t)aodCent->GetCentralityPercentile("V0M");
867  if (fCentAxis->GetXmin() > fCent || fCent >= fCentAxis->GetXmax()) {
868  fHistEventSel->Fill(kInvCent);
869  return kFALSE;
870  }
871  }
872  else {
873  fCent = 97.5;
874  fHistEventSel->Fill(kNoCent);
875  }
876  return kTRUE;
877  }
878 }
879 //_____________________________________________________________________
881 {
882  //
883  // Function to look for vertex determination in the event.
884  // First check for info in forward mult object, if not there, use the AOD header
885  //
886  // Parameters:
887  // AliAODForwardMult: forward mult object with trigger and vertex info
888  //
889  // Return: true if vertex is determined
890  //
891  if (aodfm) {
892  if (aodfm->HasIpZ()) {
893  fVtx = aodfm->GetIpZ();
894  if (fVtx < fVtxAxis->GetXmin() || fVtx >= fVtxAxis->GetXmax()) {
895  fHistEventSel->Fill(kInvVtx);
896  return kFALSE;
897  }
898  } else {
899  fVtx = 9999;
900  fHistEventSel->Fill(kNoVtx);
901  return kFALSE;
902  }
903  return kTRUE;
904  } else {
905  AliAODVertex* aodVtx = fAOD->GetPrimaryVertex();
906  if (aodVtx) {
907  fVtx = aodVtx->GetZ();
908  if (fVtx < fVtxAxis->GetXmin() || fVtx >= fVtxAxis->GetXmax()) {
909  fHistEventSel->Fill(kInvVtx);
910  return kFALSE;
911  }
912  } else {
913  fVtx = 9999;
914  fHistEventSel->Fill(kNoVtx);
915  return kFALSE;
916  }
917  return kTRUE;
918  }
919 }
920 // _____________________________________________________________________
922 {
923  //
924  // Get VZERO object from ESD or AOD
925  //
926  // Return: VZERO data object
927  //
928  AliVVZERO* vzero = 0;
929  // Get input type
931  switch (input) {
932  // If AOD input, simply get the track array from the event
933  case 1: vzero = (AliVVZERO*)fAOD->GetVZEROData();
934  break;
935  case 2: {
936  // If ESD input get event, apply track cuts
937  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
938  if (!esd) return 0;
939  vzero = (AliVVZERO*)esd->GetVZEROData();
940  break;
941  }
942  default: AliFatal("Neither ESD or AOD input. This should never happen");
943  break;
944  }
945  return vzero;
946 }
947 // _____________________________________________________________________
949 {
950  //
951  // Loops over VZERO data object and fill up d^2N/detadphi histogram for flow analysis
952  //
953  // Parameters:
954  // vzero: VZERO AOD data object
955  //
956  Int_t ring = 0;
957  Int_t bin = 0;
958  Double_t eta = 0;
959  // Sort of right for 2010 data, do not use for 2011!
960  Double_t eq[64] = { 1.43536, 1.45727, 1.44993, 1.30051, 1.17425, 1.2335, 1.22247, 1.14362,
961  1.14647, 1.25208, 1.17681, 1.21642, 1.16604, 1.05532, 1.03212, 1.1032,
962  1.22941, 1.36986, 1.14652, 1.20056, 0.927086, 1.10809, 1.03343, 1.29472,
963  1.21204, 1.29217, 1.2003, 2.10382, 1.28513, 1.40558, 1.25784, 1.21848,
964  0.475162, 0.50421, 0.503617, 0.512471, 0.515276, 0.39831, 0.415199, 0.444664,
965  0.521922, 0.785915, 0.703658, 0.832479, 0.77461, 0.73129, 0.778697, 0.710265,
966  0.89686, 0.967688, 0.974225, 0.873445, 0.811096, 0.828493, 0.889609, 0.586056,
967  1.15877, 0.954656, 0.914557, 0.979028, 1.04907, 0.748518, 0.928043, 0.98175 };
968 
969  for (Int_t i = 0; i < 64; i++) {
970  if (i % 8 == 0) {
971  ring++;
972  bin = (ring < 5 ? 11-ring : ring-3);
973  eta = fHistdNdedpV0.GetXaxis()->GetBinCenter(bin);
974  fHistdNdedpV0.SetBinContent(bin, 0, 1);
975  }
976  Float_t amp = vzero->GetMultiplicity(i);
977  amp /= eq[i];
978  Double_t phi = TMath::Pi()/8.+TMath::TwoPi()*i/8.;
979  while (phi > TMath::TwoPi()) phi -= TMath::TwoPi();
980  fHistdNdedpV0.Fill(eta, phi, amp);
981  }
982 }
983 //_____________________________________________________________________
985  : TNamed(),
986  fMaxMoment(0), // Max flow moment for this vertexbin
987  fVzMin(0), // Vertex z-coordinate min [cm]
988  fVzMax(0), // Vertex z-coordinate max [cm]
989  fType(), // Data type name e.g., FMD/SPD/FMDTR/SPDTR/MC
990  fFlags(0), // Flow flags
991  fSigmaCut(-1), // Sigma cut to remove outlier events
992  fEtaGap(-1), // Eta gap value
993  fEtaLims(), // Limits for binning in 3Cor method
994  fCumuRef(), // Histogram for reference flow
995  fCumuDiff(), // Histogram for differential flow
996  fCumuHists(0,0), // CumuHists object for keeping track of results
997  fCumuNUARef(), // Histogram for ref NUA terms
998  fCumuNUADiff(), // Histogram for diff NUA terms
999  fdNdedpRefAcc(), // Diagnostics histogram for acc. maps
1000  fdNdedpDiffAcc(),// Diagnostics histogram for acc. maps
1001  fOutliers(), // Histogram for sigma distribution
1002  fDebug() // Debug level
1003 {
1004  //
1005  // Default constructor
1006  //
1007 }
1008 //_____________________________________________________________________
1010  UShort_t moment, TString name,
1011  UShort_t flags, Double_t cut,
1012  Double_t etaGap)
1013  : TNamed("", ""),
1014  fMaxMoment(moment), // Max flow moment for this vertexbin
1015  fVzMin(vLow), // Vertex z-coordinate min [cm]
1016  fVzMax(vHigh), // Vertex z-coordinate max [cm]
1017  fType(name), // Data type name e.g., FMD/SPD/FMDTR/SPDTR/MC
1018  fFlags(flags), // Flow flags
1019  fSigmaCut(cut), // Sigma cut to remove outlier events
1020  fEtaGap(etaGap), // Eta gap value
1021  fEtaLims(), // Limits for binning in 3Cor method
1022  fCumuRef(), // Histogram for reference flow
1023  fCumuDiff(), // Histogram for differential flow
1024  fCumuHists(moment,0),// CumuHists object for keeping track of results
1025  fCumuNUARef(), // Histogram for ref NUA terms
1026  fCumuNUADiff(), // Histogram for diff NUA terms
1027  fdNdedpRefAcc(), // Diagnostics histogram for acc. maps
1028  fdNdedpDiffAcc(), // Diagnostics histogram for acc. maps
1029  fOutliers(), // Histogram for sigma distribution
1030  fDebug(0) // Debug level
1031 {
1032  //
1033  // Constructor
1034  //
1035  // Parameters
1036  // vLow: min z-coordinate
1037  // vHigh: max z-coordinate
1038  // moment: max flow moment
1039  // name: data type name (FMD/SPD/FMDTR/SPDTR/MC)
1040  // flags: flow flags
1041  // cut: sigma cut
1042  // etaGap: eta-gap value
1043  //
1044  fType.ToUpper();
1045 
1046  SetName(Form("%svertexBin%d_%d_%d%s", fType.Data(), moment, vLow, vHigh, GetQCType(fFlags)));
1047  SetTitle(Form("%svertexBin%d_%d_%d%s", fType.Data(), moment, vLow, vHigh, GetQCType(fFlags)));
1048 
1049  fDebug = AliAnalysisManager::GetAnalysisManager()->GetDebugLevel();
1050  if (fDebug > 0) Printf("AliForwardFlowTaskQC::VertexBin()\tDebugMode: %d", fDebug);
1051 
1052  // Set limits for 3 correlator method
1053  if ((fFlags & kFMD)) {
1054  fEtaLims[0] = -6.;
1055  fEtaLims[1] = -1.5;
1056  fEtaLims[2] = -0.5;
1057  fEtaLims[3] = 2.;
1058  fEtaLims[4] = 3.;
1059  fEtaLims[5] = 6.;
1060  } else if ((fFlags & kVZERO)) {
1061  fEtaLims[0] = -6;
1062  fEtaLims[1] = -2.7;
1063  fEtaLims[2] = -2.0;
1064  fEtaLims[3] = 2.0;
1065  fEtaLims[4] = 3.9;
1066  fEtaLims[5] = 6;
1067  }
1068 }
1069 //_____________________________________________________________________
1072 {
1073  //
1074  // Assignment operator
1075  //
1076  // Parameters
1077  // o: AliForwardFlowTaskQC::VertexBin
1078  //
1079  if (&o == this) return *this;
1080  fMaxMoment = o.fMaxMoment;
1081  fVzMin = o.fVzMin;
1082  fVzMax = o.fVzMax;
1083  fType = o.fType;
1084  fFlags = o.fFlags;
1085  fSigmaCut = o.fSigmaCut;
1086  fEtaGap = o.fEtaGap;
1087  fCumuRef = o.fCumuRef;
1088  fCumuDiff = o.fCumuDiff;
1089  fCumuHists = o.fCumuHists;
1094  fOutliers = o.fOutliers;
1095  fDebug = o.fDebug;
1096  for (UInt_t i = 0; i < sizeof(fEtaLims)/sizeof(Double_t); i++) fEtaLims[i] = o.fEtaLims[i];
1097 
1098  return *this;
1099 }
1100 //_____________________________________________________________________
1102 {
1103  //
1104  // Add histograms to outputlist
1105  //
1106  // Parameters
1107  // outputlist: list of histograms
1108  // centAxis: centrality axis
1109  //
1110 
1111  // First we try to find an outputlist for this vertexbin
1112  TList* list = (TList*)outputlist->FindObject(Form("%svertex_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1113  // If it doesn't exist we make one
1114  if (!list) {
1115  list = new TList();
1116  list->SetName(Form("%svertex_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1117  outputlist->Add(list);
1118  }
1119 
1120  // Get bin numbers and binning defined
1121  Int_t nHBins = GetBinNumberSin();
1122  Int_t nEtaBins = 24;
1123  if ((fFlags & k3Cor)) {
1124  if ((fFlags & kFMD)) nEtaBins = 24;
1125  else if ((fFlags & kVZERO)) nEtaBins = 19;
1126  }
1127  else if ((fFlags & kVZERO) && (fFlags & kEtaGap)) nEtaBins = 19;
1128  else if ((fFlags & kVZERO)) nEtaBins = 11;
1129 
1130  Double_t vzeroBins[12] = { -6, -3.7, -3.2, -2.7, -2.2, -1.7,
1131  2.8, 3.4, 3.9, 4.5, 5.1, 6 };
1132  Double_t vzeroBins2[20] = { -6, -3.7, -3.2, -2.7, -2.2, // VZERO
1133  -2.0, -1.5, -1.0, -0.5 , 0., 0.5, 1.0, 1.5, 2.0, // SPD
1134  2.8, 3.4, 3.9, 4.5, 5.1, 6 }; // VZERO
1135 
1136  Int_t nRefBins = nEtaBins; // needs to be something as default
1137  if ((fFlags & kStdQC)) {
1138  if ((fFlags & kSymEta) && !((fFlags & kTracks) && (fFlags & kSPD))) nRefBins = 1;
1139  else nRefBins = 2;
1140  } else if ((fFlags & kEtaGap )) {
1141  nRefBins = 2;
1142  } else if ((fFlags & k3Cor)) {
1143  nRefBins = 24;
1144  }
1145 
1146  // We initiate the reference histogram
1147  fCumuRef = new TH2D(Form("%s_%d_%d%s_ref", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1148  Form("%s_%d_%d%s_ref", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1149  nRefBins, -6., 6.,
1150  nHBins, 0.5, nHBins+0.5);
1151  if ((fFlags & kVZERO) && (fFlags & k3Cor)) fCumuRef->GetXaxis()->Set(nEtaBins, vzeroBins2);
1152  SetupNUALabels(fCumuRef->GetYaxis());
1153  list->Add(fCumuRef);
1154  // We initiate the differential histogram
1155  fCumuDiff = new TH2D(Form("%s_%d_%d%s_diff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1156  Form("%s_%d_%d%s_diff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1157  nEtaBins, -6., 6., nHBins, 0.5, nHBins+0.5);
1158  if ((fFlags & kVZERO)) {
1159  if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1160  fCumuDiff->GetXaxis()->Set(nEtaBins, vzeroBins2);
1161  else fCumuDiff->GetXaxis()->Set(nEtaBins, vzeroBins);
1162  }
1163  SetupNUALabels(fCumuDiff->GetYaxis());
1164  list->Add(fCumuDiff);
1165 
1166  // Cumulants sum hists
1167  Int_t cBins = centAxis->GetNbins();
1168  fCumuHists.ConnectList(Form("%sCumu_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)), list);
1169  TH3D* cumuHist = 0;
1170  Int_t nC2Bins = ((fFlags & kEtaGap) || (fFlags & k3Cor) ? kW2 : k3pWeight);
1171  Int_t nC4Bins = ((fFlags & kEtaGap) ? kW2 : ((fFlags & k3Cor) ? kW4 : kSinphi1phi2phi3p));
1172  for (Int_t n = 2; n <= fMaxMoment; n++) {
1173  // Initiate the ref cumulant sum histogram
1174  cumuHist = new TH3D(Form("%sv%d_vertex_%d_%d%s_cumuRef", fType.Data(), n, fVzMin, fVzMax, GetQCType(fFlags)),
1175  Form("%sv%d_vertex_%d_%d%s_cumuRef", fType.Data(), n, fVzMin, fVzMax, GetQCType(fFlags)),
1176  nRefBins, -6., 6.,
1177  cBins, 0., 100., nC2Bins, 0.5, nC2Bins+0.5);
1178  if ((fFlags & kVZERO) && (fFlags & k3Cor)) cumuHist->GetXaxis()->Set(nEtaBins, vzeroBins2);
1179  cumuHist->GetYaxis()->Set(cBins, centAxis->GetXbins()->GetArray());
1180  fCumuHists.Add(cumuHist);
1181  // Initiate the diff cumulant sum histogram
1182  cumuHist = new TH3D(Form("%sv%d_vertex_%d_%d%s_cumuDiff", fType.Data(), n, fVzMin, fVzMax, GetQCType(fFlags)),
1183  Form("%sv%d_vertex_%d_%d%s_cumuDiff", fType.Data(), n, fVzMin, fVzMax, GetQCType(fFlags)),
1184  nEtaBins, -6., 6., cBins, 0., 100., nC4Bins, 0.5, nC4Bins+0.5);
1185  if ((fFlags & kVZERO)) {
1186  if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1187  cumuHist->GetXaxis()->Set(nEtaBins, vzeroBins2);
1188  else cumuHist->GetXaxis()->Set(nEtaBins, vzeroBins);
1189  }
1190  cumuHist->GetYaxis()->Set(cBins, centAxis->GetXbins()->GetArray());
1191  fCumuHists.Add(cumuHist);
1192  }
1193 
1194  // Common NUA histograms
1195  fCumuNUARef = new TH3D(Form("%s_vertex_%d_%d%s_cumuNUARef", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1196  Form("%s_vertex_%d_%d%s_cumuNUARef", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1197  nRefBins, -6., 6.,
1198  cBins, 0., 100., nHBins, 0.5, nHBins+0.5);
1199  if ((fFlags & kVZERO) && (fFlags & k3Cor)) fCumuNUARef->GetXaxis()->Set(nEtaBins, vzeroBins2);
1200  fCumuNUARef->GetYaxis()->Set(cBins, centAxis->GetXbins()->GetArray());
1201  SetupNUALabels(fCumuNUARef->GetZaxis());
1202  fCumuNUARef->Sumw2();
1203  list->Add(fCumuNUARef);
1204 
1205  fCumuNUADiff = new TH3D(Form("%s_vertex_%d_%d%s_cumuNUADiff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1206  Form("%s_vertex_%d_%d%s_cumuNUADiff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1207  nEtaBins, -6., 6., cBins, 0., 100., nHBins, 0.5, nHBins+0.5);
1208  if ((fFlags & kVZERO)) {
1209  if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1210  fCumuNUADiff->GetXaxis()->Set(nEtaBins, vzeroBins2);
1211  else fCumuNUADiff->GetXaxis()->Set(nEtaBins, vzeroBins);
1212  }
1213  fCumuNUADiff->GetYaxis()->Set(cBins, centAxis->GetXbins()->GetArray());
1214  SetupNUALabels(fCumuNUADiff->GetZaxis());
1215  fCumuNUADiff->Sumw2();
1216  list->Add(fCumuNUADiff);
1217 
1218  // We create diagnostic histograms.
1219  TList* dList = (TList*)outputlist->FindObject("Diagnostics");
1220  if (!dList) AliFatal("No diagnostics list found");
1221 
1222  // Acceptance hist
1223  Int_t nPhiBins = ((fFlags & kFMD) ? 20 : 8);
1224  fdNdedpRefAcc = new TH2F(Form("h%sdNdedpRefAcc_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1225  Form("%s reference flow acceptance map for %d cm < v_{z} < %d cm", fType.Data(), fVzMin, fVzMax),
1226  nEtaBins, -6., 6., nPhiBins, 0, TMath::TwoPi());
1227  if ((fFlags & kVZERO)) {
1228  if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1229  fdNdedpRefAcc->GetXaxis()->Set(nEtaBins, vzeroBins2);
1230  else fdNdedpRefAcc->GetXaxis()->Set(nEtaBins, vzeroBins);
1231  }
1232  dList->Add(fdNdedpRefAcc);
1233 
1234  fdNdedpDiffAcc = new TH2F(Form("h%sdNdedpDiffAcc_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1235  Form("%s differential flow acceptance map for %d cm < v_{z} < %d cm", fType.Data(), fVzMin, fVzMax),
1236  nEtaBins, -6., 6., nPhiBins, 0, TMath::TwoPi());
1237  if ((fFlags & kVZERO)) {
1238  if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1239  fdNdedpDiffAcc->GetXaxis()->Set(nEtaBins, vzeroBins2);
1240  else fdNdedpDiffAcc->GetXaxis()->Set(nEtaBins, vzeroBins);
1241  }
1242  dList->Add(fdNdedpDiffAcc);
1243 
1244  fOutliers = new TH2F(Form("hOutliers_%s_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1245  Form("Maximum #sigma from mean N_{ch} pr. bin - %s, %d < v_{z} < %d",
1246  fType.Data(), fVzMin, fVzMax),
1247  20, 0., 100., 500, 0., ((fFlags & kMC) ? 15. : 5.));
1248  dList->Add(fOutliers);
1249 
1250  return;
1251 }
1252 //_____________________________________________________________________
1254 {
1255  //
1256  // Fill reference and differential eta-histograms
1257  //
1258  // Parameters:
1259  // dNdetadphi: 2D histogram with input data
1260  // cent: centrality
1261  // mode: filling mode: kFillRef/kFillDiff/kFillBoth
1262  //
1263  if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1264  Bool_t useEvent = kTRUE;
1265 
1266  // Fist we reset histograms
1267  if ((mode & kReset)) {
1268  if ((mode & kFillRef)) fCumuRef->Reset();
1269  if ((mode & kFillDiff)) fCumuDiff->Reset();
1270  }
1271  // Then we loop over the input and calculate sum cos(k*n*phi)
1272  // and fill it in the reference and differential histograms
1273  Int_t nBadBins = 0;
1274  Double_t limit = 9999.;
1275  for (Int_t etaBin = 1; etaBin <= dNdetadphi.GetNbinsX(); etaBin++) {
1276  Double_t eta = dNdetadphi.GetXaxis()->GetBinCenter(etaBin);
1277  // Numbers to cut away bad events
1278  Double_t runAvg = 0;
1279  Double_t max = 0;
1280  Int_t nInAvg = 0;
1281  Double_t avgSqr = 0;
1282  for (Int_t phiBin = 0; phiBin <= dNdetadphi.GetNbinsY(); phiBin++) {
1283  // Check for acceptance
1284  if (phiBin == 0) {
1285  if (dNdetadphi.GetBinContent(etaBin, 0) == 0) break;
1286  // Central limit for eta gap break for reference flow
1287  if ((fFlags & kEtaGap) && (mode & kFillRef) &&
1288  TMath::Abs(eta) < fEtaGap) break;
1289  // Backward and forward eta gap break for reference flow
1290  if ((fFlags & kEtaGap) && (mode & kFillRef) && TMath::Abs(eta) > TMath::Abs(limit)) break;
1291  if ((fFlags & kStdQC) && (fFlags & kMC) && !(fFlags & kTracks)) {
1292  if (!(fFlags & kSPD) && TMath::Abs(eta) < 1.75) break;
1293  if ((fFlags & kSPD) && TMath::Abs(eta) > 2.00) break;
1294  }
1295  if (limit > 1e3) limit = dNdetadphi.GetXaxis()->GetBinLowEdge(etaBin);
1296  continue;
1297  } // End of phiBin == 0
1298  Double_t phi = dNdetadphi.GetYaxis()->GetBinCenter(phiBin);
1299  Double_t weight = dNdetadphi.GetBinContent(etaBin, phiBin);
1300 
1301  // We calculate the average Nch per. bin
1302  avgSqr += weight*weight;
1303  runAvg += weight;
1304  nInAvg++;
1305  if (weight == 0) continue;
1306  if (weight > max) max = weight;
1307  // Fill into Cos() and Sin() hists
1308  if ((mode & kFillRef) && !((fFlags & kTracks) && (fFlags & kMC) && TMath::Abs(eta) > 0.75)) {
1309  fCumuRef->Fill(eta, 0., weight);// mult goes in underflowbin - no visual, but not needed?
1310  fdNdedpRefAcc->Fill(eta, phi, weight);
1311  }
1312  if ((mode & kFillDiff)) {
1313  fCumuDiff->Fill(eta, 0., weight);
1314  fdNdedpDiffAcc->Fill(eta, phi, weight);
1315  }
1316  for (Int_t n = 1; n <= 2*fMaxMoment; n++) {
1317  Double_t cosBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberCos(n));
1318  Double_t sinBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberSin(n));
1319  Double_t cosnPhi = weight*TMath::Cos(n*phi);
1320  Double_t sinnPhi = weight*TMath::Sin(n*phi);
1321  // fill ref
1322  if ((mode & kFillRef) && !((fFlags & kTracks) && (fFlags & kMC) && TMath::Abs(eta) > 0.75)) {
1323  fCumuRef->Fill(eta, cosBin, cosnPhi);
1324  fCumuRef->Fill(eta, sinBin, sinnPhi);
1325  }
1326  // fill diff
1327  if ((mode & kFillDiff)) {
1328  fCumuDiff->Fill(eta, cosBin, cosnPhi);
1329  fCumuDiff->Fill(eta, sinBin, sinnPhi);
1330  }
1331  } // End of NUA loop
1332  } // End of phi loop
1333  // Outlier cut calculations
1334  if (nInAvg > 0) {
1335  runAvg /= nInAvg;
1336  avgSqr /= nInAvg;
1337  Double_t stdev = (nInAvg > 1 ? TMath::Sqrt(nInAvg/(nInAvg-1))*TMath::Sqrt(avgSqr - runAvg*runAvg) : 0);
1338  Double_t nSigma = (stdev == 0 ? 0 : (max-runAvg)/stdev);
1339  if (fSigmaCut > 0. && nSigma >= fSigmaCut && cent < 60) nBadBins++;
1340  else nBadBins = 0;
1341  fOutliers->Fill(cent, nSigma);
1342  // We still finish the loop, for fOutliers to make sense,
1343  // but we do no keep the event for analysis
1344  if (nBadBins > 3) useEvent = kFALSE;
1345  }
1346  } // End of eta bin
1347 
1348  return useEvent;
1349 }
1350 //_____________________________________________________________________
1352  AliAnalysisFilter* trFilter, UShort_t mode)
1353 {
1354  //
1355  // Fill reference and differential eta-histograms
1356  //
1357  // Parameters:
1358  // trList: list of tracks
1359  // mode: filling mode: kFillRef/kFillDiff/kFillBoth
1360  //
1361  if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1362  if (!trList && !esd) {
1363  AliError("FillTracks: No AOD track list or ESD event - something might be wrong!");
1364  return kFALSE;
1365  }
1366 
1367  // Fist we reset histograms
1368  if ((mode & kReset)) {
1369  if ((mode & kFillRef)) fCumuRef->Reset();
1370  if ((mode & kFillDiff)) fCumuDiff->Reset();
1371  }
1372 
1373  // Then we loop over the input and calculate sum cos(k*n*phi)
1374  // and fill it in the reference and differential histograms
1375  Int_t nTr = 0;
1376  if (trList) nTr = trList->GetEntries();
1377  if (esd) nTr = esd->GetNumberOfTracks();
1378  if (nTr == 0) return kFALSE;
1379  AliVTrack* tr = 0;
1380  // Cuts for AOD tracks (have already been applied to ESD tracks) - except dEdx
1381  // const tpcdEdx = 10;
1382  for (Int_t i = 0; i < nTr; i++) { // track loop
1383  tr = (trList ? (AliVTrack*)trList->At(i) : (AliVTrack*)esd->GetTrack(i));
1384  if (!tr) continue;
1385  if (esd) {
1386  AliESDtrack* esdTr = (AliESDtrack*)tr;
1387  if (!trFilter->IsSelected(esdTr)) continue;
1388  }
1389  else if (trList) { // If AOD input
1390  Double_t pTMin = 0, pTMax = 0, etaMin = 0, etaMax = 0, minNCl = 0;
1391  UInt_t bit = 0;
1392  if ((fFlags & kTPC) == kTPC) pTMin = 0.2, pTMax = 5., etaMin = -0.8, etaMax = 0.8, minNCl = 70, bit = 128;
1393  if ((fFlags & kHybrid) == kHybrid) pTMin = 0.2, pTMax = 5., etaMin = -0.8, etaMax = 0.8, minNCl = 70, bit = 272;
1394 
1395  AliAODTrack* aodTr = (AliAODTrack*)tr;
1396  if (aodTr->GetID() > -1) continue;
1397  if (!aodTr->TestFilterBit(bit) || aodTr->Pt() > pTMax || aodTr->Pt() < pTMin ||
1398  aodTr->Eta() > etaMax || aodTr->Eta() < etaMin || aodTr->GetTPCNcls() < minNCl) continue;
1399  }
1400 
1401  // if (tr->GetTPCsignal() < tpcdEdx) continue;
1402  // Track accepted
1403  Double_t eta = tr->Eta();
1404  if (((fFlags & kSPD) || (fFlags & kEtaGap)) && TMath::Abs(eta) < fEtaGap) continue;
1405  Double_t phi = tr->Phi();
1406  Double_t weight = 1.;
1407 
1408  if ((mode & kFillRef)) {
1409  fCumuRef->Fill(eta, 0., weight);// mult goes in underflowbin - no visual, but not needed?
1410  fdNdedpRefAcc->Fill(eta, phi, weight);
1411  }
1412  if ((mode & kFillDiff)) {
1413  fCumuDiff->Fill(eta, 0., weight);
1414  fdNdedpDiffAcc->Fill(eta, phi, weight);
1415  }
1416  for (Int_t n = 1; n <= 2*fMaxMoment; n++) {
1417  Double_t cosBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberCos(n));
1418  Double_t sinBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberSin(n));
1419  Double_t cosnPhi = weight*TMath::Cos(n*phi);
1420  Double_t sinnPhi = weight*TMath::Sin(n*phi);
1421  // fill ref
1422  if ((mode & kFillRef)) {
1423  fCumuRef->Fill(eta, cosBin, cosnPhi);
1424  fCumuRef->Fill(eta, sinBin, sinnPhi);
1425  }
1426  // fill diff
1427  if ((mode & kFillDiff)) {
1428  fCumuDiff->Fill(eta, cosBin, cosnPhi);
1429  fCumuDiff->Fill(eta, sinBin, sinnPhi);
1430  }
1431  } // End of NUA loop
1432  } // End of track loop
1433  return kTRUE;
1434 }
1435 //_____________________________________________________________________
1437 {
1438  //
1439  // Calculate the Q cumulant up to order fMaxMoment
1440  //
1441  // Parameters:
1442  // cent: centrality of event
1443  //
1444  if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1445 
1446  // Fill out NUA hists
1447  for (Int_t etaBin = 1; etaBin <= fCumuRef->GetNbinsX(); etaBin++) {
1448  Double_t eta = fCumuRef->GetXaxis()->GetBinCenter(etaBin);
1449  if (fCumuRef->GetBinContent(etaBin, 0) <= 3) continue;
1450  if ((fFlags & kTracks) && (fFlags & kSPD) && !(fFlags & kEtaGap)) eta = -eta;
1451  for (Int_t qBin = 0; qBin <= fCumuRef->GetNbinsY(); qBin++) {
1452  fCumuNUARef->Fill(eta, cent, Double_t(qBin), fCumuRef->GetBinContent(etaBin, qBin));
1453  }
1454  }
1455  for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1456  Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1457  Int_t refetaBin = fCumuRef->GetXaxis()->FindBin(eta);
1458  if (fCumuRef->GetBinContent(refetaBin, 0) <= 3) continue;
1459  if (fCumuDiff->GetBinContent(etaBin, 0) == 0) continue;
1460  for (Int_t qBin = 0; qBin <= fCumuDiff->GetNbinsY(); qBin++) {
1461  fCumuNUADiff->Fill(eta, cent, Double_t(qBin), fCumuDiff->GetBinContent(etaBin, qBin));
1462  }
1463  }
1464 
1465  // We create the objects needed for the analysis
1466  TH3D* cumuRef = 0;
1467  TH3D* cumuDiff = 0;
1468  // For each n we loop over the hists
1469  for (Int_t n = 2; n <= fMaxMoment; n++) {
1470  cumuRef = (TH3D*)fCumuHists.Get('r',n);
1471  cumuDiff = (TH3D*)fCumuHists.Get('d',n);
1472  Int_t prevRefEtaBin = 0;
1473 
1474  // Per mom. quantities
1475  Double_t dQnReA = 0, dQnImA = 0, multA = 0;
1476  Double_t dQnReB = 0, dQnImB = 0, multB = 0;
1477  Double_t dQ2nReA = 0, dQ2nImA = 0;
1478  Double_t two = 0, w2 = 0, four = 0, w4 = 0;
1479  Double_t cosPhi1Phi2 = 0, cosPhi1Phi2Phi3m = 0;
1480  Double_t sinPhi1Phi2 = 0, sinPhi1Phi2Phi3m = 0;
1481  for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1482  Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1483  Double_t refEta = eta;
1484  if ((fFlags & kTracks) && (fFlags & kSPD) && !(fFlags & kEtaGap)) refEta = -eta;
1485  Int_t refEtaBinA = fCumuRef->GetXaxis()->FindBin(refEta);
1486  if ((fFlags & kEtaGap)) refEta = -eta;
1487  Int_t refEtaBinB = fCumuRef->GetXaxis()->FindBin(refEta);
1488  if (refEtaBinA != prevRefEtaBin) {
1489  // Reference flow
1490  multA = fCumuRef->GetBinContent(refEtaBinA, 0);
1491  dQnReA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberCos(n));
1492  dQnImA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberSin(n));
1493  dQ2nReA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberCos(n*2));
1494  dQ2nImA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberSin(n*2));
1495 
1496  multB = fCumuRef->GetBinContent(refEtaBinB, 0);
1497  dQnReB = fCumuRef->GetBinContent(refEtaBinB, GetBinNumberCos(n));
1498  dQnImB = fCumuRef->GetBinContent(refEtaBinB, GetBinNumberSin(n));
1499 
1500  if (multA <= 3 || multB <= 3) continue;
1501  // The reference flow is calculated
1502  // 2-particle
1503  if ((fFlags & kStdQC)) {
1504  w2 = multA * (multA - 1.);
1505  two = dQnReA*dQnReA + dQnImA*dQnImA - multA;
1506  } else {
1507  w2 = multA * multB;
1508  two = dQnReA*dQnReB + dQnImA*dQnImB;
1509  }
1510  cumuRef->Fill(eta, cent, kW2Two, two);
1511  cumuRef->Fill(eta, cent, kW2, w2);
1512 
1513  // The reference flow is calculated
1514  // 4-particle
1515  if ((fFlags & kStdQC)) {
1516  w4 = multA * (multA - 1.) * (multA - 2.) * (multA - 3.);
1517 
1518  four = 2.*multA*(multA-3.) + TMath::Power((TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.)),2.)
1519  -4.*(multA-2.)*(TMath::Power(dQnReA,2.) + TMath::Power(dQnImA,2.))
1520  -2.*(TMath::Power(dQnReA,2.)*dQ2nReA+2.*dQnReA*dQnImA*dQ2nImA-TMath::Power(dQnImA,2.)*dQ2nReA)
1521  +(TMath::Power(dQ2nReA,2.)+TMath::Power(dQ2nImA,2.));
1522 
1523  cumuRef->Fill(eta, cent, kW4Four, four);
1524  cumuRef->Fill(eta, cent, kW4, w4);
1525 
1526  // NUA
1527  cosPhi1Phi2 = dQnReA*dQnReA - dQnImA*dQnImA - dQ2nReA;
1528  sinPhi1Phi2 = 2.*dQnReA*dQnImA - dQ2nImA;
1529 
1530  cosPhi1Phi2Phi3m = dQnReA*(TMath::Power(dQnReA,2)+TMath::Power(dQnImA,2))
1531  -dQnReA*dQ2nReA-dQnImA*dQ2nImA-2.*(multA-1)*dQnReA;
1532 
1533  sinPhi1Phi2Phi3m = -dQnImA*(TMath::Power(dQnReA,2)+TMath::Power(dQnImA,2))
1534  +dQnReA*dQ2nImA-dQnImA*dQ2nReA+2.*(multA-1)*dQnImA;
1535 
1536  cumuRef->Fill(eta, cent, kCosphi1phi2, cosPhi1Phi2);
1537  cumuRef->Fill(eta, cent, kSinphi1phi2, sinPhi1Phi2);
1538  cumuRef->Fill(eta, cent, kCosphi1phi2phi3m, cosPhi1Phi2Phi3m);
1539  cumuRef->Fill(eta, cent, kSinphi1phi2phi3m, sinPhi1Phi2Phi3m);
1540  cumuRef->Fill(eta, cent, k3pWeight, multA*(multA-1.)*(multA-2.));
1541  } // End of QC{4}
1542  prevRefEtaBin = refEtaBinA;
1543  } // End of reference flow
1544  // For each etaBin bin the necessary values for differential flow is calculated
1545  Double_t mp = fCumuDiff->GetBinContent(etaBin, 0);
1546  Double_t pnRe = fCumuDiff->GetBinContent(etaBin, GetBinNumberCos(n));
1547  Double_t pnIm = fCumuDiff->GetBinContent(etaBin, GetBinNumberSin(n));
1548  Double_t p2nRe = fCumuDiff->GetBinContent(etaBin, GetBinNumberCos(n*2));
1549  Double_t p2nIm = fCumuDiff->GetBinContent(etaBin, GetBinNumberSin(n*2));
1550  if (mp == 0) continue;
1551  Double_t mq = 0;
1552  Double_t qnRe = 0;
1553  Double_t qnIm = 0;
1554  Double_t q2nRe = 0;
1555  Double_t q2nIm = 0;
1556 
1557  // Differential flow calculations for each eta bin is done:
1558  // 2-particle differential flow
1559  if ((fFlags & kStdQC) && (!(fFlags & kTracks) || ((fFlags & kTracks) && (fFlags & kMC) && !(fFlags & kSPD) && TMath::Abs(eta) < 0.75))) {
1560  mq = mp;
1561  qnRe = pnRe;
1562  qnIm = pnIm;
1563  q2nRe = p2nRe;
1564  q2nIm = p2nIm;
1565  }
1566 
1567  Double_t w2p = mp * multB - mq;
1568  Double_t twoPrime = pnRe*dQnReB + pnIm*dQnImB - mq;
1569 
1570  cumuDiff->Fill(eta, cent, kW2Two, twoPrime);
1571  cumuDiff->Fill(eta, cent, kW2, w2p);
1572 
1573  if ((fFlags & kEtaGap)) continue;
1574  // Differential flow calculations for each eta bin bin is done:
1575  // 4-particle differential flow
1576  Double_t w4p = (mp * multA - 3.*mq)*(multA - 1.)*(multA - 2.);
1577 
1578  Double_t fourPrime = (TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.))*(pnRe*dQnReA+pnIm*dQnImA)
1579  - q2nRe*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))
1580  - 2.*q2nIm*dQnReA*dQnImA
1581  - pnRe*(dQnReA*dQ2nReA+dQnImA*dQ2nImA)
1582  + pnIm*(dQnImA*dQ2nReA-dQnReA*dQ2nImA)
1583  - 2.*multA*(pnRe*dQnReA+pnIm*dQnImA)
1584  - 2.*(TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.))*mq
1585  + 6.*(qnRe*dQnReA+qnIm*dQnImA)
1586  + 1.*(q2nRe*dQ2nReA+q2nIm*dQ2nImA)
1587  + 2.*(pnRe*dQnReA+pnIm*dQnImA)
1588  + 2.*mq*multA
1589  - 6.*mq;
1590 
1591  cumuDiff->Fill(eta, cent, kW4Four, fourPrime);
1592  cumuDiff->Fill(eta, cent, kW4, w4p);
1593 
1594  // NUA
1595  Double_t cosPsi1Phi2 = pnRe*dQnReA - pnIm*dQnImA - q2nRe;
1596  Double_t sinPsi1Phi2 = pnRe*dQnImA + pnIm*dQnReA - q2nIm;
1597 
1598  Double_t cosPsi1Phi2Phi3p = pnRe*(TMath::Power(dQnImA,2.)+TMath::Power(dQnReA,2.)-multA)
1599  - 1.*(q2nRe*dQnReA+q2nIm*dQnImA)
1600  - mq*dQnReA+2.*qnRe;
1601 
1602  Double_t sinPsi1Phi2Phi3p = pnIm*(TMath::Power(dQnImA,2.)+TMath::Power(dQnReA,2.)-multA)
1603  - 1.*(q2nIm*dQnReA-q2nRe*dQnImA)
1604  - mq*dQnImA+2.*qnIm;
1605 
1606  Double_t cosPsi1Phi2Phi3m = pnRe*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))+2.*pnIm*dQnReA*dQnImA
1607  - 1.*(pnRe*dQ2nReA+pnIm*dQ2nImA)
1608  - 2.*mq*dQnReA+2.*qnRe;
1609 
1610  Double_t sinPsi1Phi2Phi3m = pnIm*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))-2.*pnRe*dQnReA*dQnImA
1611  - 1.*(pnIm*dQ2nReA-pnRe*dQ2nImA)
1612  + 2.*mq*dQnImA-2.*qnIm;
1613 
1614  cumuDiff->Fill(eta, cent, kCosphi1phi2, cosPsi1Phi2);
1615  cumuDiff->Fill(eta, cent, kSinphi1phi2, sinPsi1Phi2);
1616  cumuDiff->Fill(eta, cent, kCosphi1phi2phi3m, cosPsi1Phi2Phi3m);
1617  cumuDiff->Fill(eta, cent, kSinphi1phi2phi3m, sinPsi1Phi2Phi3m);
1618  cumuDiff->Fill(eta, cent, k3pWeight, (mp*multA-2.*mq)*(multA-1.));
1619  cumuDiff->Fill(eta, cent, kCosphi1phi2phi3p, cosPsi1Phi2Phi3p);
1620  cumuDiff->Fill(eta, cent, kSinphi1phi2phi3p, sinPsi1Phi2Phi3p);
1621  } // End of eta loop
1622  // Event count
1623  cumuRef->Fill(-7., cent, -0.5, 1.);
1624  } // End of moment loop
1625  return;
1626 }
1627 //_____________________________________________________________________
1629  Int_t& bLow, Int_t& bHigh) const
1630 {
1631  //
1632  // Get the limits for the 3 correlator method
1633  //
1634  // Parameters:
1635  // bin : reference bin #
1636  // aLow : Lowest bin to be included in v_A calculations
1637  // aHigh: Highest bin to be included in v_A calculations
1638  // bLow : Lowest bin to be included in v_B calculations
1639  // bHigh: Highest bin to be included in v_B calculations
1640  //
1641  if ((fFlags & kFMD)) {
1642  switch(bin) {
1643  case 0:
1644  aLow = 14; aHigh = 15;
1645  bLow = 20; bHigh = 22;
1646  break;
1647  case 1:
1648  aLow = 16; aHigh = 16;
1649  bLow = 21; bHigh = 22;
1650  break;
1651  case 2:
1652  aLow = 6; aHigh = 7;
1653  bLow = 21; bHigh = 22;
1654  break;
1655  case 3:
1656  aLow = 6; aHigh = 7;
1657  bLow = 12; bHigh = 12;
1658  break;
1659  case 4:
1660  aLow = 6; aHigh = 8;
1661  bLow = 13; bHigh = 14;
1662  break;
1663  default:
1664  AliFatal(Form("No limits for this eta region! (%d)", bin));
1665  }
1666  }
1667  else if ((fFlags & kVZERO)) {
1668  switch(bin) {
1669  case 0:
1670  aLow = 6; aHigh = 13;
1671  bLow = 17; bHigh = 18;
1672  break;
1673  case 1:
1674  aLow = 6; aHigh = 9;
1675  bLow = 17; bHigh = 18;
1676  break;
1677  case 2:
1678  aLow = 2; aHigh = 3;
1679  bLow = 17; bHigh = 18;
1680  break;
1681  case 3:
1682  aLow = 2; aHigh = 3;
1683  bLow = 6; bHigh = 9;
1684  break;
1685  case 4:
1686  aLow = 2; aHigh = 3;
1687  bLow = 6; bHigh = 13;
1688  break;
1689  default:
1690  AliFatal(Form("No limits for this eta region! (%d)", bin));
1691  }
1692  }
1693  // Try to catch cases where fEtaLimits and these values do not correspond to each other
1694  if (aHigh > fCumuNUARef->GetNbinsX() || bHigh > fCumuNUARef->GetNbinsX())
1695  AliFatal(Form("Limits outside vtx range! (%d) - aHigh = %d, bHigh = %d, Nbins = %d",
1696  bin, aHigh, bHigh, fCumuNUARef->GetNbinsX()));
1697 }
1698 //_____________________________________________________________________
1700 {
1701  //
1702  // Calculate the Q cumulant up to order fMaxMoment
1703  //
1704  // Parameters:
1705  // cent: centrality of event
1706  //
1707  if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1708 
1709  // Fill out NUA hists
1710  for (Int_t etaBin = 1; etaBin <= fCumuRef->GetNbinsX(); etaBin++) {
1711  Double_t eta = fCumuRef->GetXaxis()->GetBinCenter(etaBin);
1712  if (fCumuRef->GetBinContent(etaBin, 0) == 0) continue;
1713  for (Int_t qBin = 0; qBin <= fCumuRef->GetNbinsY(); qBin++) {
1714  fCumuNUARef->Fill(eta, cent, Double_t(qBin), fCumuRef->GetBinContent(etaBin, qBin));
1715  }
1716  }
1717  for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1718  Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1719  if (fCumuDiff->GetBinContent(etaBin, 0) == 0) continue;
1720  for (Int_t qBin = 0; qBin <= fCumuDiff->GetNbinsY(); qBin++) {
1721  fCumuNUADiff->Fill(eta, cent, Double_t(qBin), fCumuDiff->GetBinContent(etaBin, qBin));
1722  }
1723  }
1724 
1725  // We create the objects needed for the analysis
1726  TH3D* cumuRef = 0;
1727  TH3D* cumuDiff = 0;
1728  // For each n we loop over the hists
1729  for (Int_t n = 2; n <= fMaxMoment; n++) {
1730  cumuRef = (TH3D*)fCumuHists.Get('r',n);
1731  cumuDiff = (TH3D*)fCumuHists.Get('d',n);
1732 
1733  // Per mom. quantities
1734  Int_t prevLim = 0;
1735  Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
1736  Double_t dQnReA = 0, dQnImA = 0, multA = 0;
1737  Double_t dQnReB = 0, dQnImB = 0, multB = 0;
1738  Double_t two = 0, w2 = 0;
1739  for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1740  Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1741  if (fEtaLims[prevLim] < eta) {
1742  GetLimits(prevLim, aLow, aHigh, bLow, bHigh);
1743  prevLim++;
1744  multA = 0; dQnReA = 0; dQnImA = 0;
1745  multB = 0; dQnReB = 0; dQnImB = 0;
1746  // Reference flow
1747  for (Int_t a = aLow; a <= aHigh; a++) {
1748  multA += fCumuRef->GetBinContent(a, 0);
1749  dQnReA += fCumuRef->GetBinContent(a, GetBinNumberCos(n));
1750  dQnImA += fCumuRef->GetBinContent(a, GetBinNumberSin(n));
1751  }
1752  for (Int_t b = bLow; b <= bHigh; b++) {
1753  multB += fCumuRef->GetBinContent(b, 0);
1754  dQnReB += fCumuRef->GetBinContent(b, GetBinNumberCos(n));
1755  dQnImB += fCumuRef->GetBinContent(b, GetBinNumberSin(n));
1756  }
1757  // The reference flow is calculated
1758  // 2-particle
1759  w2 = multA * multB;
1760  two = dQnReA*dQnReB + dQnImA*dQnImB;
1761  } // End of reference flow
1762  cumuRef->Fill(eta, cent, kW2Two, two);
1763  cumuRef->Fill(eta, cent, kW2, w2);
1764 
1765  // For each etaBin bin the necessary values for differential flow is calculated
1766  Double_t mp = fCumuDiff->GetBinContent(etaBin, 0);
1767  Double_t pnRe = fCumuDiff->GetBinContent(etaBin, GetBinNumberCos(n));
1768  Double_t pnIm = fCumuDiff->GetBinContent(etaBin, GetBinNumberSin(n));
1769  if (mp == 0) continue;
1770 
1771  // Differential flow calculations for each eta bin is done:
1772  // 2-particle differential flow
1773  Double_t w2pA = mp * multA;
1774  Double_t twoPrimeA = pnRe*dQnReA + pnIm*dQnImA;
1775  cumuDiff->Fill(eta, cent, kW2Two, twoPrimeA);
1776  cumuDiff->Fill(eta, cent, kW2, w2pA);
1777 
1778  Double_t w2pB = mp * multB;
1779  Double_t twoPrimeB = pnRe*dQnReB + pnIm*dQnImB;
1780  cumuDiff->Fill(eta, cent, kW4Four, twoPrimeB);
1781  cumuDiff->Fill(eta, cent, kW4, w2pB);
1782  } // End of eta loop
1783  // Event count
1784  cumuRef->Fill(-7., cent, -0.5, 1.);
1785  } // End of moment loop
1786  return;
1787 
1788 }
1789 //_____________________________________________________________________
1791 {
1792  //
1793  // Finalizes the Q cumulant calculations
1794  //
1795  // Parameters:
1796  // inlist: input sumlist
1797  // outlist: output result list
1798  //
1799 
1800  // Re-find cumulants hist if Terminate is called separately
1801  if (!fCumuHists.IsConnected()) {
1802  TList* list = (TList*)inlist->FindObject(Form("%svertex_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1803  fCumuHists.ConnectList(Form("%sCumu_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)), list);
1804 
1805  if (!fCumuNUARef)
1806  fCumuNUARef = (TH3D*)list->FindObject(Form("%s_vertex_%d_%d%s_cumuNUARef", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1807  if (!fCumuNUADiff)
1808  fCumuNUADiff = (TH3D*)list->FindObject(Form("%s_vertex_%d_%d%s_cumuNUADiff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1809  }
1810  // Clone to avoid normalization problems when redoing terminate locally
1811  fCumuNUARef = (TH3D*)fCumuNUARef->Clone(Form("%s_vertex_%d_%d%s_cumuNUARefNorm", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1812  fCumuNUADiff = (TH3D*)fCumuNUADiff->Clone(Form("%s_vertex_%d_%d%s_cumuNUADiffNorm", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1813 
1814  // Diagnostics histograms
1815  TH2I* quality = (TH2I*)outlist->FindObject(Form("hQCQuality%s%s", fType.Data(), GetQCType(fFlags)));
1816  if (!quality) {
1817  quality = MakeQualityHist(Form("hQCQuality%s%s", fType.Data(), GetQCType(fFlags)));
1818  outlist->Add(quality);
1819  }
1820  TH1D* cent = (TH1D*)outlist->FindObject(Form("%s%s_cent", fType.Data(), GetQCType(fFlags)));
1821  if (!cent) {
1822  cent = new TH1D(Form("%s%s_cent", fType.Data(), GetQCType(fFlags)),
1823  Form("%s%s_cent", fType.Data(), GetQCType(fFlags)),
1824  fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXmin(), fCumuNUARef->GetYaxis()->GetXmax());
1825  cent->GetXaxis()->Set(fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXbins()->GetArray());
1826  outlist->Add(cent);
1827  }
1828  TH2D* dNdetaRef = (TH2D*)outlist->FindObject(Form("%s%s_dNdetaRef", fType.Data(), GetQCType(fFlags)));
1829  if (!dNdetaRef) {
1830  dNdetaRef = new TH2D(Form("%s%s_dNdetaRef", fType.Data(), GetQCType(fFlags)),
1831  Form("%s%s_dNdetaRef", fType.Data(), GetQCType(fFlags)),
1832  fCumuNUARef->GetNbinsX(), fCumuNUARef->GetXaxis()->GetXmin(), fCumuNUARef->GetXaxis()->GetXmax(),
1833  fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXmin(), fCumuNUARef->GetYaxis()->GetXmax());
1834  dNdetaRef->GetYaxis()->Set(fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXbins()->GetArray());
1835  dNdetaRef->Sumw2();
1836  outlist->Add(dNdetaRef);
1837  }
1838  TH2D* dNdetaDiff = (TH2D*)outlist->FindObject(Form("%s%s_dNdetaDiff", fType.Data(), GetQCType(fFlags)));
1839  if (!dNdetaDiff) {
1840  dNdetaDiff = new TH2D(Form("%s%s_dNdetaDiff", fType.Data(), GetQCType(fFlags)),
1841  Form("%s%s_dNdetaDiff", fType.Data(), GetQCType(fFlags)),
1842  fCumuNUADiff->GetNbinsX(), fCumuNUADiff->GetXaxis()->GetXmin(), fCumuNUADiff->GetXaxis()->GetXmax(),
1843  fCumuNUADiff->GetNbinsY(), fCumuNUADiff->GetYaxis()->GetXmin(), fCumuNUADiff->GetYaxis()->GetXmax());
1844  dNdetaDiff->GetYaxis()->Set(fCumuNUADiff->GetNbinsY(), fCumuNUADiff->GetYaxis()->GetXbins()->GetArray());
1845  dNdetaDiff->Sumw2();
1846  outlist->Add(dNdetaDiff);
1847  }
1848 
1849  // Setting up outputs
1850  // Create output lists and diagnostics
1851  TList* vtxList = (TList*)outlist->FindObject("vtxList");
1852  if (!vtxList) {
1853  vtxList = new TList();
1854  vtxList->SetName("vtxList");
1855  outlist->Add(vtxList);
1856  }
1857  vtxList->Add(fCumuNUARef);
1858  vtxList->Add(fCumuNUADiff);
1859 
1860  // Setup output profiles
1861  CumuHistos cumu2(fMaxMoment, ((fFlags & kNUAcorr) ? 2 : 0));
1862  CumuHistos cumu4(fMaxMoment, ((fFlags & kNUAcorr) ? 1 : 0));
1863 
1864  cumu2.ConnectList(Form("%sQC2_Cumu%s_vtx_%d_%d", fType.Data(), GetQCType(fFlags), fVzMin, fVzMax), vtxList);
1865  if ((fFlags & kStdQC))
1866  cumu4.ConnectList(Form("%sQC4_Cumu%s_vtx_%d_%d", fType.Data(), GetQCType(fFlags), fVzMin, fVzMax), vtxList);
1867 
1868  for (Int_t n = 2; n <= fMaxMoment; n++) {
1869  // 2-particle
1870  cumu2.Add(MakeOutputHist(2, n, "Ref", CumuHistos::kNoNUA));
1871  if ((fFlags & k3Cor)){
1872  cumu2.Add(MakeOutputHist(2, n, "DiffA", CumuHistos::kNoNUA));
1873  cumu2.Add(MakeOutputHist(2, n, "DiffB", CumuHistos::kNoNUA));
1874  } else {
1875  cumu2.Add(MakeOutputHist(2, n, "Diff", CumuHistos::kNoNUA));
1876  }
1877  // 4-particle
1878  if ((fFlags & kStdQC)) {
1879  cumu4.Add(MakeOutputHist(4, n, "Ref", CumuHistos::kNoNUA));
1880  cumu4.Add(MakeOutputHist(4, n, "Diff", CumuHistos::kNoNUA));
1881  }
1882  } // End of v_n result loop
1883  // NUA corrected
1884  if ((fFlags & kNUAcorr)) {
1885  for (Int_t n = 2; n <= fMaxMoment; n++) {
1886  // 2-particle
1887  cumu2.Add(MakeOutputHist(2, n, "Ref", CumuHistos::kNUAOld));
1888  if ((fFlags & k3Cor)) {
1889  cumu2.Add(MakeOutputHist(2, n, "DiffA", CumuHistos::kNUAOld));
1890  cumu2.Add(MakeOutputHist(2, n, "DiffB", CumuHistos::kNUAOld));
1891  } else {
1892  cumu2.Add(MakeOutputHist(2, n, "Diff", CumuHistos::kNUAOld));
1893  }
1894  // 4-particle
1895  if ((fFlags & kStdQC)) {
1896  cumu4.Add(MakeOutputHist(4, n, "Ref", CumuHistos::kNUAOld));
1897  cumu4.Add(MakeOutputHist(4, n, "Diff", CumuHistos::kNUAOld));
1898  }
1899  }
1900  for (Int_t n = 2; n <= fMaxMoment; n++) {
1901  // 2-particle
1902  cumu2.Add(MakeOutputHist(2, n, "Ref", CumuHistos::kNUA));
1903  if ((fFlags & k3Cor)) {
1904  cumu2.Add(MakeOutputHist(2, n, "DiffA", CumuHistos::kNUA));
1905  cumu2.Add(MakeOutputHist(2, n, "DiffB", CumuHistos::kNUA));
1906  } else {
1907  cumu2.Add(MakeOutputHist(2, n, "Diff", CumuHistos::kNUA));
1908  }
1909  }
1910  }
1911 
1912  // Calculating the cumulants
1913  if ((fFlags & k3Cor)) {
1914  Calculate3CorFlow(cumu2, quality, cent, dNdetaRef, dNdetaDiff);
1915  } else {
1916  CalculateReferenceFlow(cumu2, cumu4, quality, cent, dNdetaRef);
1917  CalculateDifferentialFlow(cumu2, cumu4, quality, dNdetaDiff);
1918  }
1919  if ((fFlags & kNUAcorr)) {
1920  SolveCoupledFlowEquations(cumu2, 'r');
1921  if ((fFlags & k3Cor)) {
1922  SolveCoupledFlowEquations(cumu2, 'a');
1923  SolveCoupledFlowEquations(cumu2, 'b');
1924  } else {
1925  SolveCoupledFlowEquations(cumu2, 'd');
1926  }
1927  }
1928 
1929  // Add to output for immediate viewing - individual vtx bins are used for final results
1930  AddVertexBins(cumu2, outlist, ((fFlags & kNUAcorr) ? 2 : 0));
1931  if ((fFlags & kStdQC)) AddVertexBins(cumu4, outlist, ((fFlags & kNUAcorr) ? 1 : 0));
1932 
1933  // Setup NUA diagnoastics histograms
1934  TList* nualist = (TList*)outlist->FindObject("NUATerms");
1935  if (!nualist) {
1936  nualist = new TList();
1937  nualist->SetName("NUATerms");
1938  outlist->Add(nualist);
1939  }
1940  // Reference
1941  TH2D* nuaRef = (TH2D*)nualist->FindObject(Form("%sReferenceNUA%s", fType.Data(), GetQCType(fFlags)));
1942  TH2D* temp = 0;
1943  if (!nuaRef) {
1944  nuaRef = (TH2D*)fCumuNUARef->Project3D("yz");
1945  nuaRef->Scale(1./fCumuNUARef->GetNbinsX());
1946  nuaRef->SetName(Form("%sReferenceNUA%s", fType.Data(), GetQCType(fFlags)));
1947  nuaRef->SetTitle(Form("%sReferenceNUA%s", fType.Data(), GetQCType(fFlags)));
1948  nualist->Add(nuaRef);
1949  } else {
1950  temp = (TH2D*)fCumuNUARef->Project3D("yz");
1951  temp->Scale(1./fCumuNUARef->GetNbinsX());
1952  nuaRef->Add(temp);
1953  delete temp;
1954  }
1955  // Filling in underflow to make scaling possible in Terminate()
1956  nuaRef->Fill(0., -1., 1.);
1957  // Differential
1958  TH2D* nuaDiff = (TH2D*)nualist->FindObject(Form("%sDifferentialNUA%s", fType.Data(), GetQCType(fFlags)));
1959  if (!nuaDiff) {
1960  nuaDiff = (TH2D*)fCumuNUADiff->Project3D("yz");
1961  nuaDiff->SetName(Form("%sDifferentialNUA%s", fType.Data(), GetQCType(fFlags)));
1962  nuaDiff->SetTitle(Form("%sDifferentialNUA%s", fType.Data(), GetQCType(fFlags)));
1963  nualist->Add(nuaDiff);
1964  } else {
1965  temp = (TH2D*)fCumuNUADiff->Project3D("yz");
1966  nuaDiff->Add(temp);
1967  delete temp;
1968  }
1969  // Filling in underflow to make scaling possible in Terminate()
1970  nuaDiff->Fill(0., -1., 1.);
1971 
1972  return;
1973 }
1974 //_____________________________________________________________________
1976  TH1D* chist, TH2D* dNdetaRef) const
1977 {
1978  //
1979  // Calculates the reference flow
1980  //
1981  // Parameters:
1982  // cumu2h: CumuHistos object with QC{2} cumulants
1983  // cumu4h: CumuHistos object with QC{4} cumulants
1984  // quality: Histogram for success rate of cumulants
1985  // chist: Centrality histogram
1986  // dNdetaRef: dN/deta histogram for estimating multiplicity used for ref calculations
1987  //
1988 
1989  // Normalizing common NUA hists
1990  for (Int_t cBin = 1; cBin <= fCumuNUARef->GetNbinsY(); cBin++) {
1991  Double_t cent = fCumuNUARef->GetYaxis()->GetBinCenter(cBin);
1992  for (Int_t eBin = 1; eBin <= fCumuNUARef->GetNbinsX(); eBin++) {
1993  Double_t eta = fCumuNUARef->GetXaxis()->GetBinCenter(eBin);
1994  Double_t mult = fCumuNUARef->GetBinContent(eBin, cBin, 0);
1995  if (mult == 0) continue;
1996  for (Int_t qBin = 1; qBin <= fCumuNUARef->GetNbinsZ(); qBin++) {
1997  fCumuNUARef->SetBinContent(eBin, cBin, qBin, fCumuNUARef->GetBinContent(eBin, cBin, qBin)/mult);
1998  fCumuNUARef->SetBinError(eBin, cBin, qBin, fCumuNUARef->GetBinError(eBin, cBin, qBin)/mult);
1999  }
2000  // Fill dN/deta diagnostics
2001  dNdetaRef->Fill(eta, cent, mult);
2002  }
2003  }
2004 
2005  // For flow calculations
2006  TH3D* cumuRef = 0;
2007  TH2D* cumu2 = 0;
2008  TH2D* cumu2NUAold = 0;
2009  TH2D* cumu4 = 0;
2010  TH2D* cumu4NUA = 0;
2011  Int_t qualityFactor = ((fFlags & kStdQC) ? 8 : 4); // used for correctly filling in quality hists
2012  // Loop over cumulant histogram for final calculations
2013  for (Int_t n = 2; n <= fMaxMoment; n++) { // Moment loop begins
2014  cumu2 = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNoNUA);
2015  if ((fFlags & kNUAcorr)) {
2016  cumu2NUAold = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNUAOld);
2017  }
2018  if ((fFlags & kStdQC)) {
2019  cumu4 = (TH2D*)cumu4h.Get('r', n, CumuHistos::kNoNUA);
2020  if ((fFlags & kNUAcorr)) cumu4NUA = (TH2D*)cumu4h.Get('r', n, CumuHistos::kNUAOld);
2021  }
2022  cumuRef = (TH3D*)fCumuHists.Get('r', n);
2023  // Begin loops
2024  for (Int_t cBin = 1; cBin <= cumuRef->GetNbinsY(); cBin++) { // Centrality loop begins
2025  Double_t cent = cumuRef->GetYaxis()->GetBinCenter(cBin);
2026  if (n == 2) chist->Fill(cent, cumuRef->GetBinContent(0, cBin, 0));
2027  if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), n, cent));
2028  for (Int_t etaBin = 1; etaBin <= cumuRef->GetNbinsX(); etaBin++) { // Eta loop begins
2029  Double_t eta = cumuRef->GetXaxis()->GetBinCenter(etaBin);
2030  Double_t refEta = eta;
2031  Int_t refEtaBinA = fCumuNUARef->GetXaxis()->FindBin(refEta);
2032  if ((fFlags & kEtaGap)) refEta = -eta;
2033  Int_t refEtaBinB = fCumuNUARef->GetXaxis()->FindBin(refEta);
2034  // 2-particle reference flow
2035  Double_t w2Two = cumuRef->GetBinContent(refEtaBinA, cBin, kW2Two);
2036  Double_t w2 = cumuRef->GetBinContent(refEtaBinA, cBin, kW2);
2037  if (w2 == 0) continue;
2038  Double_t cosP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberCos(n));
2039  Double_t sinP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberSin(n));
2040  Double_t cosP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(n));
2041  Double_t sinP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(n));
2042  Double_t cos2nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberCos(2*n));
2043  Double_t sin2nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberSin(2*n));
2044  Double_t cos2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(2*n));
2045  Double_t sin2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(2*n));
2046  Double_t two = w2Two / w2;
2047  Double_t qc2 = two;
2048  if (qc2 >= 0) cumu2->Fill(eta, cent, TMath::Sqrt(qc2));
2049 
2050  if ((fFlags & kNUAcorr)) {
2051  // Old NUA
2052  // With no eta gap the last two terms are <<cos(phi)>>^2 and <<sin(phi)>>^2,
2053  // with eta gap the different coverage is taken into account.
2054  // The next line covers both cases.
2055  qc2 -= cosP1nPhiA*cosP1nPhiB + sinP1nPhiA*sinP1nPhiB;
2056  // Extra NUA term from 2n cosines and sines
2057  Double_t den = 1+(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB);
2058  if (den != 0) qc2 /= den;
2059  else qc2 = 0;
2060  }
2061  if (qc2 <= 0) {
2062  if (fDebug > 0)
2063  AliInfo(Form("%s: QC_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f - skipping",
2064  fType.Data(), n, qc2, eta, cent));
2065  quality->Fill((n-2)*qualityFactor+2, Int_t(cent));
2066  continue;
2067  }
2068  Double_t vnTwo = TMath::Sqrt(qc2);
2069  if (!TMath::IsNaN(vnTwo)) {
2070  quality->Fill((n-2)*qualityFactor+1, Int_t(cent));
2071  if ((fFlags & kNUAcorr)) cumu2NUAold->Fill(eta, cent, vnTwo);
2072  }
2073 
2074  if (!(fFlags & kStdQC)) continue;
2075  // 4-particle reference flow
2076  Double_t w4Four = cumuRef->GetBinContent(refEtaBinA, cBin, kW4Four);
2077  Double_t w4 = cumuRef->GetBinContent(refEtaBinA, cBin, kW4);
2078  Double_t multm1m2 = cumuRef->GetBinContent(refEtaBinA, cBin, k3pWeight);
2079  if (w4 == 0 || multm1m2 == 0) continue;
2080  Double_t cosP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2);
2081  Double_t sinP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2);
2082  Double_t cosP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2phi3m);
2083  Double_t sinP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2phi3m);
2084 
2085  cosP1nPhi1P1nPhi2 /= w2;
2086  sinP1nPhi1P1nPhi2 /= w2;
2087  cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2088  sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2089  Double_t four = w4Four / w4;
2090  Double_t qc4 = four-2.*TMath::Power(two,2.);
2091  if (qc4 < 0) cumu4->Fill(eta, cent, TMath::Power(-qc4, 0.25));
2092 
2093  if ((fFlags & kNUAcorr)) {
2094  qc4 += - 4.*cosP1nPhiA*cosP1nPhi1M1nPhi2M1nPhi3
2095  + 4.*sinP1nPhiA*sinP1nPhi1M1nPhi2M1nPhi3-TMath::Power(cosP1nPhi1P1nPhi2,2.)-TMath::Power(sinP1nPhi1P1nPhi2,2.)
2096  + 4.*cosP1nPhi1P1nPhi2*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
2097  + 8.*sinP1nPhi1P1nPhi2*sinP1nPhiA*cosP1nPhiA
2098  + 8.*two*(TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.))
2099  - 6.*TMath::Power((TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.)),2.);
2100  }
2101  if (qc4 >= 0) {
2102  if (fDebug > 0)
2103  AliInfo(Form("%s: QC_%d{4} = %1.3f for eta = %1.2f and centrality %3.1f - skipping",
2104  fType.Data(), n, qc2, eta, cent));
2105  quality->Fill((n-2)*qualityFactor+6, Int_t(cent));
2106  continue;
2107  }
2108  Double_t vnFour = TMath::Power(-qc4, 0.25);
2109  if (!TMath::IsNaN(vnFour*multm1m2)) {
2110  quality->Fill((n-2)*qualityFactor+5, Int_t(cent));
2111  if ((fFlags & kNUAcorr)) cumu4NUA->Fill(eta, cent, vnFour);
2112  }
2113  } // End of eta
2114  } // End of cent
2115  } // End of moment
2116 
2117  return;
2118 }
2119 //_____________________________________________________________________
2121  TH2I* quality, TH2D* dNdetaDiff) const
2122 {
2123  //
2124  // Calculates the differential flow
2125  //
2126  // Parameters:
2127  // cumu2h: CumuHistos object with QC{2} cumulants
2128  // cumu4h: CumuHistos object with QC{4} cumulants
2129  // quality: Histogram for success rate of cumulants
2130  // dNdetaDiff: dN/deta histogram for estimating multiplicity used for diff calculations
2131  //
2132 
2133  for (Int_t cBin = 1; cBin <= fCumuNUADiff->GetNbinsY(); cBin++) {
2134  Double_t cent = fCumuNUADiff->GetYaxis()->GetBinCenter(cBin);
2135  for (Int_t eBin = 1; eBin <= fCumuNUADiff->GetNbinsX(); eBin++) {
2136  Double_t eta = fCumuNUADiff->GetXaxis()->GetBinCenter(eBin);
2137  Double_t mult = fCumuNUADiff->GetBinContent(eBin, cBin, 0);
2138  if (mult == 0) continue;
2139  for (Int_t qBin = 1; qBin <= fCumuNUADiff->GetNbinsZ(); qBin++) {
2140  fCumuNUADiff->SetBinContent(eBin, cBin, qBin, fCumuNUADiff->GetBinContent(eBin, cBin, qBin)/mult);
2141  fCumuNUADiff->SetBinError(eBin, cBin, qBin, fCumuNUADiff->GetBinError(eBin, cBin, qBin)/mult);
2142  }
2143  dNdetaDiff->Fill(eta, cent, mult);
2144  }
2145  }
2146 
2147  // For flow calculations
2148  TH3D* cumuRef = 0;
2149  TH3D* cumuDiff = 0;
2150  TH2D* cumu2 = 0;
2151  TH2D* cumu2NUAold = 0;
2152  TH2D* cumu4 = 0;
2153  TH2D* cumu4NUA = 0;
2154  Int_t qualityFactor = ((fFlags & kStdQC) ? 8 : 4); // used for correctly filling in quality hists
2155  // Loop over cumulant histogram for final calculations
2156  for (Int_t n = 2; n <= fMaxMoment; n++) { // Moment loop begins
2157  cumu2 = (TH2D*)cumu2h.Get('d', n, CumuHistos::kNoNUA);
2158  if ((fFlags & kNUAcorr)) {
2159  cumu2NUAold = (TH2D*)cumu2h.Get('d', n, CumuHistos::kNUAOld);
2160  }
2161  if ((fFlags & kStdQC)) {
2162  cumu4 = (TH2D*)cumu4h.Get('d',n);
2163  if ((fFlags & kNUAcorr)) cumu4NUA = (TH2D*)cumu4h.Get('d', n, CumuHistos::kNUAOld);
2164  }
2165  cumuRef = (TH3D*)fCumuHists.Get('r',n);
2166  cumuDiff = (TH3D*)fCumuHists.Get('d',n);
2167  for (Int_t cBin = 1; cBin <= cumuDiff->GetNbinsY(); cBin++) { // Centrality loop begins
2168  Double_t cent = cumuDiff->GetYaxis()->GetBinCenter(cBin);
2169  if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), n, cent));
2170  for (Int_t etaBin = 1; etaBin <= cumuDiff->GetNbinsX(); etaBin++) { // Eta loop begins
2171  Double_t eta = cumuDiff->GetXaxis()->GetBinCenter(etaBin);
2172  Double_t refEta = eta;
2173  Int_t refEtaBinA = fCumuNUARef->GetXaxis()->FindBin(refEta);
2174  if ((fFlags & kEtaGap)) refEta = -eta;
2175  Int_t refEtaBinB = fCumuNUARef->GetXaxis()->FindBin(refEta);
2176 
2177  // Reference objects
2178  Double_t w2 = cumuRef->GetBinContent(refEtaBinA, cBin, kW2);
2179  if (w2 == 0) continue;
2180  Double_t two = cumuRef->GetBinContent(refEtaBinA, cBin, kW2Two);
2181  two /= w2;
2182  Double_t cosP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberCos(n));
2183  Double_t sinP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberSin(n));
2184  Double_t cosP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(n));
2185  Double_t sinP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(n));
2186  Double_t cos2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(2*n));
2187  Double_t sin2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(2*n));
2188 
2189  // 2-particle differential flow
2190  Double_t w2pTwoPrime = cumuDiff->GetBinContent(etaBin, cBin, kW2Two);
2191  Double_t w2p = cumuDiff->GetBinContent(etaBin, cBin, kW2);
2192  if (w2p == 0) continue;
2193  Double_t cosP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(n));
2194  Double_t sinP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(n));
2195  Double_t cos2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(2*n));
2196  Double_t sin2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(2*n));
2197  Double_t twoPrime = w2pTwoPrime / w2p;
2198 
2199  Double_t qc2Prime = twoPrime;
2200  cumu2->Fill(eta, cent, qc2Prime);
2201  if ((fFlags & kNUAcorr)) {
2202  // Old nua
2203  qc2Prime -= cosP1nPsi*cosP1nPhiB + sinP1nPsi*sinP1nPhiB;
2204  // Extra NUA term from 2n cosines and sines
2205  qc2Prime /= (1.+(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
2206  }
2207  if (!TMath::IsNaN(qc2Prime)) {
2208  quality->Fill((n-2)*qualityFactor+3, Int_t(cent));
2209  if ((fFlags & kNUAcorr)) cumu2NUAold->Fill(eta, cent, qc2Prime);
2210  }
2211  else
2212  quality->Fill((n-2)*qualityFactor+4, Int_t(cent));
2213  if (fDebug > 1)
2214  AliInfo(Form("%s: QC'_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f",
2215  fType.Data(), n, qc2Prime, eta, cent));
2216 
2217  if (!(fFlags & kStdQC)) continue;
2218  // Reference objects
2219  Double_t cosP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2);
2220  Double_t sinP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2);
2221  Double_t cosP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2phi3m);
2222  Double_t sinP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2phi3m);
2223  Double_t multm1m2 = cumuRef->GetBinContent(refEtaBinA, cBin, k3pWeight);
2224  cosP1nPhi1P1nPhi2 /= w2;
2225  sinP1nPhi1P1nPhi2 /= w2;
2226  cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2227  sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2228 
2229  // 4-particle differential flow
2230  Double_t w4pFourPrime = cumuDiff->GetBinContent(etaBin, cBin, kW4Four);
2231  Double_t w4p = cumuDiff->GetBinContent(etaBin, cBin, kW4);
2232  Double_t mpqMult = cumuDiff->GetBinContent(etaBin, cBin, k3pWeight);
2233  if (w4p == 0 || mpqMult == 0) continue;
2234  Double_t cosP1nPsi1P1nPhi2 = cumuDiff->GetBinContent(etaBin, cBin, kCosphi1phi2);
2235  Double_t sinP1nPsi1P1nPhi2 = cumuDiff->GetBinContent(etaBin, cBin, kSinphi1phi2);
2236  Double_t cosP1nPsi1M1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kCosphi1phi2phi3m);
2237  Double_t sinP1nPsi1M1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kSinphi1phi2phi3m);
2238  Double_t cosP1nPsi1P1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kCosphi1phi2phi3p);
2239  Double_t sinP1nPsi1P1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kSinphi1phi2phi3p);
2240 
2241  cosP1nPsi1P1nPhi2 /= w2p;
2242  sinP1nPsi1P1nPhi2 /= w2p;
2243  cosP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
2244  sinP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
2245  cosP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
2246  sinP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
2247 
2248  Double_t fourPrime = w4pFourPrime / w4p;
2249  Double_t qc4Prime = fourPrime-2.*twoPrime*two;
2250  if (cumu4) cumu4->Fill(eta, cent, qc4Prime);
2251 
2252  if ((fFlags & kNUAcorr)) {
2253  qc4Prime += - cosP1nPsi*cosP1nPhi1M1nPhi2M1nPhi3
2254  + sinP1nPsi*sinP1nPhi1M1nPhi2M1nPhi3
2255  - cosP1nPhiA*cosP1nPsi1M1nPhi2M1nPhi3
2256  + sinP1nPhiA*sinP1nPsi1M1nPhi2M1nPhi3
2257  - 2.*cosP1nPhiA*cosP1nPsi1P1nPhi2M1nPhi3
2258  - 2.*sinP1nPhiA*sinP1nPsi1P1nPhi2M1nPhi3
2259  - cosP1nPsi1P1nPhi2*cosP1nPhi1P1nPhi2
2260  - sinP1nPsi1P1nPhi2*sinP1nPhi1P1nPhi2
2261  + 2.*cosP1nPhi1P1nPhi2*(cosP1nPsi*cosP1nPhiA-sinP1nPsi*sinP1nPhiA)
2262  + 2.*sinP1nPhi1P1nPhi2*(cosP1nPsi*sinP1nPhiA+sinP1nPsi*cosP1nPhiA)
2263  + 4.*two*(cosP1nPsi*cosP1nPhiA+sinP1nPsi*sinP1nPhiA)
2264  + 2.*cosP1nPsi1P1nPhi2*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
2265  + 4.*sinP1nPsi1P1nPhi2*cosP1nPhiA*sinP1nPhiA
2266  + 4.*twoPrime*(TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.))
2267  - 6.*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
2268  * (cosP1nPsi*cosP1nPhiA-sinP1nPsi*sinP1nPhiA)
2269  - 12.*cosP1nPhiA*sinP1nPhiA
2270  * (sinP1nPsi*cosP1nPhiA+cosP1nPsi*sinP1nPhiA);
2271  }
2272  // Double_t vnFourDiff = - qc4Prime / TMath::Power(-qc4, 0.75);
2273  if (!TMath::IsNaN(qc4Prime*mpqMult)) {
2274  quality->Fill((n-2)*qualityFactor+7, Int_t(cent));
2275  if (cumu4NUA) cumu4NUA->Fill(eta, cent, qc4Prime);
2276  }
2277  else
2278  quality->Fill((n-2)*qualityFactor+8, Int_t(cent));
2279  if (fDebug > 1)
2280  AliInfo(Form("%s: v_%d{4} = %1.3f for eta = %1.2f and centrality %3.1f",
2281  fType.Data(), n, qc4Prime, eta, cent));
2282  } // End of eta loop
2283  } // End of centrality loop
2284  } // End of moment
2285 
2286  return;
2287 }
2288 //_____________________________________________________________________
2290  TH2I* quality,
2291  TH1D* chist,
2292  TH2D* dNdetaRef,
2293  TH2D* dNdetaDiff) const
2294 {
2295  //
2296  // Calculates the 3 sub flow
2297  //
2298  // Parameters:
2299  // cumu2h: CumuHistos object with QC{2} cumulants
2300  // quality: Histogram for success rate of cumulants
2301  // chist: Centrality histogram
2302  // dNdetaDiff: dN/deta histogram for estimating multiplicity used for diff calculations
2303  //
2304 
2305  // For flow calculations
2306  TH3D* cumuRef = 0;
2307  TH3D* cumuDiff = 0;
2308  TH2D* cumu2r = 0;
2309  TH2D* cumu2rNUAold = 0;
2310  TH2D* cumu2a = 0;
2311  TH2D* cumu2aNUAold = 0;
2312  TH2D* cumu2b = 0;
2313  TH2D* cumu2bNUAold = 0;
2314  // Loop over cumulant histogram for final calculations
2315  for (Int_t n = 2; n <= fMaxMoment; n++) { // Moment loop begins
2316  cumu2r = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNoNUA);
2317  cumu2a = (TH2D*)cumu2h.Get('a', n, CumuHistos::kNoNUA);
2318  cumu2b = (TH2D*)cumu2h.Get('b', n, CumuHistos::kNoNUA);
2319  if ((fFlags & kNUAcorr)) {
2320  cumu2rNUAold = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNUAOld);
2321  cumu2aNUAold = (TH2D*)cumu2h.Get('a', n, CumuHistos::kNUAOld);
2322  cumu2bNUAold = (TH2D*)cumu2h.Get('b', n, CumuHistos::kNUAOld);
2323  }
2324  cumuRef = (TH3D*)fCumuHists.Get('r',n);
2325  cumuDiff = (TH3D*)fCumuHists.Get('d',n);
2326  for (Int_t cBin = 1; cBin <= cumuRef->GetNbinsY(); cBin++) { // Centrality loop begins
2327  Double_t cent = cumuRef->GetYaxis()->GetBinCenter(cBin);
2328  if (n == 2) chist->Fill(cent, cumuRef->GetBinContent(0, cBin, 0));
2329  if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), n, cent));
2330  // Here it starts!
2331  Int_t prevLim = 0;
2332  Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2333  Double_t cosP1nPhiA = 0;
2334  Double_t sinP1nPhiA = 0;
2335  Double_t cos2nPhiA = 0;
2336  Double_t sin2nPhiA = 0;
2337  Double_t cosP1nPhiB = 0;
2338  Double_t sinP1nPhiB = 0;
2339  Double_t cos2nPhiB = 0;
2340  Double_t sin2nPhiB = 0;
2341  Double_t multA = 0;
2342  Double_t multB = 0;
2343 
2344  for (Int_t etaBin = 1; etaBin <= cumuDiff->GetNbinsX(); etaBin++) { // Eta loop begins
2345  Double_t eta = cumuDiff->GetXaxis()->GetBinCenter(etaBin);
2346  // 2-particle reference flow
2347  Double_t w2Two = cumuRef->GetBinContent(etaBin, cBin, kW2Two);
2348  Double_t w2 = cumuRef->GetBinContent(etaBin, cBin, kW2);
2349  if (w2 == 0) continue;
2350 
2351  // Update NUA for new range!
2352  if (fEtaLims[prevLim] < eta) {
2353  GetLimits(prevLim, aLow, aHigh, bLow, bHigh);
2354  prevLim++;
2355  cosP1nPhiA = 0; sinP1nPhiA = 0; cos2nPhiA = 0; sin2nPhiA = 0; multA = 0;
2356  cosP1nPhiB = 0; sinP1nPhiB = 0; cos2nPhiB = 0; sin2nPhiB = 0; multB = 0;
2357  for (Int_t a = aLow; a <= aHigh; a++) {
2358  cosP1nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n));
2359  sinP1nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n));
2360  cos2nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(2*n));
2361  sin2nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(2*n));
2362  multA += fCumuNUARef->GetBinContent(a, cBin, 0);
2363  }
2364  for (Int_t b = bLow; b <= bHigh; b++) {
2365  cosP1nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n));
2366  sinP1nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n));
2367  cos2nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(2*n));
2368  sin2nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(2*n));
2369  multB += fCumuNUARef->GetBinContent(b, cBin, 0);
2370  }
2371  if (multA == 0 || multB == 0) {
2372  AliWarning(Form("Empty NUA values for 3Cor! (%s)", cumuRef->GetName()));
2373  continue;
2374  }
2375  cosP1nPhiA /= multA;
2376  sinP1nPhiA /= multA;
2377  cos2nPhiA /= multA;
2378  sin2nPhiA /= multA;
2379  cosP1nPhiB /= multB;
2380  sinP1nPhiB /= multB;
2381  cos2nPhiB /= multB;
2382  sin2nPhiB /= multB;
2383 
2384  dNdetaRef->Fill(eta, cent, multA+multB);
2385  }
2386  Double_t two = w2Two / w2;
2387 
2388  Double_t qc2 = two;
2389  if (qc2 >= 0) cumu2r->Fill(eta, cent, TMath::Sqrt(qc2));
2390 
2391  if ((fFlags & kNUAcorr)) {
2392  // Old nua
2393  qc2 -= cosP1nPhiA*cosP1nPhiB + sinP1nPhiA*sinP1nPhiB;
2394  // Extra NUA term from 2n cosines and sines
2395  qc2 /= (1+(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB));
2396  }
2397  if (qc2 <= 0) {
2398  if (fDebug > 0)
2399  AliInfo(Form("%s: QC_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f - skipping",
2400  fType.Data(), n, qc2, eta, cent));
2401  quality->Fill((n-2)*4+2, Int_t(cent));
2402  continue;
2403  }
2404  Double_t vnTwo = TMath::Sqrt(qc2);
2405  if (!TMath::IsNaN(vnTwo)) {
2406  quality->Fill((n-2)*4+1, Int_t(cent));
2407  if ((fFlags & kNUAcorr)) cumu2rNUAold->Fill(eta, cent, vnTwo);
2408  }
2409 
2410  // 2-particle differential flow
2411  Double_t w2pTwoPrimeA = cumuDiff->GetBinContent(etaBin, cBin, kW2Two);
2412  Double_t w2pA = cumuDiff->GetBinContent(etaBin, cBin, kW2);
2413  Double_t w2pTwoPrimeB = cumuDiff->GetBinContent(etaBin, cBin, kW4Four);
2414  Double_t w2pB = cumuDiff->GetBinContent(etaBin, cBin, kW4);
2415  if (w2pA == 0 || w2pB == 0) continue;
2416  Double_t cosP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(n));
2417  Double_t sinP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(n));
2418  Double_t cos2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(2*n));
2419  Double_t sin2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(2*n));
2420  Double_t mult = fCumuNUADiff->GetBinContent(etaBin, cBin, 0);
2421  if (mult == 0) continue;
2422  cosP1nPsi /= mult;
2423  sinP1nPsi /= mult;
2424  cos2nPsi /= mult;
2425  sin2nPsi /= mult;
2426  Double_t twoPrimeA = w2pTwoPrimeA / w2pA;
2427  Double_t twoPrimeB = w2pTwoPrimeB / w2pB;
2428  dNdetaDiff->Fill(eta, cent, mult);
2429 
2430  Double_t qc2PrimeA = twoPrimeA;
2431  Double_t qc2PrimeB = twoPrimeB;
2432  if (qc2PrimeA*qc2PrimeB >= 0) {
2433  cumu2a->Fill(eta, cent, qc2PrimeA);
2434  cumu2b->Fill(eta, cent, qc2PrimeB);
2435  }
2436  if ((fFlags & kNUAcorr)) {
2437  // Old nua
2438  qc2PrimeA -= cosP1nPsi*cosP1nPhiA + sinP1nPsi*sinP1nPhiA;
2439  qc2PrimeB -= cosP1nPsi*cosP1nPhiB + sinP1nPsi*sinP1nPhiB; // Is this OK?
2440  // Extra NUA term from 2n cosines and sines
2441  if (cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA != -1.) qc2PrimeA /= (1.+(cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA));
2442  if (cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB != -1.) qc2PrimeB /= (1.+(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
2443  }
2444  if (!TMath::IsNaN(qc2PrimeA) && !TMath::IsNaN(qc2PrimeB) && qc2 != 0) {
2445  if (qc2PrimeA*qc2PrimeB >= 0) {
2446  quality->Fill((n-2)*4+3, Int_t(cent));
2447  if ((fFlags & kNUAcorr)) cumu2aNUAold->Fill(eta, cent, qc2PrimeA);
2448  if ((fFlags & kNUAcorr)) cumu2bNUAold->Fill(eta, cent, qc2PrimeB);
2449  }
2450  }
2451  else
2452  quality->Fill((n-2)*4+4, Int_t(cent));
2453  if (fDebug > 1)
2454  AliInfo(Form("%s: QC'a_%d{2} = %1.3f, QC'b_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f",
2455  fType.Data(), n, qc2PrimeA, n, qc2PrimeB, eta, cent));
2456  } // End of eta loop
2457  } // End of centrality loop
2458  } // End of moment
2459 
2460  return;
2461 }
2462 //_____________________________________________________________________
2464 {
2465  //
2466  // Function to solve the coupled flow equations
2467  // We solve it by using matrix calculations:
2468  // A*v_n = V => v_n = A^-1*V
2469  // First we set up a TMatrixD container to make ROOT
2470  // do the inversions in an efficient way, we multiply the current <<2>> estimates.
2471  // Then we fill new TH2D's if the corrected <<2>>'s (cumuNUA object).
2472  //
2473  // Parameters:
2474  // cumu: CumuHistos object - uncorrected
2475  // type: Reference ('r') or differential ('d') or ('a' or 'b') for 3 correlator
2476  //
2477 
2478  // We start by declaring Matrix and vector objects, as their constructors are quite heavy
2479  TMatrixD mNUA(fMaxMoment-1, fMaxMoment-1);
2480  TVectorD vQC2(fMaxMoment-1);
2481 
2482  for (Int_t cBin = 1; cBin <= cumu.Get(type, 2, CumuHistos::kNUAOld)->GetNbinsY(); cBin++) { // cent loop
2483  Double_t cent = cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinCenter(cBin);
2484  for (Int_t eBin = 1; eBin <= cumu.Get(type, 2, CumuHistos::kNUAOld)->GetNbinsX(); eBin++) { // eta loop
2485  Double_t eta = cumu.Get(type, 2, CumuHistos::kNUAOld)->GetXaxis()->GetBinCenter(eBin);
2486  mNUA.Zero(); // reset matrix
2487  vQC2.Zero(); // reset vector
2488  for (Int_t n = 0; n < fMaxMoment-1; n++) { // moment loop
2489  vQC2(n) = static_cast<TH2D*>(cumu.Get(type, n+2, CumuHistos::kNUAOld))->GetBinContent(eBin, cBin);
2490  if (type == 'r' || type == 'R') vQC2(n) *= vQC2(n); // going back to <<2>>
2491  for (Int_t m = 0; m < fMaxMoment-1; m++) { // cross moment
2492  mNUA(n,m) = CalculateNUAMatrixElement(n, m, type, eBin, cBin);
2493  } // End of cross moment loop
2494  } // End of moment loop
2495  // Invert matrix
2496  Double_t det = 0;
2497  mNUA.Invert(&det);
2498  // If determinant is non-zero we go with corrected results
2499  if (det != 0 ) vQC2 = mNUA*vQC2;
2500  else AliWarning(Form("Determinant == 0 - cent: %d-%d, eta: %f, type: '%c', data: %s, vtx: %d-%d%s",
2501  Int_t(cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinLowEdge(cBin)),
2502  Int_t(cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinUpEdge(cBin)),
2503  cumu.Get(type, 2, CumuHistos::kNUAOld)->GetXaxis()->GetBinCenter(eBin),
2504  type, fType.Data(), fVzMin, fVzMax,
2505  GetQCType(fFlags, kTRUE)));
2506  // Go back to v_n for ref. keep <<2'>> for diff. flow).
2507  for (Int_t n = 0; n < fMaxMoment-1; n++) {
2508  Double_t vnTwo = 0;
2509  if (type == 'r' || type == 'R')
2510  vnTwo = (vQC2(n) > 0. ? TMath::Sqrt(vQC2(n)) : 0.);
2511  else {
2512  // is really more <<2'>> in this case
2513  vnTwo = vQC2(n);
2514  }
2515  // Fill in corrected v_n
2516  if (vnTwo != 0) static_cast<TH2D*>(cumu.Get(type, n+2, CumuHistos::kNUA))->Fill(eta, cent, vnTwo);
2517  } // End of moment loop
2518  } // End of eta loop
2519  } // End of centrality loop
2520  return;
2521 }
2522 //_____________________________________________________________________
2524 {
2525  //
2526  // Calculates the (n,m)-th element in the NUA matrix: 1 if n == m, otherwise:
2527  // <<cos[(n-m)phi1]>>*<<cos[(n-m)phi2]>> + <<sin[(n-m)phi1]>>*<<sin[(n-m)phi2]>>
2528  // NUA(n,m) = -----------------------------------------------------------------------------
2529  // 1 + <<cos[2nphi1]>>*<<cos[2nphi2]>> + <<sin[2nphi1]>>*<<sin[2nphi2]>>
2530  //
2531  // <<cos[(n+m)phi1]>>*<<cos[(n+m)phi2]>> + <<sin[(n+m)phi1]>>*<<sin[(n+m)phi2]>>
2532  // + -----------------------------------------------------------------------------
2533  // 1 + <<cos[2nphi1]>>*<<cos[2nphi2]>> + <<sin[2nphi1]>>*<<sin[2nphi2]>>
2534  //
2535  // Parameters:
2536  // n: row
2537  // m: coumn
2538  // type: Reference ('r') or differential ('d') or ('a' or 'b')
2539  // binA: eta bin of phi1
2540  // cBin: centrality bin
2541  //
2542  // Return: NUA(n,m)
2543  //
2544  if (n == m) return 1.;
2545  n += 2;
2546  m += 2;
2547 
2548  Double_t cosnMmPhi1 = 0, cosnMmPhi2 = 0, sinnMmPhi1 = 0, sinnMmPhi2 = 0;
2549  Double_t cosnPmPhi1 = 0, cosnPmPhi2 = 0, sinnPmPhi1 = 0, sinnPmPhi2 = 0;
2550  Double_t cos2nPhi1 = 0, cos2nPhi2 = 0, sin2nPhi1 = 0, sin2nPhi2 = 0;
2551 
2552  // reference flow
2553  if (type == 'r' || type == 'R') {
2554  if ((fFlags & k3Cor)) {
2555  Double_t eta = fCumuNUARef->GetXaxis()->GetBinCenter(binA);
2556  Int_t i = 0;
2557  while (fEtaLims[i] < eta) i++;
2558  Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2559  GetLimits(i-1, aLow, aHigh, bLow, bHigh);
2560  Double_t multA = 0, multB = 0;
2561  for (Int_t a = aLow; a <= aHigh; a++) {
2562  cosnMmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n-m));
2563  sinnMmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n-m));
2564  cosnPmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n+m));
2565  sinnPmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n+m));
2566  cos2nPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(2*n));
2567  sin2nPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(2*n));
2568  multA += fCumuNUARef->GetBinContent(a, cBin, 0);
2569  }
2570  for (Int_t b = bLow; b <= bHigh; b++) {
2571  cosnMmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n-m));
2572  sinnMmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n-m));
2573  cosnPmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n+m));
2574  sinnPmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n+m));
2575  cos2nPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(2*n));
2576  sin2nPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(2*n));
2577  multB += fCumuNUARef->GetBinContent(b, cBin, 0);
2578  }
2579  if (multA == 0 || multB == 0) {
2580  if (fDebug > 0) AliWarning("multA or multB == 0 in matrix elements, aborting NUA");
2581  return 0.;
2582  }
2583  cosnMmPhi1 /= multA;
2584  sinnMmPhi1 /= multA;
2585  cosnPmPhi1 /= multA;
2586  sinnPmPhi1 /= multA;
2587  cos2nPhi1 /= multA;
2588  sin2nPhi1 /= multA;
2589  cosnMmPhi2 /= multB;
2590  sinnMmPhi2 /= multB;
2591  cosnPmPhi2 /= multB;
2592  sinnPmPhi2 /= multB;
2593  cos2nPhi2 /= multB;
2594  sin2nPhi2 /= multB;
2595  } else {
2596  Int_t binB = fCumuNUARef->GetXaxis()->FindBin(-fCumuNUARef->GetXaxis()->GetBinCenter(binA));
2597  cosnMmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2598  sinnMmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2599  cosnPmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2600  sinnPmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2601  cos2nPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2602  sin2nPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2603  cosnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n-m));
2604  sinnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n-m));
2605  cosnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n+m));
2606  sinnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n+m));
2607  cos2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(2*n));
2608  sin2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(2*n));
2609  }
2610  } // differential flow
2611  else if (type == 'd' || type == 'D') {
2612  Int_t binB = fCumuNUARef->GetXaxis()->FindBin(-fCumuNUADiff->GetXaxis()->GetBinCenter(binA));
2613  cosnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2614  sinnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2615  cosnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2616  sinnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2617  cos2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2618  sin2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2619  cosnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n-m));
2620  sinnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n-m));
2621  cosnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n+m));
2622  sinnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n+m));
2623  cos2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(2*n));
2624  sin2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(2*n));
2625  } // 3 correlator part a or b
2626  else if (type == 'a' || type == 'A' || type == 'b' || type == 'B') {
2627  Double_t mult1 = 0, mult2 = 0;
2628  // POIs
2629  cosnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2630  sinnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2631  cosnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2632  sinnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2633  cos2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2634  sin2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2635  mult1 = fCumuNUADiff->GetBinContent(binA, cBin, 0);
2636  // RPs
2637  Double_t eta = fCumuNUADiff->GetXaxis()->GetBinCenter(binA);
2638  Int_t i = 0;
2639  while (fEtaLims[i] < eta) i++;
2640  Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2641  GetLimits(i-1, aLow, aHigh, bLow, bHigh);
2642  Int_t lLow = ((type == 'a' || type == 'A') ? aLow : bLow);
2643  Int_t lHigh = ((type == 'a' || type == 'A') ? aHigh : bHigh);
2644  for (Int_t l = lLow; l <= lHigh; l++) {
2645  cosnMmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(n-m));
2646  sinnMmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(n-m));
2647  cosnPmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(n+m));
2648  sinnPmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(n+m));
2649  cos2nPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(2*n));
2650  sin2nPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(2*n));
2651  mult2 += fCumuNUARef->GetBinContent(l, cBin, 0);
2652  }
2653  if (mult1 == 0 || mult2 == 0) {
2654  if (fDebug > 0) AliWarning("mult1 or mult2 == 0 in matrix elements, aborting NUA");
2655  return 0.;
2656  }
2657  cosnMmPhi1 /= mult1;
2658  sinnMmPhi1 /= mult1;
2659  cosnPmPhi1 /= mult1;
2660  sinnPmPhi1 /= mult1;
2661  cos2nPhi1 /= mult1;
2662  sin2nPhi1 /= mult1;
2663  cosnMmPhi2 /= mult2;
2664  sinnMmPhi2 /= mult2;
2665  cosnPmPhi2 /= mult2;
2666  sinnPmPhi2 /= mult2;
2667  cos2nPhi2 /= mult2;
2668  sin2nPhi2 /= mult2;
2669  }
2670 
2671  // Actual calculation
2672  Double_t e = cosnMmPhi1*cosnMmPhi2 + sinnMmPhi1*sinnMmPhi2 + cosnPmPhi1*cosnPmPhi2 + sinnPmPhi1*sinnPmPhi2;
2673  Double_t den = 1 + cos2nPhi1*cos2nPhi2 + sin2nPhi1*sin2nPhi2;
2674  if (den != 0) e /= den;
2675  else return 0.;
2676 
2677  return e;
2678 }
2679 //_____________________________________________________________________
2681 {
2682  //
2683  // Add up vertex bins with flow results
2684  //
2685  // Parameters:
2686  // cumu: CumuHistos object with vtxbin results
2687  // list: Outout list with added results
2688  // nNUA: # of NUA histograms to loop over
2689  //
2690  TH2D* vtxHist = 0;
2691  TProfile2D* avgProf = 0;
2692  TString name;
2693  Int_t nT = ((fFlags & k3Cor) ? 3 : 2);
2694  Char_t ct = '\0';
2695  for (UInt_t nua = 0; nua <= nNUA; nua++) { // NUA loop
2696  for (Int_t n = 2; n <= fMaxMoment; n++) { // moment loop
2697  for (Int_t t = 0; t < nT; t++) { // type loop (r/d/a/b)
2698  // Find type
2699  switch (t) {
2700  case 0: ct = 'r'; break;
2701  case 1: ct = ((fFlags & k3Cor) ? 'a' : 'd'); break;
2702  case 2: ct = 'b'; break;
2703  default: ct = '\0'; break;
2704  }
2705  vtxHist = static_cast<TH2D*>(cumu.Get(ct,n,nua));
2706  if (!vtxHist) {
2707  AliWarning("VertexBin::AddVertexBins: vtxHist not found!");
2708  continue;
2709  }
2710  name = vtxHist->GetName();
2711  // Strip name of vtx info
2712  Ssiz_t l = name.Last('x')-3;
2713  name.Resize(l);
2714  avgProf = (TProfile2D*)list->FindObject(name.Data());
2715  // if no output profile yet, make one
2716  if (!avgProf) {
2717  avgProf = new TProfile2D(name.Data(), name.Data(),
2718  vtxHist->GetNbinsX(), vtxHist->GetXaxis()->GetXmin(), vtxHist->GetXaxis()->GetXmax(),
2719  vtxHist->GetNbinsY(), vtxHist->GetYaxis()->GetXmin(), vtxHist->GetYaxis()->GetXmax());
2720  if (vtxHist->GetXaxis()->IsVariableBinSize())
2721  avgProf->GetXaxis()->Set(vtxHist->GetNbinsX(), vtxHist->GetXaxis()->GetXbins()->GetArray());
2722  if (vtxHist->GetYaxis()->IsVariableBinSize())
2723  avgProf->GetYaxis()->Set(vtxHist->GetNbinsY(), vtxHist->GetYaxis()->GetXbins()->GetArray());
2724  list->Add(avgProf);
2725  }
2726  // Fill in, cannot be done with Add function.
2727  for (Int_t e = 1; e <= vtxHist->GetNbinsX(); e++) { // eta loop
2728  Double_t eta = vtxHist->GetXaxis()->GetBinCenter(e);
2729  for (Int_t c = 1; c <= vtxHist->GetNbinsY(); c++) { // cent loop
2730  Double_t cent = vtxHist->GetYaxis()->GetBinCenter(c);
2731  Double_t cont = vtxHist->GetBinContent(e, c);
2732  if (cont == 0) continue;
2733  avgProf->Fill(eta, cent, cont);
2734  } // End of cent loop
2735  } // End of eta loop
2736  } // End of type loop
2737  } // End of moment loop
2738  } // End of nua loop
2739 }
2740 //_____________________________________________________________________
2742 {
2743  //
2744  // Get the bin number of <<cos(nphi)>>
2745  //
2746  // Parameters:
2747  // n: moment
2748  //
2749  // Return: bin number
2750  //
2751  Int_t bin = 0;
2752  n = TMath::Abs(n);
2753 
2754  if (n == 0) bin = fMaxMoment*4-1;
2755  else bin = n*2-1;
2756 
2757  return bin;
2758 }
2759 //_____________________________________________________________________
2761 {
2762  //
2763  // Get the bin number of <<sin(nphi)>>
2764  //
2765  // Parameters:
2766  // n: moment
2767  //
2768  // Return: bin number
2769  //
2770  Int_t bin = 0;
2771  n = TMath::Abs(n);
2772 
2773  if (n == 0) bin = fMaxMoment*4;
2774  else bin = n*2;
2775 
2776  return bin;
2777 }
2778 //_____________________________________________________________________
2780 {
2781  //
2782  // Setup NUA labels on axis
2783  //
2784  // Parameters:
2785  // a: Axis to set up
2786  //
2787  if (a->GetNbins() != GetBinNumberSin()) AliFatal("SetupNUALabels: Wrong number of bins on axis");
2788 
2789  Int_t i = 1, j= 1;
2790  while (i <= a->GetNbins()) {
2791  a->SetBinLabel(i++, Form("<<cos(%d#varphi)>>", j));
2792  a->SetBinLabel(i++, Form("<<sin(%d#varphi)>>", j++));
2793  }
2794 
2795  return;
2796 }
2797 //_____________________________________________________________________
2799 {
2800  //
2801  // Add a histogram for checking the analysis quality
2802  //
2803  // Parameters:
2804  // const char*: name of data type
2805  //
2806  // Return: histogram for tracking successful calculations
2807  //
2808  Int_t nBins = ((fFlags & kStdQC) ? 8 : 4);
2809  TH2I* quality = new TH2I(name, name, (fMaxMoment-1)*nBins, 1, (fMaxMoment-1)*nBins+1, fCumuNUARef->GetNbinsY(),
2810  fCumuNUARef->GetYaxis()->GetXmin(), fCumuNUARef->GetYaxis()->GetXmax());
2811  quality->GetYaxis()->Set(fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXbins()->GetArray());
2812  for (Int_t i = 2, j = 1; i <= fMaxMoment; i++) {
2813  quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{2} > 0", i));
2814  quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{2} <= 0", i));
2815  quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{2} > 0", i));
2816  quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{2} <= 0", i));
2817  if ((fFlags & kStdQC)) {
2818  quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{4} < 0", i));
2819  quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{4} >= 0", i));
2820  quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{4} < 0", i));
2821  quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{4} >= 0", i));
2822  }
2823  }
2824 
2825  return quality;
2826 }
2827 //_____________________________________________________________________
2829 {
2830  //
2831  // Setup a TH2D for the output
2832  //
2833  // Parameters:
2834  // qc # of particle correlations
2835  // n flow moment
2836  // ref Type: r/d/a/b
2837  // nua For nua corrected hists?
2838  //
2839  // Return: 2D hist for results
2840  //
2841  Bool_t ref = ((ctype[0] == 'R') || (ctype[0] == 'r'));
2842  TAxis* xAxis = (ref ? fCumuNUARef->GetXaxis() : fCumuNUADiff->GetXaxis());
2843  TAxis* yAxis = (ref ? fCumuNUARef->GetYaxis() : fCumuNUADiff->GetYaxis());
2844  TString ntype = "";
2845  switch (nua) {
2846  case CumuHistos::kNoNUA: ntype = "Un"; break;
2847  case CumuHistos::kNUAOld: ntype = "NUAOld"; break;
2848  case CumuHistos::kNUA: ntype = "NUA"; break;
2849  default: break;
2850  }
2851  TH2D* h = new TH2D(Form("%sQC%d_v%d_%s_%sCorr%s_vtx_%d_%d",
2852  fType.Data(), qc, n, ctype, ntype.Data(),
2854  Form("%sQC%d_v%d_%s_%sCorr%s_vtx_%d_%d",
2855  fType.Data(), qc, n, ctype, ntype.Data(),
2857  xAxis->GetNbins(), xAxis->GetXmin(), xAxis->GetXmax(),
2858  yAxis->GetNbins(), yAxis->GetXmin(), yAxis->GetXmax());
2859  if (xAxis->IsVariableBinSize()) h->GetXaxis()->Set(xAxis->GetNbins(), xAxis->GetXbins()->GetArray());
2860  h->GetYaxis()->Set(yAxis->GetNbins(), yAxis->GetXbins()->GetArray());
2861 
2862  return h;
2863 }
2864 //_____________________________________________________________________
2866 {
2867  //
2868  // Print the setup of the flow task
2869  //
2870  Printf("=======================================================");
2871  Printf("%s::Print", this->IsA()->GetName());
2872  Printf("Forward detector: :\t%s", ((fFlowFlags & kFMD) ? "FMD" : "VZERO"));
2873  Printf("Number of bins in vertex axis :\t%d", fVtxAxis->GetNbins());
2874  Printf("Range of vertex axis :\t[%3.1f,%3.1f]",
2875  fVtxAxis->GetXmin(), fVtxAxis->GetXmax());
2876  Printf("Number of bins in centrality axis :\t%d", fCentAxis->GetNbins());
2877  printf("Centrality binning :\t");
2878  for (Int_t cBin = 1; cBin <= fCentAxis->GetNbins(); cBin++) {
2879  printf("%02d-%02d%% ", Int_t(fCentAxis->GetBinLowEdge(cBin)), Int_t(fCentAxis->GetBinUpEdge(cBin)));
2880  if (cBin == fCentAxis->GetNbins()) printf("\n");
2881  else if (cBin % 4 == 0) printf("\n\t\t\t\t\t");
2882  }
2883  printf("Doing flow analysis for :\t");
2884  for (Int_t n = 2; n <= fMaxMoment; n++) printf("v%d ", n);
2885  printf("\n");
2886  TString type = "Standard QC{2} and QC{4} calculations.";
2887  if ((fFlowFlags & kEtaGap)) type = "QC{2} with a rapidity gap.";
2888  if ((fFlowFlags & k3Cor)) type = "QC{2} with 3 correlators.";
2889  if ((fFlowFlags & kTPC) == kTPC) type.ReplaceAll(".", " with TPC tracks for reference.");
2890  if ((fFlowFlags & kHybrid) == kHybrid) type.ReplaceAll(".", " with hybrid tracks for reference.");
2891  Printf("QC calculation type :\t%s", type.Data());
2892  Printf("Symmetrize ref. flow wrt. eta = 0 :\t%s", ((fFlowFlags & kSymEta) ? "true" : "false"));
2893  Printf("Apply NUA correction terms :\t%s", ((fFlowFlags & kNUAcorr) ? "true" : "false"));
2894  Printf("Satellite vertex flag :\t%s", ((fFlowFlags & kSatVtx) ? "true" : "false"));
2895  Printf("FMD sigma cut: :\t%f", fFMDCut);
2896  Printf("SPD sigma cut: :\t%f", fSPDCut);
2897  if ((fFlowFlags & kEtaGap) || (fFlowFlags & kTracks))
2898  Printf("Eta gap: :\t%f", fEtaGap);
2899  Printf("=======================================================");
2900 }
2901 //_____________________________________________________________________
2903 {
2904  //
2905  // Get the type of the QC calculations
2906  // Used for naming of objects in the VertexBin class,
2907  // important to avoid memory problems when running multiple
2908  // initializations of the class with different setups
2909  //
2910  // Parameters:
2911  // flags: Flow flags for type determination
2912  // prependUS: Prepend an underscore (_) to the name
2913  //
2914  // Return: QC calculation type
2915  //
2916  static TString type = "";
2917  if ((flags & kStdQC)) type = "StdQC";
2918  else if ((flags & kEtaGap)) type = "EtaGap";
2919  else if ((flags & k3Cor)) type = "3Cor";
2920  else type = "UNKNOWN";
2921  if (prependUS) type.Prepend("_");
2922  if ((flags & kTPC) == kTPC) type.Append("TPCTr");
2923  if ((flags & kHybrid) == kHybrid) type.Append("HybTr");
2924 
2925  return type.Data();
2926 }
2927 //_____________________________________________________________________
2929 {
2930  //
2931  // Get histogram/profile for type t and moment n
2932  //
2933  // Parameters:
2934  // t: type = 'r'/'d'
2935  // n: moment
2936  // nua: NUA type
2937  //
2938  n = GetPos(n, nua);
2939  if (n < 0) AliFatal(Form("CumuHistos out of range: (%c,%d)", t, n));
2940 
2941  TH1* h = 0;
2942  Int_t pos = -1;
2943  if (t == 'r' || t == 'R') pos = n;
2944  else if (t == 'd' || t == 'D') pos = n;
2945  else if (t == 'a' || t == 'A') pos = 2*n;
2946  else if (t == 'b' || t == 'B') pos = 2*n+1;
2947  else AliFatal(Form("CumuHistos wrong type: %c", t));
2948 
2949  if (t == 'r' || t == 'R') {
2950  if (pos < fRefHists->GetEntries()) {
2951  h = (TH1*)fRefHists->At(pos);
2952  }
2953  } else if (pos < fDiffHists->GetEntries()) {
2954  h = (TH1*)fDiffHists->At(pos);
2955  }
2956  if (!h) AliFatal(Form("No hist found in list %c at pos %d", t, pos));
2957 
2958  return h;
2959 }
2960 //_____________________________________________________________________
2962 {
2963  //
2964  // Connect an output list with this object
2965  //
2966  // Parameters:
2967  // name: base name
2968  // l: list to keep outputs in
2969  //
2970  TString ref = name;
2971  ref.ReplaceAll("Cumu","CumuRef");
2972  fRefHists = (TList*)l->FindObject(ref.Data());
2973  if (!fRefHists) {
2974  fRefHists = new TList();
2975  fRefHists->SetName(ref.Data());
2976  l->Add(fRefHists);
2977  } else {
2978  // Check that the list correspond to fMaxMoments
2979  if (fRefHists->GetEntries() != (fMaxMoment-1.)*(fNUA+1))
2980  AliError("CumuHistos::ConnectList Wrong number of hists in ref list,"
2981  "you are doing something wrong!");
2982  }
2983  TString diff = name;
2984  diff.ReplaceAll("Cumu","CumuDiff");
2985  fDiffHists = (TList*)l->FindObject(diff.Data());
2986  if (!fDiffHists) {
2987  fDiffHists = new TList();
2988  fDiffHists->SetName(diff.Data());
2989  l->Add(fDiffHists);
2990  } else {
2991  // Check that the list correspond to fMaxMoment
2992  if ((fDiffHists->GetEntries() != (fMaxMoment-1.)*(fNUA+1)) &&
2993  (fDiffHists->GetEntries() != 2*(fMaxMoment-1.)*(fNUA+1)))
2994  AliError(Form("CumuHistos::ConnectList Wrong number of hists in diff list,"
2995  "you are doing something wrong! (%s)",name.Data()));
2996  }
2997 
2998 }
2999 //_____________________________________________________________________
3001 {
3002  //
3003  // Add a histogram to this objects lists
3004  //
3005  // Parameters:
3006  // h: histogram/profile to add
3007  //
3008  TString name = h->GetName();
3009  if (name.Contains("Ref")) {
3010  if (fRefHists) fRefHists->Add(h);
3011  else AliFatal("CumuHistos::Add() - fRefHists does not exist");
3012  // Check that the list correspond to fMaxMoments
3013  if (fRefHists->GetEntries() > (fNUA+1)*(fMaxMoment-1.))
3014  AliError("CumuHistos::Add wrong number of hists in ref list, "
3015  "you are doing something wrong!");
3016  }
3017  else if (name.Contains("Diff")) {
3018  if (fDiffHists) fDiffHists->Add(h);
3019  else AliFatal("CumuHistos::Add() - fDiffHists does not exist");
3020  // Check that the list correspond to fMaxMoment
3021  if (fDiffHists->GetEntries() > 2*(fNUA+1)*(fMaxMoment-1.))
3022  AliError("CumuHistos::Add wrong number of hists in diff list, "
3023  "you are doing something wrong!");
3024  }
3025  return;
3026 }
3027 //_____________________________________________________________________
3029 {
3030  //
3031  // Get position in list of moment n objects
3032  // To take care of different numbering schemes
3033  //
3034  // Parameters:
3035  // n: moment
3036  // nua: # of nua corrections applied
3037  //
3038  // Return: position #
3039  //
3040  if (n > fMaxMoment) return -1;
3041  else return (n-2)+nua*(fMaxMoment-1);
3042 }
3043 //_____________________________________________________________________
3044 //
3045 //
3046 // EOF
virtual void UserExec(Option_t *option)
Bool_t HasCentrality() const
double Double_t
Definition: External.C:58
Double_t CalculateNUAMatrixElement(Int_t n, Int_t m, Char_t type, Int_t binA, Int_t cBin) const
void Calculate3CorFlow(CumuHistos &cumu2h, TH2I *quality, TH1D *chist, TH2D *dNdetaRef, TH2D *dNdetaDiff) const
AliVVZERO * GetVZERO() const
void AddOutput(TList *list, TAxis *centAxis)
void EndVtxBinList(const TList &list) const
Definition: External.C:236
void CalculateDifferentialFlow(CumuHistos &cumu2h, CumuHistos &cumu4h, TH2I *quality, TH2D *dNdetaDiff) const
Bool_t FillTracks(VertexBin *bin, UShort_t mode) const
void AddVertexBins(CumuHistos &cumu, TList *list, UInt_t nNUA) const
char Char_t
Definition: External.C:18
virtual Bool_t GetVertex(const AliAODForwardMult *aodfm)
TH1 * Get(Char_t t, Int_t n, UInt_t nua=0) const
TCanvas * c
Definition: TestFitELoss.C:172
static AliAODEvent * GetAODEvent(AliAnalysisTaskSE *task)
Float_t GetIpZ() const
virtual void Terminate(Option_t *option)
void FillVtxBinList3Cor(const TList &list, TH2D &hcent, TH2D &hfwd, Int_t vtx, UShort_t flags=0x0)
Per-event per bin.
int Int_t
Definition: External.C:63
Definition: External.C:204
unsigned int UInt_t
Definition: External.C:33
static Bool_t IsTriggerBits(UInt_t bits, UInt_t trg)
float Float_t
Definition: External.C:68
void ConnectList(TString name, TList *l)
Various utilities used in PWGLF/FORWARD.
void CalculateReferenceFlow(CumuHistos &cumu2h, CumuHistos &cumu4h, TH2I *quality, TH1D *chist, TH2D *dNdetaRef) const
void CumulantsTerminate(TList *inlist, TList *outlist)
Definition: External.C:252
Bool_t FillTracks(TObjArray *trList, AliESDEvent *esd, AliAnalysisFilter *trFilter, UShort_t mode)
VertexBin & operator=(const VertexBin &v)
void FillVtxBinList(const TList &list, TH2D &h1, Int_t vtx, UShort_t flags=0x0) const
Definition: External.C:228
Definition: External.C:212
AliAnalysisFilter * fTrackCuts
Bool_t HasIpZ() const
void GetLimits(Int_t bin, Int_t &aLow, Int_t &aHigh, Int_t &bLow, Int_t &bHigh) const
TH2D & CombineHists(TH2D &hcent, TH2D &hfwd)
Double_t nEvents
plot quality messages
Int_t GetPos(Int_t n, UInt_t nua) const
Int_t mode
Definition: anaM.C:41
virtual Bool_t GetCentrality(const AliAODForwardMult *aodfm)
AliForwardFlowTaskQC & operator=(const AliForwardFlowTaskQC &)
virtual Bool_t CheckEvent(const AliAODForwardMult *aodfm)
virtual Bool_t CheckTrigger(const AliAODForwardMult *aodfm) const
const TH2D & GetHistogram() const
TH2D * MakeOutputHist(Int_t qc, Int_t n, const Char_t *ctype, UInt_t nua) const
static UShort_t CheckForAOD()
void SetFlowFlags(UShort_t flags)
void MakeCentralityHists(TList *list) const
Bool_t FillHists(TH2D &dNdetadphi, Double_t cent, UShort_t mode)
TH2I * MakeQualityHist(const Char_t *name) const
unsigned short UShort_t
Definition: External.C:28
Float_t GetCentrality() const
const char Option_t
Definition: External.C:48
void SolveCoupledFlowEquations(CumuHistos &cumu, Char_t type) const
virtual void UserCreateOutputObjects()
static const Char_t * GetQCType(UShort_t flags, Bool_t prependUS=kTRUE)
bool Bool_t
Definition: External.C:53
void SetCentralityAxis(TAxis *axis)
void FillVZEROHist(AliVVZERO *vzero)
Definition: External.C:196
const TH2D & GetHistogram() const
void FillVtxBinListEtaGap(const TList &list, TH2D &href, TH2D &hdiff, Int_t vtx, UShort_t flags=0x0) const