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