3 // Calculate flow in the forward and central regions using the Q cumulants method.
9 // - AnalysisResults.root or forward_flow.root
13 #include <TInterpreter.h>
19 #include <TProfile2D.h>
20 #include <TParameter.h>
25 #include "AliForwardFlowTaskQC.h"
26 #include "AliAnalysisManager.h"
27 #include "AliAODHandler.h"
28 #include "AliAODInputHandler.h"
29 #include "AliAODForwardMult.h"
30 #include "AliAODCentralMult.h"
31 #include "AliAODEvent.h"
32 #include "AliForwardUtil.h"
33 #include "AliVVZERO.h"
34 #include "AliAODVertex.h"
35 #include "AliCentrality.h"
36 #include "AliESDEvent.h"
37 #include "AliVTrack.h"
38 #include "AliESDtrack.h"
39 #include "AliAODTrack.h"
40 #include "AliAnalysisFilter.h"
41 #include "AliAODMCHeader.h"
42 #include "AliForwardFlowWeights.h"
44 ClassImp(AliForwardFlowTaskQC)
49 AliForwardFlowTaskQC::AliForwardFlowTaskQC()
50 : AliAnalysisTaskSE(),
51 fVtxAxis(), // Axis to control vertex binning
52 fCentAxis(), // Axis to control centrality/multiplicity binning
53 fFMDCut(-1), // FMD sigma cut
54 fSPDCut(-1), // SPD sigma cut
55 fFlowFlags(0), // Flow flags
56 fEtaGap(0.), // Eta gap value
57 fBinsForward(), // List with forward flow hists
58 fBinsCentral(), // List with central flow hists
59 fSumList(0), // Event sum list
60 fOutputList(0), // Result output list
61 fAOD(0), // AOD input event
62 fTrackCuts(0), // ESD track cuts
63 fMaxMoment(0), // Max flow moment
64 fVtx(1111), // Z vertex coordinate
65 fCent(-1), // Centrality
66 fHistdNdedpV0(), // Hist for v0
67 fHistdNdedp3Cor(), // Hist for combining detectors
68 fHistFMDSPDCorr(), // FMD SPD correlation
69 fHistCent(), // Hist for centrality
70 fHistVertexSel(), // Hist for selected vertices
71 fHistEventSel() // Hist for event selection
74 // Default constructor
77 //_____________________________________________________________________
78 AliForwardFlowTaskQC::AliForwardFlowTaskQC(const char* name)
79 : AliAnalysisTaskSE(name),
80 fVtxAxis(), // Axis to control vertex binning
81 fCentAxis(), // Axis to control centrality/multiplicity binning
82 fFMDCut(-1), // FMD sigma cut
83 fSPDCut(-1), // SPD sigma cut
84 fFlowFlags(0x0), // Flow flags
85 fEtaGap(0.), // Eta gap value
86 fBinsForward(), // List with forward flow hists
87 fBinsCentral(), // List with central flow hists
88 fSumList(0), // Event sum list
89 fOutputList(0), // Result output list
90 fAOD(0), // AOD input event
91 fTrackCuts(0), // ESD track cuts
92 fMaxMoment(4), // Max flow moment
93 fVtx(1111), // Z vertex coordinate
94 fCent(-1), // Centrality
95 fHistdNdedpV0(), // Histo for v0
96 fHistdNdedp3Cor(), // Histo for combining detectors
97 fHistFMDSPDCorr(), // FMD SPD correlation
98 fHistCent(), // Hist for centrality
99 fHistVertexSel(), // Hist for selected vertices
100 fHistEventSel() // Hist for event selection
106 // name: Name of task
108 DefineOutput(1, TList::Class());
109 DefineOutput(2, TList::Class());
111 //_____________________________________________________________________
112 AliForwardFlowTaskQC::AliForwardFlowTaskQC(const AliForwardFlowTaskQC& o)
113 : AliAnalysisTaskSE(o),
114 fVtxAxis(o.fVtxAxis), // Axis to control vertex binning
115 fCentAxis(o.fCentAxis), // Array to control centrality/multiplicity binning
116 fFMDCut(o.fFMDCut), // FMD sigma cut
117 fSPDCut(o.fSPDCut), // SPD sigma cut
118 fFlowFlags(o.fFlowFlags), // Flow flags
119 fEtaGap(o.fEtaGap), // Eta gap value
120 fBinsForward(), // List with forward flow hists
121 fBinsCentral(), // List with central flow hists
122 fSumList(o.fSumList), // Event sum list
123 fOutputList(o.fOutputList), // Result output list
124 fAOD(o.fAOD), // AOD input event
125 fTrackCuts(o.fTrackCuts), // ESD track cuts
126 fMaxMoment(o.fMaxMoment), // Flow moments
127 fVtx(o.fVtx), // Z vertex coordinate
128 fCent(o.fCent), // Centrality
129 fHistdNdedpV0(o.fHistdNdedpV0), // Histo for v0
130 fHistdNdedp3Cor(o.fHistdNdedp3Cor),// Histo for combining detectors
131 fHistFMDSPDCorr(o.fHistFMDSPDCorr),// FMD SPD correlation
132 fHistCent(o.fHistCent), // Hist for centrality
133 fHistVertexSel(o.fHistVertexSel), // Hist for selected vertices
134 fHistEventSel(o.fHistEventSel) // Hist for event selection
140 // o: Object to copy from
143 //_____________________________________________________________________
144 AliForwardFlowTaskQC&
145 AliForwardFlowTaskQC::operator=(const AliForwardFlowTaskQC& o)
148 // Assignment operator
150 if (&o == this) return *this;
151 fVtxAxis = o.fVtxAxis;
152 fCentAxis = o.fCentAxis;
155 fFlowFlags = o.fFlowFlags;
157 fSumList = o.fSumList;
158 fOutputList = o.fOutputList;
160 fTrackCuts = o.fTrackCuts;
161 fMaxMoment = o.fMaxMoment;
164 fHistdNdedpV0 = o.fHistdNdedpV0;
165 fHistdNdedp3Cor = o.fHistdNdedp3Cor;
166 fHistFMDSPDCorr = o.fHistFMDSPDCorr;
167 fHistCent = o.fHistCent;
168 fHistVertexSel = o.fHistVertexSel;
169 fHistEventSel = o.fHistEventSel;
173 //_____________________________________________________________________
174 void AliForwardFlowTaskQC::SetFlowFlags(UShort_t flags)
177 // Set flow flags, making sure the detector setup is right
182 if ((flags & kFMD) && (flags & kVZERO))
183 AliFatal("Cannot do analysis on more than one forward detector!");
184 else if (!(flags & kFMD) && !(flags & kVZERO))
185 AliFatal("You need to add a forward detector!");
186 else fFlowFlags = flags;
188 //_____________________________________________________________________
189 void AliForwardFlowTaskQC::UserCreateOutputObjects()
192 // Create output objects
196 if ((fFlowFlags & kTracks) && !fTrackCuts) AliFatal("No track cuts set!");
199 PostData(1, fSumList);
201 //_____________________________________________________________________
202 void AliForwardFlowTaskQC::InitVertexBins()
205 // Init vertexbin objects for forward and central detectors, and add them to the lists
207 for (Int_t v = 1; v <= fVtxAxis->GetNbins(); v++) {
208 Int_t vL = Int_t(fVtxAxis->GetBinLowEdge(v));
209 Int_t vH = Int_t(fVtxAxis->GetBinUpEdge(v));
210 if ((fFlowFlags & kFMD)) {
211 fBinsForward.Add(new VertexBin(vL, vH, fMaxMoment, "FMD", fFlowFlags, fFMDCut, fEtaGap));
212 if (!(fFlowFlags & k3Cor))
213 fBinsCentral.Add(new VertexBin(vL, vH, fMaxMoment, "SPD-FMD", fFlowFlags|kNUAcorr|kSPD, fSPDCut, fEtaGap));
215 else if ((fFlowFlags & kVZERO)) {
216 fBinsForward.Add(new VertexBin(vL, vH, fMaxMoment, "VZERO", fFlowFlags, 0, fEtaGap));
217 if ((fFlowFlags & kEtaGap) && !(fFlowFlags & kTracks))
218 fBinsCentral.Add(new VertexBin(vL, vH, fMaxMoment, "SPD-VZERO", fFlowFlags|kNUAcorr|kSPD, fSPDCut, fEtaGap));
222 //_____________________________________________________________________
223 void AliForwardFlowTaskQC::InitHists()
226 // Init histograms and add vertex bin histograms to the sum list
229 fSumList = new TList();
230 fSumList->SetName("Sums");
231 fSumList->SetOwner();
233 if (!fVtxAxis) fVtxAxis = new TAxis(20, -10, 10);
234 fVtxAxis->SetName("VtxAxis");
235 if (!fCentAxis) fCentAxis = new TAxis(20, 0, 100);
236 fVtxAxis->SetName("CentAxis");
238 fHistCent = new TH1D("hCent", "Centralities", 100, 0, 100);
239 fHistVertexSel = new TH1D("hVertexSel", "Selected vertices", fVtxAxis->GetNbins(), fVtxAxis->GetXmin(), fVtxAxis->GetXmax());
240 fHistEventSel = new TH1I("hEventSel", "Event Selection", kOK, 0.5, kOK+0.5);
241 fHistEventSel->GetXaxis()->SetBinLabel(kNoEvent, "No AOD event");
242 fHistEventSel->GetXaxis()->SetBinLabel(kNoForward, "No forward det");
243 fHistEventSel->GetXaxis()->SetBinLabel(kNoCentral, "No central det");
244 fHistEventSel->GetXaxis()->SetBinLabel(kNoTrigger, "Not triggered");
245 fHistEventSel->GetXaxis()->SetBinLabel(kNoCent, "No centrality");
246 fHistEventSel->GetXaxis()->SetBinLabel(kInvCent, "Centrality outside range");
247 fHistEventSel->GetXaxis()->SetBinLabel(kNoVtx, "No vertex");
248 fHistEventSel->GetXaxis()->SetBinLabel(kInvVtx, "Vtx outside range");
249 fHistEventSel->GetXaxis()->SetBinLabel(kOK, "OK!");
251 fHistFMDSPDCorr = new TH2D("hFMDSPDCorr", "hFMDSPCCorr", 200, 0., 20000., 200, 0, 7500);
253 TList* dList = new TList();
254 dList->SetName("Diagnostics");
255 dList->Add(fHistCent);
256 dList->Add(fHistVertexSel);
257 dList->Add(fHistEventSel);
258 dList->Add(fHistFMDSPDCorr);
259 fSumList->Add(dList);
261 fHistdNdedp3Cor = TH2D(Form("hdNdedpCombined_%s", GetQCType(fFlowFlags)), Form("hdNdedpCombined_%s", GetQCType(fFlowFlags)),
262 200, -4., 6., 20, 0., TMath::TwoPi());
263 if ((fFlowFlags & kVZERO)) {
264 Double_t bins[12] = { -6, -3.7, -3.2, -2.7, -2.2, -1.7,
265 2.8, 3.4, 3.9, 4.5, 5.1, 6 };
266 fHistdNdedpV0 = TH2D(Form("hdNdedpv0%s", GetQCType(fFlowFlags)), Form("hdNdedpv0%s", GetQCType(fFlowFlags)),
267 11, -6, 6, 8, 0, TMath::TwoPi());
268 fHistdNdedpV0.GetXaxis()->Set(11, bins);
269 if ((fFlowFlags & k3Cor)) {
270 Double_t bins2[20] = { -6, -3.7, -3.2, -2.7, -2.2, // VZERO
271 -2.0, -1.5, -1.0, -0.5 , 0., 0.5, 1.0, 1.5, 2.0, // SPD
272 2.8, 3.4, 3.9, 4.5, 5.1, 6 }; // VZERO
273 fHistdNdedp3Cor.GetXaxis()->Set(19, bins2);
274 fHistdNdedp3Cor.GetYaxis()->Set(8, 0., TMath::TwoPi());
278 TIter nextForward(&fBinsForward);
280 while ((bin = static_cast<VertexBin*>(nextForward()))) {
281 bin->AddOutput(fSumList, fCentAxis);
283 TIter nextCentral(&fBinsCentral);
284 while ((bin = static_cast<VertexBin*>(nextCentral()))) {
285 bin->AddOutput(fSumList, fCentAxis);
288 //_____________________________________________________________________
289 void AliForwardFlowTaskQC::UserExec(Option_t */*option*/)
292 // Calls the analyze function - called every event
298 // Reset data members
304 PostData(1, fSumList);
308 //_____________________________________________________________________
309 Bool_t AliForwardFlowTaskQC::Analyze()
312 // Load forward and central detector objects from aod tree and call the
313 // cumulants calculation for the correct vertex bin
315 // Return: true on success
319 fAOD = dynamic_cast<AliAODEvent*>(AliForwardUtil::GetAODEvent(this));
321 fHistEventSel->Fill(kNoEvent);
325 // Get detector objects
326 AliAODForwardMult* aodfmult = static_cast<AliAODForwardMult*>(fAOD->FindListObject("Forward"));
327 AliAODCentralMult* aodcmult = static_cast<AliAODCentralMult*>(fAOD->FindListObject("CentralClusters"));
328 AliVVZERO* vzero = GetVZERO();
329 if ((fFlowFlags & kVZERO)) {
331 fHistdNdedpV0.Reset();
332 FillVZEROHist(vzero);
336 // We make sure that the necessary forward object is there
337 if ((fFlowFlags & kFMD) && !aodfmult) {
338 fHistEventSel->Fill(kNoForward);
341 else if ((fFlowFlags & kVZERO) && !vzero) {
342 fHistEventSel->Fill(kNoForward);
345 if (!aodcmult) fHistEventSel->Fill(kNoCentral);
347 // Check event for triggers, get centrality, vtx etc.
348 if (!CheckEvent(aodfmult)) return kFALSE;
349 Int_t vtx = fVtxAxis->FindBin(fVtx)-1;
351 // Then we assign a reference to the forward histogram of interest
352 TH2D& forwarddNdedp = ((fFlowFlags & kFMD) ? aodfmult->GetHistogram() : fHistdNdedpV0);
353 TH2D& spddNdedp = aodcmult->GetHistogram();
354 if ((fFlowFlags & kStdQC)) {
355 FillVtxBinList(fBinsForward, forwarddNdedp, vtx);
356 } else if ((fFlowFlags & kEtaGap)) {
357 FillVtxBinListEtaGap(fBinsForward, forwarddNdedp, forwarddNdedp, vtx);
359 // At the moment only clusters are supported for the central region (some day add tracks?)
360 // So no extra checks necessary
362 if ((fFlowFlags & kStdQC)) {
363 FillVtxBinList(fBinsCentral, spddNdedp, vtx);
364 } else if ((fFlowFlags & kEtaGap)) {
365 FillVtxBinListEtaGap(fBinsCentral, forwarddNdedp, spddNdedp, vtx);
366 } else if ((fFlowFlags & k3Cor)) {
367 FillVtxBinList3Cor(fBinsForward, spddNdedp, forwarddNdedp, vtx);
371 Double_t totForward = forwarddNdedp.Integral(1, forwarddNdedp.GetNbinsX(), 1, forwarddNdedp.GetNbinsY());
372 Double_t totSPD = spddNdedp.Integral(1, spddNdedp.GetNbinsX(), 1, spddNdedp.GetNbinsY());
373 fHistFMDSPDCorr->Fill(totForward, totSPD);
379 //_____________________________________________________________________
380 void AliForwardFlowTaskQC::FillVtxBinList(const TList& list, TH2D& h, Int_t vtx, UShort_t flags) const
383 // Loops over list of VtxBins, fills hists of bins for current vertex
384 // and runs analysis on those bins
387 // list: list of VtxBins
388 // h: dN/detadphi histogram
389 // vtx: current vertex bin
390 // flags: extra flags to handle calculations.
392 // Note: The while loop is used in this function and the next 2 for historical reasons,
393 // as originally each moment had it's own VertexBin object.
396 Int_t nVtxBins = fVtxAxis->GetNbins();
398 while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
400 // If no tracks do things normally
401 if (!(fFlowFlags & kTracks) || (flags & kMC)) {
402 if (!bin->FillHists(h, fCent, kFillBoth|flags|kReset)) continue;
404 // if tracks things are more complicated
405 else if ((fFlowFlags & kTracks)) {
406 if (!FillTracks(bin, kFillRef|kReset|flags)) continue;
407 if (!bin->FillHists(h, fCent, kFillDiff|kReset|flags)) continue;
409 bin->CumulantsAccumulate(fCent);
414 //_____________________________________________________________________
415 void AliForwardFlowTaskQC::FillVtxBinListEtaGap(const TList& list, TH2D& href,
416 TH2D& hdiff, Int_t vtx, UShort_t flags) const
419 // Loops over list of VtxBins, fills hists of bins for current vertex
420 // and runs analysis on those bins
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.
431 Int_t nVtxBins = fVtxAxis->GetNbins();
433 while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
435 if (!(fFlowFlags & kTracks) || (flags & kMC)) {
436 if(!bin->FillHists(href, fCent, kFillRef|flags|kReset)) continue;
438 else if ((fFlowFlags & kTracks)) {
439 if (!FillTracks(bin, kFillRef|kReset|flags)) continue;
441 if (!bin->FillHists(hdiff, fCent, kFillDiff|kReset|flags)) continue;
442 bin->CumulantsAccumulate(fCent);
447 //_____________________________________________________________________
448 void AliForwardFlowTaskQC::FillVtxBinList3Cor(const TList& list, TH2D& hcent,
449 TH2D& hfwd, Int_t vtx, UShort_t flags)
452 // Loops over list of VtxBins, fills hists of bins for current vertex
453 // and runs analysis on those bins
456 // list: list of VtxBins
457 // hcent: dN/detadphi histogram for central coverage
458 // hfwd: dN/detadphi histogram for forward coverage
459 // vBin: current vertex bin
460 // flags: extra flags to handle calculations.
464 Int_t nVtxBins = fVtxAxis->GetNbins();
466 TH2D& h = CombineHists(hcent, hfwd);
468 while ((bin = static_cast<VertexBin*>(list.At(vtx+(nVtxBins*i))))) {
470 if (!bin->FillHists(h, fCent, kFillBoth|flags|kReset)) continue;
471 bin->CumulantsAccumulate3Cor(fCent);
476 //_____________________________________________________________________
477 TH2D& AliForwardFlowTaskQC::CombineHists(TH2D& hcent, TH2D& hfwd)
480 // Combines a forward and central d^2N/detadphi histogram.
481 // At some point it might need a flag to choose which histogram gets
482 // priority when there is an overlap, at the moment the average is chosen
485 // hcent: Central barrel detector
486 // hfwd: Forward detector
488 // Return: reference to combined hist
491 // If hists are the same (MC input) don't do anything
492 if (&hcent == &hfwd) return hcent;
494 fHistdNdedp3Cor.Reset();
496 if ((fFlowFlags & kFMD)) {
497 for (Int_t e = 1; e <= fHistdNdedp3Cor.GetNbinsX(); e++) {
498 Double_t eta = fHistdNdedp3Cor.GetXaxis()->GetBinCenter(e);
499 Bool_t fwdCov = (hfwd.GetBinContent(e, 0) != 0);
500 Bool_t centCov = (hcent.GetBinContent(e, 0) != 0);
501 if (!fwdCov && !centCov) continue;
502 else fHistdNdedp3Cor.SetBinContent(e, 0, 1);
503 for (Int_t p = 1; p <= fHistdNdedp3Cor.GetNbinsY(); p++) {
504 Double_t phi = fHistdNdedp3Cor.GetYaxis()->GetBinCenter(p);
508 cont += hfwd.GetBinContent(e, p);
512 cont += hcent.GetBinContent(e, p);
515 if (cont == 0 || n == 0) continue;
517 fHistdNdedp3Cor.Fill(eta, phi, cont);
520 // VZERO, SPD input, here we do not average but cut to avoid
521 // (too much) overlap.
522 } else if ((fFlowFlags & kVZERO)) {
524 for (Int_t eV = 1; eV <= hfwd.GetNbinsX(); eV++) {
525 Double_t eta = hfwd.GetXaxis()->GetBinLowEdge(eV)+0.1;
526 if (hfwd.GetBinContent(eV, 0) == 0) continue;
528 Int_t he = fHistdNdedp3Cor.GetXaxis()->FindBin(eta);
529 fHistdNdedp3Cor.SetBinContent(he, 0, 1);
531 for (Int_t p = 1; p <= hfwd.GetNbinsY(); p++) {
532 Double_t phi = hfwd.GetYaxis()->GetBinCenter(p);
533 Double_t cont = hfwd.GetBinContent(eV, p);
534 fHistdNdedp3Cor.Fill(eta, phi, cont);
538 Int_t eSs = hcent.GetXaxis()->FindBin(-1.99);
539 Int_t eSe = hcent.GetXaxis()->FindBin(1.99);
540 for (Int_t eS = eSs; eS <= eSe; eS++) {
541 Double_t eta = hcent.GetXaxis()->GetBinCenter(eS);
542 if (hcent.GetBinContent(eS, 0) == 0) continue;
544 Int_t he = fHistdNdedp3Cor.GetXaxis()->FindBin(eta);
545 fHistdNdedp3Cor.SetBinContent(he, 0, 1);
547 for (Int_t p = 1; p <= hcent.GetNbinsY(); p++) {
548 Double_t phi = hcent.GetYaxis()->GetBinCenter(p);
549 Double_t cont = hcent.GetBinContent(eS, p);
550 fHistdNdedp3Cor.Fill(eta, phi, cont);
554 return fHistdNdedp3Cor;
556 //_____________________________________________________________________
557 Bool_t AliForwardFlowTaskQC::FillTracks(VertexBin* bin, UShort_t mode) const
560 // Get TPC tracks to use for reference flow.
562 // Return: TObjArray with tracks
564 TObjArray* trList = 0;
565 AliESDEvent* esdEv = 0;
566 if (AliForwardUtil::CheckForAOD() == 1) // AOD tracks
567 trList = static_cast<TObjArray*>(fAOD->GetTracks());
569 esdEv = dynamic_cast<AliESDEvent*>(InputEvent());
571 Bool_t useEvent = bin->FillTracks(trList, esdEv, fTrackCuts, mode);
574 //_____________________________________________________________________
575 void AliForwardFlowTaskQC::Terminate(Option_t */*option*/)
578 // Calls the finalize function, done at the end of the analysis
584 // Make sure pointers are set to the correct lists
585 fSumList = dynamic_cast<TList*> (GetOutputData(1));
587 AliError("Could not retrieve TList fSumList");
591 fOutputList = new TList();
592 fOutputList->SetName("Results");
593 fOutputList->SetOwner();
595 if ((fFlowFlags & kEtaGap) || (fFlowFlags & kTracks)) {
596 TParameter<Double_t>* etaGap = new TParameter<Double_t>("EtaGap", fEtaGap);
597 fOutputList->Add(etaGap);
599 // We only add axes in terminate, as TAxis object do not merge well,
600 // and so we get a mess when running on the grid if we put them in the sum list...
601 fVtxAxis->SetName("VtxAxis");
602 fOutputList->Add(fVtxAxis);
603 fCentAxis->SetName("CentAxis");
604 fOutputList->Add(fCentAxis);
606 // Run finalize on VertexBins
609 // Loop over output to get dN/deta hists - used for diagnostics
610 TIter next(fOutputList);
615 while ((o = next())) {
617 if (name.Contains("dNdeta")) {
618 dNdeta = dynamic_cast<TH2D*>(o);
619 name.ReplaceAll("dNdeta", "cent");
620 name.ReplaceAll("Ref", "");
621 name.ReplaceAll("Diff", "");
622 cent = dynamic_cast<TH1D*>(fOutputList->FindObject(name.Data()));
623 if (!dNdeta || !cent) continue;
624 for (Int_t cBin = 1; cBin <= dNdeta->GetNbinsY(); cBin++) {
625 Double_t nEvents = cent->GetBinContent(cBin);
626 if (nEvents == 0) continue;
627 for (Int_t eBin = 1; eBin <= dNdeta->GetNbinsX(); eBin++) {
628 dNdeta->SetBinContent(eBin, cBin, dNdeta->GetBinContent(eBin, cBin)/nEvents);
629 dNdeta->SetBinError(eBin, cBin, dNdeta->GetBinError(eBin, cBin)/nEvents);
635 // Loop over output and make 1D projections for fast look at results
636 MakeCentralityHists(fOutputList);
637 TList* vtxList = (TList*)fOutputList->FindObject("vtxList");
638 if (vtxList) MakeCentralityHists(vtxList);
639 TList* nuaList = (TList*)fOutputList->FindObject("NUATerms");
640 TIter nextNua(nuaList);
643 while ((o = nextNua())) {
644 if (!(h = dynamic_cast<TH2D*>(o))) continue;
645 Double_t evts = h->GetBinContent(0, 0);
646 if (evts != 0) h->Scale(1./evts);
648 if (nuaList) MakeCentralityHists(nuaList);
650 PostData(2, fOutputList);
654 //_____________________________________________________________________
655 void AliForwardFlowTaskQC::Finalize()
658 // Finalize command, called by Terminate()
661 // Reinitiate vertex bins if Terminate is called separately!
662 if (fBinsForward.GetEntries() == 0) InitVertexBins();
664 // Iterate over all vertex bins objects and finalize cumulants
666 EndVtxBinList(fBinsForward);
667 EndVtxBinList(fBinsCentral);
671 //_____________________________________________________________________
672 void AliForwardFlowTaskQC::EndVtxBinList(const TList& list) const
675 // Loop over VertexBin list and call terminate on each
678 // list: VertexBin list
682 while ((bin = static_cast<VertexBin*>(next()))) {
683 bin->CumulantsTerminate(fSumList, fOutputList);
687 // _____________________________________________________________________
688 void AliForwardFlowTaskQC::MakeCentralityHists(TList* list) const
691 // Loop over a list containing a TH2D with flow results
692 // and project to TH1's in specific centrality bins
695 // list: Flow results list
701 TIter nextHist(list);
702 while ((helper = dynamic_cast<TObject*>(nextHist()))) {
703 if (!(hist2D = dynamic_cast<TH2D*>(helper))) continue;
704 for (Int_t cBin = 1; cBin <= hist2D->GetNbinsY(); cBin++) {
705 Int_t cMin = Int_t(hist2D->GetYaxis()->GetBinLowEdge(cBin));
706 Int_t cMax = Int_t(hist2D->GetYaxis()->GetBinUpEdge(cBin));
707 TString name = Form("cent_%d-%d", cMin, cMax);
708 centList = (TList*)list->FindObject(name.Data());
710 centList = new TList();
711 centList->SetName(name.Data());
714 hist1D = hist2D->ProjectionX(Form("%s_%s", hist2D->GetName(), name.Data()),
716 hist1D->SetTitle(hist1D->GetName());
717 centList->Add(hist1D);
721 // _____________________________________________________________________
722 Bool_t AliForwardFlowTaskQC::CheckEvent(const AliAODForwardMult* aodfm)
725 // Function to check that an AOD event meets the cuts
728 // AliAODForwardMult: forward mult object with trigger and vertex info
730 // Return: false if there is no trigger or if the centrality or vertex
731 // is out of range. Otherwise true.
734 // First check for trigger
735 if (!CheckTrigger(aodfm)) {
736 fHistEventSel->Fill(kNoTrigger);
740 // Then check for centrality
741 if (!GetCentrality(aodfm)) {
745 // And finally check for vertex
746 if (!GetVertex(aodfm)) {
750 // Ev. accepted - filling diag. hists
751 fHistCent->Fill(fCent);
752 fHistVertexSel->Fill(fVtx);
753 fHistEventSel->Fill(kOK);
757 // _____________________________________________________________________
758 Bool_t AliForwardFlowTaskQC::CheckTrigger(const AliAODForwardMult* aodfm) const
761 // Function to look for a trigger string in the event.
762 // First check for info in forward mult object, if not there, use the AOD header
765 // AliAODForwardMult: forward mult object with trigger and vertex info
767 // Return: true if offline trigger is present
769 if (aodfm) return aodfm->IsTriggerBits(AliAODForwardMult::kOffline);
770 // this may need to be changed for 2011 data to handle kCentral and so on...
771 else return (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
772 ->IsEventSelected() & AliVEvent::kMB);
774 // _____________________________________________________________________
775 Bool_t AliForwardFlowTaskQC::GetCentrality(const AliAODForwardMult* aodfm)
778 // Function to look get centrality of the event.
779 // First check for info in forward mult object, if not there, use the AOD header
782 // AliAODForwardMult: forward mult object with trigger and vertex info
784 // Return: true if centrality determination is present
787 if (aodfm->HasCentrality()) {
788 fCent = (Double_t)aodfm->GetCentrality();
789 if (fCentAxis->GetXmin() > fCent || fCent >= fCentAxis->GetXmax()) {
790 fHistEventSel->Fill(kInvCent);
796 fHistEventSel->Fill(kNoCent);
800 AliCentrality* aodCent = fAOD->GetCentrality();
802 fCent = (Double_t)aodCent->GetCentralityPercentile("V0M");
803 if (fCentAxis->GetXmin() > fCent || fCent >= fCentAxis->GetXmax()) {
804 fHistEventSel->Fill(kInvCent);
810 fHistEventSel->Fill(kNoCent);
815 //_____________________________________________________________________
816 Bool_t AliForwardFlowTaskQC::GetVertex(const AliAODForwardMult* aodfm)
819 // Function to look for vertex determination in the event.
820 // First check for info in forward mult object, if not there, use the AOD header
823 // AliAODForwardMult: forward mult object with trigger and vertex info
825 // Return: true if vertex is determined
828 if (aodfm->HasIpZ()) {
829 fVtx = aodfm->GetIpZ();
830 if (fVtx < fVtxAxis->GetXmin() || fVtx >= fVtxAxis->GetXmax()) {
831 fHistEventSel->Fill(kInvVtx);
836 fHistEventSel->Fill(kNoVtx);
841 AliAODVertex* aodVtx = fAOD->GetPrimaryVertex();
843 fVtx = aodVtx->GetZ();
844 if (fVtx < fVtxAxis->GetXmin() || fVtx >= fVtxAxis->GetXmax()) {
845 fHistEventSel->Fill(kInvVtx);
850 fHistEventSel->Fill(kNoVtx);
856 // _____________________________________________________________________
857 AliVVZERO* AliForwardFlowTaskQC::GetVZERO() const
860 // Get VZERO object from ESD or AOD
862 // Return: VZERO data object
864 AliVVZERO* vzero = 0;
866 UShort_t input = AliForwardUtil::CheckForAOD();
868 // If AOD input, simply get the track array from the event
869 case 1: vzero = (AliVVZERO*)fAOD->GetVZEROData();
872 // If ESD input get event, apply track cuts
873 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
875 vzero = (AliVVZERO*)esd->GetVZEROData();
878 default: AliFatal("Neither ESD or AOD input. This should never happen");
883 // _____________________________________________________________________
884 void AliForwardFlowTaskQC::FillVZEROHist(AliVVZERO* vzero)
887 // Loops over VZERO data object and fill up d^2N/detadphi histogram for flow analysis
890 // vzero: VZERO AOD data object
895 // Sort of right for 2010 data, do not use for 2011!
896 Double_t eq[64] = { 1.43536, 1.45727, 1.44993, 1.30051, 1.17425, 1.2335, 1.22247, 1.14362,
897 1.14647, 1.25208, 1.17681, 1.21642, 1.16604, 1.05532, 1.03212, 1.1032,
898 1.22941, 1.36986, 1.14652, 1.20056, 0.927086, 1.10809, 1.03343, 1.29472,
899 1.21204, 1.29217, 1.2003, 2.10382, 1.28513, 1.40558, 1.25784, 1.21848,
900 0.475162, 0.50421, 0.503617, 0.512471, 0.515276, 0.39831, 0.415199, 0.444664,
901 0.521922, 0.785915, 0.703658, 0.832479, 0.77461, 0.73129, 0.778697, 0.710265,
902 0.89686, 0.967688, 0.974225, 0.873445, 0.811096, 0.828493, 0.889609, 0.586056,
903 1.15877, 0.954656, 0.914557, 0.979028, 1.04907, 0.748518, 0.928043, 0.98175 };
904 for (Int_t i = 0; i < 64; i++) {
907 bin = (ring < 5 ? ring+1 : 15-ring);
908 eta = fHistdNdedpV0.GetXaxis()->GetBinCenter(bin);
909 fHistdNdedpV0.SetBinContent(bin, 0, 1);
911 Float_t amp = vzero->GetMultiplicity(i);
913 Double_t phi = TMath::Pi()/8.+TMath::TwoPi()*i/8.;
914 while (phi > TMath::TwoPi()) phi -= TMath::TwoPi();
915 fHistdNdedpV0.Fill(eta, phi, amp);
918 //_____________________________________________________________________
919 AliForwardFlowTaskQC::VertexBin::VertexBin()
921 fMaxMoment(0), // Max flow moment for this vertexbin
922 fVzMin(0), // Vertex z-coordinate min [cm]
923 fVzMax(0), // Vertex z-coordinate max [cm]
924 fType(), // Data type name e.g., FMD/SPD/FMDTR/SPDTR/MC
925 fFlags(0), // Flow flags
926 fSigmaCut(-1), // Sigma cut to remove outlier events
927 fEtaGap(-1), // Eta gap value
928 fEtaLims(), // Limits for binning in 3Cor method
929 fCumuRef(), // Histogram for reference flow
930 fCumuDiff(), // Histogram for differential flow
931 fCumuHists(0,0), // CumuHists object for keeping track of results
932 fCumuNUARef(), // Histogram for ref NUA terms
933 fCumuNUADiff(), // Histogram for diff NUA terms
934 fdNdedpRefAcc(), // Diagnostics histogram for acc. maps
935 fdNdedpDiffAcc(),// Diagnostics histogram for acc. maps
936 fOutliers(), // Histogram for sigma distribution
937 fDebug() // Debug level
940 // Default constructor
943 //_____________________________________________________________________
944 AliForwardFlowTaskQC::VertexBin::VertexBin(Int_t vLow, Int_t vHigh,
945 UShort_t moment, TString name,
946 UShort_t flags, Double_t cut,
949 fMaxMoment(moment), // Max flow moment for this vertexbin
950 fVzMin(vLow), // Vertex z-coordinate min [cm]
951 fVzMax(vHigh), // Vertex z-coordinate max [cm]
952 fType(name), // Data type name e.g., FMD/SPD/FMDTR/SPDTR/MC
953 fFlags(flags), // Flow flags
954 fSigmaCut(cut), // Sigma cut to remove outlier events
955 fEtaGap(etaGap), // Eta gap value
956 fEtaLims(), // Limits for binning in 3Cor method
957 fCumuRef(), // Histogram for reference flow
958 fCumuDiff(), // Histogram for differential flow
959 fCumuHists(moment,0),// CumuHists object for keeping track of results
960 fCumuNUARef(), // Histogram for ref NUA terms
961 fCumuNUADiff(), // Histogram for diff NUA terms
962 fdNdedpRefAcc(), // Diagnostics histogram for acc. maps
963 fdNdedpDiffAcc(), // Diagnostics histogram for acc. maps
964 fOutliers(), // Histogram for sigma distribution
965 fDebug(0) // Debug level
971 // vLow: min z-coordinate
972 // vHigh: max z-coordinate
973 // moment: max flow moment
974 // name: data type name (FMD/SPD/FMDTR/SPDTR/MC)
977 // etaGap: eta-gap value
981 SetName(Form("%svertexBin%d_%d_%d%s", fType.Data(), moment, vLow, vHigh, GetQCType(fFlags)));
982 SetTitle(Form("%svertexBin%d_%d_%d%s", fType.Data(), moment, vLow, vHigh, GetQCType(fFlags)));
984 fDebug = AliAnalysisManager::GetAnalysisManager()->GetDebugLevel();
985 if (fDebug > 0) Printf("AliForwardFlowTaskQC::VertexBin()\tDebugMode: %d", fDebug);
987 // Set limits for 3 correlator method
988 if ((fFlags & kFMD)) {
995 } else if ((fFlags & kVZERO)) {
1004 //_____________________________________________________________________
1005 AliForwardFlowTaskQC::VertexBin&
1006 AliForwardFlowTaskQC::VertexBin::operator=(const AliForwardFlowTaskQC::VertexBin& o)
1009 // Assignment operator
1012 // o: AliForwardFlowTaskQC::VertexBin
1014 if (&o == this) return *this;
1015 fMaxMoment = o.fMaxMoment;
1020 fSigmaCut = o.fSigmaCut;
1021 fEtaGap = o.fEtaGap;
1022 fCumuRef = o.fCumuRef;
1023 fCumuDiff = o.fCumuDiff;
1024 fCumuHists = o.fCumuHists;
1025 fCumuNUARef = o.fCumuNUARef;
1026 fCumuNUADiff = o.fCumuNUADiff;
1027 fdNdedpRefAcc = o.fdNdedpRefAcc;
1028 fdNdedpDiffAcc = o.fdNdedpDiffAcc;
1029 fOutliers = o.fOutliers;
1031 for (UInt_t i = 0; i < sizeof(fEtaLims)/sizeof(Double_t); i++) fEtaLims[i] = o.fEtaLims[i];
1035 //_____________________________________________________________________
1036 void AliForwardFlowTaskQC::VertexBin::AddOutput(TList* outputlist, TAxis* centAxis)
1039 // Add histograms to outputlist
1042 // outputlist: list of histograms
1043 // centAxis: centrality axis
1046 // First we try to find an outputlist for this vertexbin
1047 TList* list = (TList*)outputlist->FindObject(Form("%svertex_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1048 // If it doesn't exist we make one
1051 list->SetName(Form("%svertex_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1052 outputlist->Add(list);
1055 // Get bin numbers and binning defined
1056 Int_t nHBins = GetBinNumberSin();
1057 Int_t nEtaBins = 48;
1058 if ((fFlags & k3Cor)) {
1059 if ((fFlags & kFMD)) nEtaBins = 24;
1060 else if ((fFlags & kVZERO)) nEtaBins = 19;
1062 else if ((fFlags & kVZERO) && (fFlags & kEtaGap)) nEtaBins = 19;
1063 else if ((fFlags & kVZERO)) nEtaBins = 11;
1065 Double_t vzeroBins[12] = { -6, -3.7, -3.2, -2.7, -2.2, -1.7,
1066 2.8, 3.4, 3.9, 4.5, 5.1, 6 };
1067 Double_t vzeroBins2[20] = { -6, -3.7, -3.2, -2.7, -2.2, // VZERO
1068 -2.0, -1.5, -1.0, -0.5 , 0., 0.5, 1.0, 1.5, 2.0, // SPD
1069 2.8, 3.4, 3.9, 4.5, 5.1, 6 }; // VZERO
1071 Int_t nRefBins = nEtaBins; // needs to be something as default
1072 if ((fFlags & kStdQC)) {
1073 if ((fFlags & kSymEta) && !((fFlags & kTracks) && (fFlags & kSPD))) nRefBins = 1;
1075 } else if ((fFlags & kEtaGap )) {
1077 } else if ((fFlags & k3Cor)) {
1081 // We initiate the reference histogram
1082 fCumuRef = new TH2D(Form("%s_%d_%d%s_ref", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1083 Form("%s_%d_%d%s_ref", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1085 nHBins, 0.5, nHBins+0.5);
1086 if ((fFlags & kVZERO) && (fFlags & k3Cor)) fCumuRef->GetXaxis()->Set(nEtaBins, vzeroBins2);
1087 SetupNUALabels(fCumuRef->GetYaxis());
1088 list->Add(fCumuRef);
1089 // We initiate the differential histogram
1090 fCumuDiff = new TH2D(Form("%s_%d_%d%s_diff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1091 Form("%s_%d_%d%s_diff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1092 nEtaBins, -6., 6., nHBins, 0.5, nHBins+0.5);
1093 if ((fFlags & kVZERO)) {
1094 if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1095 fCumuDiff->GetXaxis()->Set(nEtaBins, vzeroBins2);
1096 else fCumuDiff->GetXaxis()->Set(nEtaBins, vzeroBins);
1098 SetupNUALabels(fCumuDiff->GetYaxis());
1099 list->Add(fCumuDiff);
1101 // Cumulants sum hists
1102 Int_t cBins = centAxis->GetNbins();
1103 fCumuHists.ConnectList(Form("%sCumu_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)), list);
1105 Int_t nC2Bins = ((fFlags & kEtaGap) || (fFlags & k3Cor) ? kW2 : k3pWeight);
1106 Int_t nC4Bins = ((fFlags & kEtaGap) ? kW2 : ((fFlags & k3Cor) ? kW4 : kSinphi1phi2phi3p));
1107 for (Int_t n = 2; n <= fMaxMoment; n++) {
1108 // Initiate the ref cumulant sum histogram
1109 cumuHist = new TH3D(Form("%sv%d_vertex_%d_%d%s_cumuRef", fType.Data(), n, fVzMin, fVzMax, GetQCType(fFlags)),
1110 Form("%sv%d_vertex_%d_%d%s_cumuRef", fType.Data(), n, fVzMin, fVzMax, GetQCType(fFlags)),
1112 cBins, 0., 100., nC2Bins, 0.5, nC2Bins+0.5);
1113 if ((fFlags & kVZERO) && (fFlags & k3Cor)) cumuHist->GetXaxis()->Set(nEtaBins, vzeroBins2);
1114 cumuHist->GetYaxis()->Set(cBins, centAxis->GetXbins()->GetArray());
1115 fCumuHists.Add(cumuHist);
1116 // Initiate the diff cumulant sum histogram
1117 cumuHist = new TH3D(Form("%sv%d_vertex_%d_%d%s_cumuDiff", fType.Data(), n, fVzMin, fVzMax, GetQCType(fFlags)),
1118 Form("%sv%d_vertex_%d_%d%s_cumuDiff", fType.Data(), n, fVzMin, fVzMax, GetQCType(fFlags)),
1119 nEtaBins, -6., 6., cBins, 0., 100., nC4Bins, 0.5, nC4Bins+0.5);
1120 if ((fFlags & kVZERO)) {
1121 if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1122 cumuHist->GetXaxis()->Set(nEtaBins, vzeroBins2);
1123 else cumuHist->GetXaxis()->Set(nEtaBins, vzeroBins);
1125 cumuHist->GetYaxis()->Set(cBins, centAxis->GetXbins()->GetArray());
1126 fCumuHists.Add(cumuHist);
1129 // Common NUA histograms
1130 fCumuNUARef = new TH3D(Form("%s_vertex_%d_%d%s_cumuNUARef", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1131 Form("%s_vertex_%d_%d%s_cumuNUARef", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1133 cBins, 0., 100., nHBins, 0.5, nHBins+0.5);
1134 if ((fFlags & kVZERO) && (fFlags & k3Cor)) fCumuNUARef->GetXaxis()->Set(nEtaBins, vzeroBins2);
1135 fCumuNUARef->GetYaxis()->Set(cBins, centAxis->GetXbins()->GetArray());
1136 SetupNUALabels(fCumuNUARef->GetZaxis());
1137 fCumuNUARef->Sumw2();
1138 list->Add(fCumuNUARef);
1140 fCumuNUADiff = new TH3D(Form("%s_vertex_%d_%d%s_cumuNUADiff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1141 Form("%s_vertex_%d_%d%s_cumuNUADiff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1142 nEtaBins, -6., 6., cBins, 0., 100., nHBins, 0.5, nHBins+0.5);
1143 if ((fFlags & kVZERO)) {
1144 if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1145 fCumuNUADiff->GetXaxis()->Set(nEtaBins, vzeroBins2);
1146 else fCumuNUADiff->GetXaxis()->Set(nEtaBins, vzeroBins);
1148 fCumuNUADiff->GetYaxis()->Set(cBins, centAxis->GetXbins()->GetArray());
1149 SetupNUALabels(fCumuNUADiff->GetZaxis());
1150 fCumuNUADiff->Sumw2();
1151 list->Add(fCumuNUADiff);
1153 // We create diagnostic histograms.
1154 TList* dList = (TList*)outputlist->FindObject("Diagnostics");
1155 if (!dList) AliFatal("No diagnostics list found");
1158 Int_t nPhiBins = ((fFlags & kFMD) ? 20 : 8);
1159 fdNdedpRefAcc = new TH2F(Form("h%sdNdedpRefAcc_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1160 Form("%s reference flow acceptance map for %d cm < v_{z} < %d cm", fType.Data(), fVzMin, fVzMax),
1161 nEtaBins, -6., 6., nPhiBins, 0, TMath::TwoPi());
1162 if ((fFlags & kVZERO)) {
1163 if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1164 fdNdedpRefAcc->GetXaxis()->Set(nEtaBins, vzeroBins2);
1165 else fdNdedpRefAcc->GetXaxis()->Set(nEtaBins, vzeroBins);
1167 dList->Add(fdNdedpRefAcc);
1169 fdNdedpDiffAcc = new TH2F(Form("h%sdNdedpDiffAcc_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1170 Form("%s differential flow acceptance map for %d cm < v_{z} < %d cm", fType.Data(), fVzMin, fVzMax),
1171 nEtaBins, -6., 6., nPhiBins, 0, TMath::TwoPi());
1172 if ((fFlags & kVZERO)) {
1173 if ((fFlags & kSPD) || (fFlags & k3Cor) || (fFlags & kEtaGap))
1174 fdNdedpDiffAcc->GetXaxis()->Set(nEtaBins, vzeroBins2);
1175 else fdNdedpDiffAcc->GetXaxis()->Set(nEtaBins, vzeroBins);
1177 dList->Add(fdNdedpDiffAcc);
1179 fOutliers = new TH2F(Form("hOutliers_%s_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)),
1180 Form("Maximum #sigma from mean N_{ch} pr. bin - %s, %d < v_{z} < %d",
1181 fType.Data(), fVzMin, fVzMax),
1182 20, 0., 100., 500, 0., ((fFlags & kMC) ? 15. : 5.));
1183 dList->Add(fOutliers);
1187 //_____________________________________________________________________
1188 Bool_t AliForwardFlowTaskQC::VertexBin::FillHists(TH2D& dNdetadphi, Double_t cent, UShort_t mode)
1191 // Fill reference and differential eta-histograms
1194 // dNdetadphi: 2D histogram with input data
1196 // mode: filling mode: kFillRef/kFillDiff/kFillBoth
1198 if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1199 Bool_t useEvent = kTRUE;
1201 // Fist we reset histograms
1202 if ((mode & kReset)) {
1203 if ((mode & kFillRef)) fCumuRef->Reset();
1204 if ((mode & kFillDiff)) fCumuDiff->Reset();
1206 // Then we loop over the input and calculate sum cos(k*n*phi)
1207 // and fill it in the reference and differential histograms
1209 Double_t limit = 9999.;
1210 for (Int_t etaBin = 1; etaBin <= dNdetadphi.GetNbinsX(); etaBin++) {
1211 Double_t eta = dNdetadphi.GetXaxis()->GetBinCenter(etaBin);
1212 // Numbers to cut away bad events
1213 Double_t runAvg = 0;
1216 Double_t avgSqr = 0;
1217 for (Int_t phiBin = 0; phiBin <= dNdetadphi.GetNbinsY(); phiBin++) {
1218 // Check for acceptance
1220 if (dNdetadphi.GetBinContent(etaBin, 0) == 0) break;
1221 // Central limit for eta gap break for reference flow
1222 if ((fFlags & kEtaGap) && (mode & kFillRef) &&
1223 TMath::Abs(eta) < fEtaGap) break;
1224 // Backward and forward eta gap break for reference flow
1225 if ((fFlags & kEtaGap) && (mode & kFillRef) && TMath::Abs(eta) > TMath::Abs(limit)) break;
1226 if ((fFlags & kStdQC) && (fFlags & kMC) && !(fFlags & kTracks)) {
1227 if (!(fFlags & kSPD) && TMath::Abs(eta) < 1.75) break;
1228 if ((fFlags & kSPD) && TMath::Abs(eta) > 2.00) break;
1230 if (limit > 1e3) limit = dNdetadphi.GetXaxis()->GetBinLowEdge(etaBin);
1232 } // End of phiBin == 0
1233 Double_t phi = dNdetadphi.GetYaxis()->GetBinCenter(phiBin);
1234 Double_t weight = dNdetadphi.GetBinContent(etaBin, phiBin);
1236 // We calculate the average Nch per. bin
1237 avgSqr += weight*weight;
1240 if (weight == 0) continue;
1241 if (weight > max) max = weight;
1242 // Fill into Cos() and Sin() hists
1243 if ((mode & kFillRef) && !((fFlags & kTracks) && (fFlags & kMC) && TMath::Abs(eta) > 0.75)) {
1244 fCumuRef->Fill(eta, 0., weight);// mult goes in underflowbin - no visual, but not needed?
1245 fdNdedpRefAcc->Fill(eta, phi, weight);
1247 if ((mode & kFillDiff)) {
1248 fCumuDiff->Fill(eta, 0., weight);
1249 fdNdedpDiffAcc->Fill(eta, phi, weight);
1251 for (Int_t n = 1; n <= 2*fMaxMoment; n++) {
1252 Double_t cosBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberCos(n));
1253 Double_t sinBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberSin(n));
1254 Double_t cosnPhi = weight*TMath::Cos(n*phi);
1255 Double_t sinnPhi = weight*TMath::Sin(n*phi);
1257 if ((mode & kFillRef) && !((fFlags & kTracks) && (fFlags & kMC) && TMath::Abs(eta) > 0.75)) {
1258 fCumuRef->Fill(eta, cosBin, cosnPhi);
1259 fCumuRef->Fill(eta, sinBin, sinnPhi);
1262 if ((mode & kFillDiff)) {
1263 fCumuDiff->Fill(eta, cosBin, cosnPhi);
1264 fCumuDiff->Fill(eta, sinBin, sinnPhi);
1266 } // End of NUA loop
1267 } // End of phi loop
1268 // Outlier cut calculations
1272 Double_t stdev = (nInAvg > 1 ? TMath::Sqrt(nInAvg/(nInAvg-1))*TMath::Sqrt(avgSqr - runAvg*runAvg) : 0);
1273 Double_t nSigma = (stdev == 0 ? 0 : (max-runAvg)/stdev);
1274 if (fSigmaCut > 0. && nSigma >= fSigmaCut && cent < 60) nBadBins++;
1276 fOutliers->Fill(cent, nSigma);
1277 // We still finish the loop, for fOutliers to make sense,
1278 // but we do no keep the event for analysis
1279 if (nBadBins > 3) useEvent = kFALSE;
1285 //_____________________________________________________________________
1286 Bool_t AliForwardFlowTaskQC::VertexBin::FillTracks(TObjArray* trList, AliESDEvent* esd,
1287 AliAnalysisFilter* trFilter, UShort_t mode)
1290 // Fill reference and differential eta-histograms
1293 // trList: list of tracks
1294 // mode: filling mode: kFillRef/kFillDiff/kFillBoth
1296 if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1297 if (!trList && !esd) {
1298 AliError("FillTracks: No AOD track list or ESD event - something might be wrong!");
1302 // Fist we reset histograms
1303 if ((mode & kReset)) {
1304 if ((mode & kFillRef)) fCumuRef->Reset();
1305 if ((mode & kFillDiff)) fCumuDiff->Reset();
1308 // Then we loop over the input and calculate sum cos(k*n*phi)
1309 // and fill it in the reference and differential histograms
1311 if (trList) nTr = trList->GetEntries();
1312 if (esd) nTr = esd->GetNumberOfTracks();
1313 if (nTr == 0) return kFALSE;
1315 // Cuts for AOD tracks (have already been applied to ESD tracks) - except dEdx
1316 // const tpcdEdx = 10;
1317 for (Int_t i = 0; i < nTr; i++) { // track loop
1318 tr = (trList ? (AliVTrack*)trList->At(i) : (AliVTrack*)esd->GetTrack(i));
1321 AliESDtrack* esdTr = (AliESDtrack*)tr;
1322 if (!trFilter->IsSelected(esdTr)) continue;
1324 else if (trList) { // If AOD input
1325 Double_t pTMin = 0, pTMax = 0, etaMin = 0, etaMax = 0, minNCl = 0, bit = 0;
1326 if ((fFlags & kTPC) == kTPC) pTMin = 0.2, pTMax = 5., etaMin = -0.8, etaMax = 0.8, minNCl = 70, bit = 128;
1327 if ((fFlags & kHybrid) == kHybrid) pTMin = 0.2, pTMax = 5., etaMin = -0.8, etaMax = 0.8, minNCl = 70, bit = 272;
1329 AliAODTrack* aodTr = (AliAODTrack*)tr;
1330 if (aodTr->GetID() > -1) continue;
1331 if (!aodTr->TestFilterBit(bit) || aodTr->Pt() > pTMax || aodTr->Pt() < pTMin ||
1332 aodTr->Eta() > etaMax || aodTr->Eta() < etaMin || aodTr->GetTPCNcls() < minNCl) continue;
1335 // if (tr->GetTPCsignal() < tpcdEdx) continue;
1337 Double_t eta = tr->Eta();
1338 if (((fFlags & kSPD) || (fFlags & kEtaGap)) && TMath::Abs(eta) < fEtaGap) continue;
1339 Double_t phi = tr->Phi();
1340 // Double_t pT = tr->Pt();
1341 // AliAODMCHeader* head = static_cast<AliAODMCHeader*>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
1342 // Double_t rp = head->GetReactionPlaneAngle();
1343 // Double_t b = head->GetImpactParameter();
1344 Double_t weight = 1.;//AliForwardFlowWeights::Instance().CalcWeight(eta, pT, phi, 0, rp, b);
1346 if ((mode & kFillRef)) {
1347 fCumuRef->Fill(eta, 0., weight);// mult goes in underflowbin - no visual, but not needed?
1348 fdNdedpRefAcc->Fill(eta, phi, weight);
1350 if ((mode & kFillDiff)) {
1351 fCumuDiff->Fill(eta, 0., weight);
1352 fdNdedpDiffAcc->Fill(eta, phi, weight);
1354 for (Int_t n = 1; n <= 2*fMaxMoment; n++) {
1355 Double_t cosBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberCos(n));
1356 Double_t sinBin = fCumuDiff->GetYaxis()->GetBinCenter(GetBinNumberSin(n));
1357 Double_t cosnPhi = weight*TMath::Cos(n*phi);
1358 Double_t sinnPhi = weight*TMath::Sin(n*phi);
1360 if ((mode & kFillRef)) {
1361 fCumuRef->Fill(eta, cosBin, cosnPhi);
1362 fCumuRef->Fill(eta, sinBin, sinnPhi);
1365 if ((mode & kFillDiff)) {
1366 fCumuDiff->Fill(eta, cosBin, cosnPhi);
1367 fCumuDiff->Fill(eta, sinBin, sinnPhi);
1369 } // End of NUA loop
1370 } // End of track loop
1373 //_____________________________________________________________________
1374 void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate(Double_t cent)
1377 // Calculate the Q cumulant up to order fMaxMoment
1380 // cent: centrality of event
1382 if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1384 // Fill out NUA hists
1385 for (Int_t etaBin = 1; etaBin <= fCumuRef->GetNbinsX(); etaBin++) {
1386 Double_t eta = fCumuRef->GetXaxis()->GetBinCenter(etaBin);
1387 if (fCumuRef->GetBinContent(etaBin, 0) == 0) continue;
1388 if ((fFlags & kTracks) && (fFlags && kSPD) && !(fFlags & kEtaGap)) eta = -eta;
1389 for (Int_t qBin = 0; qBin <= fCumuRef->GetNbinsY(); qBin++) {
1390 fCumuNUARef->Fill(eta, cent, Double_t(qBin), fCumuRef->GetBinContent(etaBin, qBin));
1393 for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1394 Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1395 if (fCumuDiff->GetBinContent(etaBin, 0) == 0) continue;
1396 for (Int_t qBin = 0; qBin <= fCumuDiff->GetNbinsY(); qBin++) {
1397 fCumuNUADiff->Fill(eta, cent, Double_t(qBin), fCumuDiff->GetBinContent(etaBin, qBin));
1401 // We create the objects needed for the analysis
1404 // For each n we loop over the hists
1405 for (Int_t n = 2; n <= fMaxMoment; n++) {
1406 cumuRef = (TH3D*)fCumuHists.Get('r',n);
1407 cumuDiff = (TH3D*)fCumuHists.Get('d',n);
1408 Int_t prevRefEtaBin = 0;
1410 // Per mom. quantities
1411 Double_t dQnReA = 0, dQnImA = 0, multA = 0;
1412 Double_t dQnReB = 0, dQnImB = 0, multB = 0;
1413 Double_t dQ2nReA = 0, dQ2nImA = 0;
1414 Double_t two = 0, w2 = 0, four = 0, w4 = 0;
1415 Double_t cosPhi1Phi2 = 0, cosPhi1Phi2Phi3m = 0;
1416 Double_t sinPhi1Phi2 = 0, sinPhi1Phi2Phi3m = 0;
1417 for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1418 Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1419 Double_t refEta = eta;
1420 if ((fFlags & kTracks) && (fFlags && kSPD) && !(fFlags & kEtaGap)) refEta = -eta;
1421 Int_t refEtaBinA = fCumuRef->GetXaxis()->FindBin(refEta);
1422 if ((fFlags & kEtaGap)) refEta = -eta;
1423 Int_t refEtaBinB = fCumuRef->GetXaxis()->FindBin(refEta);
1424 if (refEtaBinA != prevRefEtaBin) {
1425 prevRefEtaBin = refEtaBinA;
1427 multA = fCumuRef->GetBinContent(refEtaBinA, 0);
1428 dQnReA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberCos(n));
1429 dQnImA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberSin(n));
1430 dQ2nReA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberCos(n*2));
1431 dQ2nImA = fCumuRef->GetBinContent(refEtaBinA, GetBinNumberSin(n*2));
1433 multB = fCumuRef->GetBinContent(refEtaBinB, 0);
1434 dQnReB = fCumuRef->GetBinContent(refEtaBinB, GetBinNumberCos(n));
1435 dQnImB = fCumuRef->GetBinContent(refEtaBinB, GetBinNumberSin(n));
1437 if (multA <= 3 || multB <= 3) return;
1438 // The reference flow is calculated
1440 if ((fFlags & kStdQC)) {
1441 w2 = multA * (multA - 1.);
1442 two = dQnReA*dQnReA + dQnImA*dQnImA - multA;
1445 two = dQnReA*dQnReB + dQnImA*dQnImB;
1447 cumuRef->Fill(eta, cent, kW2Two, two);
1448 cumuRef->Fill(eta, cent, kW2, w2);
1450 // The reference flow is calculated
1452 if ((fFlags & kStdQC)) {
1453 w4 = multA * (multA - 1.) * (multA - 2.) * (multA - 3.);
1455 four = 2.*multA*(multA-3.) + TMath::Power((TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.)),2.)
1456 -4.*(multA-2.)*(TMath::Power(dQnReA,2.) + TMath::Power(dQnImA,2.))
1457 -2.*(TMath::Power(dQnReA,2.)*dQ2nReA+2.*dQnReA*dQnImA*dQ2nImA-TMath::Power(dQnImA,2.)*dQ2nReA)
1458 +(TMath::Power(dQ2nReA,2.)+TMath::Power(dQ2nImA,2.));
1460 cumuRef->Fill(eta, cent, kW4Four, four);
1461 cumuRef->Fill(eta, cent, kW4, w4);
1464 cosPhi1Phi2 = dQnReA*dQnReA - dQnImA*dQnImA - dQ2nReA;
1465 sinPhi1Phi2 = 2.*dQnReA*dQnImA - dQ2nImA;
1467 cosPhi1Phi2Phi3m = dQnReA*(TMath::Power(dQnReA,2)+TMath::Power(dQnImA,2))
1468 -dQnReA*dQ2nReA-dQnImA*dQ2nImA-2.*(multA-1)*dQnReA;
1470 sinPhi1Phi2Phi3m = -dQnImA*(TMath::Power(dQnReA,2)+TMath::Power(dQnImA,2))
1471 +dQnReA*dQ2nImA-dQnImA*dQ2nReA+2.*(multA-1)*dQnImA;
1473 cumuRef->Fill(eta, cent, kCosphi1phi2, cosPhi1Phi2);
1474 cumuRef->Fill(eta, cent, kSinphi1phi2, sinPhi1Phi2);
1475 cumuRef->Fill(eta, cent, kCosphi1phi2phi3m, cosPhi1Phi2Phi3m);
1476 cumuRef->Fill(eta, cent, kSinphi1phi2phi3m, sinPhi1Phi2Phi3m);
1477 cumuRef->Fill(eta, cent, k3pWeight, multA*(multA-1.)*(multA-2.));
1479 } // End of reference flow
1480 // For each etaBin bin the necessary values for differential flow is calculated
1481 Double_t mp = fCumuDiff->GetBinContent(etaBin, 0);
1482 Double_t pnRe = fCumuDiff->GetBinContent(etaBin, GetBinNumberCos(n));
1483 Double_t pnIm = fCumuDiff->GetBinContent(etaBin, GetBinNumberSin(n));
1484 Double_t p2nRe = fCumuDiff->GetBinContent(etaBin, GetBinNumberCos(n*2));
1485 Double_t p2nIm = fCumuDiff->GetBinContent(etaBin, GetBinNumberSin(n*2));
1486 if (mp == 0) continue;
1493 // Differential flow calculations for each eta bin is done:
1494 // 2-particle differential flow
1495 if ((fFlags & kStdQC) && (!(fFlags & kTracks) || ((fFlags & kTracks) && (fFlags & kMC) && !(fFlags & kSPD) && TMath::Abs(eta) < 0.75))) {
1503 Double_t w2p = mp * multB - mq;
1504 Double_t twoPrime = pnRe*dQnReB + pnIm*dQnImB - mq;
1506 cumuDiff->Fill(eta, cent, kW2Two, twoPrime);
1507 cumuDiff->Fill(eta, cent, kW2, w2p);
1509 if ((fFlags & kEtaGap)) continue;
1510 // Differential flow calculations for each eta bin bin is done:
1511 // 4-particle differential flow
1512 Double_t w4p = (mp * multA - 3.*mq)*(multA - 1.)*(multA - 2.);
1514 Double_t fourPrime = (TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.))*(pnRe*dQnReA+pnIm*dQnImA)
1515 - q2nRe*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))
1516 - 2.*q2nIm*dQnReA*dQnImA
1517 - pnRe*(dQnReA*dQ2nReA+dQnImA*dQ2nImA)
1518 + pnIm*(dQnImA*dQ2nReA-dQnReA*dQ2nImA)
1519 - 2.*multA*(pnRe*dQnReA+pnIm*dQnImA)
1520 - 2.*(TMath::Power(dQnReA,2.)+TMath::Power(dQnImA,2.))*mq
1521 + 6.*(qnRe*dQnReA+qnIm*dQnImA)
1522 + 1.*(q2nRe*dQ2nReA+q2nIm*dQ2nImA)
1523 + 2.*(pnRe*dQnReA+pnIm*dQnImA)
1527 cumuDiff->Fill(eta, cent, kW4Four, fourPrime);
1528 cumuDiff->Fill(eta, cent, kW4, w4p);
1531 Double_t cosPsi1Phi2 = pnRe*dQnReA - pnIm*dQnImA - q2nRe;
1532 Double_t sinPsi1Phi2 = pnRe*dQnImA + pnIm*dQnReA - q2nIm;
1534 Double_t cosPsi1Phi2Phi3p = pnRe*(TMath::Power(dQnImA,2.)+TMath::Power(dQnReA,2.)-multA)
1535 - 1.*(q2nRe*dQnReA+q2nIm*dQnImA)
1536 - mq*dQnReA+2.*qnRe;
1538 Double_t sinPsi1Phi2Phi3p = pnIm*(TMath::Power(dQnImA,2.)+TMath::Power(dQnReA,2.)-multA)
1539 - 1.*(q2nIm*dQnReA-q2nRe*dQnImA)
1540 - mq*dQnImA+2.*qnIm;
1542 Double_t cosPsi1Phi2Phi3m = pnRe*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))+2.*pnIm*dQnReA*dQnImA
1543 - 1.*(pnRe*dQ2nReA+pnIm*dQ2nImA)
1544 - 2.*mq*dQnReA+2.*qnRe;
1546 Double_t sinPsi1Phi2Phi3m = pnIm*(TMath::Power(dQnReA,2.)-TMath::Power(dQnImA,2.))-2.*pnRe*dQnReA*dQnImA
1547 - 1.*(pnIm*dQ2nReA-pnRe*dQ2nImA)
1548 + 2.*mq*dQnImA-2.*qnIm;
1550 cumuDiff->Fill(eta, cent, kCosphi1phi2, cosPsi1Phi2);
1551 cumuDiff->Fill(eta, cent, kSinphi1phi2, sinPsi1Phi2);
1552 cumuDiff->Fill(eta, cent, kCosphi1phi2phi3m, cosPsi1Phi2Phi3m);
1553 cumuDiff->Fill(eta, cent, kSinphi1phi2phi3m, sinPsi1Phi2Phi3m);
1554 cumuDiff->Fill(eta, cent, k3pWeight, (mp*multA-2.*mq)*(multA-1.));
1555 cumuDiff->Fill(eta, cent, kCosphi1phi2phi3p, cosPsi1Phi2Phi3p);
1556 cumuDiff->Fill(eta, cent, kSinphi1phi2phi3p, sinPsi1Phi2Phi3p);
1557 } // End of eta loop
1559 cumuRef->Fill(-7., cent, -0.5, 1.);
1560 } // End of moment loop
1563 //_____________________________________________________________________
1564 void AliForwardFlowTaskQC::VertexBin::GetLimits(Int_t bin, Int_t& aLow, Int_t& aHigh,
1565 Int_t& bLow, Int_t& bHigh) const
1568 // Get the limits for the 3 correlator method
1571 // bin : reference bin #
1572 // aLow : Lowest bin to be included in v_A calculations
1573 // aHigh: Highest bin to be included in v_A calculations
1574 // bLow : Lowest bin to be included in v_B calculations
1575 // bHigh: Highest bin to be included in v_B calculations
1577 if ((fFlags & kFMD)) {
1580 aLow = 14; aHigh = 15;
1581 bLow = 20; bHigh = 22;
1584 aLow = 16; aHigh = 16;
1585 bLow = 21; bHigh = 22;
1588 aLow = 6; aHigh = 7;
1589 bLow = 21; bHigh = 22;
1592 aLow = 6; aHigh = 7;
1593 bLow = 12; bHigh = 12;
1596 aLow = 6; aHigh = 8;
1597 bLow = 13; bHigh = 14;
1600 AliFatal(Form("No limits for this eta region! (%d)", bin));
1603 else if ((fFlags & kVZERO)) {
1606 aLow = 6; aHigh = 13;
1607 bLow = 17; bHigh = 18;
1610 aLow = 6; aHigh = 9;
1611 bLow = 17; bHigh = 18;
1614 aLow = 2; aHigh = 3;
1615 bLow = 17; bHigh = 18;
1618 aLow = 2; aHigh = 3;
1619 bLow = 6; bHigh = 9;
1622 aLow = 2; aHigh = 3;
1623 bLow = 6; bHigh = 13;
1626 AliFatal(Form("No limits for this eta region! (%d)", bin));
1629 // Try to catch cases where fEtaLimits and these values do not correspond to each other
1630 if (aHigh > fCumuNUARef->GetNbinsX() || bHigh > fCumuNUARef->GetNbinsX())
1631 AliFatal(Form("Limits outside vtx range! (%d) - aHigh = %d, bHigh = %d, Nbins = %d",
1632 bin, aHigh, bHigh, fCumuNUARef->GetNbinsX()));
1634 //_____________________________________________________________________
1635 void AliForwardFlowTaskQC::VertexBin::CumulantsAccumulate3Cor(Double_t cent)
1638 // Calculate the Q cumulant up to order fMaxMoment
1641 // cent: centrality of event
1643 if (!fCumuRef) AliFatal("You have not called AddOutput() - Terminating!");
1645 // Fill out NUA hists
1646 for (Int_t etaBin = 1; etaBin <= fCumuRef->GetNbinsX(); etaBin++) {
1647 Double_t eta = fCumuRef->GetXaxis()->GetBinCenter(etaBin);
1648 if (fCumuRef->GetBinContent(etaBin, 0) == 0) continue;
1649 for (Int_t qBin = 0; qBin <= fCumuRef->GetNbinsY(); qBin++) {
1650 fCumuNUARef->Fill(eta, cent, Double_t(qBin), fCumuRef->GetBinContent(etaBin, qBin));
1653 for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1654 Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1655 if (fCumuDiff->GetBinContent(etaBin, 0) == 0) continue;
1656 for (Int_t qBin = 0; qBin <= fCumuDiff->GetNbinsY(); qBin++) {
1657 fCumuNUADiff->Fill(eta, cent, Double_t(qBin), fCumuDiff->GetBinContent(etaBin, qBin));
1661 // We create the objects needed for the analysis
1664 // For each n we loop over the hists
1665 for (Int_t n = 2; n <= fMaxMoment; n++) {
1666 cumuRef = (TH3D*)fCumuHists.Get('r',n);
1667 cumuDiff = (TH3D*)fCumuHists.Get('d',n);
1669 // Per mom. quantities
1671 Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
1672 Double_t dQnReA = 0, dQnImA = 0, multA = 0;
1673 Double_t dQnReB = 0, dQnImB = 0, multB = 0;
1674 Double_t two = 0, w2 = 0;
1675 for (Int_t etaBin = 1; etaBin <= fCumuDiff->GetNbinsX(); etaBin++) {
1676 Double_t eta = fCumuDiff->GetXaxis()->GetBinCenter(etaBin);
1677 if (fEtaLims[prevLim] < eta) {
1678 GetLimits(prevLim, aLow, aHigh, bLow, bHigh);
1680 multA = 0; dQnReA = 0; dQnImA = 0;
1681 multB = 0; dQnReB = 0; dQnImB = 0;
1683 for (Int_t a = aLow; a <= aHigh; a++) {
1684 multA += fCumuRef->GetBinContent(a, 0);
1685 dQnReA += fCumuRef->GetBinContent(a, GetBinNumberCos(n));
1686 dQnImA += fCumuRef->GetBinContent(a, GetBinNumberSin(n));
1688 for (Int_t b = bLow; b <= bHigh; b++) {
1689 multB += fCumuRef->GetBinContent(b, 0);
1690 dQnReB += fCumuRef->GetBinContent(b, GetBinNumberCos(n));
1691 dQnImB += fCumuRef->GetBinContent(b, GetBinNumberSin(n));
1693 // The reference flow is calculated
1696 two = dQnReA*dQnReB + dQnImA*dQnImB;
1697 } // End of reference flow
1698 cumuRef->Fill(eta, cent, kW2Two, two);
1699 cumuRef->Fill(eta, cent, kW2, w2);
1701 // For each etaBin bin the necessary values for differential flow is calculated
1702 Double_t mp = fCumuDiff->GetBinContent(etaBin, 0);
1703 Double_t pnRe = fCumuDiff->GetBinContent(etaBin, GetBinNumberCos(n));
1704 Double_t pnIm = fCumuDiff->GetBinContent(etaBin, GetBinNumberSin(n));
1705 if (mp == 0) continue;
1707 // Differential flow calculations for each eta bin is done:
1708 // 2-particle differential flow
1709 Double_t w2pA = mp * multA;
1710 Double_t twoPrimeA = pnRe*dQnReA + pnIm*dQnImA;
1711 cumuDiff->Fill(eta, cent, kW2Two, twoPrimeA);
1712 cumuDiff->Fill(eta, cent, kW2, w2pA);
1714 Double_t w2pB = mp * multB;
1715 Double_t twoPrimeB = pnRe*dQnReB + pnIm*dQnImB;
1716 cumuDiff->Fill(eta, cent, kW4Four, twoPrimeB);
1717 cumuDiff->Fill(eta, cent, kW4, w2pB);
1718 } // End of eta loop
1720 cumuRef->Fill(-7., cent, -0.5, 1.);
1721 } // End of moment loop
1725 //_____________________________________________________________________
1726 void AliForwardFlowTaskQC::VertexBin::CumulantsTerminate(TList* inlist, TList* outlist)
1729 // Finalizes the Q cumulant calculations
1732 // inlist: input sumlist
1733 // outlist: output result list
1736 // Re-find cumulants hist if Terminate is called separately
1737 if (!fCumuHists.IsConnected()) {
1738 TList* list = (TList*)inlist->FindObject(Form("%svertex_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1739 fCumuHists.ConnectList(Form("%sCumu_%d_%d%s", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)), list);
1742 fCumuNUARef = (TH3D*)list->FindObject(Form("%s_vertex_%d_%d%s_cumuNUARef", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1744 fCumuNUADiff = (TH3D*)list->FindObject(Form("%s_vertex_%d_%d%s_cumuNUADiff", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1746 // Clone to avoid normalization problems when redoing terminate locally
1747 fCumuNUARef = (TH3D*)fCumuNUARef->Clone(Form("%s_vertex_%d_%d%s_cumuNUARefNorm", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1748 fCumuNUADiff = (TH3D*)fCumuNUADiff->Clone(Form("%s_vertex_%d_%d%s_cumuNUADiffNorm", fType.Data(), fVzMin, fVzMax, GetQCType(fFlags)));
1750 // Diagnostics histograms
1751 TH2I* quality = (TH2I*)outlist->FindObject(Form("hQCQuality%s%s", fType.Data(), GetQCType(fFlags)));
1753 quality = MakeQualityHist(Form("hQCQuality%s%s", fType.Data(), GetQCType(fFlags)));
1754 outlist->Add(quality);
1756 TH1D* cent = (TH1D*)outlist->FindObject(Form("%s%s_cent", fType.Data(), GetQCType(fFlags)));
1758 cent = new TH1D(Form("%s%s_cent", fType.Data(), GetQCType(fFlags)),
1759 Form("%s%s_cent", fType.Data(), GetQCType(fFlags)),
1760 fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXmin(), fCumuNUARef->GetYaxis()->GetXmax());
1761 cent->GetXaxis()->Set(fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXbins()->GetArray());
1764 TH2D* dNdetaRef = (TH2D*)outlist->FindObject(Form("%s%s_dNdetaRef", fType.Data(), GetQCType(fFlags)));
1766 dNdetaRef = new TH2D(Form("%s%s_dNdetaRef", fType.Data(), GetQCType(fFlags)),
1767 Form("%s%s_dNdetaRef", fType.Data(), GetQCType(fFlags)),
1768 fCumuNUARef->GetNbinsX(), fCumuNUARef->GetXaxis()->GetXmin(), fCumuNUARef->GetXaxis()->GetXmax(),
1769 fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXmin(), fCumuNUARef->GetYaxis()->GetXmax());
1770 dNdetaRef->GetYaxis()->Set(fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXbins()->GetArray());
1772 outlist->Add(dNdetaRef);
1774 TH2D* dNdetaDiff = (TH2D*)outlist->FindObject(Form("%s%s_dNdetaDiff", fType.Data(), GetQCType(fFlags)));
1776 dNdetaDiff = new TH2D(Form("%s%s_dNdetaDiff", fType.Data(), GetQCType(fFlags)),
1777 Form("%s%s_dNdetaDiff", fType.Data(), GetQCType(fFlags)),
1778 fCumuNUADiff->GetNbinsX(), fCumuNUADiff->GetXaxis()->GetXmin(), fCumuNUADiff->GetXaxis()->GetXmax(),
1779 fCumuNUADiff->GetNbinsY(), fCumuNUADiff->GetYaxis()->GetXmin(), fCumuNUADiff->GetYaxis()->GetXmax());
1780 dNdetaDiff->GetYaxis()->Set(fCumuNUADiff->GetNbinsY(), fCumuNUADiff->GetYaxis()->GetXbins()->GetArray());
1781 dNdetaDiff->Sumw2();
1782 outlist->Add(dNdetaDiff);
1785 // Setting up outputs
1786 // Create output lists and diagnostics
1787 TList* vtxList = (TList*)outlist->FindObject("vtxList");
1789 vtxList = new TList();
1790 vtxList->SetName("vtxList");
1791 outlist->Add(vtxList);
1793 vtxList->Add(fCumuNUARef);
1794 vtxList->Add(fCumuNUADiff);
1796 // Setup output profiles
1797 CumuHistos cumu2(fMaxMoment, ((fFlags & kNUAcorr) ? 2 : 0));
1798 CumuHistos cumu4(fMaxMoment, ((fFlags & kNUAcorr) ? 1 : 0));
1800 cumu2.ConnectList(Form("%sQC2_Cumu%s_vtx_%d_%d", fType.Data(), GetQCType(fFlags), fVzMin, fVzMax), vtxList);
1801 if ((fFlags & kStdQC))
1802 cumu4.ConnectList(Form("%sQC4_Cumu%s_vtx_%d_%d", fType.Data(), GetQCType(fFlags), fVzMin, fVzMax), vtxList);
1804 for (Int_t n = 2; n <= fMaxMoment; n++) {
1806 cumu2.Add(MakeOutputHist(2, n, "Ref", CumuHistos::kNoNUA));
1807 if ((fFlags & k3Cor)){
1808 cumu2.Add(MakeOutputHist(2, n, "DiffA", CumuHistos::kNoNUA));
1809 cumu2.Add(MakeOutputHist(2, n, "DiffB", CumuHistos::kNoNUA));
1811 cumu2.Add(MakeOutputHist(2, n, "Diff", CumuHistos::kNoNUA));
1814 if ((fFlags & kStdQC)) {
1815 cumu4.Add(MakeOutputHist(4, n, "Ref", CumuHistos::kNoNUA));
1816 cumu4.Add(MakeOutputHist(4, n, "Diff", CumuHistos::kNoNUA));
1818 } // End of v_n result loop
1820 if ((fFlags & kNUAcorr)) {
1821 for (Int_t n = 2; n <= fMaxMoment; n++) {
1823 cumu2.Add(MakeOutputHist(2, n, "Ref", CumuHistos::kNUAOld));
1824 if ((fFlags & k3Cor)) {
1825 cumu2.Add(MakeOutputHist(2, n, "DiffA", CumuHistos::kNUAOld));
1826 cumu2.Add(MakeOutputHist(2, n, "DiffB", CumuHistos::kNUAOld));
1828 cumu2.Add(MakeOutputHist(2, n, "Diff", CumuHistos::kNUAOld));
1831 if ((fFlags & kStdQC)) {
1832 cumu4.Add(MakeOutputHist(4, n, "Ref", CumuHistos::kNUAOld));
1833 cumu4.Add(MakeOutputHist(4, n, "Diff", CumuHistos::kNUAOld));
1836 for (Int_t n = 2; n <= fMaxMoment; n++) {
1838 cumu2.Add(MakeOutputHist(2, n, "Ref", CumuHistos::kNUA));
1839 if ((fFlags & k3Cor)) {
1840 cumu2.Add(MakeOutputHist(2, n, "DiffA", CumuHistos::kNUA));
1841 cumu2.Add(MakeOutputHist(2, n, "DiffB", CumuHistos::kNUA));
1843 cumu2.Add(MakeOutputHist(2, n, "Diff", CumuHistos::kNUA));
1848 // Calculating the cumulants
1849 if ((fFlags & k3Cor)) {
1850 Calculate3CorFlow(cumu2, quality, cent, dNdetaRef, dNdetaDiff);
1852 CalculateReferenceFlow(cumu2, cumu4, quality, cent, dNdetaRef);
1853 CalculateDifferentialFlow(cumu2, cumu4, quality, dNdetaDiff);
1855 if ((fFlags & kNUAcorr)) {
1856 SolveCoupledFlowEquations(cumu2, 'r');
1857 if ((fFlags & k3Cor)) {
1858 SolveCoupledFlowEquations(cumu2, 'a');
1859 SolveCoupledFlowEquations(cumu2, 'b');
1861 SolveCoupledFlowEquations(cumu2, 'd');
1865 // Add to output for immediate viewing - individual vtx bins are used for final results
1866 AddVertexBins(cumu2, outlist, ((fFlags & kNUAcorr) ? 2 : 0));
1867 if ((fFlags & kStdQC)) AddVertexBins(cumu4, outlist, ((fFlags & kNUAcorr) ? 1 : 0));
1869 // Setup NUA diagnoastics histograms
1870 TList* nualist = (TList*)outlist->FindObject("NUATerms");
1872 nualist = new TList();
1873 nualist->SetName("NUATerms");
1874 outlist->Add(nualist);
1877 TH2D* nuaRef = (TH2D*)nualist->FindObject(Form("%sReferenceNUA%s", fType.Data(), GetQCType(fFlags)));
1880 nuaRef = (TH2D*)fCumuNUARef->Project3D("yz");
1881 nuaRef->Scale(1./fCumuNUARef->GetNbinsX());
1882 nuaRef->SetName(Form("%sReferenceNUA%s", fType.Data(), GetQCType(fFlags)));
1883 nuaRef->SetTitle(Form("%sReferenceNUA%s", fType.Data(), GetQCType(fFlags)));
1884 nualist->Add(nuaRef);
1886 temp = (TH2D*)fCumuNUARef->Project3D("yz");
1887 temp->Scale(1./fCumuNUARef->GetNbinsX());
1891 // Filling in underflow to make scaling possible in Terminate()
1892 nuaRef->Fill(0., -1., 1.);
1894 TH2D* nuaDiff = (TH2D*)nualist->FindObject(Form("%sDifferentialNUA%s", fType.Data(), GetQCType(fFlags)));
1896 nuaDiff = (TH2D*)fCumuNUADiff->Project3D("yz");
1897 nuaDiff->SetName(Form("%sDifferentialNUA%s", fType.Data(), GetQCType(fFlags)));
1898 nuaDiff->SetTitle(Form("%sDifferentialNUA%s", fType.Data(), GetQCType(fFlags)));
1899 nualist->Add(nuaDiff);
1901 temp = (TH2D*)fCumuNUADiff->Project3D("yz");
1905 // Filling in underflow to make scaling possible in Terminate()
1906 nuaDiff->Fill(0., -1., 1.);
1910 //_____________________________________________________________________
1911 void AliForwardFlowTaskQC::VertexBin::CalculateReferenceFlow(CumuHistos& cumu2h, CumuHistos& cumu4h, TH2I* quality,
1912 TH1D* chist, TH2D* dNdetaRef) const
1915 // Calculates the reference flow
1918 // cumu2h: CumuHistos object with QC{2} cumulants
1919 // cumu4h: CumuHistos object with QC{4} cumulants
1920 // quality: Histogram for success rate of cumulants
1921 // chist: Centrality histogram
1922 // dNdetaRef: dN/deta histogram for estimating multiplicity used for ref calculations
1925 // Normalizing common NUA hists
1926 for (Int_t cBin = 1; cBin <= fCumuNUARef->GetNbinsY(); cBin++) {
1927 Double_t cent = fCumuNUARef->GetYaxis()->GetBinCenter(cBin);
1928 for (Int_t eBin = 1; eBin <= fCumuNUARef->GetNbinsX(); eBin++) {
1929 Double_t eta = fCumuNUARef->GetXaxis()->GetBinCenter(eBin);
1930 Double_t mult = fCumuNUARef->GetBinContent(eBin, cBin, 0);
1931 if (mult == 0) continue;
1932 for (Int_t qBin = 1; qBin <= fCumuNUARef->GetNbinsZ(); qBin++) {
1933 fCumuNUARef->SetBinContent(eBin, cBin, qBin, fCumuNUARef->GetBinContent(eBin, cBin, qBin)/mult);
1934 fCumuNUARef->SetBinError(eBin, cBin, qBin, fCumuNUARef->GetBinError(eBin, cBin, qBin)/mult);
1936 // Fill dN/deta diagnostics
1937 dNdetaRef->Fill(eta, cent, mult);
1941 // For flow calculations
1944 TH2D* cumu2NUAold = 0;
1947 Int_t qualityFactor = ((fFlags & kStdQC) ? 8 : 4); // used for correctly filling in quality hists
1948 // Loop over cumulant histogram for final calculations
1949 for (Int_t n = 2; n <= fMaxMoment; n++) { // Moment loop begins
1950 cumu2 = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNoNUA);
1951 if ((fFlags & kNUAcorr)) {
1952 cumu2NUAold = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNUAOld);
1954 if ((fFlags & kStdQC)) {
1955 cumu4 = (TH2D*)cumu4h.Get('r', n, CumuHistos::kNoNUA);
1956 if ((fFlags & kNUAcorr)) cumu4NUA = (TH2D*)cumu4h.Get('r', n, CumuHistos::kNUAOld);
1958 cumuRef = (TH3D*)fCumuHists.Get('r', n);
1960 for (Int_t cBin = 1; cBin <= cumuRef->GetNbinsY(); cBin++) { // Centrality loop begins
1961 Double_t cent = cumuRef->GetYaxis()->GetBinCenter(cBin);
1962 if (n == 2) chist->Fill(cent, cumuRef->GetBinContent(0, cBin, 0));
1963 if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), n, cent));
1964 for (Int_t etaBin = 1; etaBin <= cumuRef->GetNbinsX(); etaBin++) { // Eta loop begins
1965 Double_t eta = cumuRef->GetXaxis()->GetBinCenter(etaBin);
1966 Double_t refEta = eta;
1967 Int_t refEtaBinA = fCumuNUARef->GetXaxis()->FindBin(refEta);
1968 if ((fFlags & kEtaGap)) refEta = -eta;
1969 Int_t refEtaBinB = fCumuNUARef->GetXaxis()->FindBin(refEta);
1970 // 2-particle reference flow
1971 Double_t w2Two = cumuRef->GetBinContent(refEtaBinA, cBin, kW2Two);
1972 Double_t w2 = cumuRef->GetBinContent(refEtaBinA, cBin, kW2);
1973 if (w2 == 0) continue;
1974 Double_t cosP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberCos(n));
1975 Double_t sinP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberSin(n));
1976 Double_t cosP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(n));
1977 Double_t sinP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(n));
1978 Double_t cos2nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberCos(2*n));
1979 Double_t sin2nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberSin(2*n));
1980 Double_t cos2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(2*n));
1981 Double_t sin2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(2*n));
1982 Double_t two = w2Two / w2;
1984 if (qc2 >= 0) cumu2->Fill(eta, cent, TMath::Sqrt(qc2));
1986 if ((fFlags & kNUAcorr)) {
1988 // With no eta gap the last two terms are <<cos(phi)>>^2 and <<sin(phi)>>^2,
1989 // with eta gap the different coverage is taken into account.
1990 // The next line covers both cases.
1991 qc2 -= cosP1nPhiA*cosP1nPhiB + sinP1nPhiA*sinP1nPhiB;
1992 // Extra NUA term from 2n cosines and sines
1993 Double_t den = 1-(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB);
1994 if (den != 0) qc2 /= den;
1999 AliInfo(Form("%s: QC_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f - skipping",
2000 fType.Data(), n, qc2, eta, cent));
2001 quality->Fill((n-2)*qualityFactor+2, Int_t(cent));
2004 Double_t vnTwo = TMath::Sqrt(qc2);
2005 if (!TMath::IsNaN(vnTwo)) {
2006 quality->Fill((n-2)*qualityFactor+1, Int_t(cent));
2007 if ((fFlags & kNUAcorr)) cumu2NUAold->Fill(eta, cent, vnTwo);
2010 if (!(fFlags & kStdQC)) continue;
2011 // 4-particle reference flow
2012 Double_t w4Four = cumuRef->GetBinContent(refEtaBinA, cBin, kW4Four);
2013 Double_t w4 = cumuRef->GetBinContent(refEtaBinA, cBin, kW4);
2014 Double_t multm1m2 = cumuRef->GetBinContent(refEtaBinA, cBin, k3pWeight);
2015 if (w4 == 0 || multm1m2 == 0) continue;
2016 Double_t cosP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2);
2017 Double_t sinP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2);
2018 Double_t cosP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2phi3m);
2019 Double_t sinP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2phi3m);
2021 cosP1nPhi1P1nPhi2 /= w2;
2022 sinP1nPhi1P1nPhi2 /= w2;
2023 cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2024 sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2025 Double_t four = w4Four / w4;
2026 Double_t qc4 = four-2.*TMath::Power(two,2.);
2027 if (qc4 < 0) cumu4->Fill(eta, cent, TMath::Power(-qc4, 0.25));
2029 if ((fFlags & kNUAcorr)) {
2030 qc4 += - 4.*cosP1nPhiA*cosP1nPhi1M1nPhi2M1nPhi3
2031 + 4.*sinP1nPhiA*sinP1nPhi1M1nPhi2M1nPhi3-TMath::Power(cosP1nPhi1P1nPhi2,2.)-TMath::Power(sinP1nPhi1P1nPhi2,2.)
2032 + 4.*cosP1nPhi1P1nPhi2*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
2033 + 8.*sinP1nPhi1P1nPhi2*sinP1nPhiA*cosP1nPhiA
2034 + 8.*two*(TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.))
2035 - 6.*TMath::Power((TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.)),2.);
2039 AliInfo(Form("%s: QC_%d{4} = %1.3f for eta = %1.2f and centrality %3.1f - skipping",
2040 fType.Data(), n, qc2, eta, cent));
2041 quality->Fill((n-2)*qualityFactor+6, Int_t(cent));
2044 Double_t vnFour = TMath::Power(-qc4, 0.25);
2045 if (!TMath::IsNaN(vnFour*multm1m2)) {
2046 quality->Fill((n-2)*qualityFactor+5, Int_t(cent));
2047 if ((fFlags & kNUAcorr)) cumu4NUA->Fill(eta, cent, vnFour);
2055 //_____________________________________________________________________
2056 void AliForwardFlowTaskQC::VertexBin::CalculateDifferentialFlow(CumuHistos& cumu2h, CumuHistos& cumu4h,
2057 TH2I* quality, TH2D* dNdetaDiff) const
2060 // Calculates the differential flow
2063 // cumu2h: CumuHistos object with QC{2} cumulants
2064 // cumu4h: CumuHistos object with QC{4} cumulants
2065 // quality: Histogram for success rate of cumulants
2066 // dNdetaDiff: dN/deta histogram for estimating multiplicity used for diff calculations
2069 for (Int_t cBin = 1; cBin <= fCumuNUADiff->GetNbinsY(); cBin++) {
2070 Double_t cent = fCumuNUADiff->GetYaxis()->GetBinCenter(cBin);
2071 for (Int_t eBin = 1; eBin <= fCumuNUADiff->GetNbinsX(); eBin++) {
2072 Double_t eta = fCumuNUADiff->GetXaxis()->GetBinCenter(eBin);
2073 Double_t mult = fCumuNUADiff->GetBinContent(eBin, cBin, 0);
2074 if (mult == 0) continue;
2075 for (Int_t qBin = 1; qBin <= fCumuNUADiff->GetNbinsZ(); qBin++) {
2076 fCumuNUADiff->SetBinContent(eBin, cBin, qBin, fCumuNUADiff->GetBinContent(eBin, cBin, qBin)/mult);
2077 fCumuNUADiff->SetBinError(eBin, cBin, qBin, fCumuNUADiff->GetBinError(eBin, cBin, qBin)/mult);
2079 dNdetaDiff->Fill(eta, cent, mult);
2083 // For flow calculations
2087 TH2D* cumu2NUAold = 0;
2090 Int_t qualityFactor = ((fFlags & kStdQC) ? 8 : 4); // used for correctly filling in quality hists
2091 // Loop over cumulant histogram for final calculations
2092 for (Int_t n = 2; n <= fMaxMoment; n++) { // Moment loop begins
2093 cumu2 = (TH2D*)cumu2h.Get('d', n, CumuHistos::kNoNUA);
2094 if ((fFlags & kNUAcorr)) {
2095 cumu2NUAold = (TH2D*)cumu2h.Get('d', n, CumuHistos::kNUAOld);
2097 if ((fFlags & kStdQC)) {
2098 cumu4 = (TH2D*)cumu4h.Get('d',n);
2099 if ((fFlags & kNUAcorr)) cumu4NUA = (TH2D*)cumu4h.Get('d', n, CumuHistos::kNUAOld);
2101 cumuRef = (TH3D*)fCumuHists.Get('r',n);
2102 cumuDiff = (TH3D*)fCumuHists.Get('d',n);
2103 for (Int_t cBin = 1; cBin <= cumuDiff->GetNbinsY(); cBin++) { // Centrality loop begins
2104 Double_t cent = cumuDiff->GetYaxis()->GetBinCenter(cBin);
2105 if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), n, cent));
2106 for (Int_t etaBin = 1; etaBin <= cumuDiff->GetNbinsX(); etaBin++) { // Eta loop begins
2107 Double_t eta = cumuDiff->GetXaxis()->GetBinCenter(etaBin);
2108 Double_t refEta = eta;
2109 Int_t refEtaBinA = fCumuNUARef->GetXaxis()->FindBin(refEta);
2110 if ((fFlags & kEtaGap)) refEta = -eta;
2111 Int_t refEtaBinB = fCumuNUARef->GetXaxis()->FindBin(refEta);
2113 // Reference objects
2114 Double_t w2 = cumuRef->GetBinContent(refEtaBinA, cBin, kW2);
2115 if (w2 == 0) continue;
2116 Double_t two = cumuRef->GetBinContent(refEtaBinA, cBin, kW2Two);
2118 Double_t cosP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberCos(n));
2119 Double_t sinP1nPhiA = fCumuNUARef->GetBinContent(refEtaBinA, cBin, GetBinNumberSin(n));
2120 Double_t cosP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(n));
2121 Double_t sinP1nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(n));
2122 Double_t cos2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberCos(2*n));
2123 Double_t sin2nPhiB = fCumuNUARef->GetBinContent(refEtaBinB, cBin, GetBinNumberSin(2*n));
2125 // 2-particle differential flow
2126 Double_t w2pTwoPrime = cumuDiff->GetBinContent(etaBin, cBin, kW2Two);
2127 Double_t w2p = cumuDiff->GetBinContent(etaBin, cBin, kW2);
2128 if (w2p == 0) continue;
2129 Double_t cosP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(n));
2130 Double_t sinP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(n));
2131 Double_t cos2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(2*n));
2132 Double_t sin2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(2*n));
2133 Double_t twoPrime = w2pTwoPrime / w2p;
2135 Double_t qc2Prime = twoPrime;
2136 cumu2->Fill(eta, cent, qc2Prime);
2137 if ((fFlags & kNUAcorr)) {
2139 qc2Prime -= cosP1nPsi*cosP1nPhiB + sinP1nPsi*sinP1nPhiB;
2140 // Extra NUA term from 2n cosines and sines
2141 qc2Prime /= (1.-(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
2143 if (!TMath::IsNaN(qc2Prime)) {
2144 quality->Fill((n-2)*qualityFactor+3, Int_t(cent));
2145 if ((fFlags & kNUAcorr)) cumu2NUAold->Fill(eta, cent, qc2Prime);
2148 quality->Fill((n-2)*qualityFactor+4, Int_t(cent));
2150 AliInfo(Form("%s: QC'_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f",
2151 fType.Data(), n, qc2Prime, eta, cent));
2153 if (!(fFlags & kStdQC)) continue;
2154 // Reference objects
2155 Double_t cosP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2);
2156 Double_t sinP1nPhi1P1nPhi2 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2);
2157 Double_t cosP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kCosphi1phi2phi3m);
2158 Double_t sinP1nPhi1M1nPhi2M1nPhi3 = cumuRef->GetBinContent(refEtaBinA, cBin, kSinphi1phi2phi3m);
2159 Double_t multm1m2 = cumuRef->GetBinContent(refEtaBinA, cBin, k3pWeight);
2160 cosP1nPhi1P1nPhi2 /= w2;
2161 sinP1nPhi1P1nPhi2 /= w2;
2162 cosP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2163 sinP1nPhi1M1nPhi2M1nPhi3 /= multm1m2;
2165 // 4-particle differential flow
2166 Double_t w4pFourPrime = cumuDiff->GetBinContent(etaBin, cBin, kW4Four);
2167 Double_t w4p = cumuDiff->GetBinContent(etaBin, cBin, kW4);
2168 Double_t mpqMult = cumuDiff->GetBinContent(etaBin, cBin, k3pWeight);
2169 if (w4p == 0 || mpqMult == 0) continue;
2170 Double_t cosP1nPsi1P1nPhi2 = cumuDiff->GetBinContent(etaBin, cBin, kCosphi1phi2);
2171 Double_t sinP1nPsi1P1nPhi2 = cumuDiff->GetBinContent(etaBin, cBin, kSinphi1phi2);
2172 Double_t cosP1nPsi1M1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kCosphi1phi2phi3m);
2173 Double_t sinP1nPsi1M1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kSinphi1phi2phi3m);
2174 Double_t cosP1nPsi1P1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kCosphi1phi2phi3p);
2175 Double_t sinP1nPsi1P1nPhi2M1nPhi3 = cumuDiff->GetBinContent(etaBin, cBin, kSinphi1phi2phi3p);
2177 cosP1nPsi1P1nPhi2 /= w2p;
2178 sinP1nPsi1P1nPhi2 /= w2p;
2179 cosP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
2180 sinP1nPsi1M1nPhi2M1nPhi3 /= mpqMult;
2181 cosP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
2182 sinP1nPsi1P1nPhi2M1nPhi3 /= mpqMult;
2184 Double_t fourPrime = w4pFourPrime / w4p;
2185 Double_t qc4Prime = fourPrime-2.*twoPrime*two;
2186 if (cumu4) cumu4->Fill(eta, cent, qc4Prime);
2188 if ((fFlags & kNUAcorr)) {
2189 qc4Prime += - cosP1nPsi*cosP1nPhi1M1nPhi2M1nPhi3
2190 + sinP1nPsi*sinP1nPhi1M1nPhi2M1nPhi3
2191 - cosP1nPhiA*cosP1nPsi1M1nPhi2M1nPhi3
2192 + sinP1nPhiA*sinP1nPsi1M1nPhi2M1nPhi3
2193 - 2.*cosP1nPhiA*cosP1nPsi1P1nPhi2M1nPhi3
2194 - 2.*sinP1nPhiA*sinP1nPsi1P1nPhi2M1nPhi3
2195 - cosP1nPsi1P1nPhi2*cosP1nPhi1P1nPhi2
2196 - sinP1nPsi1P1nPhi2*sinP1nPhi1P1nPhi2
2197 + 2.*cosP1nPhi1P1nPhi2*(cosP1nPsi*cosP1nPhiA-sinP1nPsi*sinP1nPhiA)
2198 + 2.*sinP1nPhi1P1nPhi2*(cosP1nPsi*sinP1nPhiA+sinP1nPsi*cosP1nPhiA)
2199 + 4.*two*(cosP1nPsi*cosP1nPhiA+sinP1nPsi*sinP1nPhiA)
2200 + 2.*cosP1nPsi1P1nPhi2*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
2201 + 4.*sinP1nPsi1P1nPhi2*cosP1nPhiA*sinP1nPhiA
2202 + 4.*twoPrime*(TMath::Power(cosP1nPhiA,2.)+TMath::Power(sinP1nPhiA,2.))
2203 - 6.*(TMath::Power(cosP1nPhiA,2.)-TMath::Power(sinP1nPhiA,2.))
2204 * (cosP1nPsi*cosP1nPhiA-sinP1nPsi*sinP1nPhiA)
2205 - 12.*cosP1nPhiA*sinP1nPhiA
2206 * (sinP1nPsi*cosP1nPhiA+cosP1nPsi*sinP1nPhiA);
2208 // Double_t vnFourDiff = - qc4Prime / TMath::Power(-qc4, 0.75);
2209 if (!TMath::IsNaN(qc4Prime*mpqMult)) {
2210 quality->Fill((n-2)*qualityFactor+7, Int_t(cent));
2211 if (cumu4NUA) cumu4NUA->Fill(eta, cent, qc4Prime);
2214 quality->Fill((n-2)*qualityFactor+8, Int_t(cent));
2216 AliInfo(Form("%s: v_%d{4} = %1.3f for eta = %1.2f and centrality %3.1f",
2217 fType.Data(), n, qc4Prime, eta, cent));
2218 } // End of eta loop
2219 } // End of centrality loop
2224 //_____________________________________________________________________
2225 void AliForwardFlowTaskQC::VertexBin::Calculate3CorFlow(CumuHistos& cumu2h, TH2I* quality, TH1D* chist,
2226 TH2D* dNdetaRef, TH2D* dNdetaDiff) const
2229 // Calculates the 3 sub flow
2232 // cumu2h: CumuHistos object with QC{2} cumulants
2233 // quality: Histogram for success rate of cumulants
2234 // chist: Centrality histogram
2235 // dNdetaDiff: dN/deta histogram for estimating multiplicity used for diff calculations
2238 // For flow calculations
2242 TH2D* cumu2rNUAold = 0;
2244 TH2D* cumu2aNUAold = 0;
2246 TH2D* cumu2bNUAold = 0;
2247 // Loop over cumulant histogram for final calculations
2248 for (Int_t n = 2; n <= fMaxMoment; n++) { // Moment loop begins
2249 cumu2r = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNoNUA);
2250 cumu2a = (TH2D*)cumu2h.Get('a', n, CumuHistos::kNoNUA);
2251 cumu2b = (TH2D*)cumu2h.Get('b', n, CumuHistos::kNoNUA);
2252 if ((fFlags & kNUAcorr)) {
2253 cumu2rNUAold = (TH2D*)cumu2h.Get('r', n, CumuHistos::kNUAOld);
2254 cumu2aNUAold = (TH2D*)cumu2h.Get('a', n, CumuHistos::kNUAOld);
2255 cumu2bNUAold = (TH2D*)cumu2h.Get('b', n, CumuHistos::kNUAOld);
2257 cumuRef = (TH3D*)fCumuHists.Get('r',n);
2258 cumuDiff = (TH3D*)fCumuHists.Get('d',n);
2259 for (Int_t cBin = 1; cBin <= cumuRef->GetNbinsY(); cBin++) { // Centrality loop begins
2260 Double_t cent = cumuRef->GetYaxis()->GetBinCenter(cBin);
2261 if (n == 2) chist->Fill(cent, cumuRef->GetBinContent(0, cBin, 0));
2262 if (fDebug > 0) AliInfo(Form("%s - v_%d: centrality %3.1f:..", fType.Data(), n, cent));
2265 Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2266 Double_t cosP1nPhiA = 0;
2267 Double_t sinP1nPhiA = 0;
2268 Double_t cos2nPhiA = 0;
2269 Double_t sin2nPhiA = 0;
2270 Double_t cosP1nPhiB = 0;
2271 Double_t sinP1nPhiB = 0;
2272 Double_t cos2nPhiB = 0;
2273 Double_t sin2nPhiB = 0;
2277 for (Int_t etaBin = 1; etaBin <= cumuDiff->GetNbinsX(); etaBin++) { // Eta loop begins
2278 Double_t eta = cumuDiff->GetXaxis()->GetBinCenter(etaBin);
2279 // 2-particle reference flow
2280 Double_t w2Two = cumuRef->GetBinContent(etaBin, cBin, kW2Two);
2281 Double_t w2 = cumuRef->GetBinContent(etaBin, cBin, kW2);
2282 if (w2 == 0) continue;
2284 // Update NUA for new range!
2285 if (fEtaLims[prevLim] < eta) {
2286 GetLimits(prevLim, aLow, aHigh, bLow, bHigh);
2288 cosP1nPhiA = 0; sinP1nPhiA = 0; cos2nPhiA = 0; sin2nPhiA = 0; multA = 0;
2289 cosP1nPhiB = 0; sinP1nPhiB = 0; cos2nPhiB = 0; sin2nPhiB = 0; multB = 0;
2290 for (Int_t a = aLow; a <= aHigh; a++) {
2291 cosP1nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n));
2292 sinP1nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n));
2293 cos2nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(2*n));
2294 sin2nPhiA += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(2*n));
2295 multA += fCumuNUARef->GetBinContent(a, cBin, 0);
2297 for (Int_t b = bLow; b <= bHigh; b++) {
2298 cosP1nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n));
2299 sinP1nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n));
2300 cos2nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(2*n));
2301 sin2nPhiB += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(2*n));
2302 multB += fCumuNUARef->GetBinContent(b, cBin, 0);
2304 if (multA == 0 || multB == 0) {
2305 AliWarning(Form("Empty NUA values for 3Cor! (%s)", cumuRef->GetName()));
2308 cosP1nPhiA /= multA;
2309 sinP1nPhiA /= multA;
2312 cosP1nPhiB /= multB;
2313 sinP1nPhiB /= multB;
2317 dNdetaRef->Fill(eta, cent, multA+multB);
2319 Double_t two = w2Two / w2;
2322 if (qc2 >= 0) cumu2r->Fill(eta, cent, TMath::Sqrt(qc2));
2324 if ((fFlags & kNUAcorr)) {
2326 qc2 -= cosP1nPhiA*cosP1nPhiB + sinP1nPhiA*sinP1nPhiB;
2327 // Extra NUA term from 2n cosines and sines
2328 qc2 /= (1-(cos2nPhiA*cos2nPhiB + sin2nPhiA*sin2nPhiB));
2332 AliInfo(Form("%s: QC_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f - skipping",
2333 fType.Data(), n, qc2, eta, cent));
2334 quality->Fill((n-2)*4+2, Int_t(cent));
2337 Double_t vnTwo = TMath::Sqrt(qc2);
2338 if (!TMath::IsNaN(vnTwo)) {
2339 quality->Fill((n-2)*4+1, Int_t(cent));
2340 if ((fFlags & kNUAcorr)) cumu2rNUAold->Fill(eta, cent, vnTwo);
2343 // 2-particle differential flow
2344 Double_t w2pTwoPrimeA = cumuDiff->GetBinContent(etaBin, cBin, kW2Two);
2345 Double_t w2pA = cumuDiff->GetBinContent(etaBin, cBin, kW2);
2346 Double_t w2pTwoPrimeB = cumuDiff->GetBinContent(etaBin, cBin, kW4Four);
2347 Double_t w2pB = cumuDiff->GetBinContent(etaBin, cBin, kW4);
2348 if (w2pA == 0 || w2pB == 0) continue;
2349 Double_t cosP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(n));
2350 Double_t sinP1nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(n));
2351 Double_t cos2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberCos(2*n));
2352 Double_t sin2nPsi = fCumuNUADiff->GetBinContent(etaBin, cBin, GetBinNumberSin(2*n));
2353 Double_t mult = fCumuNUADiff->GetBinContent(etaBin, cBin, 0);
2354 if (mult == 0) continue;
2359 Double_t twoPrimeA = w2pTwoPrimeA / w2pA;
2360 Double_t twoPrimeB = w2pTwoPrimeB / w2pB;
2361 dNdetaDiff->Fill(eta, cent, mult);
2363 Double_t qc2PrimeA = twoPrimeA;
2364 Double_t qc2PrimeB = twoPrimeB;
2365 if (qc2PrimeA*qc2PrimeB >= 0) {
2366 cumu2a->Fill(eta, cent, qc2PrimeA);
2367 cumu2b->Fill(eta, cent, qc2PrimeB);
2369 if ((fFlags & kNUAcorr)) {
2371 qc2PrimeA -= cosP1nPsi*cosP1nPhiA + sinP1nPsi*sinP1nPhiA;
2372 qc2PrimeB -= cosP1nPsi*cosP1nPhiB + sinP1nPsi*sinP1nPhiB; // Is this OK?
2373 // Extra NUA term from 2n cosines and sines
2374 if (cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA != 1.) qc2PrimeA /= (1.-(cos2nPsi*cos2nPhiA + sin2nPsi*sin2nPhiA));
2375 if (cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB != 1.) qc2PrimeB /= (1.-(cos2nPsi*cos2nPhiB + sin2nPsi*sin2nPhiB));
2377 if (!TMath::IsNaN(qc2PrimeA) && !TMath::IsNaN(qc2PrimeB) && qc2 != 0) {
2378 if (qc2PrimeA*qc2PrimeB >= 0) {
2379 quality->Fill((n-2)*4+3, Int_t(cent));
2380 if ((fFlags & kNUAcorr)) cumu2aNUAold->Fill(eta, cent, qc2PrimeA);
2381 if ((fFlags & kNUAcorr)) cumu2bNUAold->Fill(eta, cent, qc2PrimeB);
2385 quality->Fill((n-2)*4+4, Int_t(cent));
2387 AliInfo(Form("%s: QC'a_%d{2} = %1.3f, QC'b_%d{2} = %1.3f for eta = %1.2f and centrality %3.1f",
2388 fType.Data(), n, qc2PrimeA, n, qc2PrimeB, eta, cent));
2389 } // End of eta loop
2390 } // End of centrality loop
2395 //_____________________________________________________________________
2396 void AliForwardFlowTaskQC::VertexBin::SolveCoupledFlowEquations(CumuHistos& cumu, const Char_t type) const
2399 // Function to solve the coupled flow equations
2400 // We solve it by using matrix calculations:
2401 // A*v_n = V => v_n = A^-1*V
2402 // First we set up a TMatrixD container to make ROOT
2403 // do the inversions in an efficient way, we multiply the current <<2>> estimates.
2404 // Then we fill new TH2D's if the corrected <<2>>'s (cumuNUA object).
2407 // cumu: CumuHistos object - uncorrected
2408 // type: Reference ('r') or differential ('d') or ('a' or 'b') for 3 correlator
2411 // We start by declaring Matrix and vector objects, as their constructors are quite heavy
2412 TMatrixD mNUA(fMaxMoment-1, fMaxMoment-1);
2413 TVectorD vQC2(fMaxMoment-1);
2415 for (Int_t cBin = 1; cBin <= cumu.Get(type, 2, CumuHistos::kNUAOld)->GetNbinsY(); cBin++) { // cent loop
2416 Double_t cent = cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinCenter(cBin);
2417 for (Int_t eBin = 1; eBin <= cumu.Get(type, 2, CumuHistos::kNUAOld)->GetNbinsX(); eBin++) { // eta loop
2418 Double_t eta = cumu.Get(type, 2, CumuHistos::kNUAOld)->GetXaxis()->GetBinCenter(eBin);
2419 mNUA.Zero(); // reset matrix
2420 vQC2.Zero(); // reset vector
2421 for (Int_t n = 0; n < fMaxMoment-1; n++) { // moment loop
2422 vQC2(n) = static_cast<TH2D*>(cumu.Get(type, n+2, CumuHistos::kNUAOld))->GetBinContent(eBin, cBin);
2423 if (type == 'r' || type == 'R') vQC2(n) *= vQC2(n); // going back to <<2>>
2424 for (Int_t m = 0; m < fMaxMoment-1; m++) { // cross moment
2425 mNUA(n,m) = CalculateNUAMatrixElement(n, m, type, eBin, cBin);
2426 } // End of cross moment loop
2427 } // End of moment loop
2431 // If determinant is non-zero we go with corrected results
2432 if (det != 0 ) vQC2 = mNUA*vQC2;
2433 else AliWarning(Form("Determinant == 0 - cent: %d-%d, eta: %f, type: '%c', data: %s, vtx: %d-%d%s",
2434 Int_t(cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinLowEdge(cBin)),
2435 Int_t(cumu.Get(type, 2, CumuHistos::kNUAOld)->GetYaxis()->GetBinUpEdge(cBin)),
2436 cumu.Get(type, 2, CumuHistos::kNUAOld)->GetXaxis()->GetBinCenter(eBin),
2437 type, fType.Data(), fVzMin, fVzMax,
2438 GetQCType(fFlags, kTRUE)));
2439 // Go back to v_n for ref. keep <<2'>> for diff. flow).
2440 for (Int_t n = 0; n < fMaxMoment-1; n++) {
2442 if (type == 'r' || type == 'R')
2443 vnTwo = (vQC2(n) > 0. ? TMath::Sqrt(vQC2(n)) : 0.);
2445 // is really more <<2'>> in this case
2448 // Fill in corrected v_n
2449 if (vnTwo != 0) static_cast<TH2D*>(cumu.Get(type, n+2, CumuHistos::kNUA))->Fill(eta, cent, vnTwo);
2450 } // End of moment loop
2451 } // End of eta loop
2452 } // End of centrality loop
2455 //_____________________________________________________________________
2456 Double_t AliForwardFlowTaskQC::VertexBin::CalculateNUAMatrixElement(Int_t n, Int_t m, Char_t type, Int_t binA, Int_t cBin) const
2459 // Calculates the (n,m)-th element in the NUA matrix: 1 if n == m, otherwise:
2460 // <<cos[(n-m)phi1]>>*<<cos[(n-m)phi2]>> + <<sin[(n-m)phi1]>>*<<sin[(n-m)phi2]>>
2461 // NUA(n,m) = -----------------------------------------------------------------------------
2462 // 1 + <<cos[2nphi1]>>*<<cos[2nphi2]>> + <<sin[2nphi1]>>*<<sin[2nphi2]>>
2464 // <<cos[(n+m)phi1]>>*<<cos[(n+m)phi2]>> + <<sin[(n+m)phi1]>>*<<sin[(n+m)phi2]>>
2465 // + -----------------------------------------------------------------------------
2466 // 1 + <<cos[2nphi1]>>*<<cos[2nphi2]>> + <<sin[2nphi1]>>*<<sin[2nphi2]>>
2471 // type: Reference ('r') or differential ('d') or ('a' or 'b')
2472 // binA: eta bin of phi1
2473 // cBin: centrality bin
2477 if (n == m) return 1.;
2481 Double_t cosnMmPhi1 = 0, cosnMmPhi2 = 0, sinnMmPhi1 = 0, sinnMmPhi2 = 0;
2482 Double_t cosnPmPhi1 = 0, cosnPmPhi2 = 0, sinnPmPhi1 = 0, sinnPmPhi2 = 0;
2483 Double_t cos2nPhi1 = 0, cos2nPhi2 = 0, sin2nPhi1 = 0, sin2nPhi2 = 0;
2486 if (type == 'r' || type == 'R') {
2487 if ((fFlags & k3Cor)) {
2488 Double_t eta = fCumuNUARef->GetXaxis()->GetBinCenter(binA);
2490 while (fEtaLims[i] < eta) i++;
2491 Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2492 GetLimits(i-1, aLow, aHigh, bLow, bHigh);
2493 Double_t multA = 0, multB = 0;
2494 for (Int_t a = aLow; a <= aHigh; a++) {
2495 cosnMmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n-m));
2496 sinnMmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n-m));
2497 cosnPmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(n+m));
2498 sinnPmPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(n+m));
2499 cos2nPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberCos(2*n));
2500 sin2nPhi1 += fCumuNUARef->GetBinContent(a, cBin, GetBinNumberSin(2*n));
2501 multA += fCumuNUARef->GetBinContent(a, cBin, 0);
2503 for (Int_t b = bLow; b <= bHigh; b++) {
2504 cosnMmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n-m));
2505 sinnMmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n-m));
2506 cosnPmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(n+m));
2507 sinnPmPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(n+m));
2508 cos2nPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberCos(2*n));
2509 sin2nPhi2 += fCumuNUARef->GetBinContent(b, cBin, GetBinNumberSin(2*n));
2510 multB += fCumuNUARef->GetBinContent(b, cBin, 0);
2512 if (multA == 0 || multB == 0) {
2513 if (fDebug > 0) AliWarning("multA or multB == 0 in matrix elements, aborting NUA");
2516 cosnMmPhi1 /= multA;
2517 sinnMmPhi1 /= multA;
2518 cosnPmPhi1 /= multA;
2519 sinnPmPhi1 /= multA;
2522 cosnMmPhi2 /= multB;
2523 sinnMmPhi2 /= multB;
2524 cosnPmPhi2 /= multB;
2525 sinnPmPhi2 /= multB;
2529 Int_t binB = fCumuNUARef->GetXaxis()->FindBin(-fCumuNUARef->GetXaxis()->GetBinCenter(binA));
2530 cosnMmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2531 sinnMmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2532 cosnPmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2533 sinnPmPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2534 cos2nPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2535 sin2nPhi1 = fCumuNUARef->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2536 cosnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n-m));
2537 sinnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n-m));
2538 cosnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n+m));
2539 sinnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n+m));
2540 cos2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(2*n));
2541 sin2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(2*n));
2543 } // differential flow
2544 else if (type == 'd' || type == 'D') {
2545 Int_t binB = fCumuNUARef->GetXaxis()->FindBin(-fCumuNUADiff->GetXaxis()->GetBinCenter(binA));
2546 cosnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2547 sinnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2548 cosnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2549 sinnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2550 cos2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2551 sin2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2552 cosnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n-m));
2553 sinnMmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n-m));
2554 cosnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(n+m));
2555 sinnPmPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(n+m));
2556 cos2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberCos(2*n));
2557 sin2nPhi2 = fCumuNUARef->GetBinContent(binB, cBin, GetBinNumberSin(2*n));
2558 } // 3 correlator part a or b
2559 else if (type == 'a' || type == 'A' || type == 'b' || type == 'B') {
2560 Double_t mult1 = 0, mult2 = 0;
2562 cosnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n-m));
2563 sinnMmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n-m));
2564 cosnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(n+m));
2565 sinnPmPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(n+m));
2566 cos2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberCos(2*n));
2567 sin2nPhi1 = fCumuNUADiff->GetBinContent(binA, cBin, GetBinNumberSin(2*n));
2568 mult1 = fCumuNUADiff->GetBinContent(binA, cBin, 0);
2570 Double_t eta = fCumuNUADiff->GetXaxis()->GetBinCenter(binA);
2572 while (fEtaLims[i] < eta) i++;
2573 Int_t aLow = 0, aHigh = 0, bLow = 0, bHigh = 0;
2574 GetLimits(i-1, aLow, aHigh, bLow, bHigh);
2575 Int_t lLow = ((type == 'a' || type == 'A') ? aLow : bLow);
2576 Int_t lHigh = ((type == 'a' || type == 'A') ? aHigh : bHigh);
2577 for (Int_t l = lLow; l <= lHigh; l++) {
2578 cosnMmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(n-m));
2579 sinnMmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(n-m));
2580 cosnPmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(n+m));
2581 sinnPmPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(n+m));
2582 cos2nPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberCos(2*n));
2583 sin2nPhi2 += fCumuNUARef->GetBinContent(l, cBin, GetBinNumberSin(2*n));
2584 mult2 += fCumuNUARef->GetBinContent(l, cBin, 0);
2586 if (mult1 == 0 || mult2 == 0) {
2587 if (fDebug > 0) AliWarning("mult1 or mult2 == 0 in matrix elements, aborting NUA");
2590 cosnMmPhi1 /= mult1;
2591 sinnMmPhi1 /= mult1;
2592 cosnPmPhi1 /= mult1;
2593 sinnPmPhi1 /= mult1;
2596 cosnMmPhi2 /= mult2;
2597 sinnMmPhi2 /= mult2;
2598 cosnPmPhi2 /= mult2;
2599 sinnPmPhi2 /= mult2;
2604 // Actual calculation
2605 Double_t e = cosnMmPhi1*cosnMmPhi2 + sinnMmPhi1*sinnMmPhi2 + cosnPmPhi1*cosnPmPhi2 + sinnPmPhi1*sinnPmPhi2;
2606 Double_t den = 1 + cos2nPhi1*cos2nPhi2 + sin2nPhi1*sin2nPhi2;
2607 if (den != 0) e /= den;
2612 //_____________________________________________________________________
2613 void AliForwardFlowTaskQC::VertexBin::AddVertexBins(CumuHistos& cumu, TList* list, UInt_t nNUA) const
2616 // Add up vertex bins with flow results
2619 // cumu: CumuHistos object with vtxbin results
2620 // list: Outout list with added results
2621 // nNUA: # of NUA histograms to loop over
2624 TProfile2D* avgProf = 0;
2626 Int_t nT = ((fFlags & k3Cor) ? 3 : 2);
2628 for (UInt_t nua = 0; nua <= nNUA; nua++) { // NUA loop
2629 for (Int_t n = 2; n <= fMaxMoment; n++) { // moment loop
2630 for (Int_t t = 0; t < nT; t++) { // type loop (r/d/a/b)
2633 case 0: ct = 'r'; break;
2634 case 1: ct = ((fFlags & k3Cor) ? 'a' : 'd'); break;
2635 case 2: ct = 'b'; break;
2636 default: ct = '\0'; break;
2638 vtxHist = static_cast<TH2D*>(cumu.Get(ct,n,nua));
2640 AliWarning("VertexBin::AddVertexBins: vtxHist not found!");
2643 name = vtxHist->GetName();
2644 // Strip name of vtx info
2645 Ssiz_t l = name.Last('x')-3;
2647 avgProf = (TProfile2D*)list->FindObject(name.Data());
2648 // if no output profile yet, make one
2650 avgProf = new TProfile2D(name.Data(), name.Data(),
2651 vtxHist->GetNbinsX(), vtxHist->GetXaxis()->GetXmin(), vtxHist->GetXaxis()->GetXmax(),
2652 vtxHist->GetNbinsY(), vtxHist->GetYaxis()->GetXmin(), vtxHist->GetYaxis()->GetXmax());
2653 if (vtxHist->GetXaxis()->IsVariableBinSize())
2654 avgProf->GetXaxis()->Set(vtxHist->GetNbinsX(), vtxHist->GetXaxis()->GetXbins()->GetArray());
2655 if (vtxHist->GetYaxis()->IsVariableBinSize())
2656 avgProf->GetYaxis()->Set(vtxHist->GetNbinsY(), vtxHist->GetYaxis()->GetXbins()->GetArray());
2659 // Fill in, cannot be done with Add function.
2660 for (Int_t e = 1; e <= vtxHist->GetNbinsX(); e++) { // eta loop
2661 Double_t eta = vtxHist->GetXaxis()->GetBinCenter(e);
2662 for (Int_t c = 1; c <= vtxHist->GetNbinsY(); c++) { // cent loop
2663 Double_t cent = vtxHist->GetYaxis()->GetBinCenter(c);
2664 Double_t cont = vtxHist->GetBinContent(e, c);
2665 if (cont == 0) continue;
2666 avgProf->Fill(eta, cent, cont);
2667 } // End of cent loop
2668 } // End of eta loop
2669 } // End of type loop
2670 } // End of moment loop
2671 } // End of nua loop
2673 //_____________________________________________________________________
2674 Int_t AliForwardFlowTaskQC::VertexBin::GetBinNumberCos(Int_t n) const
2677 // Get the bin number of <<cos(nphi)>>
2682 // Return: bin number
2687 if (n == 0) bin = fMaxMoment*4-1;
2692 //_____________________________________________________________________
2693 Int_t AliForwardFlowTaskQC::VertexBin::GetBinNumberSin(Int_t n) const
2696 // Get the bin number of <<sin(nphi)>>
2701 // Return: bin number
2706 if (n == 0) bin = fMaxMoment*4;
2711 //_____________________________________________________________________
2712 void AliForwardFlowTaskQC::VertexBin::SetupNUALabels(TAxis* a) const
2715 // Setup NUA labels on axis
2718 // a: Axis to set up
2720 if (a->GetNbins() != GetBinNumberSin()) AliFatal("SetupNUALabels: Wrong number of bins on axis");
2723 while (i <= a->GetNbins()) {
2724 a->SetBinLabel(i++, Form("<<cos(%d#varphi)>>", j));
2725 a->SetBinLabel(i++, Form("<<sin(%d#varphi)>>", j++));
2730 //_____________________________________________________________________
2731 TH2I* AliForwardFlowTaskQC::VertexBin::MakeQualityHist(const Char_t* name) const
2734 // Add a histogram for checking the analysis quality
2737 // const char*: name of data type
2739 // Return: histogram for tracking successful calculations
2741 Int_t nBins = ((fFlags & kStdQC) ? 8 : 4);
2742 TH2I* quality = new TH2I(name, name, (fMaxMoment-1)*nBins, 1, (fMaxMoment-1)*nBins+1, fCumuNUARef->GetNbinsY(),
2743 fCumuNUARef->GetYaxis()->GetXmin(), fCumuNUARef->GetYaxis()->GetXmax());
2744 quality->GetYaxis()->Set(fCumuNUARef->GetNbinsY(), fCumuNUARef->GetYaxis()->GetXbins()->GetArray());
2745 for (Int_t i = 2, j = 1; i <= fMaxMoment; i++) {
2746 quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{2} > 0", i));
2747 quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{2} <= 0", i));
2748 quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{2} > 0", i));
2749 quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{2} <= 0", i));
2750 if ((fFlags & kStdQC)) {
2751 quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{4} < 0", i));
2752 quality->GetXaxis()->SetBinLabel(j++, Form("QC_{%d}{4} >= 0", i));
2753 quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{4} < 0", i));
2754 quality->GetXaxis()->SetBinLabel(j++, Form("QC'_{%d}{4} >= 0", i));
2760 //_____________________________________________________________________
2761 TH2D* AliForwardFlowTaskQC::VertexBin::MakeOutputHist(Int_t qc, Int_t n, const Char_t* ctype, UInt_t nua) const
2764 // Setup a TH2D for the output
2767 // qc # of particle correlations
2769 // ref Type: r/d/a/b
2770 // nua For nua corrected hists?
2772 // Return: 2D hist for results
2774 Bool_t ref = ((ctype[0] == 'R') || (ctype[0] == 'r'));
2775 TAxis* xAxis = (ref ? fCumuNUARef->GetXaxis() : fCumuNUADiff->GetXaxis());
2776 TAxis* yAxis = (ref ? fCumuNUARef->GetYaxis() : fCumuNUADiff->GetYaxis());
2779 case CumuHistos::kNoNUA: ntype = "Un"; break;
2780 case CumuHistos::kNUAOld: ntype = "NUAOld"; break;
2781 case CumuHistos::kNUA: ntype = "NUA"; break;
2784 TH2D* h = new TH2D(Form("%sQC%d_v%d_%s_%sCorr%s_vtx_%d_%d",
2785 fType.Data(), qc, n, ctype, ntype.Data(),
2786 GetQCType(fFlags), fVzMin, fVzMax),
2787 Form("%sQC%d_v%d_%s_%sCorr%s_vtx_%d_%d",
2788 fType.Data(), qc, n, ctype, ntype.Data(),
2789 GetQCType(fFlags), fVzMin, fVzMax),
2790 xAxis->GetNbins(), xAxis->GetXmin(), xAxis->GetXmax(),
2791 yAxis->GetNbins(), yAxis->GetXmin(), yAxis->GetXmax());
2792 if (xAxis->IsVariableBinSize()) h->GetXaxis()->Set(xAxis->GetNbins(), xAxis->GetXbins()->GetArray());
2793 h->GetYaxis()->Set(yAxis->GetNbins(), yAxis->GetXbins()->GetArray());
2797 //_____________________________________________________________________
2798 void AliForwardFlowTaskQC::PrintFlowSetup() const
2801 // Print the setup of the flow task
2803 Printf("=======================================================");
2804 Printf("%s::Print", this->IsA()->GetName());
2805 Printf("Forward detector: :\t%s", ((fFlowFlags & kFMD) ? "FMD" : "VZERO"));
2806 Printf("Number of bins in vertex axis :\t%d", fVtxAxis->GetNbins());
2807 Printf("Range of vertex axis :\t[%3.1f,%3.1f]",
2808 fVtxAxis->GetXmin(), fVtxAxis->GetXmax());
2809 Printf("Number of bins in centrality axis :\t%d", fCentAxis->GetNbins());
2810 printf("Centrality binning :\t");
2811 for (Int_t cBin = 1; cBin <= fCentAxis->GetNbins(); cBin++) {
2812 printf("%02d-%02d%% ", Int_t(fCentAxis->GetBinLowEdge(cBin)), Int_t(fCentAxis->GetBinUpEdge(cBin)));
2813 if (cBin == fCentAxis->GetNbins()) printf("\n");
2814 else if (cBin % 4 == 0) printf("\n\t\t\t\t\t");
2816 printf("Doing flow analysis for :\t");
2817 for (Int_t n = 2; n <= fMaxMoment; n++) printf("v%d ", n);
2819 TString type = "Standard QC{2} and QC{4} calculations.";
2820 if ((fFlowFlags & kEtaGap)) type = "QC{2} with a rapidity gap.";
2821 if ((fFlowFlags & k3Cor)) type = "QC{2} with 3 correlators.";
2822 if ((fFlowFlags & kTPC) == kTPC) type.ReplaceAll(".", " with TPC tracks for reference.");
2823 if ((fFlowFlags & kHybrid) == kHybrid) type.ReplaceAll(".", " with hybrid tracks for reference.");
2824 Printf("QC calculation type :\t%s", type.Data());
2825 Printf("Symmetrize ref. flow wrt. eta = 0 :\t%s", ((fFlowFlags & kSymEta) ? "true" : "false"));
2826 Printf("Apply NUA correction terms :\t%s", ((fFlowFlags & kNUAcorr) ? "true" : "false"));
2827 Printf("Satellite vertex flag :\t%s", ((fFlowFlags & kSatVtx) ? "true" : "false"));
2828 Printf("FMD sigma cut: :\t%f", fFMDCut);
2829 Printf("SPD sigma cut: :\t%f", fSPDCut);
2830 if ((fFlowFlags & kEtaGap) || (fFlowFlags & kTracks))
2831 Printf("Eta gap: :\t%f", fEtaGap);
2832 Printf("=======================================================");
2834 //_____________________________________________________________________
2835 const Char_t* AliForwardFlowTaskQC::GetQCType(UShort_t flags, Bool_t prependUS)
2838 // Get the type of the QC calculations
2839 // Used for naming of objects in the VertexBin class,
2840 // important to avoid memory problems when running multiple
2841 // initializations of the class with different setups
2844 // flags: Flow flags for type determination
2845 // prependUS: Prepend an underscore (_) to the name
2847 // Return: QC calculation type
2850 if ((flags & kStdQC)) type = "StdQC";
2851 else if ((flags & kEtaGap)) type = "EtaGap";
2852 else if ((flags & k3Cor)) type = "3Cor";
2853 else type = "UNKNOWN";
2854 if (prependUS) type.Prepend("_");
2855 if ((flags & kTPC) == kTPC) type.Append("TPCTr");
2856 if ((flags & kHybrid) == kHybrid) type.Append("HybTr");
2860 //_____________________________________________________________________
2861 TH1* AliForwardFlowTaskQC::CumuHistos::Get(Char_t t, Int_t n, UInt_t nua) const
2864 // Get histogram/profile for type t and moment n
2867 // t: type = 'r'/'d'
2872 if (n < 0) AliFatal(Form("CumuHistos out of range: (%c,%d)", t, n));
2876 if (t == 'r' || t == 'R') pos = n;
2877 else if (t == 'd' || t == 'D') pos = n;
2878 else if (t == 'a' || t == 'A') pos = 2*n;
2879 else if (t == 'b' || t == 'B') pos = 2*n+1;
2880 else AliFatal(Form("CumuHistos wrong type: %c", t));
2882 if (t == 'r' || t == 'R') {
2883 if (pos < fRefHists->GetEntries()) {
2884 h = (TH1*)fRefHists->At(pos);
2886 } else if (pos < fDiffHists->GetEntries()) {
2887 h = (TH1*)fDiffHists->At(pos);
2889 if (!h) AliFatal(Form("No hist found in list %c at pos %d", t, pos));
2893 //_____________________________________________________________________
2894 void AliForwardFlowTaskQC::CumuHistos::ConnectList(TString name, TList* l)
2897 // Connect an output list with this object
2901 // l: list to keep outputs in
2904 ref.ReplaceAll("Cumu","CumuRef");
2905 fRefHists = (TList*)l->FindObject(ref.Data());
2907 fRefHists = new TList();
2908 fRefHists->SetName(ref.Data());
2911 // Check that the list correspond to fMaxMoments
2912 if (fRefHists->GetEntries() != (fMaxMoment-1.)*(fNUA+1))
2913 AliError("CumuHistos::ConnectList Wrong number of hists in ref list,"
2914 "you are doing something wrong!");
2916 TString diff = name;
2917 diff.ReplaceAll("Cumu","CumuDiff");
2918 fDiffHists = (TList*)l->FindObject(diff.Data());
2920 fDiffHists = new TList();
2921 fDiffHists->SetName(diff.Data());
2924 // Check that the list correspond to fMaxMoment
2925 if ((fDiffHists->GetEntries() != (fMaxMoment-1.)*(fNUA+1)) &&
2926 (fDiffHists->GetEntries() != 2*(fMaxMoment-1.)*(fNUA+1)))
2927 AliError(Form("CumuHistos::ConnectList Wrong number of hists in diff list,"
2928 "you are doing something wrong! (%s)",name.Data()));
2932 //_____________________________________________________________________
2933 void AliForwardFlowTaskQC::CumuHistos::Add(TH1* h) const
2936 // Add a histogram to this objects lists
2939 // h: histogram/profile to add
2941 TString name = h->GetName();
2942 if (name.Contains("Ref")) {
2943 if (fRefHists) fRefHists->Add(h);
2944 else AliFatal("CumuHistos::Add() - fRefHists does not exist");
2945 // Check that the list correspond to fMaxMoments
2946 if (fRefHists->GetEntries() > (fNUA+1)*(fMaxMoment-1.))
2947 AliError("CumuHistos::Add wrong number of hists in ref list, "
2948 "you are doing something wrong!");
2950 else if (name.Contains("Diff")) {
2951 if (fDiffHists) fDiffHists->Add(h);
2952 else AliFatal("CumuHistos::Add() - fDiffHists does not exist");
2953 // Check that the list correspond to fMaxMoment
2954 if (fDiffHists->GetEntries() > 2*(fNUA+1)*(fMaxMoment-1.))
2955 AliError("CumuHistos::Add wrong number of hists in diff list, "
2956 "you are doing something wrong!");
2960 //_____________________________________________________________________
2961 Int_t AliForwardFlowTaskQC::CumuHistos::GetPos(Int_t n, UInt_t nua) const
2964 // Get position in list of moment n objects
2965 // To take care of different numbering schemes
2969 // nua: # of nua corrections applied
2971 // Return: position #
2973 if (n > fMaxMoment) return -1;
2974 else return (n-2)+nua*(fMaxMoment-1);
2976 //_____________________________________________________________________