AliPhysics  cdeda5a (cdeda5a)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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 //_____________________________________________________________________
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;
1090  fCumuNUARef = o.fCumuNUARef;
1091  fCumuNUADiff = o.fCumuNUADiff;
1092  fdNdedpRefAcc = o.fdNdedpRefAcc;
1093  fdNdedpDiffAcc = o.fdNdedpDiffAcc;
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  TH2D* dNdetaRef, TH2D* dNdetaDiff) const
2291 {
2292  //
2293  // Calculates the 3 sub flow
2294  //
2295  // Parameters:
2296  // cumu2h: CumuHistos object with QC{2} cumulants
2297  // quality: Histogram for success rate of cumulants
2298  // chist: Centrality histogram
2299  // dNdetaDiff: dN/deta histogram for estimating multiplicity used for diff calculations
2300  //
2301 
2302  // For flow calculations
2303  TH3D* cumuRef = 0;
2304  TH3D* cumuDiff = 0;
2305  TH2D* cumu2r = 0;
2306  TH2D* cumu2rNUAold = 0;
2307  TH2D* cumu2a = 0;
2308  TH2D* cumu2aNUAold = 0;
2309  TH2D* cumu2b = 0;
2310  TH2D* cumu2bNUAold = 0;
2311  // Loop over cumulant histogram for final calculations
2312  for (Int_t n = 2; n <= fMaxMoment; n++) { // Moment loop begins
2313  cumu2r = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNoNUA);
2314  cumu2a = (TH2D*)cumu2h.Get('a', n, CumuHistos::kNoNUA);
2315  cumu2b = (TH2D*)cumu2h.Get('b', n, CumuHistos::kNoNUA);
2316  if ((fFlags & kNUAcorr)) {
2317  cumu2rNUAold = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNUAOld);
2318  cumu2aNUAold = (TH2D*)cumu2h.Get('a', n, CumuHistos::kNUAOld);
2319  cumu2bNUAold = (TH2D*)cumu2h.Get('b', n, CumuHistos::kNUAOld);
2320  }
2321  cumuRef = (TH3D*)fCumuHists.Get('r',n);
2322  cumuDiff = (TH3D*)fCumuHists.Get('d',n);
2323  for (Int_t cBin = 1; cBin <= cumuRef->GetNbinsY(); cBin++) { // Centrality loop begins
2324  Double_t cent = cumuRef->GetYaxis()->GetBinCenter(cBin);
2325  if (n == 2) chist->Fill(cent, cumuRef->GetBinContent(0, cBin, 0));
2326  if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), n, cent));
2327  // Here it starts!
2328  Int_t prevLim = 0;
2329  Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2330  Double_t cosP1nPhiA = 0;
2331  Double_t sinP1nPhiA = 0;
2332  Double_t cos2nPhiA = 0;
2333  Double_t sin2nPhiA = 0;
2334  Double_t cosP1nPhiB = 0;
2335  Double_t sinP1nPhiB = 0;
2336  Double_t cos2nPhiB = 0;
2337  Double_t sin2nPhiB = 0;
2338  Double_t multA = 0;
2339  Double_t multB = 0;
2340 
2341  for (Int_t etaBin = 1; etaBin <= cumuDiff->GetNbinsX(); etaBin++) { // Eta loop begins
2342  Double_t eta = cumuDiff->GetXaxis()->GetBinCenter(etaBin);
2343  // 2-particle reference flow
2344  Double_t w2Two = cumuRef->GetBinContent(etaBin, cBin, kW2Two);
2345  Double_t w2 = cumuRef->GetBinContent(etaBin, cBin, kW2);
2346  if (w2 == 0) continue;
2347 
2348  // Update NUA for new range!
2349  if (fEtaLims[prevLim] < eta) {
2350  GetLimits(prevLim, aLow, aHigh, bLow, bHigh);
2351  prevLim++;
2352  cosP1nPhiA = 0; sinP1nPhiA = 0; cos2nPhiA = 0; sin2nPhiA = 0; multA = 0;
2353  cosP1nPhiB = 0; sinP1nPhiB = 0; cos2nPhiB = 0; sin2nPhiB = 0; multB = 0;
2354  for (Int_t a = aLow; a <= aHigh; a++) {
2355  cosP1nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n));
2356  sinP1nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n));
2357  cos2nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(2*n));
2358  sin2nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(2*n));
2359  multA += fCumuNUARef->GetBinContent(a, cBin, 0);
2360  }
2361  for (Int_t b = bLow; b <= bHigh; b++) {
2362  cosP1nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n));
2363  sinP1nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n));
2364  cos2nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(2*n));
2365  sin2nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(2*n));
2366  multB += fCumuNUARef->GetBinContent(b, cBin, 0);
2367  }
2368  if (multA == 0 || multB == 0) {
2369  AliWarning(Form("Empty NUA values for 3Cor! (%s)", cumuRef->GetName()));
2370  continue;
2371  }
2372  cosP1nPhiA /= multA;
2373  sinP1nPhiA /= multA;
2374  cos2nPhiA /= multA;
2375  sin2nPhiA /= multA;
2376  cosP1nPhiB /= multB;
2377  sinP1nPhiB /= multB;
2378  cos2nPhiB /= multB;
2379  sin2nPhiB /= multB;
2380 
2381  dNdetaRef->Fill(eta, cent, multA+multB);
2382  }
2383  Double_t two = w2Two / w2;
2384 
2385  Double_t qc2 = two;
2386  if (qc2 >= 0) cumu2r->Fill(eta, cent, TMath::Sqrt(qc2));
2387 
2388  if ((fFlags & kNUAcorr)) {
2389  // Old nua
2390  qc2 -= cosP1nPhiA*cosP1nPhiB + sinP1nPhiA*sinP1nPhiB;
2391  // Extra NUA term from 2n cosines and sines
2392  qc2 /= (1+(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB));
2393  }
2394  if (qc2 <= 0) {
2395  if (fDebug > 0)
2396  AliInfo(Form("%s: QC_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f - skipping",
2397  fType.Data(), n, qc2, eta, cent));
2398  quality->Fill((n-2)*4+2, Int_t(cent));
2399  continue;
2400  }
2401  Double_t vnTwo = TMath::Sqrt(qc2);
2402  if (!TMath::IsNaN(vnTwo)) {
2403  quality->Fill((n-2)*4+1, Int_t(cent));
2404  if ((fFlags & kNUAcorr)) cumu2rNUAold->Fill(eta, cent, vnTwo);
2405  }
2406 
2407  // 2-particle differential flow
2408  Double_t w2pTwoPrimeA = cumuDiff->GetBinContent(etaBin, cBin, kW2Two);
2409  Double_t w2pA = cumuDiff->GetBinContent(etaBin, cBin, kW2);
2410  Double_t w2pTwoPrimeB = cumuDiff->GetBinContent(etaBin, cBin, kW4Four);
2411  Double_t w2pB = cumuDiff->GetBinContent(etaBin, cBin, kW4);
2412  if (w2pA == 0 || w2pB == 0) continue;
2413  Double_t cosP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(n));
2414  Double_t sinP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(n));
2415  Double_t cos2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(2*n));
2416  Double_t sin2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(2*n));
2417  Double_t mult = fCumuNUADiff->GetBinContent(etaBin, cBin, 0);
2418  if (mult == 0) continue;
2419  cosP1nPsi /= mult;
2420  sinP1nPsi /= mult;
2421  cos2nPsi /= mult;
2422  sin2nPsi /= mult;
2423  Double_t twoPrimeA = w2pTwoPrimeA / w2pA;
2424  Double_t twoPrimeB = w2pTwoPrimeB / w2pB;
2425  dNdetaDiff->Fill(eta, cent, mult);
2426 
2427  Double_t qc2PrimeA = twoPrimeA;
2428  Double_t qc2PrimeB = twoPrimeB;
2429  if (qc2PrimeA*qc2PrimeB >= 0) {
2430  cumu2a->Fill(eta, cent, qc2PrimeA);
2431  cumu2b->Fill(eta, cent, qc2PrimeB);
2432  }
2433  if ((fFlags & kNUAcorr)) {
2434  // Old nua
2435  qc2PrimeA -= cosP1nPsi*cosP1nPhiA + sinP1nPsi*sinP1nPhiA;
2436  qc2PrimeB -= cosP1nPsi*cosP1nPhiB + sinP1nPsi*sinP1nPhiB; // Is this OK?
2437  // Extra NUA term from 2n cosines and sines
2438  if (cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA != -1.) qc2PrimeA /= (1.+(cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA));
2439  if (cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB != -1.) qc2PrimeB /= (1.+(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
2440  }
2441  if (!TMath::IsNaN(qc2PrimeA) && !TMath::IsNaN(qc2PrimeB) && qc2 != 0) {
2442  if (qc2PrimeA*qc2PrimeB >= 0) {
2443  quality->Fill((n-2)*4+3, Int_t(cent));
2444  if ((fFlags & kNUAcorr)) cumu2aNUAold->Fill(eta, cent, qc2PrimeA);
2445  if ((fFlags & kNUAcorr)) cumu2bNUAold->Fill(eta, cent, qc2PrimeB);
2446  }
2447  }
2448  else
2449  quality->Fill((n-2)*4+4, Int_t(cent));
2450  if (fDebug > 1)
2451  AliInfo(Form("%s: QC'a_%d{2} = %1.3f, QC'b_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f",
2452  fType.Data(), n, qc2PrimeA, n, qc2PrimeB, eta, cent));
2453  } // End of eta loop
2454  } // End of centrality loop
2455  } // End of moment
2456 
2457  return;
2458 }
2459 //_____________________________________________________________________
2461 {
2462  //
2463  // Function to solve the coupled flow equations
2464  // We solve it by using matrix calculations:
2465  // A*v_n = V => v_n = A^-1*V
2466  // First we set up a TMatrixD container to make ROOT
2467  // do the inversions in an efficient way, we multiply the current <<2>> estimates.
2468  // Then we fill new TH2D's if the corrected <<2>>'s (cumuNUA object).
2469  //
2470  // Parameters:
2471  // cumu: CumuHistos object - uncorrected
2472  // type: Reference ('r') or differential ('d') or ('a' or 'b') for 3 correlator
2473  //
2474 
2475  // We start by declaring Matrix and vector objects, as their constructors are quite heavy
2476  TMatrixD mNUA(fMaxMoment-1, fMaxMoment-1);
2477  TVectorD vQC2(fMaxMoment-1);
2478 
2479  for (Int_t cBin = 1; cBin <= cumu.Get(type, 2, CumuHistos::kNUAOld)->GetNbinsY(); cBin++) { // cent loop
2480  Double_t cent = cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinCenter(cBin);
2481  for (Int_t eBin = 1; eBin <= cumu.Get(type, 2, CumuHistos::kNUAOld)->GetNbinsX(); eBin++) { // eta loop
2482  Double_t eta = cumu.Get(type, 2, CumuHistos::kNUAOld)->GetXaxis()->GetBinCenter(eBin);
2483  mNUA.Zero(); // reset matrix
2484  vQC2.Zero(); // reset vector
2485  for (Int_t n = 0; n < fMaxMoment-1; n++) { // moment loop
2486  vQC2(n) = static_cast<TH2D*>(cumu.Get(type, n+2, CumuHistos::kNUAOld))->GetBinContent(eBin, cBin);
2487  if (type == 'r' || type == 'R') vQC2(n) *= vQC2(n); // going back to <<2>>
2488  for (Int_t m = 0; m < fMaxMoment-1; m++) { // cross moment
2489  mNUA(n,m) = CalculateNUAMatrixElement(n, m, type, eBin, cBin);
2490  } // End of cross moment loop
2491  } // End of moment loop
2492  // Invert matrix
2493  Double_t det = 0;
2494  mNUA.Invert(&det);
2495  // If determinant is non-zero we go with corrected results
2496  if (det != 0 ) vQC2 = mNUA*vQC2;
2497  else AliWarning(Form("Determinant == 0 - cent: %d-%d, eta: %f, type: '%c', data: %s, vtx: %d-%d%s",
2498  Int_t(cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinLowEdge(cBin)),
2499  Int_t(cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinUpEdge(cBin)),
2500  cumu.Get(type, 2, CumuHistos::kNUAOld)->GetXaxis()->GetBinCenter(eBin),
2501  type, fType.Data(), fVzMin, fVzMax,
2502  GetQCType(fFlags, kTRUE)));
2503  // Go back to v_n for ref. keep <<2'>> for diff. flow).
2504  for (Int_t n = 0; n < fMaxMoment-1; n++) {
2505  Double_t vnTwo = 0;
2506  if (type == 'r' || type == 'R')
2507  vnTwo = (vQC2(n) > 0. ? TMath::Sqrt(vQC2(n)) : 0.);
2508  else {
2509  // is really more <<2'>> in this case
2510  vnTwo = vQC2(n);
2511  }
2512  // Fill in corrected v_n
2513  if (vnTwo != 0) static_cast<TH2D*>(cumu.Get(type, n+2, CumuHistos::kNUA))->Fill(eta, cent, vnTwo);
2514  } // End of moment loop
2515  } // End of eta loop
2516  } // End of centrality loop
2517  return;
2518 }
2519 //_____________________________________________________________________
2521 {
2522  //
2523  // Calculates the (n,m)-th element in the NUA matrix: 1 if n == m, otherwise:
2524  // <<cos[(n-m)phi1]>>*<<cos[(n-m)phi2]>> + <<sin[(n-m)phi1]>>*<<sin[(n-m)phi2]>>
2525  // NUA(n,m) = -----------------------------------------------------------------------------
2526  // 1 + <<cos[2nphi1]>>*<<cos[2nphi2]>> + <<sin[2nphi1]>>*<<sin[2nphi2]>>
2527  //
2528  // <<cos[(n+m)phi1]>>*<<cos[(n+m)phi2]>> + <<sin[(n+m)phi1]>>*<<sin[(n+m)phi2]>>
2529  // + -----------------------------------------------------------------------------
2530  // 1 + <<cos[2nphi1]>>*<<cos[2nphi2]>> + <<sin[2nphi1]>>*<<sin[2nphi2]>>
2531  //
2532  // Parameters:
2533  // n: row
2534  // m: coumn
2535  // type: Reference ('r') or differential ('d') or ('a' or 'b')
2536  // binA: eta bin of phi1
2537  // cBin: centrality bin
2538  //
2539  // Return: NUA(n,m)
2540  //
2541  if (n == m) return 1.;
2542  n += 2;
2543  m += 2;
2544 
2545  Double_t cosnMmPhi1 = 0, cosnMmPhi2 = 0, sinnMmPhi1 = 0, sinnMmPhi2 = 0;
2546  Double_t cosnPmPhi1 = 0, cosnPmPhi2 = 0, sinnPmPhi1 = 0, sinnPmPhi2 = 0;
2547  Double_t cos2nPhi1 = 0, cos2nPhi2 = 0, sin2nPhi1 = 0, sin2nPhi2 = 0;
2548 
2549  // reference flow
2550  if (type == 'r' || type == 'R') {
2551  if ((fFlags & k3Cor)) {
2552  Double_t eta = fCumuNUARef->GetXaxis()->GetBinCenter(binA);
2553  Int_t i = 0;
2554  while (fEtaLims[i] < eta) i++;
2555  Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2556  GetLimits(i-1, aLow, aHigh, bLow, bHigh);
2557  Double_t multA = 0, multB = 0;
2558  for (Int_t a = aLow; a <= aHigh; a++) {
2559  cosnMmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n-m));
2560  sinnMmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n-m));
2561  cosnPmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n+m));
2562  sinnPmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n+m));
2563  cos2nPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(2*n));
2564  sin2nPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(2*n));
2565  multA += fCumuNUARef->GetBinContent(a, cBin, 0);
2566  }
2567  for (Int_t b = bLow; b <= bHigh; b++) {
2568  cosnMmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n-m));
2569  sinnMmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n-m));
2570  cosnPmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n+m));
2571  sinnPmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n+m));
2572  cos2nPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(2*n));
2573  sin2nPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(2*n));
2574  multB += fCumuNUARef->GetBinContent(b, cBin, 0);
2575  }
2576  if (multA == 0 || multB == 0) {
2577  if (fDebug > 0) AliWarning("multA or multB == 0 in matrix elements, aborting NUA");
2578  return 0.;
2579  }
2580  cosnMmPhi1 /= multA;
2581  sinnMmPhi1 /= multA;
2582  cosnPmPhi1 /= multA;
2583  sinnPmPhi1 /= multA;
2584  cos2nPhi1 /= multA;
2585  sin2nPhi1 /= multA;
2586  cosnMmPhi2 /= multB;
2587  sinnMmPhi2 /= multB;
2588  cosnPmPhi2 /= multB;
2589  sinnPmPhi2 /= multB;
2590  cos2nPhi2 /= multB;
2591  sin2nPhi2 /= multB;
2592  } else {
2593  Int_t binB = fCumuNUARef->GetXaxis()->FindBin(-fCumuNUARef->GetXaxis()->GetBinCenter(binA));
2594  cosnMmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2595  sinnMmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2596  cosnPmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2597  sinnPmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2598  cos2nPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2599  sin2nPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2600  cosnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n-m));
2601  sinnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n-m));
2602  cosnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n+m));
2603  sinnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n+m));
2604  cos2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(2*n));
2605  sin2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(2*n));
2606  }
2607  } // differential flow
2608  else if (type == 'd' || type == 'D') {
2609  Int_t binB = fCumuNUARef->GetXaxis()->FindBin(-fCumuNUADiff->GetXaxis()->GetBinCenter(binA));
2610  cosnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2611  sinnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2612  cosnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2613  sinnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2614  cos2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2615  sin2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2616  cosnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n-m));
2617  sinnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n-m));
2618  cosnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n+m));
2619  sinnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n+m));
2620  cos2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(2*n));
2621  sin2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(2*n));
2622  } // 3 correlator part a or b
2623  else if (type == 'a' || type == 'A' || type == 'b' || type == 'B') {
2624  Double_t mult1 = 0, mult2 = 0;
2625  // POIs
2626  cosnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2627  sinnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2628  cosnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2629  sinnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2630  cos2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2631  sin2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2632  mult1 = fCumuNUADiff->GetBinContent(binA, cBin, 0);
2633  // RPs
2634  Double_t eta = fCumuNUADiff->GetXaxis()->GetBinCenter(binA);
2635  Int_t i = 0;
2636  while (fEtaLims[i] < eta) i++;
2637  Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2638  GetLimits(i-1, aLow, aHigh, bLow, bHigh);
2639  Int_t lLow = ((type == 'a' || type == 'A') ? aLow : bLow);
2640  Int_t lHigh = ((type == 'a' || type == 'A') ? aHigh : bHigh);
2641  for (Int_t l = lLow; l <= lHigh; l++) {
2642  cosnMmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(n-m));
2643  sinnMmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(n-m));
2644  cosnPmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(n+m));
2645  sinnPmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(n+m));
2646  cos2nPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(2*n));
2647  sin2nPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(2*n));
2648  mult2 += fCumuNUARef->GetBinContent(l, cBin, 0);
2649  }
2650  if (mult1 == 0 || mult2 == 0) {
2651  if (fDebug > 0) AliWarning("mult1 or mult2 == 0 in matrix elements, aborting NUA");
2652  return 0.;
2653  }
2654  cosnMmPhi1 /= mult1;
2655  sinnMmPhi1 /= mult1;
2656  cosnPmPhi1 /= mult1;
2657  sinnPmPhi1 /= mult1;
2658  cos2nPhi1 /= mult1;
2659  sin2nPhi1 /= mult1;
2660  cosnMmPhi2 /= mult2;
2661  sinnMmPhi2 /= mult2;
2662  cosnPmPhi2 /= mult2;
2663  sinnPmPhi2 /= mult2;
2664  cos2nPhi2 /= mult2;
2665  sin2nPhi2 /= mult2;
2666  }
2667 
2668  // Actual calculation
2669  Double_t e = cosnMmPhi1*cosnMmPhi2 + sinnMmPhi1*sinnMmPhi2 + cosnPmPhi1*cosnPmPhi2 + sinnPmPhi1*sinnPmPhi2;
2670  Double_t den = 1 + cos2nPhi1*cos2nPhi2 + sin2nPhi1*sin2nPhi2;
2671  if (den != 0) e /= den;
2672  else return 0.;
2673 
2674  return e;
2675 }
2676 //_____________________________________________________________________
2678 {
2679  //
2680  // Add up vertex bins with flow results
2681  //
2682  // Parameters:
2683  // cumu: CumuHistos object with vtxbin results
2684  // list: Outout list with added results
2685  // nNUA: # of NUA histograms to loop over
2686  //
2687  TH2D* vtxHist = 0;
2688  TProfile2D* avgProf = 0;
2689  TString name;
2690  Int_t nT = ((fFlags & k3Cor) ? 3 : 2);
2691  Char_t ct = '\0';
2692  for (UInt_t nua = 0; nua <= nNUA; nua++) { // NUA loop
2693  for (Int_t n = 2; n <= fMaxMoment; n++) { // moment loop
2694  for (Int_t t = 0; t < nT; t++) { // type loop (r/d/a/b)
2695  // Find type
2696  switch (t) {
2697  case 0: ct = 'r'; break;
2698  case 1: ct = ((fFlags & k3Cor) ? 'a' : 'd'); break;
2699  case 2: ct = 'b'; break;
2700  default: ct = '\0'; break;
2701  }
2702  vtxHist = static_cast<TH2D*>(cumu.Get(ct,n,nua));
2703  if (!vtxHist) {
2704  AliWarning("VertexBin::AddVertexBins: vtxHist not found!");
2705  continue;
2706  }
2707  name = vtxHist->GetName();
2708  // Strip name of vtx info
2709  Ssiz_t l = name.Last('x')-3;
2710  name.Resize(l);
2711  avgProf = (TProfile2D*)list->FindObject(name.Data());
2712  // if no output profile yet, make one
2713  if (!avgProf) {
2714  avgProf = new TProfile2D(name.Data(), name.Data(),
2715  vtxHist->GetNbinsX(), vtxHist->GetXaxis()->GetXmin(), vtxHist->GetXaxis()->GetXmax(),
2716  vtxHist->GetNbinsY(), vtxHist->GetYaxis()->GetXmin(), vtxHist->GetYaxis()->GetXmax());
2717  if (vtxHist->GetXaxis()->IsVariableBinSize())
2718  avgProf->GetXaxis()->Set(vtxHist->GetNbinsX(), vtxHist->GetXaxis()->GetXbins()->GetArray());
2719  if (vtxHist->GetYaxis()->IsVariableBinSize())
2720  avgProf->GetYaxis()->Set(vtxHist->GetNbinsY(), vtxHist->GetYaxis()->GetXbins()->GetArray());
2721  list->Add(avgProf);
2722  }
2723  // Fill in, cannot be done with Add function.
2724  for (Int_t e = 1; e <= vtxHist->GetNbinsX(); e++) { // eta loop
2725  Double_t eta = vtxHist->GetXaxis()->GetBinCenter(e);
2726  for (Int_t c = 1; c <= vtxHist->GetNbinsY(); c++) { // cent loop
2727  Double_t cent = vtxHist->GetYaxis()->GetBinCenter(c);
2728  Double_t cont = vtxHist->GetBinContent(e, c);
2729  if (cont == 0) continue;
2730  avgProf->Fill(eta, cent, cont);
2731  } // End of cent loop
2732  } // End of eta loop
2733  } // End of type loop
2734  } // End of moment loop
2735  } // End of nua loop
2736 }
2737 //_____________________________________________________________________
2739 {
2740  //
2741  // Get the bin number of <<cos(nphi)>>
2742  //
2743  // Parameters:
2744  // n: moment
2745  //
2746  // Return: bin number
2747  //
2748  Int_t bin = 0;
2749  n = TMath::Abs(n);
2750 
2751  if (n == 0) bin = fMaxMoment*4-1;
2752  else bin = n*2-1;
2753 
2754  return bin;
2755 }
2756 //_____________________________________________________________________
2758 {
2759  //
2760  // Get the bin number of <<sin(nphi)>>
2761  //
2762  // Parameters:
2763  // n: moment
2764  //
2765  // Return: bin number
2766  //
2767  Int_t bin = 0;
2768  n = TMath::Abs(n);
2769 
2770  if (n == 0) bin = fMaxMoment*4;
2771  else bin = n*2;
2772 
2773  return bin;
2774 }
2775 //_____________________________________________________________________
2777 {
2778  //
2779  // Setup NUA labels on axis
2780  //
2781  // Parameters:
2782  // a: Axis to set up
2783  //
2784  if (a->GetNbins() != GetBinNumberSin()) AliFatal("SetupNUALabels: Wrong number of bins on axis");
2785 
2786  Int_t i = 1, j= 1;
2787  while (i <= a->GetNbins()) {
2788  a->SetBinLabel(i++, Form("<<cos(%d#varphi)>>", j));
2789  a->SetBinLabel(i++, Form("<<sin(%d#varphi)>>", j++));
2790  }
2791 
2792  return;
2793 }
2794 //_____________________________________________________________________
2796 {
2797  //
2798  // Add a histogram for checking the analysis quality
2799  //
2800  // Parameters:
2801  // const char*: name of data type
2802  //
2803  // Return: histogram for tracking successful calculations
2804  //
2805  Int_t nBins = ((fFlags & kStdQC) ? 8 : 4);
2806  TH2I* quality = new TH2I(name, name, (fMaxMoment-1)*nBins, 1, (fMaxMoment-1)*nBins+1, fCumuNUARef->GetNbinsY(),
2807  fCumuNUARef->GetYaxis()->GetXmin(), fCumuNUARef->GetYaxis()->GetXmax());
2808  quality->GetYaxis()->Set(fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXbins()->GetArray());
2809  for (Int_t i = 2, j = 1; i <= fMaxMoment; i++) {
2810  quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{2} > 0", i));
2811  quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{2} <= 0", i));
2812  quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{2} > 0", i));
2813  quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{2} <= 0", i));
2814  if ((fFlags & kStdQC)) {
2815  quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{4} < 0", i));
2816  quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{4} >= 0", i));
2817  quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{4} < 0", i));
2818  quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{4} >= 0", i));
2819  }
2820  }
2821 
2822  return quality;
2823 }
2824 //_____________________________________________________________________
2826 {
2827  //
2828  // Setup a TH2D for the output
2829  //
2830  // Parameters:
2831  // qc # of particle correlations
2832  // n flow moment
2833  // ref Type: r/d/a/b
2834  // nua For nua corrected hists?
2835  //
2836  // Return: 2D hist for results
2837  //
2838  Bool_t ref = ((ctype[0] == 'R') || (ctype[0] == 'r'));
2839  TAxis* xAxis = (ref ? fCumuNUARef->GetXaxis() : fCumuNUADiff->GetXaxis());
2840  TAxis* yAxis = (ref ? fCumuNUARef->GetYaxis() : fCumuNUADiff->GetYaxis());
2841  TString ntype = "";
2842  switch (nua) {
2843  case CumuHistos::kNoNUA: ntype = "Un"; break;
2844  case CumuHistos::kNUAOld: ntype = "NUAOld"; break;
2845  case CumuHistos::kNUA: ntype = "NUA"; break;
2846  default: break;
2847  }
2848  TH2D* h = new TH2D(Form("%sQC%d_v%d_%s_%sCorr%s_vtx_%d_%d",
2849  fType.Data(), qc, n, ctype, ntype.Data(),
2850  GetQCType(fFlags), fVzMin, fVzMax),
2851  Form("%sQC%d_v%d_%s_%sCorr%s_vtx_%d_%d",
2852  fType.Data(), qc, n, ctype, ntype.Data(),
2853  GetQCType(fFlags), fVzMin, fVzMax),
2854  xAxis->GetNbins(), xAxis->GetXmin(), xAxis->GetXmax(),
2855  yAxis->GetNbins(), yAxis->GetXmin(), yAxis->GetXmax());
2856  if (xAxis->IsVariableBinSize()) h->GetXaxis()->Set(xAxis->GetNbins(), xAxis->GetXbins()->GetArray());
2857  h->GetYaxis()->Set(yAxis->GetNbins(), yAxis->GetXbins()->GetArray());
2858 
2859  return h;
2860 }
2861 //_____________________________________________________________________
2863 {
2864  //
2865  // Print the setup of the flow task
2866  //
2867  Printf("=======================================================");
2868  Printf("%s::Print", this->IsA()->GetName());
2869  Printf("Forward detector: :\t%s", ((fFlowFlags & kFMD) ? "FMD" : "VZERO"));
2870  Printf("Number of bins in vertex axis :\t%d", fVtxAxis->GetNbins());
2871  Printf("Range of vertex axis :\t[%3.1f,%3.1f]",
2872  fVtxAxis->GetXmin(), fVtxAxis->GetXmax());
2873  Printf("Number of bins in centrality axis :\t%d", fCentAxis->GetNbins());
2874  printf("Centrality binning :\t");
2875  for (Int_t cBin = 1; cBin <= fCentAxis->GetNbins(); cBin++) {
2876  printf("%02d-%02d%% ", Int_t(fCentAxis->GetBinLowEdge(cBin)), Int_t(fCentAxis->GetBinUpEdge(cBin)));
2877  if (cBin == fCentAxis->GetNbins()) printf("\n");
2878  else if (cBin % 4 == 0) printf("\n\t\t\t\t\t");
2879  }
2880  printf("Doing flow analysis for :\t");
2881  for (Int_t n = 2; n <= fMaxMoment; n++) printf("v%d ", n);
2882  printf("\n");
2883  TString type = "Standard QC{2} and QC{4} calculations.";
2884  if ((fFlowFlags & kEtaGap)) type = "QC{2} with a rapidity gap.";
2885  if ((fFlowFlags & k3Cor)) type = "QC{2} with 3 correlators.";
2886  if ((fFlowFlags & kTPC) == kTPC) type.ReplaceAll(".", " with TPC tracks for reference.");
2887  if ((fFlowFlags & kHybrid) == kHybrid) type.ReplaceAll(".", " with hybrid tracks for reference.");
2888  Printf("QC calculation type :\t%s", type.Data());
2889  Printf("Symmetrize ref. flow wrt. eta = 0 :\t%s", ((fFlowFlags & kSymEta) ? "true" : "false"));
2890  Printf("Apply NUA correction terms :\t%s", ((fFlowFlags & kNUAcorr) ? "true" : "false"));
2891  Printf("Satellite vertex flag :\t%s", ((fFlowFlags & kSatVtx) ? "true" : "false"));
2892  Printf("FMD sigma cut: :\t%f", fFMDCut);
2893  Printf("SPD sigma cut: :\t%f", fSPDCut);
2894  if ((fFlowFlags & kEtaGap) || (fFlowFlags & kTracks))
2895  Printf("Eta gap: :\t%f", fEtaGap);
2896  Printf("=======================================================");
2897 }
2898 //_____________________________________________________________________
2900 {
2901  //
2902  // Get the type of the QC calculations
2903  // Used for naming of objects in the VertexBin class,
2904  // important to avoid memory problems when running multiple
2905  // initializations of the class with different setups
2906  //
2907  // Parameters:
2908  // flags: Flow flags for type determination
2909  // prependUS: Prepend an underscore (_) to the name
2910  //
2911  // Return: QC calculation type
2912  //
2913  static TString type = "";
2914  if ((flags & kStdQC)) type = "StdQC";
2915  else if ((flags & kEtaGap)) type = "EtaGap";
2916  else if ((flags & k3Cor)) type = "3Cor";
2917  else type = "UNKNOWN";
2918  if (prependUS) type.Prepend("_");
2919  if ((flags & kTPC) == kTPC) type.Append("TPCTr");
2920  if ((flags & kHybrid) == kHybrid) type.Append("HybTr");
2921 
2922  return type.Data();
2923 }
2924 //_____________________________________________________________________
2926 {
2927  //
2928  // Get histogram/profile for type t and moment n
2929  //
2930  // Parameters:
2931  // t: type = 'r'/'d'
2932  // n: moment
2933  // nua: NUA type
2934  //
2935  n = GetPos(n, nua);
2936  if (n < 0) AliFatal(Form("CumuHistos out of range: (%c,%d)", t, n));
2937 
2938  TH1* h = 0;
2939  Int_t pos = -1;
2940  if (t == 'r' || t == 'R') pos = n;
2941  else if (t == 'd' || t == 'D') pos = n;
2942  else if (t == 'a' || t == 'A') pos = 2*n;
2943  else if (t == 'b' || t == 'B') pos = 2*n+1;
2944  else AliFatal(Form("CumuHistos wrong type: %c", t));
2945 
2946  if (t == 'r' || t == 'R') {
2947  if (pos < fRefHists->GetEntries()) {
2948  h = (TH1*)fRefHists->At(pos);
2949  }
2950  } else if (pos < fDiffHists->GetEntries()) {
2951  h = (TH1*)fDiffHists->At(pos);
2952  }
2953  if (!h) AliFatal(Form("No hist found in list %c at pos %d", t, pos));
2954 
2955  return h;
2956 }
2957 //_____________________________________________________________________
2959 {
2960  //
2961  // Connect an output list with this object
2962  //
2963  // Parameters:
2964  // name: base name
2965  // l: list to keep outputs in
2966  //
2967  TString ref = name;
2968  ref.ReplaceAll("Cumu","CumuRef");
2969  fRefHists = (TList*)l->FindObject(ref.Data());
2970  if (!fRefHists) {
2971  fRefHists = new TList();
2972  fRefHists->SetName(ref.Data());
2973  l->Add(fRefHists);
2974  } else {
2975  // Check that the list correspond to fMaxMoments
2976  if (fRefHists->GetEntries() != (fMaxMoment-1.)*(fNUA+1))
2977  AliError("CumuHistos::ConnectList Wrong number of hists in ref list,"
2978  "you are doing something wrong!");
2979  }
2980  TString diff = name;
2981  diff.ReplaceAll("Cumu","CumuDiff");
2982  fDiffHists = (TList*)l->FindObject(diff.Data());
2983  if (!fDiffHists) {
2984  fDiffHists = new TList();
2985  fDiffHists->SetName(diff.Data());
2986  l->Add(fDiffHists);
2987  } else {
2988  // Check that the list correspond to fMaxMoment
2989  if ((fDiffHists->GetEntries() != (fMaxMoment-1.)*(fNUA+1)) &&
2990  (fDiffHists->GetEntries() != 2*(fMaxMoment-1.)*(fNUA+1)))
2991  AliError(Form("CumuHistos::ConnectList Wrong number of hists in diff list,"
2992  "you are doing something wrong! (%s)",name.Data()));
2993  }
2994 
2995 }
2996 //_____________________________________________________________________
2998 {
2999  //
3000  // Add a histogram to this objects lists
3001  //
3002  // Parameters:
3003  // h: histogram/profile to add
3004  //
3005  TString name = h->GetName();
3006  if (name.Contains("Ref")) {
3007  if (fRefHists) fRefHists->Add(h);
3008  else AliFatal("CumuHistos::Add() - fRefHists does not exist");
3009  // Check that the list correspond to fMaxMoments
3010  if (fRefHists->GetEntries() > (fNUA+1)*(fMaxMoment-1.))
3011  AliError("CumuHistos::Add wrong number of hists in ref list, "
3012  "you are doing something wrong!");
3013  }
3014  else if (name.Contains("Diff")) {
3015  if (fDiffHists) fDiffHists->Add(h);
3016  else AliFatal("CumuHistos::Add() - fDiffHists does not exist");
3017  // Check that the list correspond to fMaxMoment
3018  if (fDiffHists->GetEntries() > 2*(fNUA+1)*(fMaxMoment-1.))
3019  AliError("CumuHistos::Add wrong number of hists in diff list, "
3020  "you are doing something wrong!");
3021  }
3022  return;
3023 }
3024 //_____________________________________________________________________
3026 {
3027  //
3028  // Get position in list of moment n objects
3029  // To take care of different numbering schemes
3030  //
3031  // Parameters:
3032  // n: moment
3033  // nua: # of nua corrections applied
3034  //
3035  // Return: position #
3036  //
3037  if (n > fMaxMoment) return -1;
3038  else return (n-2)+nua*(fMaxMoment-1);
3039 }
3040 //_____________________________________________________________________
3041 //
3042 //
3043 // EOF
return jsonbuilder str().c_str()
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
TList * list
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)
Int_t GetPos(Int_t n, UInt_t nua) const
Int_t mode
Definition: anaM.C:40
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
ClassImp(AliForwardFlowTaskQC) AliForwardFlowTaskQC
TH2D * MakeOutputHist(Int_t qc, Int_t n, const Char_t *ctype, UInt_t nua) const
Float_t nEvents[nProd]
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