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