Cluster improvements + drawing of digits like TBox.
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowAnalyser.cxx
CommitLineData
9777bfcb 1//////////////////////////////////////////////////////////////////////
2//
3// $Id$
4//
5// Author: Emanuele Simili
6//
7//////////////////////////////////////////////////////////////////////
8//_____________________________________________________________
9//
10// Description:
11// the AliFlowAnalyser provides the methods to perform an Event
12// plane flow analysis over AliFlowEvents.
13// - The method Init() generates the histograms.
14// - The method Analyse(AliFlowEvent*) calculates/extracts flow observables,
15// fills some histograms and performs the E.P. analysis.
16// - The method Resolution() calculates the resolution correction and
17// applies it to observed v_n values.
18// - The method Finish() saves the histograms and closes the file.
19// Weights calculation has been moved to the class AliFlowWeighter.
20//
21// The most important features of the Flow Analysis are:
22// - Reaction Plane calculation: for different harmonics and selections
23// with user defined cuts. (The calculation of the R.P. is a method in
24// the AliFlowEvent class).
25// - Correlation Analysis: particles are related to the R.P. and Vn
26// coefficient are extracted. The AliFlowAnalysisMaker also removes
27// auto-correlations of each particle with respect to the event plane,
28// by subtracting track by track from the Q vector before the correlation
29// study. Observable Flow is plotted.
30// - Differential Flow: for each harmonic and selection, flow values
31// (Vn) are plotted versus rapidity/pseudorapidity and pT.
32// - Resulution Estimate: event plane resolution is estimated via the
33// sub-event method (res ~ sqrt(cos(psi1-psi2))), and flow values are
34// then corrected for it. "True" Flow is plotted.
35// - Integrated Flow: for the doubly-integrated (pT & eta) Vn values
36// the error on the event plane resolution is folded in to the error.
37// - Phi and Baiesian weights calculation (from event observables).
38//
39// The AliFlowAnalysisMaker Class is adapted from the original
40// StFlowAnalysisMaker, succesfully used to study flow in the
41// STAR experiment at RICH.
42// Original Authors: Raimond Snellings & Art Poskanzer
43//
44//////////////////////////////////////////////////////////////////////
45
46#ifndef ALIFLOWANALYSER_CXX
47#define ALIFLOWANALYSER_CXX
48
49
50// ROOT things
51#include <TROOT.h>
52#include <TMath.h>
53#include <TFile.h>
54#include <TTree.h>
55#include <TChain.h>
56#include <TString.h>
57#include <TObject.h>
58#include <TObjArray.h>
59#include <TOrdCollection.h>
60#include <TParticle.h>
61#include <TParticlePDG.h>
62#include <TDatabasePDG.h>
63#include <TH1.h>
64#include <TH2.h>
65#include <TH3.h>
66#include <TProfile.h>
67#include <TProfile2D.h>
68#include <TVector2.h>
69
70// Flow things
71#include "AliFlowEvent.h"
72#include "AliFlowTrack.h"
73#include "AliFlowV0.h"
74#include "AliFlowConstants.h"
75#include "AliFlowSelection.h"
76#include "AliFlowAnalyser.h"
77
78// ANSI things
79#include <math.h>
80#include <stdio.h>
81#include <stdlib.h>
82#include <iostream>
83using namespace std; //required for resolving the 'cout' symbol
84
da4d7f17 85ClassImp(AliFlowAnalyser)
9777bfcb 86//-----------------------------------------------------------------------
87AliFlowAnalyser::AliFlowAnalyser(const AliFlowSelection* flowSelect)
88{
89 // default constructor (selection given or default selection)
90
91 if(flowSelect) { fFlowSelect = new AliFlowSelection(*flowSelect) ; }
92 else { fFlowSelect = new AliFlowSelection() ; }
93
94 // output file (histograms)
95 fHistFileName = "flowAnalysPlot.root" ;
96 fHistFile = 0 ;
97
98 // for phi weights
99 fPhiWgtFile = 0 ;
100 fPhiBins = Flow::nPhiBins ;
101 fPhiMin = 0.;
102 fPhiMax = 2*TMath::Pi() ;
103
104 // flags
105 fTrackLoop = kTRUE ; // main loop for charged tracks
106 fV0loop = kTRUE ; // correlation analysis is done also for neutral secundary vertex
107 fShuffle = kFALSE ; // randomly reshuffles tracks
108 fV1Ep1Ep2 = kFALSE ; // disabled by default
109 fEtaSub = kFALSE ; // eta subevents
110 fReadPhiWgt = kFALSE ; // if kTRUE & if flowPhiWgt file is there -> Phi Weights are used
111 fBayWgt = kFALSE ; // Bayesian Weights for P.Id. used
112 fRePid = kFALSE ; // recalculates the P.Id. (becomes kTRUE if new bayesian weights are plugged in)
113 fPtWgt = kFALSE ; // pT as a weight
114 fEtaWgt = kFALSE ; // eta as a weight
115 fOnePhiWgt = kTRUE ; // one/three phi-wgt histogram(s)
116// fMaxLabel = 1000 ; // for labels histogram
117
118 fPhiWgtHistList = 0 ;
119 fVnResHistList = 0 ;
120}
121//-----------------------------------------------------------------------
122AliFlowAnalyser::~AliFlowAnalyser()
123{
124 // default distructor (no actions)
125}
126//-----------------------------------------------------------------------
127Bool_t AliFlowAnalyser::Init()
128{
129 // sets some defaults for the analysis
130
131 cout << "* FlowAnalysis * - Init()" << endl ; cout << endl ;
132
133 // Open output files (->plots)
134 fHistFile = new TFile(fHistFileName.Data(), "RECREATE");
135 fHistFile->cd() ;
136
137 // counters and pointers to 0
138 fEventNumber = 0 ;
139 fTrackNumber = 0 ;
140 fV0Number = 0 ;
141 fPidId = -1 ;
142 //fNumberOfEvents = 0 ;
143 fNumberOfTracks = 0 ;
144 fNumberOfV0s = 0 ;
145 fFlowEvent = 0 ;
146 fFlowTrack = 0 ;
147 fFlowV0 = 0 ;
148 fFlowTracks = 0 ;
149 fFlowV0s = 0 ;
150 for(Int_t ii=0;ii<3;ii++) { fVertex[ii] = 0 ; }
151
152 // Check for Reading weights: if weight file is not there, wgt arrays are not used (phiwgt & bayesian)
153 if(!fPhiWgtFile)
154 {
155 fReadPhiWgt = kFALSE ;
156 fBayWgt = kFALSE ;
157 cout << " No wgt file. Phi & Bayesian weights will not be used ! " << endl ;
158 }
159
160 // Histogram settings
161 fLabel = "Pseudorapidity" ; if(strlen(fFlowSelect->PidPart()) != 0) { fLabel = "Rapidity"; }
162 fPtMaxPart = Flow::fPtMaxPart ; if(fFlowSelect->PtMaxPart()) { fPtMaxPart = fFlowSelect->PtMaxPart() ; }
163 fPtBinsPart = Flow::nPtBinsPart ; if(fFlowSelect->PtBinsPart()) { fPtBinsPart = fFlowSelect->PtBinsPart() ; }
164 fPtMin = Flow::fPtMin ;
165 fPtMax = Flow::fPtMax ;
166 fEtaMin = Flow::fEtaMin ;
167 fEtaMax = Flow::fEtaMax ;
168 fPtBins = Flow::nPtBins ;
169 fEtaBins = Flow::nEtaBins ;
170 // -
171 SetPtRangevEta(1.,0.); // (Flow::fPtMin , Flow::fPtMax) ;
172 SetEtaRangevPt(1.,0.); // (Flow::fEtaMin , Flow::fEtaMax) ;
173 // -
174 const float triggerMin = -1.5;
175 const float triggerMax = 10.5;
176 const float chargeMin = -2.5;
177 const float chargeMax = 2.5;
178 const float dcaMin = 0.0;
179 const float dcaMax = 10.0; // 0.5;
180 const float glDcaMax = 50.0;
181 const float chi2Min = 0.;
182 const float chi2Max = 500.;
183 const float chi2normMax = 50.;
184 const float chi2MaxC = 30.;
185 const float chi2MaxI = 60.;
186 const float fitPtsMin = -0.5;
187 const float fitPtsMaxTpc = 160.5;
188 const float fitPtsMaxIts = 6.5;
189 const float fitPtsMaxTrd = 130.5;
190 const float fitPtsMaxTof = 1.5;
191 const float fitOverMaxMin = 0.;
192 const float fitOverMaxMax = 1.1;
193 const float multOverOrigMin = 0.;
194 const float multOverOrigMax = 1.;
195 const float vertexZMin = -10.;
196 const float vertexZMax = 10.;
197 const float vertexXYMin = -0.1;
198 const float vertexXYMax = 0.1;
199 const float etaSymMinPart = -1.;
200 const float etaSymMaxPart = 1.;
201 const float etaSymMin = -2.;
202 const float etaSymMax = 2.;
203 const float psiMin = 0.;
204 const float psiMax = 2*TMath::Pi() ;
205 const float multMin = 0. ;
206 const float multMax =25000. ;
207 const float multV0 = 5000. ;
208 const float qMin = 0. ;
209 const float qMax = 20. ;
210 const float pidMin = 0. ;
211 const float pidMax = 1. ;
212 const float centMin = -0.5 ;
213 const float centMax = 9.5 ;
214 const float logpMin = -2.5 ;
215 const float logpMax = 2.5 ;
216 const float pMin = 0. ;
217 const float pMax = 100. ;
218 const float dEdxMax = 1000. ;
219 const float dEdxMaxTPC = 1500. ;
220 const float TOFmin = 10000. ;
221 const float TOFmax = 50000. ;
222 const float TRDmax = 2000. ;
223 const float lgMin = -1000. ;
224 const float lgMax = 1000. ;
225 const float lgMinV0 = 0. ;
226 const float lgMaxV0 = 50. ;
227 const float massMin = -0.1 ;
228 const float massMax = 2.1 ;
229 const float ZDCpartMin = -0.5 ;
230 const float ZDCpartMax = 1000.5 ;
231 const float ZDCeMin = 0. ;
232 const float ZDCeMax = 10000. ;
233 // -
234 enum { nTriggerBins = 12,
235 nChargeBins = 5,
236 nDcaBins = 1000,
237 nChi2Bins = 150,
238 nFitPtsBinsTpc = 161,
239 nFitPtsBinsIts = 7,
240 nFitPtsBinsTrd = 131,
241 nFitPtsBinsTof = 2,
242 nFitOverMaxBins = 55,
243 nMultOverOrigBins =100,
244 nVertexZBins = 200,
245 nVertexXYBins = 1000,
246 nEtaSymBins = 200,
247 nPhi3DBins = 60,
248 nPsiBins = 36,
249 nMultBins = 250,
250 nPidBins = 100,
251 nCentBins = 10,
252 nDedxBins = 1000,
253 nTofBins = 1000,
254 nTrdBins = 100,
255 nMomenBins = 500,
256 nQbins = 50,
257 nLgBins = 500,
258 nMassBins = 2200,
259 nZDCpartBins = 1001,
260 nZDCeBins = 200
261 } ;
262
263 // Histograms Booking ...
264 // Trigger
265 fHistTrigger = new TH1F("Flow_Trigger", "Flow_Trigger", nTriggerBins, triggerMin, triggerMax);
266 fHistTrigger->SetXTitle("Trigger: 1 minbias, 2 central, 3 laser, 10 other");
267 fHistTrigger->SetYTitle("Counts");
268
269 // Charge
270 fHistCharge = new TH1F("Flow_Charge", "Flow_Charge", nChargeBins, chargeMin, chargeMax);
271 fHistCharge->SetXTitle("Charge");
272 fHistCharge->SetYTitle("Counts");
273
274 // Reconstructed Momentum (constrained, if so)
275 fHistPtot = new TH1F("Flow_totalP","Flow_totalP", fPtBins, fPtMin, fPtMax);
276 fHistPtot->Sumw2();
277 fHistPtot->SetXTitle("P (GeV/c)");
278 fHistPtot->SetYTitle("Counts");
279
280 // Reconstructed pT (constrained, if so)
281 fHistPt = new TH1F("Flow_Pt","Flow_Pt", nMomenBins, pMin, pMax);
282 fHistPt->Sumw2();
283 fHistPt->SetXTitle("Pt (GeV/c)");
284 fHistPt->SetYTitle("Counts");
285
286 // Reconstructed P.id vs pT
287 fHistPidPt = new TH2F("Flow_PidPt","Flow_PidPt", 12, 0., 12., fPtBins, fPtMin, fPtMax);
288 fHistPidPt->Sumw2();
289 fHistPidPt->SetXTitle("e-, e+, mu-, mu+, pi-, pi+, K-, K+, pr-, pr+, d-, d+");
290 fHistPidPt->SetYTitle("Pt (GeV/c)");
291
292 // Distance of closest approach - Transverse, Unsigned
293 fHistDca = new TH1F("Flow_TransDCA", "Flow_TransverseDCA", nDcaBins, dcaMin, dcaMax);
294 fHistDca->Sumw2();
295 fHistDca->SetXTitle("|Track's Signed dca (cm)|");
296 fHistDca->SetYTitle("Counts");
297
298 // Distance of closest approach - Transverse
299 fHistTransDca = new TH1F("Flow_TransDCA_Signed", "Flow_TransverseDCA_Signed", nDcaBins, -dcaMax, dcaMax);
300 fHistTransDca->Sumw2();
301 fHistTransDca->SetXTitle("Track's Transverse dca (cm)");
302 fHistTransDca->SetYTitle("Counts");
303
304 // Distance of closest approach - 3d
305 fHistDcaGlobal = new TH1F("Flow_3dDcaGlobal", "Flow_3dDcaGlobal", nDcaBins, dcaMin, glDcaMax);
306 fHistDcaGlobal->Sumw2();
307 fHistDcaGlobal->SetXTitle("|3d Global Track's dca to Vertex (cm)|");
308 fHistDcaGlobal->SetYTitle("Counts");
309
310 // Distance of closest approach for particles correlated with the event plane - 3d
311 fHistDcaGlobalPart = new TH1F("Flow_3dDcaGlobalPart", "Flow_3dDcaGlobalPart", nDcaBins, dcaMin, glDcaMax);
312 fHistDcaGlobalPart->Sumw2();
313 fHistDcaGlobalPart->SetXTitle("|3d Global Track's dca to Vertex (cm)|");
314 fHistDcaGlobalPart->SetYTitle("Counts");
315
316 // Distance of closest approach for particles NOT SELECTED for correlation with the event plane - 3d
317 fHistDcaGlobalOut = new TH1F("Flow_3dDcaGlobalOut", "Flow_3dDcaGlobalOut", nDcaBins, dcaMin, glDcaMax);
318 fHistDcaGlobalOut->Sumw2();
319 fHistDcaGlobalOut->SetXTitle("|3d Global Track's dca to Vertex (cm)|");
320 fHistDcaGlobalOut->SetYTitle("Counts");
321
322 // Yield Pt-Phi for constrained tracks
323 fHistPhiPtCon = new TH2D("Flow_PtPhi_Con", "Flow_PtPhi_Con", fPhiBins, fPhiMin, fPhiMax, fPtBins, fPtMin, fPtMax);
324 fHistPhiPtCon->Sumw2();
325 fHistPhiPtCon->SetXTitle("Phi");
326 fHistPhiPtCon->SetYTitle("Pt (GeV/c)");
327
328 // Yield Pt-Phi for UNconstrained tracks
329 fHistPhiPtUnc = new TH2D("Flow_PtPhi_Unc", "Flow_PtPhi_Unc", fPhiBins, fPhiMin, fPhiMax, fPtBins, fPtMin, fPtMax);
330 fHistPhiPtUnc->Sumw2();
331 fHistPhiPtUnc->SetXTitle("Phi");
332 fHistPhiPtUnc->SetYTitle("Pt (GeV/c)");
333
334 // Chi2
335 // at the main vertex
336 fHistChi2 = new TH1F("Flow_Chi2", "Flow_Chi2", nChi2Bins, chi2Min, chi2MaxC);
337 fHistChi2->SetXTitle("Chi square at Main Vertex");
338 fHistChi2->SetYTitle("Counts");
339 // TPC
340 fHistChi2TPC = new TH1F("Flow_Chi2_TPC", "Flow_Chi2_TPC", nChi2Bins, chi2Min, chi2Max);
341 fHistChi2TPC->SetXTitle("Chi square for TPC");
342 fHistChi2TPC->SetYTitle("Counts");
343 // ITS
344 fHistChi2ITS = new TH1F("Flow_Chi2_ITS", "Flow_Chi2_ITS", nChi2Bins, chi2Min, chi2MaxI);
345 fHistChi2ITS->SetXTitle("Chi square for ITS");
346 fHistChi2ITS->SetYTitle("Counts");
347 // TRD
348 fHistChi2TRD = new TH1F("Flow_Chi2_TRD", "Flow_Chi2_TRD", nChi2Bins, chi2Min, chi2Max);
349 fHistChi2TRD->SetXTitle("Chi square for TRD");
350 fHistChi2TRD->SetYTitle("Counts");
351 // TOF
352 fHistChi2TOF = new TH1F("Flow_Chi2_TOF", "Flow_Chi2_TOF", nChi2Bins, chi2Min, chi2Max);
353 fHistChi2TOF->SetXTitle("Chi square for TOF");
354 fHistChi2TOF->SetYTitle("Counts");
355
356 // Chi2 normalized
357 // TPC
358 fHistChi2normTPC = new TH1F("Flow_Chi2norm_TPC", "Flow_Chi2norm_TPC", nChi2Bins, chi2Min, chi2normMax);
359 fHistChi2normTPC->SetXTitle("Normalized #Chi^{2} for TPC");
360 fHistChi2normTPC->SetYTitle("Counts");
361 // ITS
362 fHistChi2normITS = new TH1F("Flow_Chi2norm_ITS", "Flow_Chi2norm_ITS", nChi2Bins, chi2Min, chi2normMax);
363 fHistChi2normITS->SetXTitle("Normalized #Chi^{2} for ITS");
364 fHistChi2normITS->SetYTitle("Counts");
365 // TRD
366 fHistChi2normTRD = new TH1F("Flow_Chi2norm_TRD", "Flow_Chi2norm_TRD", nChi2Bins, chi2Min, chi2normMax);
367 fHistChi2normTRD->SetXTitle("Normalized #Chi^{2} for TRD");
368 fHistChi2normTRD->SetYTitle("Counts");
369 // TOF
370 fHistChi2normTOF = new TH1F("Flow_Chi2norm_TOF", "Flow_Chi2norm_TOF", nChi2Bins, chi2Min, chi2normMax);
371 fHistChi2normTOF->SetXTitle("Normalized #Chi^{2} for TOF");
372 fHistChi2normTOF->SetYTitle("Counts");
373
374 // FitPts
375 // TPC
376 fHistFitPtsTPC = new TH1F("Flow_FitPts_TPC", "Flow_FitPts_TPC", nFitPtsBinsTpc, fitPtsMin, fitPtsMaxTpc);
377 fHistFitPtsTPC->SetXTitle("Fit Points");
378 fHistFitPtsTPC->SetYTitle("Counts");
379 // ITS
380 fHistFitPtsITS = new TH1F("Flow_HitPts_ITS", "Flow_HitPts_ITS", nFitPtsBinsIts, fitPtsMin, fitPtsMaxIts);
381 fHistFitPtsITS->SetXTitle("Fit Points");
382 fHistFitPtsITS->SetYTitle("Counts");
383 // TRD
384 fHistFitPtsTRD = new TH1F("Flow_FitPts_TRD", "Flow_FitPts_TRD", nFitPtsBinsTrd, fitPtsMin, fitPtsMaxTrd);
385 fHistFitPtsTRD->SetXTitle("Fit Points");
386 fHistFitPtsTRD->SetYTitle("Counts");
387 // TOF
388 fHistFitPtsTOF = new TH1F("Flow_HitPts_TOF", "Flow_HitPts_TOF", nFitPtsBinsTof, fitPtsMin, fitPtsMaxTof);
389 fHistFitPtsTOF->SetXTitle("Fit Points");
390 fHistFitPtsTOF->SetYTitle("Counts");
391
392 // MaxPts
393 // TPC
394 fHistMaxPtsTPC = new TH1F("Flow_MaxPts_TPC", "Flow_MaxPts_TPC", nFitPtsBinsTpc, fitPtsMin, fitPtsMaxTpc);
395 fHistMaxPtsTPC->SetXTitle("Max Points");
396 fHistMaxPtsTPC->SetYTitle("Counts");
397 // ITS
398 fHistMaxPtsITS = new TH1F("Flow_MaxPts_ITS", "Flow_MaxPts_ITS", nFitPtsBinsIts, fitPtsMin, fitPtsMaxIts);
399 fHistMaxPtsITS->SetXTitle("Max Points");
400 fHistMaxPtsITS->SetYTitle("Counts");
401 // TRD
402 fHistMaxPtsTRD = new TH1F("Flow_MaxPts_TRD", "Flow_MaxPts_TRD", nFitPtsBinsTrd, fitPtsMin, fitPtsMaxTrd);
403 fHistMaxPtsTRD->SetXTitle("Max Points");
404 fHistMaxPtsTRD->SetYTitle("Counts");
405 // TOF
406 fHistMaxPtsTOF = new TH1F("Flow_MaxPts_TOF", "Flow_MaxPts_TOF", nFitPtsBinsTof, fitPtsMin, fitPtsMaxTof);
407 fHistMaxPtsTOF->SetXTitle("Max Points");
408 fHistMaxPtsTOF->SetYTitle("Counts");
409
410 // FitOverMax
411 // Tpc
412 fHistFitOverMaxTPC = new TH1F("Flow_FitOverMax_TPC", "Flow_FitOverMax_TPC", nFitOverMaxBins, fitOverMaxMin, fitOverMaxMax);
413 fHistFitOverMaxTPC->SetXTitle("(Fit Points - 1) / Max Points");
414 fHistFitOverMaxTPC->SetYTitle("Counts");
415 // All
416 fHistFitOverMax = new TH1F("Flow_FitOverMax", "Flow_FitOverMax", nFitOverMaxBins, fitOverMaxMin, fitOverMaxMax);
417 fHistFitOverMax->SetXTitle("(Fit Points - 1) / Max Points");
418 fHistFitOverMax->SetYTitle("Counts");
419
420 // lenght
421 fHistLenght = new TH1F("Flow_TrackLenght", "Flow_TrackLenght", nLgBins, lgMin, lgMax);
422 fHistLenght->SetXTitle("Lenght of the Track (cm)");
423 fHistLenght->SetYTitle("Counts");
424
425 // OrigMult
426 fHistOrigMult = new TH1F("Flow_OrigMult", "Flow_OrigMult", nMultBins, multMin, multMax);
427 fHistOrigMult->SetXTitle("Original Mult");
428 fHistOrigMult->SetYTitle("Counts");
429
430 // MultEta
431 fHistMultEta = new TH1F("Flow_MultEta", "Flow_MultEta", nMultBins, multMin, multMax);
432 fHistMultEta->SetXTitle("Mult for Centrality");
433 fHistMultEta->SetYTitle("Counts");
434
435 // Mult
436 fHistMult = new TH1F("Flow_Mult", "Flow_Mult", nMultBins, multMin, multMax);
437 fHistMult->SetXTitle("Mult");
438 fHistMult->SetYTitle("Counts");
439
440 // V0s multiplicity
441 fHistV0Mult = new TH1F("FlowV0_Mult","FlowV0_Mult", nMultBins, multMin, multV0);
442 fHistV0Mult->SetXTitle("V0s Multiplicity");
443 fHistV0Mult->SetYTitle("Counts");
444
445 // MultOverOrig
446 fHistMultOverOrig = new TH1F("Flow_MultOverOrig", "Flow_MultOverOrig", nMultOverOrigBins, multOverOrigMin, multOverOrigMax);
447 fHistMultOverOrig->SetXTitle("Mult / Orig. Mult");
448 fHistMultOverOrig->SetYTitle("Counts");
449
450 // Mult correlated with the event planes
451 fHistMultPart = new TH1F("Flow_MultPart", "Flow_MultPart", 2*nMultBins, multMin, multMax);
452 fHistMultPart->SetXTitle("Mult of Correlated Particles");
453 fHistMultPart->SetYTitle("Counts");
454
455 // Mult correlated with the event planes in 1 unit rapidity
456 fHistMultPartUnit = new TH1F("Flow_MultPartUnit", "Flow_MultPartUnit", 2*nMultBins, multMin, multMax/2);
457 fHistMultPartUnit->SetXTitle("Mult of Correlated Particles (-0.5 < eta < 0.5)");
458 fHistMultPartUnit->SetYTitle("Counts");
459
460 // Mult of V0s correlated with the event planes
461 fHistV0MultPart = new TH1F("FlowV0_MultPart", "FlowV0_MultPart", nMultBins, multMin, multV0);
462 fHistV0MultPart->SetXTitle("Mult of Correlated V0s");
463 fHistV0MultPart->SetYTitle("Counts");
464
465 // VertexZ
466 fHistVertexZ = new TH1F("Flow_VertexZ", "Flow_VertexZ", nVertexZBins, vertexZMin, vertexZMax);
467 fHistVertexZ->SetXTitle("Vertex Z (cm)");
468 fHistVertexZ->SetYTitle("Counts");
469
470 // VertexXY
471 fHistVertexXY2D = new TH2F("Flow_VertexXY2D", "Flow_VertexXY2D", nVertexXYBins, vertexXYMin, vertexXYMax, nVertexXYBins, vertexXYMin, vertexXYMax);
472 fHistVertexXY2D->SetXTitle("Vertex X (cm)");
473 fHistVertexXY2D->SetYTitle("Vertex Y (cm)");
474
475 // EtaSym vs. Vertex Z Tpc
476 fHistEtaSymVerZ2D = new TH2F("Flow_EtaSymVerZ2D", "Flow_EtaSymVerZ2D", nVertexZBins, vertexZMin, vertexZMax, nEtaSymBins, etaSymMin, etaSymMax);
477 fHistEtaSymVerZ2D->SetXTitle("Vertex Z (cm)");
478 fHistEtaSymVerZ2D->SetYTitle("Eta Symmetry TPC");
479
480 // EtaSym Tpc
481 fHistEtaSym = new TH1F("Flow_EtaSym_TPC", "Flow_EtaSym_TPC", nEtaSymBins, etaSymMin, etaSymMax);
482 fHistEtaSym->SetXTitle("Eta Symmetry Ratio TPC");
483 fHistEtaSym->SetYTitle("Counts");
484
485 // EtaSym vs. Vertex Z Tpc - correlation analysis
486 fHistEtaSymVerZ2DPart = new TH2F("Flow_EtaSymVerZ2D_part", "Flow_EtaSymVerZ2D_part", nVertexZBins, vertexZMin, vertexZMax, nEtaSymBins, etaSymMinPart, etaSymMaxPart);
487 fHistEtaSymVerZ2DPart->SetXTitle("Vertex Z (cm)");
488 fHistEtaSymVerZ2DPart->SetYTitle("Eta Symmetry TPC");
489
490 // EtaSym Tpc - correlation analysis
491 fHistEtaSymPart = new TH1F("Flow_EtaSym_TPC_part", "Flow_EtaSym_TPC_part", nEtaSymBins, etaSymMinPart, etaSymMaxPart);
492 fHistEtaSymPart->SetXTitle("Eta Symmetry Ratio TPC");
493 fHistEtaSymPart->SetYTitle("Counts");
494
495 // phi , whatever
496 fHistPhi = new TH1F("Flow_xPhi", "Flow_xPhi (all)", fPhiBins, fPhiMin, fPhiMax);
497 fHistPhi->SetXTitle("Phi (rad)");
498 fHistPhi->SetYTitle("Counts");
499
500 // phi constrained
501 fHistPhiCons = new TH1F("Flow_cPhi", "Flow_cPhi", fPhiBins, fPhiMin, fPhiMax);
502 fHistPhiCons->SetXTitle("cPhi (rad)");
503 fHistPhiCons->SetYTitle("Counts");
504
505 // EtaPtPhi , whatever
506 fHistAllEtaPtPhi3D = new TH3F("Flow_EtaPtPhi3Dall", "Flow_EtaPtPhi3Dall (whatever)", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, nPhi3DBins, fPhiMin, fPhiMax);
507 fHistAllEtaPtPhi3D->SetXTitle("Eta");
508 fHistAllEtaPtPhi3D->SetYTitle("Pt (GeV/c)");
509 fHistAllEtaPtPhi3D->SetZTitle("Phi (rad)");
510
511 // Constrained EtaPtPhi
512 fHistConsEtaPtPhi3D = new TH3F("Flow_consEtaPtPhi3D", "Flow_consEtaPtPhi3D (constrainable)", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, nPhi3DBins, fPhiMin, fPhiMax);
513 fHistConsEtaPtPhi3D->SetXTitle("cEta");
514 fHistConsEtaPtPhi3D->SetYTitle("cPt (GeV/c)");
515 fHistConsEtaPtPhi3D->SetZTitle("cPhi (rad)");
516
517 // Global EtaPtPhi
518 fHistGlobEtaPtPhi3D = new TH3F("Flow_globEtaPtPhi3D", "Flow_globEtaPtPhi3D (constrainable)", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, nPhi3DBins, fPhiMin, fPhiMax);
519 fHistGlobEtaPtPhi3D->SetXTitle("gEta");
520 fHistGlobEtaPtPhi3D->SetYTitle("gPt (GeV/c)");
521 fHistGlobEtaPtPhi3D->SetZTitle("gPhi (rad)");
522
523 // UnConstrained EtaPtPhi
524 fHistUncEtaPtPhi3D = new TH3F("Flow_uncEtaPtPhi3D", "Flow_uncEtaPtPhi3D (un-constrainable)", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, nPhi3DBins, fPhiMin, fPhiMax);
525 fHistUncEtaPtPhi3D->SetXTitle("gEta");
526 fHistUncEtaPtPhi3D->SetYTitle("gPt (GeV/c)");
527 fHistUncEtaPtPhi3D->SetZTitle("gPhi (rad)");
528
529 // EtaPtPhi for particles correlated with the event plane
530 fHistEtaPtPhi3DPart = new TH3F("Flow_EtaPtPhi3Dpart", "Flow_EtaPtPhi3Dpart (selected part)", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, nPhi3DBins, fPhiMin, fPhiMax);
531 fHistEtaPtPhi3DPart->SetXTitle("Eta");
532 fHistEtaPtPhi3DPart->SetYTitle("Pt (GeV/c)");
533 fHistEtaPtPhi3DPart->SetZTitle("Phi (rad)");
534
535 // EtaPtPhi for particles NOT SELECTED for correlation with the event plane
536 fHistEtaPtPhi3DOut = new TH3F("Flow_EtaPtPhi3Dout", "Flow_EtaPtPhi3Dout (NOT selected part)", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, nPhi3DBins, fPhiMin, fPhiMax);
537 fHistEtaPtPhi3DOut->SetXTitle("Eta");
538 fHistEtaPtPhi3DOut->SetYTitle("Pt (GeV/c)");
539 fHistEtaPtPhi3DOut->SetZTitle("Phi (rad)");
540
541 // Yield Pt-Phi for all positive
542 fHistPtPhiPos = new TH2D("Flow_PtPhi_Plus", "Flow_PtPhi_Plus", fPhiBins, fPhiMin, fPhiMax, 50, fPtMin, fPtMax);
543 fHistPtPhiPos->Sumw2();
544 fHistPtPhiPos->SetXTitle("Phi");
545 fHistPtPhiPos->SetYTitle("Pt (GeV/c)");
546
547 // Yield Pt-Phi for all negative
548 fHistPtPhiNeg = new TH2D("Flow_PtPhi_Minus", "Flow_PtPhi_Minus", fPhiBins, fPhiMin, fPhiMax, 50, fPtMin, fPtMax);
549 fHistPtPhiNeg->Sumw2();
550 fHistPtPhiNeg->SetXTitle("Phi");
551 fHistPtPhiNeg->SetYTitle("Pt (GeV/c)");
552
553 // Yield for all particles
554 fHistYieldAll2D = new TH2D("Flow_YieldAll2D", "Flow_YieldAll2D", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
555 fHistYieldAll2D->Sumw2();
556 fHistYieldAll2D->SetXTitle("Pseudorapidty");
557 fHistYieldAll2D->SetYTitle("Pt (GeV/c)");
558
559 // Yield for constrainable tracks
560 fHistYieldCon2D = new TH2D("Flow_YieldCons2D", "Flow_YieldCons2D", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
561 fHistYieldCon2D->Sumw2();
562 fHistYieldCon2D->SetXTitle("Pseudorapidty");
563 fHistYieldCon2D->SetYTitle("Pt (GeV/c)");
564
565 // Yield for un-constrainable tracks
566 fHistYieldUnc2D = new TH2D("Flow_YieldUnc2D", "Flow_YieldUnc2D", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
567 fHistYieldUnc2D->Sumw2();
568 fHistYieldUnc2D->SetXTitle("Pseudorapidty");
569 fHistYieldUnc2D->SetYTitle("Pt (GeV/c)");
570
571 // Yield for particles correlated with the event plane
572 fHistYieldPart2D = new TH2D("Flow_YieldPart2D", "Flow_YieldPart2D (selected part)", fEtaBins, fEtaMin, fEtaMax, fPtBinsPart, fPtMin, fPtMaxPart);
573 fHistYieldPart2D->Sumw2();
574 fHistYieldPart2D->SetXTitle((char*)fLabel.Data());
575 fHistYieldPart2D->SetYTitle("Pt (GeV/c)");
576
577 // Yield for particles NOT SELECTED for correlation with the event plane
578 fHistYieldOut2D = new TH2D("Flow_YieldOut2D", "Flow_YieldOut2D (NOT selected part)", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
579 fHistYieldOut2D->Sumw2();
580 fHistYieldOut2D->SetXTitle("Pseudorapidty");
581 fHistYieldOut2D->SetYTitle("Pt (GeV/c)");
582
583 // invariant Mass for all particles (from TOF)
584 fHistInvMass = new TH1F("Flow_InvMass", "Flow_InvMass (tof)", nMassBins, massMin, massMax);
585 fHistInvMass->SetXTitle("Invariant Mass (GeV)");
586 fHistInvMass->SetYTitle("Counts");
587
588 // invariant Mass for particles correlated with the event plane (from TOF)
589 fHistInvMassPart = new TH1F("Flow_InvMassPart", "Flow_InvMassPart (tof)", nMassBins, massMin, massMax);
590 fHistInvMassPart->SetXTitle("Invariant Mass (GeV)");
591 fHistInvMassPart->SetYTitle("Counts");
592
593 // invariant Mass for particles NOT SELECTED for correlation with the event plane (from TOF)
594 fHistInvMassOut = new TH1F("Flow_InvMassOut", "Flow_InvMassOut (tof)", nMassBins, massMin, massMax);
595 fHistInvMassOut->SetXTitle("Invariant Mass (GeV)");
596 fHistInvMassOut->SetYTitle("Counts");
597
598 // Mean Eta in each bin for particles correlated with the event plane
599 fHistBinEta = new TProfile("Flow_Bin_Eta", "Flow_Bin_Eta_part (selected part)", fEtaBins, fEtaMin, fEtaMax, fEtaMin, fEtaMax, "");
600 fHistBinEta->SetXTitle((char*)fLabel.Data());
601 fHistBinEta->SetYTitle((char*)fLabel.Data());
602
603 // Mean Pt in each bin for particles correlated with the event plane
604 fHistBinPt = new TProfile("Flow_Bin_Pt", "Flow_Bin_Pt_part (selected part)", fPtBinsPart, fPtMin, fPtMaxPart, fPtMin, fPtMaxPart, "");
605 fHistBinPt->SetXTitle("Pt (GeV/c)");
606 fHistBinPt->SetYTitle("<Pt> (GeV/c)");
607
608 // cos(n*phiLab)
609 fHistCosPhi = new TProfile("Flow_CosPhiLab", "Flow_CosPhiLab", Flow::nHars, 0.5, (float)(Flow::nHars) + 0.5, -100., 100., "");
610 fHistCosPhi->SetXTitle("Harmonic");
611 fHistCosPhi->SetYTitle("<cos(n*PhiLab)> (%)");
612
613 // PID pi+
614 fHistPidPiPlus = new TH1F("Flow_PidPiPlus", "Flow_PidPiPlus", nPidBins, pidMin, pidMax);
615 fHistPidPiPlus->SetXTitle("ALICE P.Id.");
616 fHistPidPiPlus->SetYTitle("Counts");
617
618 // PID pi-
619 fHistPidPiMinus = new TH1F("Flow_PidPiMinus", "Flow_PidPiMinus", nPidBins, pidMin, pidMax);
620 fHistPidPiMinus->SetXTitle("ALICE P.Id.");
621 fHistPidPiMinus->SetYTitle("Counts");
622
623 // PID proton
624 fHistPidProton = new TH1F("Flow_PidProton", "Flow_PidProton", nPidBins, pidMin, pidMax);
625 fHistPidProton->SetXTitle("ALICE P.Id.");
626 fHistPidProton->SetYTitle("Counts");
627
628 // PID anti proton
629 fHistPidAntiProton = new TH1F("Flow_PidAntiProton", "Flow_PidAntiProton", nPidBins, pidMin, pidMax);
630 fHistPidAntiProton->SetXTitle("ALICE P.Id.");
631 fHistPidAntiProton->SetYTitle("Counts");
632
633 // PID Kplus
634 fHistPidKplus = new TH1F("Flow_PidKplus", "Flow_PidKplus", nPidBins, pidMin, pidMax);
635 fHistPidKplus->SetXTitle("ALICE P.Id.");
636 fHistPidKplus->SetYTitle("Counts");
637
638 // PID Kminus
639 fHistPidKminus = new TH1F("Flow_PidKminus", "Flow_PidKminus", nPidBins, pidMin, pidMax);
640 fHistPidKminus->SetXTitle("ALICE P.Id.");
641 fHistPidKminus->SetYTitle("Counts");
642
643 // PID deuteron
644 fHistPidDeuteron = new TH1F("Flow_PidDeuteron", "Flow_PidDeuteron", nPidBins, pidMin, pidMax);
645 fHistPidDeuteron->SetXTitle("ALICE P.Id.");
646 fHistPidDeuteron->SetYTitle("Counts");
647
648 // PID anti deuteron
649 fHistPidAntiDeuteron = new TH1F("Flow_PidAntiDeuteron", "Flow_PidAntiDeuteron", nPidBins, pidMin, pidMax);
650 fHistPidAntiDeuteron->SetXTitle("ALICE P.Id.");
651 fHistPidAntiDeuteron->SetYTitle("Counts");
652
653 // PID electron
654 fHistPidElectron = new TH1F("Flow_PidElectron", "Flow_PidElectron", nPidBins, pidMin, pidMax);
655 fHistPidElectron->SetXTitle("ALICE P.Id.");
656 fHistPidElectron->SetYTitle("Counts");
657
658 // PID positron
659 fHistPidPositron = new TH1F("Flow_PidPositron", "Flow_PidPositron", nPidBins, pidMin, pidMax);
660 fHistPidPositron->SetXTitle("ALICE P.Id.");
661 fHistPidPositron->SetYTitle("Counts");
662
663 // PID Muon+
664 fHistPidMuonPlus = new TH1F("Flow_PidMuonPlus", "Flow_PidMuonPlus", nPidBins, pidMin, pidMax);
665 fHistPidMuonPlus->SetXTitle("ALICE P.Id.");
666 fHistPidMuonPlus->SetYTitle("Counts");
667
668 // PID Muon-
669 fHistPidMuonMinus = new TH1F("Flow_PidMuonMinus", "Flow_PidMuonMinus", nPidBins, pidMin, pidMax);
670 fHistPidMuonMinus->SetXTitle("ALICE P.Id.");
671 fHistPidMuonMinus->SetYTitle("Counts");
672
673 // PID pi+ selected
674 fHistPidPiPlusPart = new TH1F("Flow_PidPiPlusPart", "Flow_PidPiPlusPart", nPidBins, pidMin, pidMax);
675 fHistPidPiPlusPart->SetXTitle("ALICE P.Id. (pid = pi+)");
676 fHistPidPiPlusPart->SetYTitle("Counts");
677
678 // PID pi- selected
679 fHistPidPiMinusPart = new TH1F("Flow_PidPiMinusPart", "Flow_PidPiMinusPart", nPidBins, pidMin, pidMax);
680 fHistPidPiMinusPart->SetXTitle("ALICE P.Id. (pid = pi-)");
681 fHistPidPiMinusPart->SetYTitle("Counts");
682
683 // PID proton selected
684 fHistPidProtonPart = new TH1F("Flow_PidProtonPart", "Flow_PidProtonPart", nPidBins, pidMin, pidMax);
685 fHistPidProtonPart->SetXTitle("ALICE P.Id. (pid = p)");
686 fHistPidProtonPart->SetYTitle("Counts");
687
688 // PID anti proton selected
689 fHistPidAntiProtonPart = new TH1F("Flow_PidAntiProtonPart", "Flow_PidAntiProtonPart", nPidBins, pidMin, pidMax);
690 fHistPidAntiProtonPart->SetXTitle("ALICE P.Id. (pid = p-)");
691 fHistPidAntiProtonPart->SetYTitle("Counts");
692
693 // PID Kplus selected
694 fHistPidKplusPart = new TH1F("Flow_PidKplusPart", "Flow_PidKplusPart", nPidBins, pidMin, pidMax);
695 fHistPidKplusPart->SetXTitle("ALICE P.Id. (pid = K+)");
696 fHistPidKplusPart->SetYTitle("Counts");
697
698 // PID Kminus selected
699 fHistPidKminusPart = new TH1F("Flow_PidKminusPart", "Flow_PidKminusPart", nPidBins, pidMin, pidMax);
700 fHistPidKminusPart->SetXTitle("ALICE P.Id. (pid = K-)");
701 fHistPidKminusPart->SetYTitle("Counts");
702
703 // PID deuteron selected
704 fHistPidDeuteronPart = new TH1F("Flow_PidDeuteronPart", "Flow_PidDeuteronPart", nPidBins, pidMin, pidMax);
705 fHistPidDeuteronPart->SetXTitle("ALICE P.Id. (pid = d)");
706 fHistPidDeuteronPart->SetYTitle("Counts");
707
708 // PID anti deuteron selected
709 fHistPidAntiDeuteronPart = new TH1F("Flow_PidAntiDeuteronPart", "Flow_PidAntiDeuteronPart", nPidBins, pidMin, pidMax);
710 fHistPidAntiDeuteronPart->SetXTitle("ALICE P.Id. (pid = d--)");
711 fHistPidAntiDeuteronPart->SetYTitle("Counts");
712
713 // PID electron selected
714 fHistPidElectronPart = new TH1F("Flow_PidElectronPart", "Flow_PidElectronPart", nPidBins, pidMin, pidMax);
715 fHistPidElectronPart->SetXTitle("ALICE P.Id. (pid = e-)");
716 fHistPidElectronPart->SetYTitle("Counts");
717
718 // PID positron selected
719 fHistPidPositronPart = new TH1F("Flow_PidPositronPart", "Flow_PidPositronPart", nPidBins, pidMin, pidMax);
720 fHistPidPositronPart->SetXTitle("ALICE P.Id. (pid = e+)");
721 fHistPidPositronPart->SetYTitle("Counts");
722
723 // PID Muon+ selected
724 fHistPidMuonPlusPart = new TH1F("Flow_PidMuonPlusPart", "Flow_PidMuonPlusPart", nPidBins, pidMin, pidMax);
725 fHistPidMuonPlusPart->SetXTitle("ALICE P.Id. (pid = mu+)");
726 fHistPidMuonPlusPart->SetYTitle("Counts");
727
728 // PID Muon- selected
729 fHistPidMuonMinusPart = new TH1F("Flow_PidMuonMinusPart", "Flow_PidMuonMinusPart", nPidBins, pidMin, pidMax);
730 fHistPidMuonMinusPart->SetXTitle("ALICE P.Id. (pid = mu-)");
731 fHistPidMuonMinusPart->SetYTitle("Counts");
732
733 // PID multiplicities (all)
734 fHistPidMult = new TProfile("Flow_PidMult", "Flow_PidMult", 15, 0.5, 15.5, "");
735 fHistPidMult->SetXTitle("all, h+, h-, pi+, pi-, pr+, pr-, K+, K-, d+, d-, e-, e+, mu-, mu+");
736 fHistPidMult->SetYTitle("Multiplicity");
737
738 // PID for all tracks
739 fHistBayPidMult = new TH1F("Flow_BayPidMult","Flow_BayPidMult",Flow::nPid,-0.5,((float)Flow::nPid-0.5));
740 fHistBayPidMult->Sumw2() ;
741 fHistBayPidMult->SetXTitle("e+/- , mu+/- , pi+/- , K+/- , p+/- , d+/- ");
742 fHistBayPidMult->SetYTitle("Counts");
743
744 // PID for particles correlated with the event plane
745 fHistBayPidMultPart = new TH1F("Flow_BayPidMult_Part","Flow_BayPidMult_Part",Flow::nPid,-0.5,((float)Flow::nPid-0.5));
746 fHistBayPidMultPart->Sumw2() ;
747 fHistBayPidMultPart->SetXTitle("e+/- , mu+/- , pi+/- , K+/- , p+/- , d+/- ");
748 fHistBayPidMultPart->SetYTitle("Counts");
749
750 // Centrality
751 fHistCent = new TH1F("Flow_Cent", "Flow_Cent", nCentBins, centMin, centMax);
752 fHistCent->SetXTitle("Centrality Bin");
753 fHistCent->SetYTitle("Counts");
754
755 // Deposited Energy in ZDC
756 fHistEnergyZDC = new TH2F("Flow_ZDC_E", "Flow_ZDC_E", 3, 0., 3., nZDCeBins, ZDCeMin, ZDCeMax);
757 fHistEnergyZDC->SetXTitle("neutron , proton , e.m.");
758 fHistEnergyZDC->SetYTitle("ZDC energy");
759
760 // Part. versus Energy in ZDC
761 fHistPartZDC = new TH1F("Flow_ZDC_Participants", "Flow_ZDC_Participants", nZDCpartBins, ZDCpartMin, ZDCpartMax);
762 fHistPartZDC->SetXTitle("ZDC part");
763 fHistPartZDC->SetYTitle("Counts");
764
765 // MeanDedxPos TPC
766 fHistMeanDedxPos2D = new TH2F("Flow_MeanDedxPos2D_TPC","Flow_MeanDedxPos2D_TPC", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMaxTPC);
767 fHistMeanDedxPos2D->SetXTitle("log(momentum) (GeV/c)");
768 fHistMeanDedxPos2D->SetYTitle("mean dEdx (+)");
769
770 // MeanDedxNeg TPC
771 fHistMeanDedxNeg2D = new TH2F("Flow_MeanDedxNeg2D_TPC","Flow_MeanDedxNeg2D_TPC", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMaxTPC);
772 fHistMeanDedxNeg2D->SetXTitle("log(momentum) (GeV/c)");
773 fHistMeanDedxNeg2D->SetYTitle("mean dEdx (-)");
774
775 // MeanDedxPos ITS
776 fHistMeanDedxPos2DITS = new TH2F("Flow_MeanDedxPos2D_ITS","Flow_MeanDedxPos2D_ITS", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMax);
777 fHistMeanDedxPos2DITS->SetXTitle("log(momentum) (GeV/c)");
778 fHistMeanDedxPos2DITS->SetYTitle("mean dEdx (+)");
779
780 // MeanDedxNeg ITS
781 fHistMeanDedxNeg2DITS = new TH2F("Flow_MeanDedxNeg2D_ITS","Flow_MeanDedxNeg2D_ITS", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMax);
782 fHistMeanDedxNeg2DITS->SetXTitle("log(momentum) (GeV/c)");
783 fHistMeanDedxNeg2DITS->SetYTitle("mean dEdx (-)");
784
785 // MeanDedxPos TPC 3D selected Part
786 fHistMeanDedxPos3DPart = new TH3F("Flow_MeanDedxPos3DPart_TPC","Flow_MeanDedxPos3DPart_TPC", nMomenBins, logpMin, logpMax, nDedxBins, 0., dEdxMaxTPC, Flow::nPid, 0., Flow::nPid);
787 fHistMeanDedxPos3DPart->SetXTitle("log(momentum) (GeV/c)");
788 fHistMeanDedxPos3DPart->SetYTitle("mean dEdx (+)");
789 fHistMeanDedxPos3DPart->SetZTitle("e , mu , pi , k , p , d");
790
791 // MeanDedxNeg TPC 3D selected Part
792 fHistMeanDedxNeg3DPart = new TH3F("Flow_MeanDedxNeg3DPart_TPC","Flow_MeanDedxNeg3DPart_TPC", nMomenBins, logpMin, logpMax, nDedxBins, 0., dEdxMaxTPC, Flow::nPid, 0., Flow::nPid);
793 fHistMeanDedxNeg3DPart->SetXTitle("log(momentum) (GeV/c)");
794 fHistMeanDedxNeg3DPart->SetYTitle("mean dEdx (-)");
795 fHistMeanDedxNeg3DPart->SetZTitle("e , mu , pi , k , p , d");
796
797 // MeanDedxPos ITS 3D selected Part
798 fHistMeanDedxPos3DPartITS = new TH3F("Flow_MeanDedxPos3DPart_ITS","Flow_MeanDedxPos3DPart_ITS", nMomenBins, logpMin, logpMax, nDedxBins, 0., dEdxMax, Flow::nPid, 0., Flow::nPid);
799 fHistMeanDedxPos3DPartITS->SetXTitle("log(momentum) (GeV/c)");
800 fHistMeanDedxPos3DPartITS->SetYTitle("mean dEdx (+)");
801 fHistMeanDedxPos3DPartITS->SetZTitle("e , mu , pi , k , p , d");
802
803 // MeanDedxNeg ITS 3D selected Part
804 fHistMeanDedxNeg3DPartITS = new TH3F("Flow_MeanDedxNeg3DPart_ITS","Flow_MeanDedxNeg3DPart_ITS", nMomenBins, logpMin, logpMax, nDedxBins, 0., dEdxMax, Flow::nPid, 0., Flow::nPid);
805 fHistMeanDedxNeg3DPartITS->SetXTitle("log(momentum) (GeV/c)");
806 fHistMeanDedxNeg3DPartITS->SetYTitle("mean dEdx (-)");
807 fHistMeanDedxNeg3DPartITS->SetZTitle("e , mu , pi , k , p , d");
808
809 // MeanSignalPos TRD
810 fHistMeanDedxPos2DTRD = new TH2F("Flow_MeanSignalPos2D_TRD","Flow_MeanSignalPos2D_TRD", nMomenBins, logpMin, logpMax, nTrdBins, 0, TRDmax);
811 fHistMeanDedxPos2DTRD->SetXTitle("log(momentum) (GeV/c)");
812 fHistMeanDedxPos2DTRD->SetYTitle("mean TRD (+)");
813
814 // MeanSignalNeg TRD
815 fHistMeanDedxNeg2DTRD = new TH2F("Flow_MeanSignalNeg2D_TRD","Flow_MeanSignalNeg2D_TRD", nMomenBins, logpMin, logpMax, nTrdBins, 0, TRDmax);
816 fHistMeanDedxNeg2DTRD->SetXTitle("log(momentum) (GeV/c)");
817 fHistMeanDedxNeg2DTRD->SetYTitle("mean TRD (-)");
818
819 // MeanTimePos TOF
820 fHistMeanDedxPos2DTOF = new TH2F("Flow_MeanTimePos2D_TOF","Flow_MeanTimePos2D_TOF", nMomenBins, logpMin, logpMax, nTofBins, TOFmin, TOFmax);
821 fHistMeanDedxPos2DTOF->SetXTitle("log(momentum) (GeV/c)");
822 fHistMeanDedxPos2DTOF->SetYTitle("mean Time of Flight (+)");
823
824 // MeanTimeNeg TOF
825 fHistMeanDedxNeg2DTOF = new TH2F("Flow_MeanTimeNeg2D_TOF","Flow_MeanTimeNeg2D_TOF", nMomenBins, logpMin, logpMax, nTofBins, TOFmin, TOFmax);
826 fHistMeanDedxNeg2DTOF->SetXTitle("log(momentum) (GeV/c)");
827 fHistMeanDedxNeg2DTOF->SetYTitle("mean Time of Flight (-)");
828
829 // Detector response for each particle { ...
830 // MeanDedx PiPlus - TPC
831 fHistMeanTPCPiPlus = new TH2F("Flow_MeanDedxPiPlusTPC", "Flow_MeanDedxPiPlusTPC", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMaxTPC);
832 fHistMeanTPCPiPlus->SetXTitle("Log10(p) (GeV/c)");
833 fHistMeanTPCPiPlus->SetYTitle("mean dEdx (pi+)");
834 // MeanDedx PiMinus
835 fHistMeanTPCPiMinus = new TH2F("Flow_MeanDedxPiMinusTPC", "Flow_MeanDedxPiMinusTPC", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMaxTPC);
836 fHistMeanTPCPiMinus->SetXTitle("Log10(p) (GeV/c)");
837 fHistMeanTPCPiMinus->SetYTitle("mean dEdx (pi-)");
838 // MeanDedx Proton
839 fHistMeanTPCProton = new TH2F("Flow_MeanDedxProtonTPC", "Flow_MeanDedxProtonTPC", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMaxTPC);
840 fHistMeanTPCProton->SetXTitle("Log10(p) (GeV/c)");
841 fHistMeanTPCProton->SetYTitle("mean dEdx (pr+)");
842 // MeanDedx Pbar
843 fHistMeanTPCPbar = new TH2F("Flow_MeanDedxPbarTPC", "Flow_MeanDedxPbarTPC", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMaxTPC);
844 fHistMeanTPCPbar->SetXTitle("Log10(p) (GeV/c)");
845 fHistMeanTPCPbar->SetYTitle("mean dEdx (pr-)");
846 // MeanDedx Kplus
847 fHistMeanTPCKplus = new TH2F("Flow_MeanDedxKplusTPC", "Flow_MeanDedxKplusTPC", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMaxTPC);
848 fHistMeanTPCKplus->SetXTitle("Log10(p) (GeV/c)");
849 fHistMeanTPCKplus->SetYTitle("mean dEdx (K+)");
850 // MeanDedx Kminus
851 fHistMeanTPCKminus = new TH2F("Flow_MeanDedxKminusTPC", "Flow_MeanDedxKminusTPC", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMaxTPC);
852 fHistMeanTPCKminus->SetXTitle("Log10(p) (GeV/c)");
853 fHistMeanTPCKminus->SetYTitle("mean dEdx (K-)");
854 // MeanDedx Deuteron
855 fHistMeanTPCDeuteron = new TH2F("Flow_MeanDedxDeuteronTPC", "Flow_MeanDedxDeuteronTPC", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMaxTPC);
856 fHistMeanTPCDeuteron->SetXTitle("Log10(p) (GeV/c)");
857 fHistMeanTPCDeuteron->SetYTitle("mean dEdx (d+)");
858 // MeanDedx AntiDeuteron
859 fHistMeanTPCAntiDeuteron = new TH2F("Flow_MeanDedxAntiDeuteronTPC", "Flow_MeanDedxAntiDeuteronTPC", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMaxTPC);
860 fHistMeanTPCAntiDeuteron->SetXTitle("Log10(p) (GeV/c)");
861 fHistMeanTPCAntiDeuteron->SetYTitle("mean dEdx (d-)");
862 // MeanDedxElectron
863 fHistMeanTPCElectron = new TH2F("Flow_MeanDedxElectronTPC", "Flow_MeanDedxElectronTPC", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMaxTPC);
864 fHistMeanTPCElectron->SetXTitle("Log10(p) (GeV/c)");
865 fHistMeanTPCElectron->SetYTitle("mean dEdx (e-)");
866 // MeanDedx Positron
867 fHistMeanTPCPositron = new TH2F("Flow_MeanDedxPositronTPC", "Flow_MeanDedxPositronTPC", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMaxTPC);
868 fHistMeanTPCPositron->SetXTitle("Log10(p) (GeV/c)");
869 fHistMeanTPCPositron->SetYTitle("mean dEdx (e+)");
870 // MeanDedx MuonPlus
871 fHistMeanTPCMuonPlus = new TH2F("Flow_MeanDedxMuonPlusTPC", "Flow_MeanDedxMuonPlusTPC", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMaxTPC);
872 fHistMeanTPCMuonPlus->SetXTitle("Log10(p) (GeV/c)");
873 fHistMeanTPCMuonPlus->SetYTitle("mean dEdx (mu+)");
874 // MeanDedx MuonMinus
875 fHistMeanTPCMuonMinus = new TH2F("Flow_MeanDedxMuonMinusTPC", "Flow_MeanDedxMuonMinusTPC", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMaxTPC);
876 fHistMeanTPCMuonMinus->SetXTitle("Log10(p) (GeV/c)");
877 fHistMeanTPCMuonMinus->SetYTitle("mean dEdx (mu-)");
878
879 // MeanDedx PiPlus - ITS
880 fHistMeanITSPiPlus = new TH2F("Flow_MeanDedxPiPlusITS", "Flow_MeanDedxPiPlusITS", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMax);
881 fHistMeanITSPiPlus->SetXTitle("Log10(p) (GeV/c)");
882 fHistMeanITSPiPlus->SetYTitle("mean dEdx (pi+)");
883 // MeanDedx PiMinus
884 fHistMeanITSPiMinus = new TH2F("Flow_MeanDedxPiMinusITS", "Flow_MeanDedxPiMinusITS", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMax);
885 fHistMeanITSPiMinus->SetXTitle("Log10(p) (GeV/c)");
886 fHistMeanITSPiMinus->SetYTitle("mean dEdx (pi-)");
887 // MeanDedx Proton
888 fHistMeanITSProton = new TH2F("Flow_MeanDedxProtonITS", "Flow_MeanDedxProtonITS", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMax);
889 fHistMeanITSProton->SetXTitle("Log10(p) (GeV/c)");
890 fHistMeanITSProton->SetYTitle("mean dEdx (pr+)");
891 // MeanDedx Pbar
892 fHistMeanITSPbar = new TH2F("Flow_MeanDedxPbarITS", "Flow_MeanDedxPbarITS", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMax);
893 fHistMeanITSPbar->SetXTitle("Log10(p) (GeV/c)");
894 fHistMeanITSPbar->SetYTitle("mean dEdx (pr-)");
895 // MeanDedx Kplus
896 fHistMeanITSKplus = new TH2F("Flow_MeanDedxKplusITS", "Flow_MeanDedxKplusITS", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMax);
897 fHistMeanITSKplus->SetXTitle("Log10(p) (GeV/c)");
898 fHistMeanITSKplus->SetYTitle("mean dEdx (K+)");
899 // MeanDedx Kminus
900 fHistMeanITSKminus = new TH2F("Flow_MeanDedxKminusITS", "Flow_MeanDedxKminusITS", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMax);
901 fHistMeanITSKminus->SetXTitle("Log10(p) (GeV/c)");
902 fHistMeanITSKminus->SetYTitle("mean dEdx (K-)");
903 // MeanDedx Deuteron
904 fHistMeanITSDeuteron = new TH2F("Flow_MeanDedxDeuteronITS", "Flow_MeanDedxDeuteronITS", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMax);
905 fHistMeanITSDeuteron->SetXTitle("Log10(p) (GeV/c)");
906 fHistMeanITSDeuteron->SetYTitle("mean dEdx (d+)");
907 // MeanDedx AntiDeuteron
908 fHistMeanITSAntiDeuteron = new TH2F("Flow_MeanDedxAntiDeuteronITS", "Flow_MeanDedxAntiDeuteronITS", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMax);
909 fHistMeanITSAntiDeuteron->SetXTitle("Log10(p) (GeV/c)");
910 fHistMeanITSAntiDeuteron->SetYTitle("mean dEdx (d-)");
911 // MeanDedx Electron
912 fHistMeanITSElectron = new TH2F("Flow_MeanDedxElectronITS", "Flow_MeanDedxElectronITS", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMax);
913 fHistMeanITSElectron->SetXTitle("Log10(p) (GeV/c)");
914 fHistMeanITSElectron->SetYTitle("mean dEdx (e-)");
915 // MeanDedx Positron
916 fHistMeanITSPositron = new TH2F("Flow_MeanDedxPositronITS", "Flow_MeanDedxPositronITS", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMax);
917 fHistMeanITSPositron->SetXTitle("Log10(p) (GeV/c)");
918 fHistMeanITSPositron->SetYTitle("mean dEdx (e+)");
919 // MeanDedx MuonPlus
920 fHistMeanITSMuonPlus = new TH2F("Flow_MeanDedxMuonPlusITS", "Flow_MeanDedxMuonPlusITS", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMax);
921 fHistMeanITSMuonPlus->SetXTitle("Log10(p) (GeV/c)");
922 fHistMeanITSMuonPlus->SetYTitle("mean dEdx (mu+)");
923 // MeanDedx MuonMinus
924 fHistMeanITSMuonMinus = new TH2F("Flow_MeanDedxMuonMinusITS", "Flow_MeanDedxMuonMinusITS", nMomenBins, logpMin, logpMax, nDedxBins, 0, dEdxMax);
925 fHistMeanITSMuonMinus->SetXTitle("Log10(p) (GeV/c)");
926 fHistMeanITSMuonMinus->SetYTitle("mean dEdx (mu-)");
927
928 // MeanTrd PiPlus - TRD
929 fHistMeanTRDPiPlus = new TH2F("Flow_MeanTrdPiPlusTRD", "Flow_MeanTrdPiPlusTRD", nMomenBins, logpMin, logpMax, nTrdBins, 0, TRDmax);
930 fHistMeanTRDPiPlus->SetXTitle("Log10(p) (GeV/c)");
931 fHistMeanTRDPiPlus->SetYTitle("mean t.r.[] (pi+)");
932 // MeanTrd PiMinus
933 fHistMeanTRDPiMinus = new TH2F("Flow_MeanTrdPiMinusTRD", "Flow_MeanTrdPiMinusTRD", nMomenBins, logpMin, logpMax, nTrdBins, 0, TRDmax);
934 fHistMeanTRDPiMinus->SetXTitle("Log10(p) (GeV/c)");
935 fHistMeanTRDPiMinus->SetYTitle("mean t.r.[] (pi-)");
936 // MeanTrd Proton
937 fHistMeanTRDProton = new TH2F("Flow_MeanTrdProtonTRD", "Flow_MeanTrdProtonTRD", nMomenBins, logpMin, logpMax, nTrdBins, 0, TRDmax);
938 fHistMeanTRDProton->SetXTitle("Log10(p) (GeV/c)");
939 fHistMeanTRDProton->SetYTitle("mean t.r.[] (pr+)");
940 // MeanTrd Pbar
941 fHistMeanTRDPbar = new TH2F("Flow_MeanTrdPbarTRD", "Flow_MeanTrdPbarTRD", nMomenBins, logpMin, logpMax, nTrdBins, 0, TRDmax);
942 fHistMeanTRDPbar->SetXTitle("Log10(p) (GeV/c)");
943 fHistMeanTRDPbar->SetYTitle("mean t.r.[] (pr-)");
944 // MeanTrd Kplus
945 fHistMeanTRDKplus = new TH2F("Flow_MeanTrdKplusTRD", "Flow_MeanTrdKplusTRD", nMomenBins, logpMin, logpMax, nTrdBins, 0, TRDmax);
946 fHistMeanTRDKplus->SetXTitle("Log10(p) (GeV/c)");
947 fHistMeanTRDKplus->SetYTitle("mean t.r.[] (K+)");
948 // MeanTrd Kminus
949 fHistMeanTRDKminus = new TH2F("Flow_MeanTrdKminusTRD", "Flow_MeanTrdKminusTRD", nMomenBins, logpMin, logpMax, nTrdBins, 0, TRDmax);
950 fHistMeanTRDKminus->SetXTitle("Log10(p) (GeV/c)");
951 fHistMeanTRDKminus->SetYTitle("mean t.r.[] (K-)");
952 // MeanTrd Deuteron
953 fHistMeanTRDDeuteron = new TH2F("Flow_MeanTrdDeuteronTRD", "Flow_MeanTrdDeuteronTRD", nMomenBins, logpMin, logpMax, nTrdBins, 0, TRDmax);
954 fHistMeanTRDDeuteron->SetXTitle("Log10(p) (GeV/c)");
955 fHistMeanTRDDeuteron->SetYTitle("mean t.r.[] (d+)");
956 // MeanTrd AntiDeuteron
957 fHistMeanTRDAntiDeuteron = new TH2F("Flow_MeanTrdAntiDeuteronTRD", "Flow_MeanTrdAntiDeuteronTRD", nMomenBins, logpMin, logpMax, nTrdBins, 0, TRDmax);
958 fHistMeanTRDAntiDeuteron->SetXTitle("Log10(p) (GeV/c)");
959 fHistMeanTRDAntiDeuteron->SetYTitle("mean t.r.[] (d-)");
960 // MeanTrd Electron
961 fHistMeanTRDElectron = new TH2F("Flow_MeanTrdElectronTRD", "Flow_MeanTrdElectronTRD", nMomenBins, logpMin, logpMax, nTrdBins, 0, TRDmax);
962 fHistMeanTRDElectron->SetXTitle("Log10(p) (GeV/c)");
963 fHistMeanTRDElectron->SetYTitle("mean t.r.[] (e-)");
964 // MeanTrd Positron
965 fHistMeanTRDPositron = new TH2F("Flow_MeanTrdPositronTRD", "Flow_MeanTrdPositronTRD", nMomenBins, logpMin, logpMax, nTrdBins, 0, TRDmax);
966 fHistMeanTRDPositron->SetXTitle("Log10(p) (GeV/c)");
967 fHistMeanTRDPositron->SetYTitle("mean t.r.[] (e+)");
968 // MeanTrd MuonPlus
969 fHistMeanTRDMuonPlus = new TH2F("Flow_MeanTrdMuonPlusTRD", "Flow_MeanTrdMuonPlusTRD", nMomenBins, logpMin, logpMax, nTrdBins, 0, TRDmax);
970 fHistMeanTRDMuonPlus->SetXTitle("Log10(p) (GeV/c)");
971 fHistMeanTRDMuonPlus->SetYTitle("mean t.r.[] (mu+)");
972 // MeanTrd MuonMinus
973 fHistMeanTRDMuonMinus = new TH2F("Flow_MeanTrdMuonMinusTRD", "Flow_MeanTrdMuonMinusTRD", nMomenBins, logpMin, logpMax, nTrdBins, 0, TRDmax);
974 fHistMeanTRDMuonMinus->SetXTitle("Log10(p) (GeV/c)");
975 fHistMeanTRDMuonMinus->SetYTitle("mean t.r.[] (mu-)");
976
977 // T.O.F. PiPlus - TOF
978 fHistMeanTOFPiPlus = new TH2F("Flow_MeanTofPiPlusTOF", "Flow_MeanTofPiPlusTOF", nMassBins, massMin, massMax, nTofBins, TOFmin, TOFmax);
979 fHistMeanTOFPiPlus->SetXTitle("invariant mass (GeV)");
980 fHistMeanTOFPiPlus->SetYTitle("mean t.o.f.[psec] (pi+)");
981 // MeanTof PiMinus
982 fHistMeanTOFPiMinus = new TH2F("Flow_MeanTofPiMinusTOF", "Flow_MeanTofPiMinusTOF", nMassBins, massMin, massMax, nTofBins, TOFmin, TOFmax);
983 fHistMeanTOFPiMinus->SetXTitle("invariant mass (GeV)");
984 fHistMeanTOFPiMinus->SetYTitle("mean t.o.f.[psec] (pi-)");
985 // MeanTof Proton
986 fHistMeanTOFProton = new TH2F("Flow_MeanTofProtonTOF", "Flow_MeanTofProtonTOF", nMassBins, massMin, massMax, nTofBins, TOFmin, TOFmax);
987 fHistMeanTOFProton->SetXTitle("invariant mass (GeV)");
988 fHistMeanTOFProton->SetYTitle("mean t.o.f.[psec] (pr+)");
989 // Mean TofPbar
990 fHistMeanTOFPbar = new TH2F("Flow_MeanTofPbarTOF", "Flow_MeanTofPbarTOF", nMassBins, massMin, massMax, nTofBins, TOFmin, TOFmax);
991 fHistMeanTOFPbar->SetXTitle("invariant mass (GeV)");
992 fHistMeanTOFPbar->SetYTitle("mean t.o.f.[psec] (pr-)");
993 // mean t.o.f.[psec]Kplus
994 fHistMeanTOFKplus = new TH2F("Flow_MeanTofKplusTOF", "Flow_MeanTofKplusTOF", nMassBins, massMin, massMax, nTofBins, TOFmin, TOFmax);
995 fHistMeanTOFKplus->SetXTitle("invariant mass (GeV)");
996 fHistMeanTOFKplus->SetYTitle("mean t.o.f.[psec] (K+)");
997 // mean t.o.f.[psec]Kminus
998 fHistMeanTOFKminus = new TH2F("Flow_MeanTofKminusTOF", "Flow_MeanTofKminusTOF", nMassBins, massMin, massMax, nTofBins, TOFmin, TOFmax);
999 fHistMeanTOFKminus->SetXTitle("invariant mass (GeV)");
1000 fHistMeanTOFKminus->SetYTitle("mean t.o.f.[psec] (K-)");
1001 // MeanTof Deuteron
1002 fHistMeanTOFDeuteron = new TH2F("Flow_MeanTofDeuteronTOF", "Flow_MeanTofDeuteronTOF", nMassBins, massMin, massMax, nTofBins, TOFmin, TOFmax);
1003 fHistMeanTOFDeuteron->SetXTitle("invariant mass (GeV)");
1004 fHistMeanTOFDeuteron->SetYTitle("mean t.o.f.[psec] (d+)");
1005 // MeanTof AntiDeuteron
1006 fHistMeanTOFAntiDeuteron = new TH2F("Flow_MeanTofAntiDeuteronTOF", "Flow_MeanTofAntiDeuteronTOF", nMassBins, massMin, massMax, nTofBins, TOFmin, TOFmax);
1007 fHistMeanTOFAntiDeuteron->SetXTitle("invariant mass (GeV)");
1008 fHistMeanTOFAntiDeuteron->SetYTitle("mean t.o.f.[psec] (d-)");
1009 // MeanTof Electron
1010 fHistMeanTOFElectron = new TH2F("Flow_MeanTofElectronTOF", "Flow_MeanTofElectronTOF", nMassBins, massMin, massMax, nTofBins, TOFmin, TOFmax);
1011 fHistMeanTOFElectron->SetXTitle("invariant mass (GeV)");
1012 fHistMeanTOFElectron->SetYTitle("mean t.o.f.[psec] (e-)");
1013 // MeanTof Positron
1014 fHistMeanTOFPositron = new TH2F("Flow_MeanTofPositronTOF", "Flow_MeanTofPositronTOF", nMassBins, massMin, massMax, nTofBins, TOFmin, TOFmax);
1015 fHistMeanTOFPositron->SetXTitle("invariant mass (GeV)");
1016 fHistMeanTOFPositron->SetYTitle("mean t.o.f.[psec] (e+)");
1017 // MeanTof MuonPlus
1018 fHistMeanTOFMuonPlus = new TH2F("Flow_MeanTofMuonPlusTOF", "Flow_MeanTofMuonPlusTOF", nMassBins, massMin, massMax, nTofBins, TOFmin, TOFmax);
1019 fHistMeanTOFMuonPlus->SetXTitle("invariant mass (GeV)");
1020 fHistMeanTOFMuonPlus->SetYTitle("mean t.o.f.[psec] (mu+)");
1021 // MeanTof MuonMinus
1022 fHistMeanTOFMuonMinus = new TH2F("Flow_MeanTofMuonMinusTOF", "Flow_MeanTofMuonMinusTOF", nMassBins, massMin, massMax, nTofBins, TOFmin, TOFmax);
1023 fHistMeanTOFMuonMinus->SetXTitle("invariant mass (GeV)");
1024 fHistMeanTOFMuonMinus->SetYTitle("mean t.o.f.[psec] (mu-)");
1025 // ... }
1026
1027 TString* histTitle;
1028 for (int n = 0; n < Flow::nSubs; n++) // for sub-events
1029 {
1030 for (int k = 0; k < Flow::nSels; k++)
1031 {
1032 for (int j = 0; j < Flow::nHars; j++)
1033 {
1034 float order = (float)(j + 1);
1035 int i = Flow::nSubs * k + n ;
1036
1037 // event planes
1038 histTitle = new TString("Flow_Psi_Sub");
1039 *histTitle += n+1;
1040 histTitle->Append("_Sel");
1041 *histTitle += k+1;
1042 histTitle->Append("_Har");
1043 *histTitle += j+1;
1044 fHistSub[i].fHistSubHar[j].fHistPsiSubs = new TH1F(histTitle->Data(),histTitle->Data(), nPsiBins, psiMin, (psiMax / order));
1045 fHistSub[i].fHistSubHar[j].fHistPsiSubs->SetXTitle("Event Plane Angle (rad)");
1046 fHistSub[i].fHistSubHar[j].fHistPsiSubs->SetYTitle("Counts");
1047 delete histTitle;
1048 }
1049 }
1050 }
1051
1052 if(fV0loop) // All V0s (if there, if flag on)
1053 {
1054 // Mass
1055 fHistV0Mass = new TH1F("FlowV0_InvMass", "FlowV0_InvMass", nMassBins, massMin, massMax);
1056 fHistV0Mass->SetXTitle("Invariant Mass (GeV)");
1057 fHistV0Mass->SetYTitle("Counts");
1058 // Distance of closest approach
1059 fHistV0Dca = new TH1F("FlowV0_Dca", "FlowV0_Dca", nDcaBins, dcaMin, glDcaMax);
1060 fHistV0Dca->SetXTitle("dca between tracks (cm)");
1061 fHistV0Dca->SetYTitle("Counts");
1062 // lenght
1063 fHistV0Lenght = new TH1F("FlowV0_Lenght", "FlowV0_Lenght", nLgBins, lgMinV0, lgMaxV0);
1064 fHistV0Lenght->SetXTitle("Distance of V0s (cm)");
1065 fHistV0Lenght->SetYTitle("Counts");
1066 // Sigma for all particles
1067 fHistV0Sigma = new TH1F("FlowV0_Sigma", "FlowV0_Sigma", nLgBins, lgMinV0, lgMaxV0 );
1068 fHistV0Sigma->SetXTitle("Sigma");
1069 fHistV0Sigma->SetYTitle("Counts");
1070 // Chi2
1071 fHistV0Chi2 = new TH1F("FlowV0_Chi2", "FlowV0_Chi2", nChi2Bins, chi2Min, chi2MaxC);
1072 fHistV0Chi2->SetXTitle("Chi square at Main Vertex");
1073 fHistV0Chi2->SetYTitle("Counts");
1074 // EtaPtPhi
1075 fHistV0EtaPtPhi3D = new TH3F("FlowV0_EtaPtPhi3D", "FlowV0_EtaPtPhi3D", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, nPhi3DBins, fPhiMin, fPhiMax);
1076 fHistV0EtaPtPhi3D->SetXTitle("Eta");
1077 fHistV0EtaPtPhi3D->SetYTitle("Pt (GeV/c)");
1078 fHistV0EtaPtPhi3D->SetZTitle("Phi (rad)");
1079 // Yield for all v0s
1080 fHistV0YieldAll2D = new TH2D("FlowV0_YieldAll2D", "FlowV0_YieldAll2D", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
1081 fHistV0YieldAll2D->Sumw2();
1082 fHistV0YieldAll2D->SetXTitle("Pseudorapidty");
1083 fHistV0YieldAll2D->SetYTitle("Pt (GeV/c)");
1084 // Mass slices on pT
1085 fHistV0MassPtSlices = new TH2D("FlowV0_MassPtSlices", "FlowV0_MassPtSlices", nMassBins, massMin, massMax, fPtBins, fPtMin, fPtMax);
1086 fHistV0MassPtSlices->Sumw2();
1087 fHistV0MassPtSlices->SetXTitle("Invariant Mass (GeV)");
1088 fHistV0MassPtSlices->SetYTitle("Pt (GeV/c)");
1089 // Yield v0s (total P and Rapidity)
1090 fHistV0PYall2D = new TH2D("FlowV0_PYall2D", "FlowV0_PYall2D", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
1091 fHistV0PYall2D->Sumw2();
1092 fHistV0PYall2D->SetXTitle("Rapidty");
1093 fHistV0PYall2D->SetYTitle("P (GeV/c)");
1094
1095 // Selected V0s ...
1096 // Yield
1097 fHistV0YieldPart2D = new TH2D("FlowV0_Yield2Dsel", "FlowV0_Yield2Dsel", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
1098 fHistV0YieldPart2D->Sumw2();
1099 fHistV0YieldPart2D->SetXTitle("Pseudorapidty");
1100 fHistV0YieldPart2D->SetYTitle("Pt (GeV/c)");
1101 // Mass Window
1102 fHistV0MassWin = new TH1F("FlowV0_MassWinPart", "FlowV0_MassWinPart", nMassBins, massMin, massMax);
1103 fHistV0MassWin->SetXTitle("Invariant Mass (GeV)");
1104 fHistV0MassWin->SetYTitle("Counts");
1105 // EtaPtPhi
1106 fHistV0EtaPtPhi3DPart = new TH3F("FlowV0_EtaPtPhi3Dpart", "FlowV0_EtaPtPhi3Dpart", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, nPhi3DBins, fPhiMin, fPhiMax);
1107 fHistV0EtaPtPhi3DPart->SetXTitle("Eta");
1108 fHistV0EtaPtPhi3DPart->SetYTitle("Pt (GeV/c)");
1109 fHistV0EtaPtPhi3DPart->SetZTitle("Phi (rad)");
1110 // Distance of closest approach
1111 fHistV0DcaPart = new TH1F("FlowV0_DcaPart", "FlowV0_DcaPart", nDcaBins, dcaMin, dcaMax);
1112 fHistV0DcaPart->Sumw2();
1113 fHistV0DcaPart->SetXTitle("dca between tracks (cm)");
1114 fHistV0DcaPart->SetYTitle("Counts");
1115 // lenght
1116 fHistV0LenghtPart = new TH1F("FlowV0_LenghtPart", "FlowV0_LenghtPart", nLgBins, lgMinV0, lgMaxV0);
1117 fHistV0LenghtPart->SetXTitle("Distance of V0s (cm)");
1118 fHistV0LenghtPart->SetYTitle("Counts");
1119 // SideBand Mass (sidebands)
1120 fHistV0sbMassSide = new TH1F("FlowV0sb_MassWinSideBands", "FlowV0sb_MassWinSideBands", nMassBins, massMin, massMax);
1121 fHistV0sbMassSide->SetXTitle("Invariant Mass (GeV)");
1122 fHistV0sbMassSide->SetYTitle("Counts");
1123 // EtaPtPhi (sidebands)
1124 fHistV0sbEtaPtPhi3DPart = new TH3F("FlowV0sb_EtaPtPhi3D", "FlowV0sb_EtaPtPhi3D", fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, nPhi3DBins, fPhiMin, fPhiMax);
1125 fHistV0sbEtaPtPhi3DPart->SetXTitle("Eta");
1126 fHistV0sbEtaPtPhi3DPart->SetYTitle("Pt (GeV/c)");
1127 fHistV0sbEtaPtPhi3DPart->SetZTitle("Phi (rad)");
1128 // Distance of closest approach (sidebands)
1129 fHistV0sbDcaPart = new TH1F("FlowV0sb_Dca", "FlowV0sb_Dca", nDcaBins, dcaMin, dcaMax);
1130 fHistV0sbDcaPart->Sumw2();
1131 fHistV0sbDcaPart->SetXTitle("dca between tracks (cm)");
1132 fHistV0sbDcaPart->SetYTitle("Counts");
1133 // lenght (sidebands)
1134 fHistV0sbLenghtPart = new TH1F("FlowV0sb_Lenght", "FlowV0sb_Lenght", nLgBins, lgMinV0, lgMaxV0);
1135 fHistV0sbLenghtPart->SetXTitle("Distance of V0s (cm)");
1136 fHistV0sbLenghtPart->SetYTitle("Counts");
1137
1138 // Mean Eta in each bin
1139 fHistV0BinEta = new TProfile("FlowV0_Bin_Eta", "FlowV0_Bin_Eta", fEtaBins, fEtaMin, fEtaMax, fEtaMin, fEtaMax, "");
1140 fHistV0BinEta->SetXTitle((char*)fLabel.Data());
1141 fHistV0BinEta->SetYTitle("<Eta>");
1142 // Mean Pt in each bin
1143 fHistV0BinPt = new TProfile("FlowV0_Bin_Pt", "FlowV0_Bin_Pt", fPtBinsPart, fPtMin, fPtMaxPart, fPtMin, fPtMaxPart, "");
1144 fHistV0BinPt->SetXTitle("Pt (GeV/c)");
1145 fHistV0BinPt->SetYTitle("<Pt> (GeV/c)");
1146 // Mean Eta in each bin (sidebands)
1147 fHistV0sbBinEta = new TProfile("FlowV0sb_Bin_Eta", "FlowV0sb_Bin_Eta", fEtaBins, fEtaMin, fEtaMax, fEtaMin, fEtaMax, "");
1148 fHistV0sbBinEta->SetXTitle((char*)fLabel.Data());
1149 fHistV0sbBinEta->SetYTitle("<Eta>");
1150 // Mean Pt in each bin (sidebands)
1151 fHistV0sbBinPt = new TProfile("FlowV0sb_Bin_Pt", "FlowV0sb_Bin_Pt", fPtBinsPart, fPtMin, fPtMaxPart, fPtMin, fPtMaxPart, "");
1152 fHistV0sbBinPt->SetXTitle("Pt (GeV/c)");
1153 fHistV0sbBinPt->SetYTitle("<Pt> (GeV/c)");
1154 }
1155
1156 for (int k = 0; k < Flow::nSels; k++) // for each selection
1157 {
1158 // cos(n*delta_Psi)
1159 histTitle = new TString("Flow_Cos_Sel");
1160 *histTitle += k+1;
1161 fHistFull[k].fHistCos = new TProfile(histTitle->Data(), histTitle->Data(), Flow::nHars, 0.5, (float)(Flow::nHars) + 0.5, -1., 1., "");
1162 fHistFull[k].fHistCos->SetXTitle("Harmonic");
1163 fHistFull[k].fHistCos->SetYTitle("<cos(n*delta_Psi)>");
1164 delete histTitle;
1165
1166 // resolution
1167 histTitle = new TString("Flow_Res_Sel");
1168 *histTitle += k+1;
1169 fHistFull[k].fHistRes = new TH1F(histTitle->Data(), histTitle->Data(), Flow::nHars, 0.5, (float)(Flow::nHars) + 0.5);
1170 fHistFull[k].fHistRes->SetXTitle("Harmonic");
1171 fHistFull[k].fHistRes->SetYTitle("Resolution");
1172 delete histTitle;
1173
1174 // vObs
1175 histTitle = new TString("Flow_vObs_Sel");
1176 *histTitle += k+1;
1177 fHistFull[k].fHistvObs = new TProfile(histTitle->Data(), histTitle->Data(), Flow::nHars, 0.5, (float)(Flow::nHars) + 0.5, -100., 100., "");
1178 fHistFull[k].fHistvObs->SetXTitle("Harmonic");
1179 fHistFull[k].fHistvObs->SetYTitle("vObs (%)");
1180 delete histTitle;
1181
1182 // vObs V0
1183 histTitle = new TString("FlowV0_vObs_Sel");
1184 *histTitle += k+1;
1185 fHistFull[k].fHistV0vObs = new TProfile(histTitle->Data(), histTitle->Data(), Flow::nHars, 0.5, (float)(Flow::nHars) + 0.5, -100., 100., "");
1186 fHistFull[k].fHistV0vObs->SetXTitle("Harmonic");
1187 fHistFull[k].fHistV0vObs->SetYTitle("vObs (%)");
1188 delete histTitle;
1189
1190 // vObs V0 sideband SX
1191 histTitle = new TString("FlowV0sb_vObs_sx_Sel");
1192 *histTitle += k+1;
1193 fHistFull[k].fHistV0sbvObsSx = new TProfile(histTitle->Data(), histTitle->Data(), Flow::nHars, 0.5, (float)(Flow::nHars) + 0.5, -100., 100., "");
1194 fHistFull[k].fHistV0sbvObsSx->SetXTitle("Harmonic");
1195 fHistFull[k].fHistV0sbvObsSx->SetYTitle("vObs (%)");
1196 delete histTitle;
1197
1198 // vObs V0 sideband DX
1199 histTitle = new TString("FlowV0sb_vObs_dx_Sel");
1200 *histTitle += k+1;
1201 fHistFull[k].fHistV0sbvObsDx = new TProfile(histTitle->Data(), histTitle->Data(), Flow::nHars, 0.5, (float)(Flow::nHars) + 0.5, -100., 100., "");
1202 fHistFull[k].fHistV0sbvObsDx->SetXTitle("Harmonic");
1203 fHistFull[k].fHistV0sbvObsDx->SetYTitle("vObs (%)");
1204 delete histTitle;
1205
1206 // PID for tracks used in R.P.
1207 histTitle = new TString("Flow_BayPidMult_Sel");
1208 *histTitle += k+1;
1209 fHistFull[k].fHistBayPidMult = new TH1F(histTitle->Data(), histTitle->Data(),Flow::nPid,-0.5,((float)Flow::nPid-0.5));
1210 fHistFull[k].fHistBayPidMult->Sumw2() ;
1211 fHistFull[k].fHistBayPidMult->SetXTitle("e+/- , mu+/- , pi+/- , K+/- , p+/- , d+/- ");
1212 fHistFull[k].fHistBayPidMult->SetYTitle("Counts");
1213 delete histTitle;
1214
1215 for (int j = 0; j < Flow::nHars; j++) // for each harmonic
1216 {
1217 float order = (float)(j+1);
1218
1219 // multiplicity
1220 histTitle = new TString("Flow_Mul_Sel");
1221 *histTitle += k+1;
1222 histTitle->Append("_Har");
1223 *histTitle += j+1;
1224 fHistFull[k].fHistFullHar[j].fHistMult = new TH1F(histTitle->Data(),histTitle->Data(), nMultBins, multMin, multMax);
1225 fHistFull[k].fHistFullHar[j].fHistMult->SetXTitle("Multiplicity");
1226 fHistFull[k].fHistFullHar[j].fHistMult->SetYTitle("Counts");
1227 delete histTitle;
1228
1229 // event plane
1230 histTitle = new TString("Flow_Psi_Sel");
1231 *histTitle += k+1;
1232 histTitle->Append("_Har");
1233 *histTitle += j+1;
1234 fHistFull[k].fHistFullHar[j].fHistPsi = new TH1F(histTitle->Data(), histTitle->Data(), nPsiBins, psiMin, psiMax / order);
1235 fHistFull[k].fHistFullHar[j].fHistPsi->SetXTitle("Event Plane Angle (rad)");
1236 fHistFull[k].fHistFullHar[j].fHistPsi->SetYTitle("Counts");
1237 delete histTitle;
1238
1239 // event plane difference of two selections
1240 histTitle = new TString("Flow_Psi_Diff_Sel");
1241 *histTitle += k+1;
1242 histTitle->Append("_Har");
1243 *histTitle += j+1;
1244 if (k == 0 )
1245 {
1246 Int_t my_order = 1;
1247 if (j == 1) { my_order = 2 ; }
1248 fHistFull[k].fHistFullHar[j].fHistPsiDiff = new TH1F(histTitle->Data(), histTitle->Data(), nPsiBins, -psiMax/my_order/2., psiMax/my_order/2.);
1249 }
1250 else
1251 {
1252 fHistFull[k].fHistFullHar[j].fHistPsiDiff = new TH1F(histTitle->Data(), histTitle->Data(), nPsiBins, -psiMax/2., psiMax/2.);
1253 }
1254 if (k == 0)
1255 {
1256 if (j == 0)
1257 {
1258 fHistFull[k].fHistFullHar[j].fHistPsiDiff->SetXTitle("#Psi_{1,Sel1} - #Psi_{1,Sel2}(rad)");
1259 }
1260 else if (j == 1)
1261 {
1262 fHistFull[k].fHistFullHar[j].fHistPsiDiff->SetXTitle("#Psi_{2,Sel1} - #Psi_{2,Sel2}(rad)");
1263 }
1264 }
1265 else if (k == 1)
1266 {
1267 if (j == 0)
1268 {
1269 fHistFull[k].fHistFullHar[j].fHistPsiDiff->SetXTitle("#Psi_{1,Sel1} - #Psi_{2,Sel2}(rad)");
1270 }
1271 else if (j == 1)
1272 {
1273 fHistFull[k].fHistFullHar[j].fHistPsiDiff->SetXTitle("#Psi_{1,Sel1} - #Psi_{2,Sel1}(rad)");
1274 }
1275 }
1276 fHistFull[k].fHistFullHar[j].fHistPsiDiff->SetYTitle("Counts");
1277 delete histTitle;
1278
1279 // correlation of sub-event planes
1280 histTitle = new TString("Flow_Psi_Sub_Corr_Sel");
1281 *histTitle += k+1;
1282 histTitle->Append("_Har");
1283 *histTitle += j+1;
1284 fHistFull[k].fHistFullHar[j].fHistPsiSubCorr = new TH1F(histTitle->Data(), histTitle->Data(), nPsiBins, psiMin, psiMax / order);
1285 fHistFull[k].fHistFullHar[j].fHistPsiSubCorr->Sumw2();
1286 fHistFull[k].fHistFullHar[j].fHistPsiSubCorr->SetXTitle("Sub-Event Correlation (rad)");
1287 fHistFull[k].fHistFullHar[j].fHistPsiSubCorr->SetYTitle("Counts");
1288 delete histTitle;
1289
1290 // correlation of sub-event planes of different order
1291 histTitle = new TString("Flow_Psi_Sub_Corr_Diff_Sel");
1292 *histTitle += k+1;
1293 histTitle->Append("_Har");
1294 *histTitle += j+1;
1295 fHistFull[k].fHistFullHar[j].fHistPsiSubCorrDiff = new TH1F(histTitle->Data(), histTitle->Data(), nPsiBins, psiMin, psiMax / (order+1.));
1296 fHistFull[k].fHistFullHar[j].fHistPsiSubCorrDiff->Sumw2();
1297 fHistFull[k].fHistFullHar[j].fHistPsiSubCorrDiff->SetXTitle("Sub-Event Correlation (rad)");
1298 fHistFull[k].fHistFullHar[j].fHistPsiSubCorrDiff->SetYTitle("Counts");
1299 delete histTitle;
1300
1301 // q
1302 histTitle = new TString("Flow_NormQ_Sel");
1303 *histTitle += k+1;
1304 histTitle->Append("_Har");
1305 *histTitle += j+1;
1306 fHistFull[k].fHistFullHar[j].fHistQnorm = new TH1F(histTitle->Data(), histTitle->Data(), nQbins, qMin, qMax);
1307 fHistFull[k].fHistFullHar[j].fHistQnorm->Sumw2();
1308 fHistFull[k].fHistFullHar[j].fHistQnorm->SetXTitle("q = |Q|/sqrt(Mult)");
1309 fHistFull[k].fHistFullHar[j].fHistQnorm->SetYTitle("Counts");
1310 delete histTitle;
1311
1312 // particle-plane azimuthal correlation
1313 histTitle = new TString("Flow_Phi_Corr_Sel");
1314 *histTitle += k+1;
1315 histTitle->Append("_Har");
1316 *histTitle += j+1;
1317 fHistFull[k].fHistFullHar[j].fHistPhiCorr = new TH1F(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax / order);
1318 fHistFull[k].fHistFullHar[j].fHistPhiCorr->Sumw2();
1319 fHistFull[k].fHistFullHar[j].fHistPhiCorr->SetXTitle("Particle-Plane Correlation (rad)");
1320 fHistFull[k].fHistFullHar[j].fHistPhiCorr->SetYTitle("Counts");
1321 delete histTitle;
1322
1323 // neutral particle-plane azimuthal correlation
1324 histTitle = new TString("FlowV0_Phi_Corr_Sel");
1325 *histTitle += k+1;
1326 histTitle->Append("_Har");
1327 *histTitle += j+1;
1328 fHistFull[k].fHistFullHar[j].fHistV0PhiCorr = new TH1F(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax / order);
1329 fHistFull[k].fHistFullHar[j].fHistV0PhiCorr->Sumw2();
1330 fHistFull[k].fHistFullHar[j].fHistV0PhiCorr->SetXTitle("V0-Plane Correlation (rad)");
1331 fHistFull[k].fHistFullHar[j].fHistV0PhiCorr->SetYTitle("Counts");
1332 delete histTitle;
1333
1334 // neutral sidebands-plane azimuthal correlation
1335 histTitle = new TString("FlowV0sb_Phi_Corr_Sel");
1336 *histTitle += k+1;
1337 histTitle->Append("_Har");
1338 *histTitle += j+1;
1339 fHistFull[k].fHistFullHar[j].fHistV0sbPhiCorr = new TH1F(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax / order);
1340 fHistFull[k].fHistFullHar[j].fHistV0sbPhiCorr->Sumw2();
1341 fHistFull[k].fHistFullHar[j].fHistV0sbPhiCorr->SetXTitle("V0sideBands-Plane Correlation (rad)");
1342 fHistFull[k].fHistFullHar[j].fHistV0sbPhiCorr->SetYTitle("Counts");
1343 delete histTitle;
1344
1345 // Yield(pt)
1346 histTitle = new TString("Flow_YieldPt_Sel");
1347 *histTitle += k+1;
1348 histTitle->Append("_Har");
1349 *histTitle += j+1;
1350 fHistFull[k].fHistFullHar[j].fHistYieldPt = new TH1F(histTitle->Data(), histTitle->Data(), fPtBins, fPtMin, fPtMax);
1351 fHistFull[k].fHistFullHar[j].fHistYieldPt->Sumw2();
1352 fHistFull[k].fHistFullHar[j].fHistYieldPt->SetXTitle("Pt (GeV/c)");
1353 fHistFull[k].fHistFullHar[j].fHistYieldPt->SetYTitle("Yield");
1354 delete histTitle;
1355
1356 // EtaPtPhi
1357 histTitle = new TString("Flow_EtaPtPhi3D_Sel");
1358 *histTitle += k+1;
1359 histTitle->Append("_Har");
1360 *histTitle += j+1;
1361 fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3D = new TH3F(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, nPhi3DBins, fPhiMin, fPhiMax);
1362 fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3D->SetXTitle("Eta");
1363 fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3D->SetYTitle("Pt (GeV/c)");
1364 fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3D->SetZTitle("Phi (rad)");
1365 delete histTitle;
1366
1367 // Yield(eta,pt)
1368 histTitle = new TString("Flow_Yield2D_Sel");
1369 *histTitle += k+1;
1370 histTitle->Append("_Har");
1371 *histTitle += j+1;
1372 fHistFull[k].fHistFullHar[j].fHistYield2D = new TH2D(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
1373 fHistFull[k].fHistFullHar[j].fHistYield2D->Sumw2();
1374 fHistFull[k].fHistFullHar[j].fHistYield2D->SetXTitle("Pseudorapidty");
1375 fHistFull[k].fHistFullHar[j].fHistYield2D->SetYTitle("Pt (GeV/c)");
1376 delete histTitle;
1377
1378 // Dca - 3D
1379 histTitle = new TString("Flow_3dDca_Sel") ;
1380 *histTitle += k+1;
1381 histTitle->Append("_Har");
1382 *histTitle += j+1;
1383 fHistFull[k].fHistFullHar[j].fHistDcaGlob = new TH1F(histTitle->Data(), histTitle->Data(), nDcaBins, dcaMin, glDcaMax);
1384 fHistFull[k].fHistFullHar[j].fHistDcaGlob->Sumw2();
1385 fHistFull[k].fHistFullHar[j].fHistDcaGlob->SetXTitle("|3d Global Track's dca to Vertex (cm)|");
1386 delete histTitle;
1387
1388 // Yield(pt) - excluded from R.P.
1389 histTitle = new TString("Flow_YieldPt_outSel");
1390 *histTitle += k+1;
1391 histTitle->Append("_Har");
1392 *histTitle += j+1;
1393 fHistFull[k].fHistFullHar[j].fHistYieldPtout = new TH1F(histTitle->Data(), histTitle->Data(), fPtBins, fPtMin, fPtMax);
1394 fHistFull[k].fHistFullHar[j].fHistYieldPtout->Sumw2();
1395 fHistFull[k].fHistFullHar[j].fHistYieldPtout->SetXTitle("Pt (GeV/c)");
1396 fHistFull[k].fHistFullHar[j].fHistYieldPtout->SetYTitle("Yield");
1397 delete histTitle;
1398
1399 // EtaPtPhi - excluded from R.P.
1400 histTitle = new TString("Flow_EtaPtPhi3D_outSel");
1401 *histTitle += k+1;
1402 histTitle->Append("_Har");
1403 *histTitle += j+1;
1404 fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3Dout = new TH3F(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax, nPhi3DBins, fPhiMin, fPhiMax);
1405 fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3Dout->SetXTitle("Eta");
1406 fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3Dout->SetYTitle("Pt (GeV/c)");
1407 fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3Dout->SetZTitle("Phi (rad)");
1408 delete histTitle;
1409
1410 // Yield(eta,pt) - excluded from R.P.
1411 histTitle = new TString("Flow_Yield2D_outSel");
1412 *histTitle += k+1;
1413 histTitle->Append("_Har");
1414 *histTitle += j+1;
1415 fHistFull[k].fHistFullHar[j].fHistYield2Dout = new TH2D(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, fPtBins, fPtMin, fPtMax);
1416 fHistFull[k].fHistFullHar[j].fHistYield2Dout->Sumw2();
1417 fHistFull[k].fHistFullHar[j].fHistYield2Dout->SetXTitle("Pseudorapidty");
1418 fHistFull[k].fHistFullHar[j].fHistYield2Dout->SetYTitle("Pt (GeV/c)");
1419 delete histTitle;
1420
1421 // Dca - 3D - excluded from R.P.
1422 histTitle = new TString("Flow_3dDca_outSel") ;
1423 *histTitle += k+1;
1424 histTitle->Append("_Har");
1425 *histTitle += j+1;
1426 fHistFull[k].fHistFullHar[j].fHistDcaGlobout = new TH1F(histTitle->Data(), histTitle->Data(), nDcaBins, dcaMin, glDcaMax);
1427 fHistFull[k].fHistFullHar[j].fHistDcaGlobout->Sumw2();
1428 fHistFull[k].fHistFullHar[j].fHistDcaGlobout->SetXTitle("|3d Global Track's dca to Vertex (cm)|");
1429 delete histTitle;
1430
1431 // Flow observed - v_obs,Pt,Eta
1432 histTitle = new TString("Flow_vObs2D_Sel");
1433 *histTitle += k+1;
1434 histTitle->Append("_Har");
1435 *histTitle += j+1;
1436 fHistFull[k].fHistFullHar[j].fHistvObs2D = new TProfile2D(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, fPtBinsPart, fPtMin, fPtMaxPart, -100., 100., "");
1437 fHistFull[k].fHistFullHar[j].fHistvObs2D->SetXTitle((char*)fLabel.Data());
1438 fHistFull[k].fHistFullHar[j].fHistvObs2D->SetYTitle("Pt (GeV/c)");
1439 delete histTitle;
1440
1441 // v_obs,Eta
1442 histTitle = new TString("Flow_vObsEta_Sel");
1443 *histTitle += k+1;
1444 histTitle->Append("_Har");
1445 *histTitle += j+1;
1446 fHistFull[k].fHistFullHar[j].fHistvObsEta = new TProfile(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, -100., 100., "");
1447 fHistFull[k].fHistFullHar[j].fHistvObsEta->SetXTitle((char*)fLabel.Data());
1448 fHistFull[k].fHistFullHar[j].fHistvObsEta->SetYTitle("v (%)");
1449 delete histTitle;
1450
1451 // v_obs,Pt
1452 histTitle = new TString("Flow_vObsPt_Sel");
1453 *histTitle += k+1;
1454 histTitle->Append("_Har");
1455 *histTitle += j+1;
1456 fHistFull[k].fHistFullHar[j].fHistvObsPt = new TProfile(histTitle->Data(), histTitle->Data(), fPtBinsPart, fPtMin, fPtMaxPart, -100., 100., "");
1457 fHistFull[k].fHistFullHar[j].fHistvObsPt->SetXTitle("Pt (GeV/c)");
1458 fHistFull[k].fHistFullHar[j].fHistvObsPt->SetYTitle("v (%)");
1459 delete histTitle;
1460
1461 // neutral Flow observed - Pt,Eta
1462 histTitle = new TString("FlowV0_vObs2D_Sel");
1463 *histTitle += k+1;
1464 histTitle->Append("_Har");
1465 *histTitle += j+1;
1466 fHistFull[k].fHistFullHar[j].fHistV0vObs2D = new TProfile2D(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, fPtBinsPart, fPtMin, fPtMaxPart, -100., 100., "");
1467 fHistFull[k].fHistFullHar[j].fHistV0vObs2D->SetXTitle((char*)fLabel.Data());
1468 fHistFull[k].fHistFullHar[j].fHistV0vObs2D->SetYTitle("Pt (GeV/c)");
1469 delete histTitle;
1470
1471 // neutral Flow observed - Eta
1472 histTitle = new TString("FlowV0_vObsEta_Sel");
1473 *histTitle += k+1;
1474 histTitle->Append("_Har");
1475 *histTitle += j+1;
1476 fHistFull[k].fHistFullHar[j].fHistV0vObsEta = new TProfile(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, -100., 100., "");
1477 fHistFull[k].fHistFullHar[j].fHistV0vObsEta->SetXTitle((char*)fLabel.Data());
1478 fHistFull[k].fHistFullHar[j].fHistV0vObsEta->SetYTitle("v (%)");
1479 delete histTitle;
1480
1481 // neutral Flow observed - Pt
1482 histTitle = new TString("FlowV0_vObsPt_Sel");
1483 *histTitle += k+1;
1484 histTitle->Append("_Har");
1485 *histTitle += j+1;
1486 fHistFull[k].fHistFullHar[j].fHistV0vObsPt = new TProfile(histTitle->Data(), histTitle->Data(), fPtBinsPart, fPtMin, fPtMaxPart, -100., 100., "");
1487 fHistFull[k].fHistFullHar[j].fHistV0vObsPt->SetXTitle("Pt (GeV/c)");
1488 fHistFull[k].fHistFullHar[j].fHistV0vObsPt->SetYTitle("v (%)");
1489 delete histTitle;
1490
1491 // neutral sidebands Flow observed - Pt,Eta
1492 histTitle = new TString("FlowV0sb_vObs2D_Sel");
1493 *histTitle += k+1;
1494 histTitle->Append("_Har");
1495 *histTitle += j+1;
1496 fHistFull[k].fHistFullHar[j].fHistV0sbvObs2D = new TProfile2D(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, fPtBinsPart, fPtMin, fPtMaxPart, -100., 100., "");
1497 fHistFull[k].fHistFullHar[j].fHistV0sbvObs2D->SetXTitle((char*)fLabel.Data());
1498 fHistFull[k].fHistFullHar[j].fHistV0sbvObs2D->SetYTitle("Pt (GeV/c)");
1499 delete histTitle;
1500
1501 // neutral sidebands Flow observed - Eta
1502 histTitle = new TString("FlowV0sb_vObsEta_Sel");
1503 *histTitle += k+1;
1504 histTitle->Append("_Har");
1505 *histTitle += j+1;
1506 fHistFull[k].fHistFullHar[j].fHistV0sbvObsEta = new TProfile(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, -100., 100., "");
1507 fHistFull[k].fHistFullHar[j].fHistV0sbvObsEta->SetXTitle((char*)fLabel.Data());
1508 fHistFull[k].fHistFullHar[j].fHistV0sbvObsEta->SetYTitle("v (%)");
1509 delete histTitle;
1510
1511 // neutral sidebands Flow observed - Pt
1512 histTitle = new TString("FlowV0sb_vObsPt_Sel");
1513 *histTitle += k+1;
1514 histTitle->Append("_Har");
1515 *histTitle += j+1;
1516 fHistFull[k].fHistFullHar[j].fHistV0sbvObsPt = new TProfile(histTitle->Data(), histTitle->Data(), fPtBinsPart, fPtMin, fPtMaxPart, -100., 100., "");
1517 fHistFull[k].fHistFullHar[j].fHistV0sbvObsPt->SetXTitle("Pt (GeV/c)");
1518 fHistFull[k].fHistFullHar[j].fHistV0sbvObsPt->SetYTitle("v (%)");
1519 delete histTitle;
1520
1521 // SX neutral sidebands Flow observed - Eta
1522 histTitle = new TString("FlowV0sb_vObsEta_sx_Sel");
1523 *histTitle += k+1;
1524 histTitle->Append("_Har");
1525 *histTitle += j+1;
1526 fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaSx = new TProfile(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, -100., 100., "");
1527 fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaSx->SetXTitle((char*)fLabel.Data());
1528 fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaSx->SetYTitle("v (%)");
1529 delete histTitle;
1530
1531 // SX neutral sidebands Flow observed - Pt
1532 histTitle = new TString("FlowV0sb_vObsPt_sx_Sel");
1533 *histTitle += k+1;
1534 histTitle->Append("_Har");
1535 *histTitle += j+1;
1536 fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtSx = new TProfile(histTitle->Data(), histTitle->Data(), fPtBinsPart, fPtMin, fPtMaxPart, -100., 100., "");
1537 fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtSx->SetXTitle("Pt (GeV/c)");
1538 fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtSx->SetYTitle("v (%)");
1539 delete histTitle;
1540
1541 // DX neutral sidebands Flow observed - Eta
1542 histTitle = new TString("FlowV0sb_vObsEta_dx_Sel");
1543 *histTitle += k+1;
1544 histTitle->Append("_Har");
1545 *histTitle += j+1;
1546 fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaDx = new TProfile(histTitle->Data(), histTitle->Data(), fEtaBins, fEtaMin, fEtaMax, -100., 100., "");
1547 fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaDx->SetXTitle((char*)fLabel.Data());
1548 fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaDx->SetYTitle("v (%)");
1549 delete histTitle;
1550
1551 // DX neutral sidebands Flow observed - Pt
1552 histTitle = new TString("FlowV0sb_vObsPt_dx_Sel");
1553 *histTitle += k+1;
1554 histTitle->Append("_Har");
1555 *histTitle += j+1;
1556 fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtDx = new TProfile(histTitle->Data(), histTitle->Data(), fPtBinsPart, fPtMin, fPtMaxPart, -100., 100., "");
1557 fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtDx->SetXTitle("Pt (GeV/c)");
1558 fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtDx->SetYTitle("v (%)");
1559 delete histTitle;
1560
1561 // Phi lab
1562 // Tpc (plus)
1563 histTitle = new TString("Flow_Phi_TPCplus_Sel");
1564 *histTitle += k+1;
1565 histTitle->Append("_Har");
1566 *histTitle += j+1;
1567 fHistFull[k].fHistFullHar[j].fHistPhiPlus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1568 fHistFull[k].fHistFullHar[j].fHistPhiPlus->SetXTitle("Azimuthal Angles (rad)");
1569 fHistFull[k].fHistFullHar[j].fHistPhiPlus->SetYTitle("Counts");
1570 delete histTitle;
1571
1572 // Tpc (minus)
1573 histTitle = new TString("Flow_Phi_TPCminus_Sel");
1574 *histTitle += k+1;
1575 histTitle->Append("_Har");
1576 *histTitle += j+1;
1577 fHistFull[k].fHistFullHar[j].fHistPhiMinus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1578 fHistFull[k].fHistFullHar[j].fHistPhiMinus->SetXTitle("Azimuthal Angles (rad)");
1579 fHistFull[k].fHistFullHar[j].fHistPhiMinus->SetYTitle("Counts");
1580 delete histTitle;
1581
1582 // Tpc (cross)
1583 histTitle = new TString("Flow_Phi_TPCcross_Sel");
1584 *histTitle += k+1;
1585 histTitle->Append("_Har");
1586 *histTitle += j+1;
1587 fHistFull[k].fHistFullHar[j].fHistPhiAll = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1588 fHistFull[k].fHistFullHar[j].fHistPhiAll->SetXTitle("Azimuthal Angles (rad)");
1589 fHistFull[k].fHistFullHar[j].fHistPhiAll->SetYTitle("Counts");
1590 delete histTitle;
1591
1592 // Tpc
1593 histTitle = new TString("Flow_Phi_TPC_Sel");
1594 *histTitle += k+1;
1595 histTitle->Append("_Har");
1596 *histTitle += j+1;
1597 fHistFull[k].fHistFullHar[j].fHistPhi = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1598 fHistFull[k].fHistFullHar[j].fHistPhi->SetXTitle("Azimuthal Angles (rad)");
1599 fHistFull[k].fHistFullHar[j].fHistPhi->SetYTitle("Counts");
1600 delete histTitle;
1601
1602 // Phi lab flattened
1603 // Tpc (Plus)
1604 histTitle = new TString("Flow_Phi_Flat_TPCplus_Sel");
1605 *histTitle += k+1;
1606 histTitle->Append("_Har");
1607 *histTitle += j+1;
1608 fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1609 fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus->SetXTitle("Azimuthal Angles (rad)");
1610 fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus->SetYTitle("Counts");
1611 delete histTitle;
1612
1613 // Tpc (Minus)
1614 histTitle = new TString("Flow_Phi_Flat_TPCminus_Sel");
1615 *histTitle += k+1;
1616 histTitle->Append("_Har");
1617 *histTitle += j+1;
1618 fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1619 fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus->SetXTitle("Azimuthal Angles (rad)");
1620 fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus->SetYTitle("Counts");
1621 delete histTitle;
1622
1623 // Tpc (cross)
1624 histTitle = new TString("Flow_Phi_Flat_TPCcross_Sel");
1625 *histTitle += k+1;
1626 histTitle->Append("_Har");
1627 *histTitle += j+1;
1628 fHistFull[k].fHistFullHar[j].fHistPhiFlatAll = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1629 fHistFull[k].fHistFullHar[j].fHistPhiFlatAll->SetXTitle("Azimuthal Angles (rad)");
1630 fHistFull[k].fHistFullHar[j].fHistPhiFlatAll->SetYTitle("Counts");
1631 delete histTitle;
1632
1633 // Tpc
1634 histTitle = new TString("Flow_Phi_Flat_TPC_Sel");
1635 *histTitle += k+1;
1636 histTitle->Append("_Har");
1637 *histTitle += j+1;
1638 fHistFull[k].fHistFullHar[j].fHistPhiFlat = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1639 fHistFull[k].fHistFullHar[j].fHistPhiFlat->SetXTitle("Azimuthal Angles (rad)");
1640 fHistFull[k].fHistFullHar[j].fHistPhiFlat->SetYTitle("Counts");
1641 delete histTitle;
1642 }
1643 }
1644 cout << " Init() - Histograms booked" << endl ;
1645
1646 return kTRUE ;
1647}
1648//-----------------------------------------------------------------------
1649
1650//-----------------------------------------------------------------------
1651Bool_t AliFlowAnalyser::Finish()
1652{
1653 // Close the analysis and saves the histograms on the histFile .
1654
1655 cout << "* FlowAnalysis * - Finish()" << endl ; cout << endl ;
1656
1657 // Write all histograms
1658 fHistFile->cd() ;
1659 fHistFile->Write() ;
1660
1661 // Write Resolution corrected histograms
1662 if(fVnResHistList)
1663 {
1664 fVnResHistList->Write();
1665 delete fVnResHistList ;
1666 }
1667 else { cout << " E.P. resolution has not been calculated. No v_n histograms!" << endl ; }
1668
1669 // Write PhiWgt histograms
1670 if(fPhiWgtHistList)
1671 {
1672 fPhiWgtHistList->Write();
1673 delete fPhiWgtHistList ;
1674 }
1675
1676 fFlowSelect->Write();
1677 // delete fFlowSelect ;
1678
1679 fHistFile->Close() ;
1680
1681 cout << " Finish() - Histograms saved : " << fHistFileName.Data() << endl ; cout << endl ;
1682
1683 return kTRUE ;
1684}
1685//-----------------------------------------------------------------------
1686// ###
1687//----------------------------------------------------------------------
1688Float_t AliFlowAnalyser::GetRunBayesian(Int_t nPid, Int_t selN)
1689{
1690 // Returns the normalized particle abundance of "e","mu","pi","k","p","d"
1691 // in all the analysed events (in selection selN).
1692 // Call at the end of the analysis.
1693
1694 if(selN>Flow::nSels) { selN = 0 ; }
1695 Double_t totCount = (fHistFull[selN].fHistBayPidMult)->GetSumOfWeights() ;
1696 if(totCount) { return (fHistFull[selN].fHistBayPidMult->GetBinContent(nPid+1) / totCount) ; }
1697 else { return 1. ; }
1698}
1699//-----------------------------------------------------------------------
1700void AliFlowAnalyser::PrintRunBayesian(Int_t selN)
1701{
1702 // Prints the normalized particle abundance of all the analysed events
1703 // (in selection selN).
1704
1705 if(selN>Flow::nSels) { selN = 0 ; }
1706 Char_t* names[Flow::nPid] = {"e","mu","pi","k","p","d"} ;
1707 Double_t bayes = 0. ;
1708 cout << " selN = " << selN << " particles normalized abundance : " ;
1709 for(int i=0;i<Flow::nPid;i++)
1710 {
1711 bayes = GetRunBayesian(i, selN) ;
1712 cout << bayes << "_" << names[i] << " ; " ;
1713 }
1714 cout << endl ;
1715
1716 return ;
1717}
1718//-----------------------------------------------------------------------
1719void AliFlowAnalyser::FillWgtArrays(TFile* wgtFile)
1720{
1721 // Loads PhiWeights & Bayesian particles' abundance from file (default: flowPhiWgt.hist.root).
1722 // Weights are stored in a static TH1D* data member, ready to be plugged into the AliFlowEvent.
1723 // The plugging is done by the method ::FillEvtPhiWgt() (if wgt file is there).
1724
1725 fPhiWgtFile = wgtFile ;
1726
1727 TString* histTitle ;
1728 TH1D* TPC_all ; TH1D* TPC_plus ; TH1D* TPC_minus ; TH1D* TPC_cross ;
1729 TH1D* PID_bay ;
1730 for(int k=0;k<Flow::nSels;k++)
1731 {
1732 for(int j=0;j<Flow::nHars;j++)
1733 {
1734 // Tpc (plus)
1735 histTitle = new TString("Flow_Phi_Weight_TPCplus_Sel");
1736 *histTitle += k+1;
1737 histTitle->Append("_Har");
1738 *histTitle += j+1;
1739 TPC_plus = (TH1D*)fPhiWgtFile->Get(histTitle->Data());
1740 delete histTitle;
1741 // Tpc (minus)
1742 histTitle = new TString("Flow_Phi_Weight_TPCminus_Sel");
1743 *histTitle += k+1;
1744 histTitle->Append("_Har");
1745 *histTitle += j+1;
1746 TPC_minus = (TH1D*)fPhiWgtFile->Get(histTitle->Data());
1747 delete histTitle;
1748 // Tpc (cross)
1749 histTitle = new TString("Flow_Phi_Weight_TPCcross_Sel");
1750 *histTitle += k+1;
1751 histTitle->Append("_Har");
1752 *histTitle += j+1;
1753 TPC_cross = (TH1D*)fPhiWgtFile->Get(histTitle->Data());
1754 delete histTitle;
1755
1756 // Tpc
1757 histTitle = new TString("Flow_Phi_Weight_TPC_Sel");
1758 *histTitle += k+1;
1759 histTitle->Append("_Har");
1760 *histTitle += j+1;
1761 TPC_all = (TH1D*)fPhiWgtFile->Get(histTitle->Data());
1762 delete histTitle;
1763
1764 for(int n=0;n<fPhiBins;n++)
1765 {
1766 fPhiWgtPlus[k][j][n] = TPC_plus->GetBinContent(n+1) ;
1767 fPhiWgtMinus[k][j][n] = TPC_minus->GetBinContent(n+1) ;
1768 fPhiWgtCross[k][j][n] = TPC_cross->GetBinContent(n+1) ;
1769 fPhiWgt[k][j][n] = TPC_all->GetBinContent(n+1) ;
1770 // cout << " Weights: " << fPhiWgt[k][j][n] << " ; " << fPhiWgtPlus[k][j][n] << " | " << fPhiWgtMinus[k][j][n] << " | " << fPhiWgtCross[k][j][n] << endl ;
1771 }
1772 }
1773
1774 // Bayesian weights
1775 histTitle = new TString("Flow_BayPidMult_Sel");
1776 *histTitle += k+1;
1777 PID_bay = (TH1D*)fPhiWgtFile->Get(histTitle->Data());
1778 delete histTitle;
1779 Double_t totCount = PID_bay->GetSumOfWeights() ;
1780 for (int n=0;n<Flow::nPid;n++)
1781 {
1782 if(totCount) { fBayesianWgt[k][n] = PID_bay->GetBinContent(n+1) / totCount ; }
1783 else { fBayesianWgt[k][n] = 1. ; }
1784 // cout << " Bayesian Weights (" << n << ") : " << fBayesianWgt[k][n] << endl ;
1785 }
1786 }
1787
1788 delete TPC_all ; delete TPC_plus ; delete TPC_minus ; delete TPC_cross ;
1789 delete PID_bay ;
1790
1791 return ;
1792}
1793//-----------------------------------------------------------------------
1794void AliFlowAnalyser::FillEvtPhiWgt(AliFlowEvent* fFlowEvent)
1795{
1796 // Plugs phi weights into the static dwgt data member of the AliFlowEvent class.
1797 // Weights are given in special Flow::PhiWgt_t arrays (see AliFlowConstants),
1798 // which are read from the wgt histograms by the method FillWgtArrays(...).
1799
1800 fFlowEvent->SetPhiWeight(fPhiWgt);
1801 fFlowEvent->SetPhiWeightPlus(fPhiWgtPlus);
1802 fFlowEvent->SetPhiWeightMinus(fPhiWgtMinus);
1803 fFlowEvent->SetPhiWeightCross(fPhiWgtCross);
1804
1805 return ;
1806}
1807//-----------------------------------------------------------------------
1808void AliFlowAnalyser::FillBayesianWgt(AliFlowEvent* fFlowEvent)
1809{
1810 // Plugs Bayesian particle abundance into the current AliFlowEvent.
1811 // A bayesian vector should be used for the PId of any different selection
1812 // (different sets of cuts -> different particle abundance), but for now
1813 // just the Selection n.0 (with no cuts) is used .
1814 // (AliFlowEvent::mBayesianCs[6] is a 1-dimensional array, change that first!).
1815
1816 Double_t bayes[Flow::nPid] ;
1817 Double_t bayCheck = 0. ;
1818 for (int n=0;n<Flow::nPid;n++)
1819 {
1820 bayes[n] = fBayesianWgt[0][n] ;
1821 bayCheck += bayes[n] ;
1822 // cout << "Bayesian V[" << n << "] = " << fBayesianWgt[0][n] << endl ;
1823 }
1824 if(bayCheck) { fFlowEvent->SetBayesian(bayes) ; fRePid = kTRUE ; }
1825 else { cout << "An empty bayesian vector is stored !!! - Bayesian weights = {1,1,1,1,1,1} " << endl ; }
1826
1827 return ;
1828}
1829//-----------------------------------------------------------------------
1830void AliFlowAnalyser::Weightening()
1831{
1832 // Calculates weights, and fills PhiWgt histograms .
1833 // This is called at the end of the event analysis.
1834
1835 cout << " AliFlowAnalyser::Weightening() " << endl ; cout << endl ;
1836
1837 // PhiWgt histogram collection
1838 fPhiWgtHistList = new TOrdCollection(4*Flow::nSels*Flow::nHars) ;
1839
1840 // Creates PhiWgt Histograms
1841 TString* histTitle ;
1842 for(int k = 0; k < Flow::nSels; k++)
1843 {
1844 for(int j = 0; j < Flow::nHars; j++)
1845 {
1846 // Tpc (plus)
1847 histTitle = new TString("Flow_Phi_Weight_TPCplus_Sel");
1848 *histTitle += k+1;
1849 histTitle->Append("_Har");
1850 *histTitle += j+1;
1851 fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1852 fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->Sumw2();
1853 fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetXTitle("Azimuthal Angles (rad)");
1854 fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetYTitle("PhiWgt");
1855 delete histTitle;
1856 // Tpc (minus)
1857 histTitle = new TString("Flow_Phi_Weight_TPCminus_Sel");
1858 *histTitle += k+1;
1859 histTitle->Append("_Har");
1860 *histTitle += j+1;
1861 fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1862 fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->Sumw2();
1863 fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetXTitle("Azimuthal Angles (rad)");
1864 fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetYTitle("PhiWgt");
1865 delete histTitle;
1866 // Tpc (cross)
1867 histTitle = new TString("Flow_Phi_Weight_TPCcross_Sel");
1868 *histTitle += k+1;
1869 histTitle->Append("_Har");
1870 *histTitle += j+1;
1871 fHistFull[k].fHistFullHar[j].fHistPhiWgtAll = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1872 fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->Sumw2();
1873 fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetXTitle("Azimuthal Angles (rad)");
1874 fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetYTitle("PhiWgt");
1875 delete histTitle;
1876 // Tpc
1877 histTitle = new TString("Flow_Phi_Weight_TPC_Sel");
1878 *histTitle += k+1;
1879 histTitle->Append("_Har");
1880 *histTitle += j+1;
1881 fHistFull[k].fHistFullHar[j].fHistPhiWgt = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
1882 fHistFull[k].fHistFullHar[j].fHistPhiWgt->Sumw2();
1883 fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetXTitle("Azimuthal Angles (rad)");
1884 fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetYTitle("PhiWgt");
1885 delete histTitle;
1886
1887 // Calculate PhiWgt
1888 double meanPlus = fHistFull[k].fHistFullHar[j].fHistPhiPlus->Integral() / (double)fPhiBins ;
1889 double meanMinus = fHistFull[k].fHistFullHar[j].fHistPhiMinus->Integral() / (double)fPhiBins ;
1890 double meanCross = fHistFull[k].fHistFullHar[j].fHistPhiAll->Integral() / (double)fPhiBins ;
1891 double meanTPC = fHistFull[k].fHistFullHar[j].fHistPhi->Integral() / (double)fPhiBins ;
1892
1893 // Tpc
1894 for (int i=0;i<fPhiBins;i++)
1895 {
1896 fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetBinContent(i+1,meanPlus);
1897 fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetBinError(i+1, 0.);
1898 fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetBinContent(i+1,meanMinus);
1899 fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetBinError(i+1, 0.);
1900 fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetBinContent(i+1,meanCross);
1901 fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetBinError(i+1, 0.);
1902 fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetBinContent(i+1,meanTPC);
1903 fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetBinError(i+1, 0.);
1904 }
1905
1906 if(meanTPC==0) { cout << " Sel." << k << " , Har." << j << " : empty phi histogram ! " << endl ; }
1907 else
1908 {
1909 fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->Divide(fHistFull[k].fHistFullHar[j].fHistPhiPlus);
1910 fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->Divide(fHistFull[k].fHistFullHar[j].fHistPhiMinus);
1911 fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->Divide(fHistFull[k].fHistFullHar[j].fHistPhiAll);
1912 fHistFull[k].fHistFullHar[j].fHistPhiWgt->Divide(fHistFull[k].fHistFullHar[j].fHistPhi);
1913 }
1914
1915 fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus);
1916 fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus);
1917 fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtAll);
1918 fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgt);
1919 }
1920 }
1921
1922 return ;
1923}
1924//-----------------------------------------------------------------------
1925// ###
1926//-----------------------------------------------------------------------
1927Bool_t AliFlowAnalyser::Analyse(AliFlowEvent* flowEvent)
1928{
1929 // Runs the analysis on the AliFlowEvent (* fFlowEvent).
1930 // This method can be inserted in a loop over a collection of
1931 // AliFlowEvents or for on-fly analysis on AliESDs if used toghether
1932 // with the AliFlowMaker.
1933
1934 cout << " AliFlowAnalyser::Analyze(" << fFlowEvent << " ) - " << fEventNumber << endl ;
1935 if(!flowEvent) { return kFALSE ; }
1936 else { fFlowEvent = flowEvent ; }
1937
1938 if(fFlowSelect->Select(fFlowEvent)) // event selected - here below the ANALYSIS FLAGS are setted -
1939 {
1940 cout << " * 1 . Load event (track & v0s) and set flags . " << endl ;
1941 fFlowTracks = fFlowEvent->TrackCollection() ;
1942 fNumberOfTracks = fFlowTracks->GetEntries() ;
1943 fFlowV0s = fFlowEvent->V0Collection() ;
1944 fNumberOfV0s = fFlowV0s->GetEntries() ;
1945 cout << " event ID = " << fFlowEvent->EventID() << " : found " << fNumberOfTracks << " AliFlowTracks, and " << fNumberOfV0s << " AliFlowV0s . " << endl ;
1946
1947 if(fReadPhiWgt) { FillEvtPhiWgt(fFlowEvent) ; } // phi and bayesian weights are filled previous to the loop (FillWgtArrays(TFile*))
1948 else { fFlowEvent->SetNoWgt() ; } // phi weights can be used or not , this plugs them into the event
1949
1950 if(fBayWgt) { FillBayesianWgt(fFlowEvent) ; } // bayesian weights can be used or not , this plugs them into the event
1951 if(fRePid) { fFlowEvent->SetPids() ; } // re-calculate all p.id. hypotesis with the (new) bayesian array
1952
1953 if(fOnePhiWgt) { fFlowEvent->SetOnePhiWgt() ; } // one phi-wgt histogram
1954 else { fFlowEvent->SetFirstLastPhiWgt() ; } // three phi-wgt histogram
1955 if(fPtWgt) { fFlowEvent->SetPtWgt(); ; } // pT as a weight
1956 if(fEtaWgt) { fFlowEvent->SetEtaWgt() ; } // eta as a weight
1957
1958 if(fShuffle) { fFlowEvent->RandomShuffle() ; } // tracks re-shuffling
1959
1960 fFlowEvent->SetSelections(fFlowSelect) ; // does the selection of tracks for r.p. calculation (sets flags in AliFlowTrack)
1961 fFlowEvent->SetEtaSubs(fEtaSub) ; // setting for the subevents (eta or random)
1962 fFlowEvent->MakeSubEvents() ; // makes the subevent, eta or random basing on the previous flag
1963
1964 cout << " * 2 . Calculating event quantities all in one shoot . " << endl ;
1965 fFlowEvent->MakeAll() ;
1966
1967 if(FillFromFlowEvent(fFlowEvent)) // calculates event quantities
1968 {
1969 cout << " * 3 . Event Histograms and Particles loop . " << endl ;
1970 FillEventHistograms(fFlowEvent); // fill histograms from AliFlowEvents
1971 if(fTrackLoop) { FillParticleHistograms(fFlowTracks) ; } // fill histograms from AliFlowTracks
1972 if(fV0loop) { FillV0Histograms(fFlowV0s) ; } // fill histograms from AliFlowV0s
1973 //FillLabels() ; // fill the histogram of MC labels (from the simulation)
1974 }
1975 else
1976 {
1977 cout << " * 3 . Event psi = 0 - Skipping! " << endl ;
1978 return kFALSE ;
1979 }
1980 }
1981 else
1982 {
1983 cout << " * 0 . Event " << fEventNumber << " (event ID = " << fFlowEvent->EventID() << ") discarded . " << endl ;
1984 delete fFlowEvent ; fFlowEvent = 0 ;
1985 return kFALSE ;
1986 }
1987 fEventNumber++ ;
1988
1989 return kTRUE ;
1990}
1991//-----------------------------------------------------------------------
1992Bool_t AliFlowAnalyser::FillFromFlowEvent(AliFlowEvent* fFlowEvent)
1993{
1994 // gets event quantities, returns kFALSE if Q vector is always 0
1995
1996 Int_t selCheck = 0 ;
1997 for(Int_t k = 0; k < Flow::nSels; k++)
1998 {
1999 fFlowSelect->SetSelection(k) ;
2000 for(Int_t j = 0; j < Flow::nHars; j++)
2001 {
2002 fFlowSelect->SetHarmonic(j) ;
2003 for(Int_t n = 0; n < Flow::nSubs; n++)
2004 {
2005 fFlowSelect->SetSubevent(n) ;
2006 fPsiSub[n][k][j] = fFlowEvent->Psi(fFlowSelect) ; // sub-event quantities
2007 fMultSub[n][k][j] = fFlowEvent->Mult(fFlowSelect) ;
2008 }
2009 fFlowSelect->SetSubevent(-1);
2010 fQ[k][j] = fFlowEvent->Q(fFlowSelect) ; // full event quantities
2011 fPsi[k][j] = fFlowEvent->Psi(fFlowSelect) ;
2012 fQnorm[k][j] = fFlowEvent->NormQ(fFlowSelect).Mod() ; // was: fFlowEvent->q(fFlowSelect) ; // but the normalization was bad (no pT,eta weight)
2013 fMult[k][j] = fFlowEvent->Mult(fFlowSelect) ;
2014 selCheck += fMult[k][j] ;
2015 }
2016 }
2017 if(!selCheck) { return kFALSE ; } // if there are no particles in the selection -> skip the event
2018
2019 return kTRUE ;
2020}
2021//-----------------------------------------------------------------------
2022void AliFlowAnalyser::FillEventHistograms(AliFlowEvent* fFlowEvent)
2023{
2024 // fill event histograms
2025
2026 cout << " Fill Event Histograms ... " << endl ;
2027
2028 float trigger = (float)fFlowEvent->L0TriggerWord() ;
2029 fHistTrigger->Fill(trigger);
2030
2031 // no selections: OrigMult, Centrality, Mult, MultOverOrig, VertexZ, VertexXY
2032 int origMult = fFlowEvent->OrigMult();
2033 fHistOrigMult->Fill((float)origMult);
2034 fHistMultEta->Fill((float)fFlowEvent->MultEta());
2035
2036 int cent = fFlowEvent->Centrality();
2037 fHistCent->Fill((float)cent);
2038
2039 fHistMult->Fill((float)fNumberOfTracks) ;
2040 fHistV0Mult->Fill((float)fNumberOfV0s) ;
2041 if(origMult) { fHistMultOverOrig->Fill((float)fNumberOfTracks/(float)origMult) ; }
2042
2043 fFlowEvent->VertexPos(fVertex) ;
2044 fHistVertexZ->Fill(fVertex[2]) ;
2045 fHistVertexXY2D->Fill(fVertex[0],fVertex[1]) ;
2046
2047 // ZDC info
2048 fHistPartZDC->Fill(fFlowEvent->ZDCpart()) ;
2049 for(int ii=0;ii<3;ii++) { fHistEnergyZDC->Fill(ii,fFlowEvent->ZDCenergy(ii)) ; }
2050
2051 // sub-event Psi_Subs
2052 for(int k = 0; k < Flow::nSels; k++)
2053 {
2054 for(int j = 0; j < Flow::nHars; j++)
2055 {
2056 for(int n = 0; n < Flow::nSubs; n++)
2057 {
2058 int iii = Flow::nSubs * k + n ; //cout << " " << k << j << n << " , " << iii << endl ;
2059 fHistSub[iii].fHistSubHar[j].fHistPsiSubs->Fill(fPsiSub[n][k][j]) ;
2060 }
2061 }
2062 }
2063
2064 // full event Psi, PsiSubCorr, PsiSubCorrDiff, cos, mult, q
2065 for(int k = 0; k < Flow::nSels; k++)
2066 {
2067 for(int j = 0; j < Flow::nHars; j++)
2068 {
2069 float order = (float)(j+1);
2070 fHistFull[k].fHistFullHar[j].fHistPsi->Fill(fPsi[k][j]);
2071 if(k<2 && j<2)
2072 {
2073 if(k==0)
2074 {
2075 float psi1 = fPsi[0][j] ;
2076 float psi2 = fPsi[1][j] ;
2077 float diff = psi1 - psi2 ;
2078 if(diff < -TMath::Pi()/(j+1)) { diff += 2*TMath::Pi()/(j+1) ; }
2079 else if(diff > +TMath::Pi()/(j+1)) { diff -= 2*TMath::Pi()/(j+1) ; }
2080 fHistFull[k].fHistFullHar[j].fHistPsiDiff->Fill(diff) ; // k=0
2081 }
2082 else if(k==1)
2083 {
51a9b7d8 2084 float psi1 = 0. ; float psi2 = 0. ;
9777bfcb 2085 if (j==0) { psi1 = fPsi[0][0] ; psi2 = fPsi[1][1] ; }
2086 else if(j==1) { psi1 = fPsi[0][0] ; psi2 = fPsi[0][1] ; }
2087 float diff = psi1 - psi2 ;
2088 diff = (TMath::Abs(diff) > TMath::Pi()) ? ((diff > 0.) ? -(2*TMath::Pi()-diff) : -(diff+2*TMath::Pi())) : diff ;
2089 fHistFull[k].fHistFullHar[j].fHistPsiDiff->Fill(diff) ; // k=1
2090 }
2091 }
2092
2093 if(fPsiSub[0][k][j] != 0. && fPsiSub[1][k][j] != 0.)
2094 {
2095 float psiSubCorr; // this is: delta_Psi
2096 if(fV1Ep1Ep2 == kFALSE || order != 1)
2097 {
2098 psiSubCorr = fPsiSub[0][k][j] - fPsiSub[1][k][j];
2099 }
2100 else // i.e. (fV1Ep1Ep2 == kTRUE && order == 1)
2101 {
2102 psiSubCorr = fPsiSub[0][k][0] + fPsiSub[1][k][0] - 2*fPsi[k][1];
2103 }
2104 fHistFull[k].fHistCos->Fill(order, (float)cos(order * psiSubCorr)) ;
2105 if(psiSubCorr < 0.) { psiSubCorr += 2*TMath::Pi()/order ; }
2106 if(psiSubCorr > 2*TMath::Pi()/order) { psiSubCorr -= 2*TMath::Pi()/order ; } // for v1Ep1Ep2 which gives -2*TMath::Pi() < psiSubCorr < 2*2*TMath::Pi()
2107 fHistFull[k].fHistFullHar[j].fHistPsiSubCorr->Fill(psiSubCorr);
2108 }
2109
2110 if(j < Flow::nHars - 1) // subevents of different harmonics
2111 {
51a9b7d8 2112 int j1 = 0 ; int j2 = 0 ;
9777bfcb 2113 float psiSubCorrDiff;
2114 if(j==0) { j1 = 1, j2 = 2 ; }
2115 else if(j==1) { j1 = 1, j2 = 3 ; }
2116 else if(j==2) { j1 = 2, j2 = 4 ; }
2117 psiSubCorrDiff = fmod((double)fPsiSub[0][k][j1-1],2*TMath::Pi()/(double)j2)-fmod((double)fPsiSub[1][k][j2-1],2*TMath::Pi()/(double)j2) ;
2118 if(psiSubCorrDiff < 0.) { psiSubCorrDiff += 2*TMath::Pi()/(float)j2 ; }
2119 fHistFull[k].fHistFullHar[j].fHistPsiSubCorrDiff->Fill(psiSubCorrDiff) ;
2120 psiSubCorrDiff = fmod((double)fPsiSub[0][k][j2-1],2*TMath::Pi()/(double)j2)-fmod((double)fPsiSub[1][k][j1-1],2*TMath::Pi()/(double)j2) ;
2121 if(psiSubCorrDiff < 0.) { psiSubCorrDiff += 2*TMath::Pi()/(float)j2 ; }
2122 fHistFull[k].fHistFullHar[j].fHistPsiSubCorrDiff->Fill(psiSubCorrDiff) ;
2123 }
2124
2125 fHistFull[k].fHistFullHar[j].fHistMult->Fill((float)fMult[k][j]) ;
2126 fHistFull[k].fHistFullHar[j].fHistQnorm->Fill(fQnorm[k][j]) ;
2127 }
2128 }
2129
2130 return ;
2131}
2132//-----------------------------------------------------------------------
2133void AliFlowAnalyser::FillParticleHistograms(TObjArray* fFlowTracks)
2134{
2135 // fills tracks histograms
2136
2137 cout << " Tracks Loop . " << endl ;
2138
2139 float corrMultUnit = 0. ;
2140 float corrMultN = 0. ;
2141 float etaSymPosTpcN = 0. ;
2142 float etaSymNegTpcN = 0. ;
2143 float etaSymPosTpcNpart = 0. ;
2144 float etaSymNegTpcNpart = 0. ;
2145 float hPlusN = 0. ;
2146 float hMinusN = 0. ;
2147 float piPlusN = 0. ;
2148 float piMinusN = 0. ;
2149 float protonN = 0. ;
2150 float pbarN = 0. ;
2151 float kMinusN = 0. ;
2152 float kPlusN = 0. ;
2153 float deuteronN = 0. ;
2154 float dbarN = 0. ;
2155 float electronN = 0. ;
2156 float positronN = 0. ;
2157 float muonMinusN = 0. ;
2158 float muonPlusN = 0. ;
2159
2160 for(fTrackNumber=0;fTrackNumber<fNumberOfTracks;fTrackNumber++)
2161 {
2162 fFlowTrack = (AliFlowTrack*)fFlowTracks->At(fTrackNumber) ;
2163 //cout << "Track n. " << fTrackNumber << endl ; fFlowTrack->Dump() ;
2164
2165 bool constrainable = fFlowTrack->IsConstrainable() ;
2166 // int label = fFlowTrack->Label() ;
2167 float dcaGlobal = TMath::Abs(fFlowTrack->Dca()) ;
2168 float dcaSigned = fFlowTrack->TransDcaSigned() ;
2169 float dcaTrans = TMath::Abs(dcaSigned) ;
2170 float eta = fFlowTrack->Eta() ;
2171 float phi = fFlowTrack->Phi() ;
2172 float pt = fFlowTrack->Pt() ;
2173 float etaGlob = fFlowTrack->EtaGlobal() ;
2174 float phiGlob = fFlowTrack->PhiGlobal() ;
2175 float ptGlob = fFlowTrack->PtGlobal() ;
2176 float totalp = fFlowTrack->P() ;
2177 // float logp = TMath::Log10(totalp) ;
2178 // float zFirstPoint = fFlowTrack->ZFirstPoint() ;
2179 // float zLastPoint = fFlowTrack->ZLastPoint() ;
2180 float lenght = fFlowTrack->TrackLength() ;
2181 int charge = fFlowTrack->Charge() ;
2182 float chi2 = fFlowTrack->Chi2() ;
2183 int fitPtsTPC = (int)((float)fFlowTrack->FitPtsTPC()) ;
2184 int maxPtsTPC = fFlowTrack->MaxPtsTPC() ;
2185 float chi2TPC = fFlowTrack->Chi2TPC() ;
2186 int fitPtsITS = fFlowTrack->FitPtsITS() ;
2187 int maxPtsITS = fFlowTrack->MaxPtsITS() ;
2188 float chi2ITS = fFlowTrack->Chi2ITS() ;
2189 int fitPtsTRD = fFlowTrack->NhitsTRD() ;
2190 int maxPtsTRD = fFlowTrack->MaxPtsTRD() ;
2191 float chi2TRD = fFlowTrack->Chi2TRD() ;
2192 int fitPtsTOF = fFlowTrack->NhitsTOF() ;
2193 int maxPtsTOF = fFlowTrack->MaxPtsTOF() ;
2194 float chi2TOF = fFlowTrack->Chi2TOF() ;
2195 float dEdx = fFlowTrack->DedxTPC() ;
2196 float its = fFlowTrack->DedxITS() ;
2197 float trd = fFlowTrack->SigTRD() ;
2198 float tof = fFlowTrack->TofTOF() ;
2199 float lpTPC = 0 ; if(fFlowTrack->PatTPC()>0) { lpTPC = TMath::Log10(fFlowTrack->PatTPC()) ; }
2200 float lpITS = 0 ; if(fFlowTrack->PatITS()>0) { lpITS = TMath::Log10(fFlowTrack->PatITS()) ; }
2201 float lpTRD = 0 ; if(fFlowTrack->PatTRD()>0) { lpTRD = TMath::Log10(fFlowTrack->PatTRD()) ; }
2202 float lpTOF = 0 ; if(fFlowTrack->PatTOF()>0) { lpTOF = TMath::Log10(fFlowTrack->PatTOF()) ; }
2203 float invMass = fFlowTrack->InvMass() ;
2204 Char_t pid[10]="0" ; strcpy(pid,fFlowTrack->Pid()) ;
2205 fPidId = -1 ; // assigned later
2206
2207 // no selections: Charge, Dca, Chi2, FitPts, MaxPts, FitOverMax, PID
2208 fHistCharge->Fill((float)charge);
2209 fHistDcaGlobal->Fill(dcaGlobal);
2210 fHistDca->Fill(dcaTrans) ;
2211 fHistTransDca->Fill(dcaSigned);
2212 fHistChi2->Fill(chi2);
2213
2214 // - here ITS (chi2 & nHits)
2215 fHistChi2ITS->Fill(chi2ITS);
2216 if (fitPtsITS>0)
2217 fHistChi2normITS->Fill(chi2ITS/((float)fitPtsITS));
2218 fHistFitPtsITS->Fill((float)fitPtsITS);
2219 fHistMaxPtsITS->Fill((float)maxPtsITS);
2220
2221 // - here TPC (chi2 & nHits)
2222 fHistChi2TPC->Fill(chi2TPC);
2223 if (fitPtsTPC>0)
2224 fHistChi2normTPC->Fill(chi2TPC/((float)fitPtsTPC));
2225 fHistFitPtsTPC->Fill((float)fitPtsTPC);
2226 fHistMaxPtsTPC->Fill((float)maxPtsTPC);
2227 if(maxPtsTPC>0)
2228 fHistFitOverMaxTPC->Fill((float)(fitPtsTPC)/(float)maxPtsTPC);
2229
2230 // - here TRD (chi2 & nHits)
2231 fHistChi2TRD->Fill(chi2TRD);
2232 if (fitPtsTRD>0)
2233 fHistChi2normTRD->Fill(chi2TRD/((float)fitPtsTRD));
2234 fHistFitPtsTRD->Fill((float)fitPtsTRD);
2235 fHistMaxPtsTRD->Fill((float)maxPtsTRD);
2236
2237 // - here TOF (chi2 & nHits)
2238 fHistChi2TOF->Fill(chi2TOF);
2239 if (fitPtsTOF>0)
2240 fHistChi2normTOF->Fill(chi2TOF/((float)fitPtsTOF));
2241 fHistFitPtsTOF->Fill((float)fitPtsTOF);
2242 fHistMaxPtsTOF->Fill((float)maxPtsTOF);
2243
2244 // fit over max (all)
2245 int maxPts = maxPtsITS + maxPtsTPC + maxPtsTRD + maxPtsTOF ;
2246 int fitPts = fitPtsITS + fitPtsTPC + fitPtsTRD + fitPtsTOF ;
2247 if(maxPts>0)
2248 fHistFitOverMax->Fill((float)(fitPts)/(float)maxPts) ;
2249
2250 // lenght
2251 fHistLenght->Fill(lenght) ;
2252 fHistInvMass->Fill(invMass) ;
2253
2254 // PID histograms & multiplicity count (for bayesian histogram)
2255 if(charge == 1)
2256 {
2257 hPlusN++ ;
2258 fHistMeanDedxPos2D->Fill(lpTPC, dEdx) ;
2259 fHistMeanDedxPos2DITS->Fill(lpITS, its) ;
2260 fHistMeanDedxPos2DTRD->Fill(lpTRD, trd) ;
2261 fHistMeanDedxPos2DTOF->Fill(lpTOF, tof) ;
2262 //
2263 float positron = fFlowTrack->ElectronPositronProb() ;
2264 fHistPidPositron->Fill(positron) ;
2265 if(strcmp(pid, "e+") == 0)
2266 {
2267 fPidId = 0 ; positronN++ ; fHistPidPt->Fill(1.5,pt) ;
2268 fHistMeanTPCPositron->Fill(lpTPC, dEdx) ;
2269 fHistMeanITSPositron->Fill(lpITS, its);
2270 fHistMeanTRDPositron->Fill(lpTRD, trd);
2271 fHistMeanTOFPositron->Fill(invMass, tof);
2272 fHistPidPositronPart->Fill(positron) ;
2273 }
2274 float muonPlus = fFlowTrack->MuonPlusMinusProb() ;
2275 fHistPidMuonPlus->Fill(muonPlus) ;
2276 if(strcmp(pid, "mu+") == 0)
2277 {
2278 fPidId = 1 ; muonPlusN++ ; fHistPidPt->Fill(3.5,pt) ;
2279 fHistMeanTPCMuonPlus->Fill(lpTPC, dEdx) ;
2280 fHistMeanITSMuonPlus->Fill(lpITS, its);
2281 fHistMeanTRDMuonPlus->Fill(lpTRD, trd);
2282 fHistMeanTOFMuonPlus->Fill(invMass, tof);
2283 fHistPidMuonPlusPart->Fill(muonPlus) ;
2284 }
2285 float piPlus = fFlowTrack->PionPlusMinusProb() ;
2286 fHistPidPiPlus->Fill(piPlus) ;
2287 if(strcmp(pid, "pi+") == 0)
2288 {
2289 fPidId = 2 ; piPlusN++ ; fHistPidPt->Fill(5.5,pt) ;
2290 fHistMeanTPCPiPlus->Fill(lpTPC, dEdx) ;
2291 fHistMeanITSPiPlus->Fill(lpITS, its);
2292 fHistMeanTRDPiPlus->Fill(lpTRD, trd);
2293 fHistMeanTOFPiPlus->Fill(invMass, tof);
2294 fHistPidPiPlusPart->Fill(piPlus) ;
2295 }
2296 float kplus = fFlowTrack->KaonPlusMinusProb() ;
2297 fHistPidKplus->Fill(kplus) ;
2298 if(strcmp(pid, "k+") == 0)
2299 {
2300 fPidId = 3 ; kPlusN++ ; fHistPidPt->Fill(7.5,pt) ;
2301 fHistMeanTPCKplus->Fill(lpTPC, dEdx) ;
2302 fHistMeanITSKplus->Fill(lpITS, its);
2303 fHistMeanTRDKplus->Fill(lpTRD, trd);
2304 fHistMeanTOFKplus->Fill(invMass, tof);
2305 fHistPidKplusPart->Fill(kplus) ;
2306 }
2307 float proton = fFlowTrack->ProtonPbarProb() ;
2308 fHistPidProton->Fill(proton) ;
2309 if(strcmp(pid, "pr+") == 0)
2310 {
2311 fPidId = 4 ; protonN++ ; fHistPidPt->Fill(9.5,pt) ;
2312 fHistMeanTPCProton->Fill(lpTPC, dEdx) ;
2313 fHistMeanITSProton->Fill(lpITS, its);
2314 fHistMeanTRDProton->Fill(lpTRD, trd);
2315 fHistMeanTOFProton->Fill(invMass, tof);
2316 fHistPidProtonPart->Fill(proton) ;
2317 }
2318 float deuteron = fFlowTrack->DeuteriumAntiDeuteriumProb() ;
2319 fHistPidDeuteron->Fill(deuteron) ;
2320 if(strcmp(pid, "d+") == 0)
2321 {
2322 fPidId = 6 ; deuteronN++ ; fHistPidPt->Fill(11.5,pt) ;
2323 fHistMeanTPCDeuteron->Fill(lpTPC, dEdx) ;
2324 fHistMeanITSDeuteron->Fill(lpITS, its);
2325 fHistMeanTRDDeuteron->Fill(lpTRD, trd);
2326 fHistMeanTOFDeuteron->Fill(invMass, tof);
2327 fHistPidDeuteronPart->Fill(deuteron) ;
2328 }
2329 }
2330 else if(charge == -1)
2331 {
2332 hMinusN++ ;
2333 fHistMeanDedxNeg2D->Fill(lpTPC, dEdx) ;
2334 fHistMeanDedxNeg2DITS->Fill(lpITS, its) ;
2335 fHistMeanDedxNeg2DTRD->Fill(lpTRD, trd) ;
2336 fHistMeanDedxNeg2DTOF->Fill(lpTOF, tof) ;
2337 //
2338 float electron = fFlowTrack->ElectronPositronProb() ;
2339 fHistPidElectron->Fill(electron);
2340 if(strcmp(pid, "e-") == 0)
2341 {
2342 fPidId = 0 ; electronN++ ; fHistPidPt->Fill(0.5,pt) ;
2343 fHistMeanTPCElectron->Fill(lpTPC, dEdx);
2344 fHistMeanITSElectron->Fill(lpITS, its);
2345 fHistMeanTRDElectron->Fill(lpTRD, trd);
2346 fHistMeanTOFElectron->Fill(invMass, tof);
2347 fHistPidElectronPart->Fill(electron);
2348 }
2349 float muonMinus = fFlowTrack->MuonPlusMinusProb() ;
2350 fHistPidMuonMinus->Fill(muonMinus) ;
2351 if(strcmp(pid, "mu-") == 0)
2352 {
2353 fPidId = 1 ; muonMinusN++ ; fHistPidPt->Fill(2.5,pt) ;
2354 fHistMeanTPCMuonMinus->Fill(lpTPC, dEdx) ;
2355 fHistMeanITSMuonMinus->Fill(lpITS, its);
2356 fHistMeanTRDMuonMinus->Fill(lpTRD, trd);
2357 fHistMeanTOFMuonMinus->Fill(invMass, tof);
2358 fHistPidMuonMinusPart->Fill(muonMinus) ;
2359 }
2360 float piMinus = fFlowTrack->PionPlusMinusProb() ;
2361 fHistPidPiMinus->Fill(piMinus) ;
2362 if(strcmp(pid, "pi-") == 0)
2363 {
2364 fPidId = 2 ; piMinusN++ ; fHistPidPt->Fill(4.5,pt) ;
2365 fHistMeanTPCPiMinus->Fill(lpTPC, dEdx);
2366 fHistMeanITSPiMinus->Fill(lpITS, its);
2367 fHistMeanTRDPiMinus->Fill(lpTRD, trd);
2368 fHistMeanTOFPiMinus->Fill(invMass, tof);
2369 fHistPidPiMinusPart->Fill(piMinus);
2370 }
2371 float kminus = fFlowTrack->KaonPlusMinusProb() ;
2372 fHistPidKminus->Fill(kminus);
2373 if(strcmp(pid, "k-") == 0)
2374 {
2375 fPidId = 3 ; kMinusN++ ; fHistPidPt->Fill(6.5,pt) ;
2376 fHistMeanTPCKminus->Fill(lpTPC, dEdx);
2377 fHistMeanITSKminus->Fill(lpITS, its);
2378 fHistMeanTRDKminus->Fill(lpTRD, trd);
2379 fHistMeanTOFKminus->Fill(invMass, tof);
2380 fHistPidKminusPart->Fill(kminus);
2381 }
2382 float antiproton = fFlowTrack->ProtonPbarProb() ;
2383 fHistPidAntiProton->Fill(antiproton);
2384 if(strcmp(pid, "pr-") == 0)
2385 {
2386 fPidId = 4 ; pbarN++ ; fHistPidPt->Fill(8.5,pt) ;
2387 fHistMeanTPCPbar->Fill(lpTPC, dEdx);
2388 fHistMeanITSPbar->Fill(lpITS, its);
2389 fHistMeanTRDPbar->Fill(lpTRD, trd);
2390 fHistMeanTOFPbar->Fill(invMass, tof);
2391 fHistPidAntiProtonPart->Fill(antiproton);
2392 }
2393 float antideuteron = fFlowTrack->DeuteriumAntiDeuteriumProb() ;
2394 fHistPidAntiDeuteron->Fill(antideuteron);
2395 if(strcmp(pid, "d-") == 0)
2396 {
2397 fPidId = 6 ; dbarN++ ; fHistPidPt->Fill(10.5,pt) ;
2398 fHistMeanTPCAntiDeuteron->Fill(lpTPC, dEdx);
2399 fHistMeanITSAntiDeuteron->Fill(lpITS, its);
2400 fHistMeanTRDAntiDeuteron->Fill(lpTRD, trd);
2401 fHistMeanTOFAntiDeuteron->Fill(invMass, tof);
2402 fHistPidAntiDeuteronPart->Fill(antideuteron);
2403 }
2404 }
2405
2406 // Yield3D, Yield2D, Eta, Pt, Phi, bayP.Id.
2407 fHistPtot->Fill(totalp) ;
2408 fHistPt->Fill(pt) ;
2409 fHistPhi->Fill(phi);
2410 fHistAllEtaPtPhi3D->Fill(eta, pt, phi) ;
2411 fHistYieldAll2D->Fill(eta, pt) ;
2412 fHistBayPidMult->Fill(fPidId) ;
2413 if(constrainable)
2414 {
2415 fHistPhiCons->Fill(phi);
2416 fHistPhiPtCon->Fill(phi, pt);
2417 fHistYieldCon2D->Fill(eta, pt) ;
2418 fHistConsEtaPtPhi3D->Fill(eta, pt, phi) ;
2419 fHistGlobEtaPtPhi3D->Fill(etaGlob, ptGlob, phiGlob) ;
2420 }
2421 else
2422 {
2423 fHistYieldUnc2D->Fill(etaGlob, ptGlob) ;
2424 fHistUncEtaPtPhi3D->Fill(etaGlob, ptGlob, phiGlob) ;
2425 fHistPhiPtUnc->Fill(phiGlob, ptGlob) ;
2426 }
2427 if(fFlowTrack->Charge()>0) { fHistPtPhiPos->Fill(phi, pt); }
2428 else if(fFlowTrack->Charge()<0) { fHistPtPhiNeg->Fill(phi, pt); }
2429
2430 // fills selected part histograms
2431 if(fFlowSelect->SelectPart(fFlowTrack))
2432 {
2433 if(strlen(fFlowSelect->PidPart()) != 0)
2434 {
2435 float rapidity = fFlowTrack->Y();
2436 fHistBinEta->Fill(rapidity, rapidity);
2437 fHistYieldPart2D->Fill(rapidity, pt);
2438 }
2439 else
2440 {
2441 fHistBinEta->Fill(eta, eta) ;
2442 fHistYieldPart2D->Fill(eta, pt) ;
2443 }
2444 fHistBayPidMultPart->Fill(fPidId) ;
2445 fHistBinPt->Fill(pt, pt) ;
2446 fHistEtaPtPhi3DPart->Fill(eta,pt,phi) ;
2447 fHistDcaGlobalPart->Fill(dcaGlobal) ;
2448 fHistInvMassPart->Fill(invMass) ;
2449 if(charge == 1)
2450 {
2451 fHistMeanDedxPos3DPart->Fill(lpITS, dEdx, fPidId) ;
2452 fHistMeanDedxPos3DPartITS->Fill(lpITS, its, fPidId) ;
2453 }
2454 else if(charge == -1)
2455 {
2456 fHistMeanDedxNeg3DPart->Fill(lpITS, dEdx, fPidId) ;
2457 fHistMeanDedxNeg3DPartITS->Fill(lpITS, its, fPidId) ;
2458 }
2459 if(eta > 0.) { etaSymPosTpcNpart++ ; } // for fHistEtaSymPart
2460 else { etaSymNegTpcNpart++ ; }
2461 }
2462 else
2463 {
2464 fHistEtaPtPhi3DOut->Fill(eta,pt,phi) ;
2465 fHistYieldOut2D->Fill(eta, pt) ;
2466 fHistDcaGlobalOut->Fill(dcaGlobal) ;
2467 fHistInvMassOut->Fill(invMass) ;
2468 }
2469
2470 // cos(n*phiLab)
2471 for(int j = 0; j < Flow::nHars; j++)
2472 {
2473 bool oddHar = (j+1) % 2 ;
2474 float order = (float)(j+1) ;
2475 float vIn = 100 * cos((double)order * phi) ;
2476 if(eta < 0 && oddHar) { vIn *= -1 ; }
2477 fHistCosPhi->Fill(order, vIn);
2478 }
2479
2480 //For Eta symmetry TPC
2481 if(fFlowTrack->FitPtsTPC())
2482 {
2483 if(eta > 0.) { etaSymPosTpcN++ ; } // for fHistEtaSym
2484 else { etaSymNegTpcN++ ; }
2485 }
2486
2487 // not to call it twice ...
2488 int nEtaS = HarmonicsLoop(fFlowTrack) ;
2489
2490 //For Correlated Multiplicity
2491 corrMultN += (float)nEtaS ;
2492
2493 //For Correlated Multiplicity in 1 unit rapidity
2494 if(TMath::Abs(eta) <= 0.5) { corrMultUnit += (float)nEtaS ; }
2495
2496 //delete pointer
2497 fFlowTrack = 0 ;
2498 } // end of tracks loop
2499
2500 // EtaSym
2501 float etaSymTpc = 0 ;
2502 if(etaSymPosTpcN || etaSymNegTpcN) { etaSymTpc = (etaSymPosTpcN - etaSymNegTpcN) / (etaSymPosTpcN + etaSymNegTpcN); }
2503 float etaSymTpcPart = 0 ;
2504 if(etaSymPosTpcNpart || etaSymNegTpcNpart) { etaSymTpcPart = (etaSymPosTpcNpart - etaSymNegTpcNpart) / (etaSymPosTpcNpart + etaSymNegTpcNpart) ; }
2505 Float_t vertexZ = fVertex[2] ;
2506 // all
2507 fHistEtaSym->Fill(etaSymTpc);
2508 fHistEtaSymVerZ2D->Fill(vertexZ , etaSymTpc);
2509 // selected
2510 fHistEtaSymPart->Fill(etaSymTpc);
2511 fHistEtaSymVerZ2DPart->Fill(vertexZ , etaSymTpcPart);
2512
2513 // PID multiplicities
2514 float totalMult = (float)fFlowTracks->GetEntries() ;
2515 fHistPidMult->Fill(1., totalMult);
2516 fHistPidMult->Fill(2., hPlusN);
2517 fHistPidMult->Fill(3., hMinusN);
2518 fHistPidMult->Fill(4., piPlusN);
2519 fHistPidMult->Fill(5., piMinusN);
2520 fHistPidMult->Fill(6., protonN);
2521 fHistPidMult->Fill(7., pbarN);
2522 fHistPidMult->Fill(8., kPlusN);
2523 fHistPidMult->Fill(9., kMinusN);
2524 fHistPidMult->Fill(10., deuteronN);
2525 fHistPidMult->Fill(11., dbarN);
2526 fHistPidMult->Fill(12., electronN);
2527 fHistPidMult->Fill(13., positronN);
2528 fHistPidMult->Fill(14., muonMinusN);
2529 fHistPidMult->Fill(15., muonPlusN);
2530
2531 // Multiplicity of particles correlated with the event planes
2532 corrMultN /= (float)(Flow::nHars * Flow::nSels) ;
2533 fHistMultPart->Fill(corrMultN) ;
2534 // ...in one unit rapidity
2535 corrMultUnit /= (float)(Flow::nHars * Flow::nSels) ;
2536 fHistMultPartUnit->Fill(corrMultUnit) ;
2537
2538 return ;
2539}
2540//-----------------------------------------------------------------------
2541Int_t AliFlowAnalyser::HarmonicsLoop(AliFlowTrack* fFlowTrack)
2542{
2543 // HarmonicsLoop, does the correlation analysis, fills histograms
2544 // (internally called by ::FillParticleHistograms(..) loop) .
2545
2546 float eta = fFlowTrack->Eta() ;
2547 float phi = fFlowTrack->Phi() ;
2548 float pt = fFlowTrack->Pt() ;
2549 float zFirstPoint = fFlowTrack->ZFirstPoint() ;
2550 // float zLastPoint = fFlowTrack->ZLastPoint() ;
2551
2552 // Looping over Selections and Harmonics
2553 int corrMultN = 0 ;
2554 for (int k = 0; k < Flow::nSels; k++)
2555 {
2556 fFlowSelect->SetSelection(k) ;
2557 for (int j = 0; j < Flow::nHars; j++)
2558 {
2559 bool oddHar = (j+1) % 2;
2560 fFlowSelect->SetHarmonic(j);
2561 double order = (double)(j+1);
51a9b7d8 2562 float psi_i = 0. ; float psi_2 = 0. ;
9777bfcb 2563 if(fFlowEvent->EtaSubs()) // particles with the opposite subevent
2564 {
2565 if(eta > 0) { psi_i = fPsiSub[1][k][j] ; } //check
2566 else { psi_i = fPsiSub[0][k][j] ; }
2567 }
2568 else if(order > 3. && !oddHar)
2569 {
2570 psi_i = fPsi[k][1]; // 2nd harmomic event plane
2571 if(psi_i > 2*TMath::Pi()/order) { psi_i -= 2*TMath::Pi()/order ; }
2572 if(psi_i > 2*TMath::Pi()/order) { psi_i -= 2*TMath::Pi()/order ; }
2573 }
2574 else // random subevents
2575 {
2576 psi_i = fPsi[k][j] ;
2577 }
2578
2579 if(fFlowSelect->Select(fFlowTrack)) // Get detID
2580 {
2581 Bool_t kTpcPlus = kFALSE ;
2582 Bool_t kTpcMinus = kFALSE ;
2583 Bool_t kTpcAll = kFALSE ;
2584
2585 fHistFull[k].fHistFullHar[j].fHistYieldPt->Fill(pt);
2586 fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3D->Fill(eta, pt, phi);
2587 fHistFull[k].fHistFullHar[j].fHistYield2D->Fill(eta, pt);
2588 fHistFull[k].fHistFullHar[j].fHistDcaGlob->Fill(TMath::Abs(fFlowTrack->Dca()));
2589
2590 // Set Tpc (+ and -)
2591 if(fFlowTrack->FitPtsTPC()) //OR*: AliESDtrack:: "TBits& GetTPCClusterMap()" or "Int_t GetTPCclusters(Int_t* idx)" ...
2592 {
2593 if(zFirstPoint >= 0. && eta > 0.) { kTpcPlus = kTRUE ; }
2594 else if(zFirstPoint <= 0. && eta < 0.) { kTpcMinus = kTRUE ; }
2595 else { kTpcAll = kTRUE ; }
2596 }
2597
2598 // PID Multiplicities (particle for R.P.) - done just one time for each selection
2599 if(j==0) { fHistFull[k].fHistBayPidMult->Fill(fPidId) ; }
2600
2601 // Calculate weights for filling histograms
2602 float wt = 1. ; // TMath::Abs(fFlowEvent->Weight(k, j, fFlowTrack)) ;
2603
2604 // Fill histograms with selections
2605 if(kTpcPlus) { fHistFull[k].fHistFullHar[j].fHistPhiPlus->Fill(phi,wt) ; }
2606 else if(kTpcMinus) { fHistFull[k].fHistFullHar[j].fHistPhiMinus->Fill(phi,wt) ; }
2607 else if(kTpcAll) { fHistFull[k].fHistFullHar[j].fHistPhiAll->Fill(phi,wt) ; }
2608 fHistFull[k].fHistFullHar[j].fHistPhi->Fill(phi,wt) ;
2609
2610 // Get phiWgt from file
2611 double phiWgt;
2612 if(order > 3. && !oddHar) { phiWgt = fFlowEvent->PhiWeight(k, 1, fFlowTrack) ; }
2613 else { phiWgt = fFlowEvent->PhiWeight(k, j, fFlowTrack) ; }
2614 if(oddHar && eta<0.) { phiWgt /= -1. ; } // only for flat hists
2615 // cout << " Weight [" << k << "][" << j << "] (" << fFlowTrack->GetName() << ") = " << phiWgt << " or " << wt << " for phi = " << phi << endl ; // chk
2616
2617 // Fill Flat histograms
2618 if(kTpcPlus) { fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus->Fill(phi, phiWgt) ; }
2619 else if(kTpcMinus) { fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus->Fill(phi, phiWgt) ; }
2620 else if(kTpcAll) { fHistFull[k].fHistFullHar[j].fHistPhiFlatAll->Fill(phi, phiWgt) ; }
2621 fHistFull[k].fHistFullHar[j].fHistPhiFlat->Fill(phi,phiWgt) ;
2622
2623 if(oddHar && eta<0.) { phiWgt *= -1. ; } // restore value
2624
2625 // Remove autocorrelations
2626 TVector2 Q_i;
2627 if(!fFlowEvent->EtaSubs()) // random subevents
2628 {
2629 if(order > 3. && !oddHar) // 2nd harmonic event plane
2630 {
2631 Q_i.Set(phiWgt * cos(phi * 2), phiWgt * sin(phi * 2));
2632 TVector2 mQ_i = fQ[k][1] - Q_i;
2633 psi_i = mQ_i.Phi() / 2;
2634 if(psi_i < 0.) { psi_i += TMath::Pi() ; }
2635 }
2636 else
2637 {
2638 Q_i.Set(phiWgt * cos(phi * order), phiWgt * sin(phi * order));
2639 TVector2 mQ_i = fQ[k][j] - Q_i;
2640 psi_i = mQ_i.Phi() / order;
2641 if(psi_i < 0.) { psi_i += 2*TMath::Pi()/order ; }
2642 }
2643 }
2644
2645 // Remove autocorrelations of the second order 'particles' which are used for v1{EP1,EP2}.
2646 if (fV1Ep1Ep2 == kTRUE && order == 1)
2647 {
2648 AliFlowSelection usedForPsi2 = *fFlowSelect ;
2649 usedForPsi2.SetHarmonic(1);
2650 if(usedForPsi2.Select(fFlowTrack)) // particle was used for Psi2
2651 {
2652 Q_i.Set(phiWgt * cos(phi * 2), phiWgt * sin(phi * 2));
2653 TVector2 mQ_i = fQ[k][1] - Q_i;
2654 psi_2 = mQ_i.Phi() / 2;
2655 if(psi_2 < 0.) { psi_2 += TMath::Pi() ; }
2656 }
2657 else // particle was not used for Psi2
2658 {
2659 psi_2 = fPsi[k][1];
2660 }
2661 }
2662 }
2663 else
2664 {
2665 fHistFull[k].fHistFullHar[j].fHistYieldPtout->Fill(pt);
2666 fHistFull[k].fHistFullHar[j].fHistEtaPtPhi3Dout->Fill(eta, pt, phi);
2667 fHistFull[k].fHistFullHar[j].fHistYield2Dout->Fill(eta, pt);
2668 fHistFull[k].fHistFullHar[j].fHistDcaGlobout->Fill(TMath::Abs(fFlowTrack->Dca()));
2669 }
2670
2671 // Caculate v for all particles selected for correlation analysis
2672 if(fFlowSelect->SelectPart(fFlowTrack))
2673 {
2674 corrMultN++;
2675
2676 float v;
2677 if (fV1Ep1Ep2 == kFALSE || order != 1)
2678 {
2679 v = 100 * cos(order * (phi - psi_i)) ;
2680 }
2681 else // i.e. (fV1Ep1Ep2 == kTRUE && order == 1)
2682 {
2683 v = 100 * cos(phi + psi_i - 2*psi_2) ;
2684 }
2685
2686 float vFlip = v;
2687 if(eta < 0 && oddHar) { vFlip *= -1 ; }
2688 if(strlen(fFlowSelect->PidPart()) != 0) // pid, fill rapidity
2689 {
2690 float rapidity = fFlowTrack->Y();
2691 fHistFull[k].fHistFullHar[j].fHistvObs2D->Fill(rapidity, pt, v);
2692
2693 if(fPtRangevEta[1] > fPtRangevEta[0]) // cut is used
2694 {
2695 if(pt < fPtRangevEta[1] && pt >= fPtRangevEta[0]) // check cut range, fill if in range
2696 {
2697 fHistFull[k].fHistFullHar[j].fHistvObsEta->Fill(rapidity, v);
2698 }
2699 }
2700 else // cut is not used, fill in any case
2701 {
2702 fHistFull[k].fHistFullHar[j].fHistvObsEta->Fill(rapidity, v);
2703 }
2704 }
2705 else // no pid, fill eta
2706 {
2707 fHistFull[k].fHistFullHar[j].fHistvObs2D->Fill(eta, pt, v);
2708
2709 if(fPtRangevEta[1] > fPtRangevEta[0]) // cut is used
2710 {
2711 if(pt < fPtRangevEta[1] && pt >= fPtRangevEta[0]) // check cut range, fill if in range
2712 {
2713 fHistFull[k].fHistFullHar[j].fHistvObsEta->Fill(eta, v);
2714 }
2715 }
2716 else // cut is not used, fill in any case
2717 {
2718 fHistFull[k].fHistFullHar[j].fHistvObsEta->Fill(eta, v);
2719 }
2720 }
2721
2722 if(fEtaRangevPt[1] > fEtaRangevPt[0]) // cut is used
2723 {
2724 if(TMath::Abs(eta) < fEtaRangevPt[1] && TMath::Abs(eta) >= fEtaRangevPt[0]) // check cut range, fill if in range
2725 {
2726 fHistFull[k].fHistFullHar[j].fHistvObsPt->Fill(pt, vFlip); // for odd harmonis /-/
2727 }
2728 }
2729 else // cut is not used, fill in any case
2730 {
2731 fHistFull[k].fHistFullHar[j].fHistvObsPt->Fill(pt, vFlip);
2732 }
2733
2734 // v_
2735 Bool_t etaPtNoCut = kTRUE;
2736 if(fPtRangevEta[1] > fPtRangevEta[0] && (pt < fPtRangevEta[0] || pt >= fPtRangevEta[1]))
2737 {
2738 etaPtNoCut = kFALSE;
2739 }
2740 if(fEtaRangevPt[1] > fEtaRangevPt[0] && (TMath::Abs(eta) < fEtaRangevPt[0] || TMath::Abs(eta) >= fEtaRangevPt[1]))
2741 {
2742 etaPtNoCut = kFALSE;
2743 }
2744 if(etaPtNoCut) { fHistFull[k].fHistvObs->Fill(order, vFlip) ; }
2745
2746 // Correlation of Phi of selected particles with Psi
2747 float phi_i = phi;
2748 if(eta < 0 && oddHar)
2749 {
2750 phi_i += TMath::Pi() ; // backward particle and odd harmonic
2751 if(phi_i > 2*TMath::Pi()) { phi_i -= 2*TMath::Pi() ; }
2752 }
2753 float dPhi = phi_i - psi_i;
2754 if(dPhi < 0.) { dPhi += 2*TMath::Pi() ; }
2755 fHistFull[k].fHistFullHar[j].fHistPhiCorr->Fill(fmod((double)dPhi, 2*TMath::Pi() / order));
2756 }
2757 }
2758 }
2759
2760 return corrMultN ;
2761}
2762//-----------------------------------------------------------------------
2763void AliFlowAnalyser::FillV0Histograms(TObjArray* fFlowV0s)
2764{
2765 // v0s histograms
2766
2767 cout << " V0s Loop . " << endl ;
2768
2769 int corrMultV0 = 0 ;
2770 for(fV0Number=0;fV0Number<fNumberOfV0s;fV0Number++)
2771 {
2772 fFlowV0 = (AliFlowV0*)fFlowV0s->At(fV0Number) ;
2773
2774 float mass = fFlowV0->Mass() ;
2775 float eta = fFlowV0->Eta() ;
2776 float rapidity = fFlowV0->Y() ;
2777 float phi = fFlowV0->Phi() ;
2778 float pt = fFlowV0->Pt() ;
2779 float totalp = fFlowV0->P() ;
2780 //int charge = fFlowV0->Charge() ;
2781 float dca = fFlowV0->Dca() ;
2782 float lenght = fFlowV0->V0Lenght() ;
2783 float sigma = fFlowV0->Sigma() ;
2784 float chi2 = fFlowV0->Chi2() ;
2785 Char_t pid[10] ; strcpy(pid, fFlowV0->Pid()) ;
2786 int labelPlus = -1 ;
2787 int labelMinus = -1 ;
2788 AliFlowTrack* daughterPlus = fFlowV0->DaughterP() ;
2789 AliFlowTrack* daughterMinus = fFlowV0->DaughterN() ;
2790 if(daughterPlus) { labelPlus = daughterPlus->Label() ; }
2791 if(daughterMinus) { labelMinus = daughterMinus->Label() ; }
2792
2793 fHistV0Mass->Fill(mass) ;
2794 fHistV0EtaPtPhi3D->Fill(eta, pt, phi) ;
2795 fHistV0YieldAll2D->Fill(eta, pt) ;
2796 fHistV0Dca->Fill(dca);
2797 fHistV0Chi2->Fill(chi2);
2798 fHistV0Lenght->Fill(lenght);
2799 fHistV0Sigma->Fill(sigma);
2800 fHistV0MassPtSlices->Fill(mass,pt);
2801 fHistV0PYall2D->Fill(rapidity, totalp) ;
2802
2803 if(fFlowSelect->SelectPart(fFlowV0))
2804 {
2805 bool inWin = fFlowSelect->SelectV0Part(fFlowV0) ;
2806 bool sx = fFlowSelect->SelectV0sxSide(fFlowV0) ;
2807 bool dx = fFlowSelect->SelectV0dxSide(fFlowV0) ;
2808 corrMultV0++ ;
2809
2810 fHistV0YieldPart2D->Fill(eta, pt) ;
2811 if(inWin)
2812 {
2813 fHistV0EtaPtPhi3DPart->Fill(eta, pt, phi) ;
2814 fHistV0LenghtPart->Fill(lenght);
2815 fHistV0DcaPart->Fill(dca);
2816 fHistV0MassWin->Fill(mass) ;
2817 fHistV0BinEta->Fill(eta, eta);
2818 fHistV0BinPt->Fill(pt, pt);
2819 }
2820 else
2821 {
2822 fHistV0sbEtaPtPhi3DPart->Fill(eta, pt, phi) ;
2823 fHistV0sbLenghtPart->Fill(lenght);
2824 fHistV0sbDcaPart->Fill(dca);
2825 fHistV0sbMassSide->Fill(mass) ;
2826 fHistV0sbBinEta->Fill(eta, eta);
2827 fHistV0sbBinPt->Fill(pt, pt);
2828 }
2829
2830 for(int k = 0; k < Flow::nSels; k++) // sort of HarmonicsLoop - selection number used
2831 {
2832 fFlowSelect->SetSelection(k) ;
2833 for(Int_t j=0;j<Flow::nHars;j++)
2834 {
2835 Bool_t oddHar = (j+1) % 2 ;
2836 Float_t order = (Float_t)(j+1) ;
2837 fFlowSelect->SetHarmonic(j);
2838
2839 // Remove autocorrelations
2840 Float_t psi_i, psi_2 ;
2841 TVector2 Q_1, Q_2 ;
2842 Float_t phiDaughter1 = 0. ; Float_t phiDaughter2 = 0. ;
2843 Double_t phiWgt1 = 0. ; Double_t phiWgt2 = 0. ;
2844 // -
2845 if(daughterPlus)
2846 { // 1
2847 fFlowTrack = daughterPlus ;
2848 phiDaughter1 = fFlowTrack->Phi() ;
2849 if(fFlowSelect->Select(fFlowTrack)) // Get phiWgt from file
2850 {
2851 if(order > 3. && !oddHar) { phiWgt1 = fFlowEvent->PhiWeight(k, 1, fFlowTrack) ; }
2852 else { phiWgt1 = fFlowEvent->PhiWeight(k, j, fFlowTrack) ; }
2853 }
2854 }
2855 if(daughterMinus)
2856 { // 2
2857 fFlowTrack = daughterMinus ;
2858 phiDaughter2 = fFlowTrack->Phi() ;
2859 if(fFlowSelect->Select(fFlowTrack)) // Get phiWgt from file
2860 {
2861 if(order > 3. && !oddHar) { phiWgt2 = fFlowEvent->PhiWeight(k, 1, fFlowTrack) ; }
2862 else { phiWgt2 = fFlowEvent->PhiWeight(k, j, fFlowTrack) ; }
2863 }
2864 }
2865
2866 // psi_2
2867 Q_1.Set(phiWgt1 * cos(phiDaughter1 * 2), phiWgt1 * sin(phiDaughter1 * 2));
2868 Q_2.Set(phiWgt2 * cos(phiDaughter2 * 2), phiWgt2 * sin(phiDaughter2 * 2));
2869 TVector2 mQ_i = fQ[k][1] ; mQ_i -= Q_1 ; mQ_i -= Q_2 ;
2870 psi_2 = mQ_i.Phi() / 2 ;
2871 if(psi_2 < 0.) { psi_2 += TMath::Pi() ; }
2872 if(psi_2 > TMath::Pi()) { psi_2 -= TMath::Pi() ; }
2873
2874 // psi_i
2875 if(order > 3. && !oddHar) { psi_i = psi_2 ; }
2876 else
2877 {
2878 Q_1.Set(phiWgt1 * cos(phiDaughter1 * order), phiWgt1 * sin(phiDaughter1 * order));
2879 Q_2.Set(phiWgt2 * cos(phiDaughter2 * order), phiWgt2 * sin(phiDaughter2 * order));
2880 TVector2 mQ_i = fQ[k][j] ; mQ_i -= Q_1 ; mQ_i -= Q_2 ;
2881 psi_i = mQ_i.Phi()/order ;
2882 if(psi_i < 0.) { psi_i += 2*TMath::Pi()/order ; }
2883 if(psi_i > 2*TMath::Pi()/order) { psi_i -= 2*TMath::Pi()/order ; }
2884 }
2885
2886 // Caculate v for all V0s selected for correlation analysis
2887 float v ;
2888 if(fV1Ep1Ep2 == kFALSE || order != 1) { v = 100 * cos(order * (phi - psi_i)) ; }
2889 else { v = 100 * cos(phi + psi_i - 2*psi_2) ; }
2890 float vFlip = v ; if(eta < 0 && oddHar) { vFlip *= -1 ; }
2891
2892 // invariant mass windows & sidebands
2893 if(inWin) { fHistFull[k].fHistFullHar[j].fHistV0vObs2D->Fill(eta, pt, v) ; }
2894 else { fHistFull[k].fHistFullHar[j].fHistV0sbvObs2D->Fill(eta, pt, v) ; }
2895
2896 if(fPtRangevEta[1] > fPtRangevEta[0]) // cut is used
2897 {
2898 if(pt < fPtRangevEta[1] && pt >= fPtRangevEta[0]) // check cut range, fill if in range
2899 {
2900 if(inWin) { fHistFull[k].fHistFullHar[j].fHistV0vObsEta->Fill(eta, v) ; }
2901 else
2902 {
2903 fHistFull[k].fHistFullHar[j].fHistV0sbvObsEta->Fill(eta, v) ;
2904 if(sx) { fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaSx->Fill(eta, v) ; }
2905 else if(dx) { fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaDx->Fill(eta, v) ; }
2906 }
2907 }
2908 }
2909 else // cut is not used, fill in any case
2910 {
2911 if(inWin) { fHistFull[k].fHistFullHar[j].fHistV0vObsEta->Fill(eta, v) ; }
2912 else
2913 {
2914 fHistFull[k].fHistFullHar[j].fHistV0sbvObsEta->Fill(eta, v) ;
2915 if(sx) { fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaSx->Fill(eta, v) ; }
2916 else if(dx) { fHistFull[k].fHistFullHar[j].fHistV0sbvObsEtaDx->Fill(eta, v) ; }
2917 }
2918 }
2919 if(fEtaRangevPt[1] > fEtaRangevPt[0]) // cut is used
2920 {
2921 if(TMath::Abs(eta) < fEtaRangevPt[1] && TMath::Abs(eta) >= fEtaRangevPt[0]) // check cut range, fill if in range
2922 {
2923 if(inWin) { fHistFull[k].fHistFullHar[j].fHistV0vObsPt->Fill(pt, vFlip) ; } // for odd harmonis /-/
2924 else
2925 {
2926 fHistFull[k].fHistFullHar[j].fHistV0sbvObsPt->Fill(pt, vFlip) ;
2927 if(sx) { fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtSx->Fill(eta, v) ; }
2928 else if(dx) { fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtDx->Fill(eta, v) ; }
2929 }
2930 }
2931 }
2932 else // cut is not used, fill in any case
2933 {
2934 if(inWin) { fHistFull[k].fHistFullHar[j].fHistV0vObsPt->Fill(pt, vFlip) ; }
2935 else
2936 {
2937 fHistFull[k].fHistFullHar[j].fHistV0sbvObsPt->Fill(pt, vFlip) ;
2938 if(sx) { fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtSx->Fill(eta, v) ; }
2939 else if(dx) { fHistFull[k].fHistFullHar[j].fHistV0sbvObsPtDx->Fill(eta, v) ; }
2940 }
2941 }
2942 // v_
2943 Bool_t etaPtNoCut = kTRUE;
2944 if(fPtRangevEta[1] > fPtRangevEta[0] && (pt < fPtRangevEta[0] || pt >= fPtRangevEta[1]))
2945 {
2946 etaPtNoCut = kFALSE;
2947 }
2948 if(fEtaRangevPt[1] > fEtaRangevPt[0] && (TMath::Abs(eta) < fEtaRangevPt[0] || TMath::Abs(eta) >= fEtaRangevPt[1]))
2949 {
2950 etaPtNoCut = kFALSE;
2951 }
2952 if(etaPtNoCut)
2953 {
2954 if(inWin) { fHistFull[k].fHistV0vObs->Fill(order, vFlip) ; }
2955 else
2956 {
2957 if(sx) { fHistFull[k].fHistV0sbvObsSx->Fill(order, vFlip) ; }
2958 else if(dx) { fHistFull[k].fHistV0sbvObsDx->Fill(order, vFlip) ; }
2959 }
2960 }
2961
2962 // Correlation of Phi of selected v0s with Psi
2963 float phi_i = phi;
2964 if(eta < 0 && oddHar)
2965 {
2966 phi_i += TMath::Pi() ; // backward particle and odd harmonic
2967 if(phi_i > 2*TMath::Pi()) { phi_i -= 2*TMath::Pi() ; }
2968 }
2969 float dPhi = phi_i - psi_i;
2970 if(dPhi < 0.) { dPhi += 2*TMath::Pi() ; }
2971
2972 if(inWin) { fHistFull[k].fHistFullHar[j].fHistV0PhiCorr->Fill(fmod((double)dPhi, 2*TMath::Pi() / order)) ; }
2973 else { fHistFull[k].fHistFullHar[j].fHistV0sbPhiCorr->Fill(fmod((double)dPhi, 2*TMath::Pi() / order)) ; }
2974 }
2975 }
2976 }
2977 //delete fFlowV0 ; fFlowV0 = 0 ;
2978 //delete daughterPlus ; daughterPlus = 0 ;
2979 //delete daughterMinus ; daughterMinus = 0 ;
2980 }
2981 fHistV0MultPart->Fill(corrMultV0) ;
2982
2983 return ;
2984}
2985//-----------------------------------------------------------------------
2986// void AliFlowAnalyser::FillLabels()
2987// {
2988// // Reads the tracks label (i.e. the link from ESD tracking to KineTree)
2989// // and ...
2990//
2991// cout << " Fill Labels . " << endl ;
2992//
2993// return ;
2994// }
2995//-----------------------------------------------------------------------
2996void AliFlowAnalyser::PrintEventQuantities()
2997{
2998 // prints event by event calculated quantities
2999
3000 cout << endl ; cout << " # " << fEventNumber << " - Event quantities : " << endl ; cout << endl ;
3001
3002 cout << "fQ[k][j] = " ;
3003 for(int k=0;k<Flow::nSels;k++)
3004 {
3005 for(int j=0;j<Flow::nHars;j++)
3006 {
3007 cout << (Float_t)fQ[k][j].X() << "," << (Float_t)fQ[k][j].Y() << " ; " ;
3008 }
3009 cout << endl ;
3010 cout << " " ;
3011 }
3012 cout << endl ; cout << "fPsi[k][j] = " ;
3013 for(int k=0;k<Flow::nSels;k++)
3014 {
3015 for(int j=0;j<Flow::nHars;j++) { Float_t aaa = (Float_t)fPsi[k][j] ; cout << aaa << " , " ; }
3016 cout << endl ;
3017 cout << " " ;
3018 }
3019 cout << endl ; cout << "fQnorm[k][j] = " ;
3020 for(int k=0;k<Flow::nSels;k++)
3021 {
3022 for(int j=0;j<Flow::nHars;j++) { Float_t aaa = (Float_t)fQnorm[k][j] ; cout << aaa << " , " ; }
3023 cout << endl ;
3024 cout << " " ;
3025 }
3026 cout << endl ; cout << "fMult[k][j] = " ;
3027 for(int k=0;k<Flow::nSels;k++)
3028 {
3029 for(int j=0;j<Flow::nHars;j++) { Float_t aaa = (Float_t)fMult[k][j] ; cout << aaa << " , " ; }
3030 cout << endl ;
3031 cout << " " ;
3032 }
3033 cout << endl ; cout << "fMultSub[n][k][j] = " ;
3034 for(int n=0;n<Flow::nSubs;n++)
3035 {
3036 for(int k=0;k<Flow::nSels;k++)
3037 {
3038 for(int j=0;j<Flow::nHars;j++) { Float_t aaa = fMultSub[n][k][j] ; cout << aaa << " , " ; }
3039 cout << endl ;
3040 cout << " " ;
3041 }
3042 }
3043 cout << endl ; cout << "fPsiSub[n][k][j] = " ;
3044 for(int n=0;n<Flow::nSubs;n++)
3045 {
3046 for(int k=0;k<Flow::nSels;k++)
3047 {
3048 for(int j=0;j<Flow::nHars;j++) { Float_t aaa = fPsiSub[n][k][j] ; cout << aaa << " , " ; }
3049 cout << endl ;
3050 cout << " " ;
3051 }
3052 }
3053 cout << endl ; cout << "Delta_PsiSub[k][j] = " ;
3054 for(int k=0;k<Flow::nSels;k++)
3055 {
3056 for(int j=0;j<Flow::nHars;j++) { Float_t aaa = fPsiSub[0][k][j]-fPsiSub[1][k][j] ; cout << aaa << " , " ; }
3057 cout << endl ;
3058 cout << " " ;
3059 }
3060 cout << endl ;
3061}
3062//----------------------------------------------------------------------
3063// ###
3064//----------------------------------------------------------------------
3065Bool_t AliFlowAnalyser::Resolution()
3066{
3067 // Calculates the resolution from cos(dPsi) and corrects the observed flow.
3068 // New histograms are then created and filled with the corrected vn
3069 // (see fVnResHistList), and saved in the otput file.
3070
3071 cout << " AliFlowAnalyser::Resolution()" << endl ;
3072
3073 // VnRes histogram collection
3074 fVnResHistList = new TOrdCollection(Flow::nSels*Flow::nHars);
3075
3076 // Calculate resolution from sqrt(fHistCos)
3077 double cosPair[Flow::nSels][Flow::nHars];
3078 double cosPairErr[Flow::nSels][Flow::nHars];
3079 double content, contentV0;
3080 double error, errorV0;
3081 double totalError;
3082 TString* histTitle;
3083
3084 for(int k = 0; k < Flow::nSels; k++)
3085 {
3086 // v for tracks (corrected for resolution)
3087 histTitle = new TString("Flow_v_Sel");
3088 *histTitle += k+1;
3089 fHistFull[k].fHistv = fHistFull[k].fHistvObs->ProjectionX(histTitle->Data());
3090 fHistFull[k].fHistv->SetTitle(histTitle->Data());
3091 fHistFull[k].fHistv->SetXTitle("Harmonic");
3092 fHistFull[k].fHistv->SetYTitle("v (%)");
3093 delete histTitle;
3094 fVnResHistList->AddLast(fHistFull[k].fHistv);
3095
3096 // v for v0s (corrected for resolution)
3097 histTitle = new TString("FlowV0_v_Sel");
3098 *histTitle += k+1;
3099 fHistFull[k].fHistV0v = fHistFull[k].fHistV0vObs->ProjectionX(histTitle->Data());
3100 fHistFull[k].fHistV0v->SetTitle(histTitle->Data());
3101 fHistFull[k].fHistV0v->SetXTitle("Harmonic");
3102 fHistFull[k].fHistV0v->SetYTitle("v (%)");
3103 delete histTitle;
3104 fVnResHistList->AddLast(fHistFull[k].fHistV0v);
3105
3106 for(int j = 0; j < Flow::nHars; j++)
3107 {
3108 double order = (double)(j+1);
3109 cosPair[k][j] = fHistFull[k].fHistCos->GetBinContent(j+1);
3110 cosPairErr[k][j] = fHistFull[k].fHistCos->GetBinError(j+1);
3111 if(cosPair[k][j] > 0.)
3112 {
3113 double resSub = 0. ;
3114 double resSubErr = 0. ;
3115 if(fV1Ep1Ep2 == kTRUE && order == 1) // calculate resolution of second order event plane first
3116 {
3117 double res2 = 0. ;
3118 double res2error = 0. ;
3119 if(fHistFull[k].fHistCos->GetBinContent(2) > 0.)
3120 {
3121 if(fEtaSub) // sub res only
3122 {
3123 res2 = TMath::Sqrt(fHistFull[k].fHistCos->GetBinContent(2));
3124 res2error = fHistFull[k].fHistCos->GetBinError(2) / (2. * res2);
3125 }
3126 else
3127 {
3128 if (fHistFull[k].fHistCos->GetBinContent(2) > 0.92) // resolution saturates
3129 {
3130 res2 = 0.99;
3131 res2error = 0.007;
3132 }
3133 else
3134 {
3135 double deltaRes2Sub = 0.005; // differential for the error propagation
3136 double res2Sub = TMath::Sqrt(fHistFull[k].fHistCos->GetBinContent(2));
3137 double res2SubErr = fHistFull[k].fHistCos->GetBinError(2) / (2. * res2Sub);
3138 double chiSub2 = Chi(res2Sub);
3139 double chiSub2Delta = Chi(res2Sub + deltaRes2Sub);
3140 res2 = ResEventPlane(TMath::Sqrt(2.) * chiSub2); // full event plane res.
3141 double mRes2Delta = ResEventPlane(TMath::Sqrt(2.) * chiSub2Delta);
3142 res2error = res2SubErr * fabs((double)res2 - mRes2Delta) / deltaRes2Sub;
3143 }
3144 }
3145 }
3146 else
3147 {
3148 res2 = 0.;
3149 res2error = 0.;
3150 }
3151 // now put everything together with first order event plane
3152 fRes[k][j] = TMath::Sqrt(cosPair[k][0]*res2);
3153 fResErr[k][j] = 1./(2.*fRes[k][j]) * TMath::Sqrt(cosPairErr[k][0]*cosPairErr[k][0] + res2error*res2error); // Gaussian error propagation
3154 }
3155 else if(fEtaSub) // sub res only
3156 {
3157 resSub = TMath::Sqrt(cosPair[k][j]);
3158 resSubErr = cosPairErr[k][j] / (2. * resSub);
3159 fRes[k][j] = resSub;
3160 fResErr[k][j] = resSubErr;
3161 }
3162 else if(order==4. || order==6.|| order==8.) // 2nd harmonic event plane
3163 {
3164 double deltaResSub = 0.005; // differential for the error propagation
3165 double resSub = TMath::Sqrt(cosPair[k][1]);
3166 double resSubErr = cosPairErr[k][1] / (2. * resSub);
3167 double chiSub = Chi(resSub);
3168 double chiSubDelta = Chi(resSub + deltaResSub);
3169 double mResDelta;
3170 if (order==4.)
3171 {
3172 fRes[k][j] = ResEventPlaneK2(TMath::Sqrt(2.) * chiSub); // full event plane res.
3173 mResDelta = ResEventPlaneK2(TMath::Sqrt(2.) * chiSubDelta);
3174 }
3175 else if(order==6.)
3176 {
3177 fRes[k][j] = ResEventPlaneK3(TMath::Sqrt(2.) * chiSub); // full event plane res.
3178 mResDelta = ResEventPlaneK3(TMath::Sqrt(2.) * chiSubDelta);
3179 }
3180 else
3181 {
3182 fRes[k][j] = ResEventPlaneK4(TMath::Sqrt(2.) * chiSub); // full event plane res.
3183 mResDelta = ResEventPlaneK4(TMath::Sqrt(2.) * chiSubDelta);
3184 }
3185 fResErr[k][j] = resSubErr * fabs((double)fRes[k][j] - mResDelta) / deltaResSub;
3186 }
3187 else
3188 {
3189 if(cosPair[k][j] > 0.92) // resolution saturates
3190 {
3191 fRes[k][j] = 0.99;
3192 fResErr[k][j] = 0.007;
3193 }
3194 else
3195 {
3196 double deltaResSub = 0.005; // differential for the error propagation
3197 double resSub = TMath::Sqrt(cosPair[k][j]);
3198 double resSubErr = cosPairErr[k][j] / (2. * resSub);
3199 double chiSub = Chi(resSub);
3200 double chiSubDelta = Chi(resSub + deltaResSub);
3201 fRes[k][j] = ResEventPlane(TMath::Sqrt(2.) * chiSub); // full event plane res.
3202 double mResDelta = ResEventPlane(TMath::Sqrt(2.) * chiSubDelta);
3203 fResErr[k][j] = resSubErr * TMath::Abs((double)fRes[k][j] - mResDelta) / deltaResSub;
3204 }
3205 }
3206 }
3207 else // subevent correlation must be positive
3208 {
3209 fRes[k][j] = 0.;
3210 fResErr[k][j] = 0.;
3211 }
3212 fHistFull[k].fHistRes->SetBinContent(j+1, fRes[k][j]);
3213 fHistFull[k].fHistRes->SetBinError(j+1, fResErr[k][j]);
3214
3215 // Create the v 2D histogram (Flow corrected for res.) - v,Pt,Eta
3216 histTitle = new TString("Flow_v2D_Sel");
3217 *histTitle += k+1;
3218 histTitle->Append("_Har");
3219 *histTitle += j+1;
3220 fHistFull[k].fHistFullHar[j].fHistv2D = fHistFull[k].fHistFullHar[j].fHistvObs2D->ProjectionXY(histTitle->Data());
3221 fHistFull[k].fHistFullHar[j].fHistv2D->SetTitle(histTitle->Data());
3222 fHistFull[k].fHistFullHar[j].fHistv2D->SetXTitle((char*)fLabel.Data());
3223 fHistFull[k].fHistFullHar[j].fHistv2D->SetYTitle("Pt (GeV/c)");
3224 fHistFull[k].fHistFullHar[j].fHistv2D->SetZTitle("v (%)");
3225 delete histTitle;
3226 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistv2D);
3227
3228 // Create the 1D v histograms - v,Eta
3229 histTitle = new TString("Flow_vEta_Sel");
3230 *histTitle += k+1;
3231 histTitle->Append("_Har");
3232 *histTitle += j+1;
3233 fHistFull[k].fHistFullHar[j].fHistvEta = fHistFull[k].fHistFullHar[j].fHistvObsEta->ProjectionX(histTitle->Data());
3234 fHistFull[k].fHistFullHar[j].fHistvEta->SetTitle(histTitle->Data());
3235 fHistFull[k].fHistFullHar[j].fHistvEta->SetXTitle((char*)fLabel.Data());
3236 fHistFull[k].fHistFullHar[j].fHistvEta->SetYTitle("v (%)");
3237 delete histTitle;
3238 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistvEta);
3239
3240 // v,Pt
3241 histTitle = new TString("Flow_vPt_Sel");
3242 *histTitle += k+1;
3243 histTitle->Append("_Har");
3244 *histTitle += j+1;
3245 fHistFull[k].fHistFullHar[j].fHistvPt = fHistFull[k].fHistFullHar[j].fHistvObsPt->ProjectionX(histTitle->Data());
3246 fHistFull[k].fHistFullHar[j].fHistvPt->SetTitle(histTitle->Data());
3247 fHistFull[k].fHistFullHar[j].fHistvPt->SetXTitle("Pt (GeV/c)");
3248 fHistFull[k].fHistFullHar[j].fHistvPt->SetYTitle("v (%)");
3249 delete histTitle;
3250 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistvPt);
3251
3252 // Create the v 2D histogram (V0s Flow corrected for res.) - v,Pt,Eta
3253 histTitle = new TString("FlowV0_v2D_Sel");
3254 *histTitle += k+1;
3255 histTitle->Append("_Har");
3256 *histTitle += j+1;
3257 fHistFull[k].fHistFullHar[j].fHistV0v2D = fHistFull[k].fHistFullHar[j].fHistV0vObs2D->ProjectionXY(histTitle->Data());
3258 fHistFull[k].fHistFullHar[j].fHistV0v2D->SetTitle(histTitle->Data());
3259 fHistFull[k].fHistFullHar[j].fHistV0v2D->SetXTitle((char*)fLabel.Data());
3260 fHistFull[k].fHistFullHar[j].fHistV0v2D->SetYTitle("Pt (GeV/c)");
3261 fHistFull[k].fHistFullHar[j].fHistV0v2D->SetZTitle("v (%)");
3262 delete histTitle;
3263 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0v2D);
3264
3265 // Create the 1D v histograms (V0s) - v,Eta
3266 histTitle = new TString("FlowV0_vEta_Sel");
3267 *histTitle += k+1;
3268 histTitle->Append("_Har");
3269 *histTitle += j+1;
3270 fHistFull[k].fHistFullHar[j].fHistV0vEta = fHistFull[k].fHistFullHar[j].fHistV0vObsEta->ProjectionX(histTitle->Data());
3271 fHistFull[k].fHistFullHar[j].fHistV0vEta->SetTitle(histTitle->Data());
3272 fHistFull[k].fHistFullHar[j].fHistV0vEta->SetXTitle((char*)fLabel.Data());
3273 fHistFull[k].fHistFullHar[j].fHistV0vEta->SetYTitle("v (%)");
3274 delete histTitle;
3275 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0vEta);
3276
3277 // (V0s) v,Pt
3278 histTitle = new TString("FlowV0_vPt_Sel");
3279 *histTitle += k+1;
3280 histTitle->Append("_Har");
3281 *histTitle += j+1;
3282 fHistFull[k].fHistFullHar[j].fHistV0vPt = fHistFull[k].fHistFullHar[j].fHistV0vObsPt->ProjectionX(histTitle->Data());
3283 fHistFull[k].fHistFullHar[j].fHistV0vPt->SetTitle(histTitle->Data());
3284 fHistFull[k].fHistFullHar[j].fHistV0vPt->SetXTitle("Pt (GeV/c)");
3285 fHistFull[k].fHistFullHar[j].fHistV0vPt->SetYTitle("v (%)");
3286 delete histTitle;
3287 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0vPt);
3288
3289 // Create the v 2D histogram (V0s sidebands Flow corrected for res.) - v,Pt,Eta
3290 histTitle = new TString("FlowV0sb_v2D_Sel");
3291 *histTitle += k+1;
3292 histTitle->Append("_Har");
3293 *histTitle += j+1;
3294 fHistFull[k].fHistFullHar[j].fHistV0sbv2D = fHistFull[k].fHistFullHar[j].fHistV0sbvObs2D->ProjectionXY(histTitle->Data());
3295 fHistFull[k].fHistFullHar[j].fHistV0sbv2D->SetTitle(histTitle->Data());
3296 fHistFull[k].fHistFullHar[j].fHistV0sbv2D->SetXTitle((char*)fLabel.Data());
3297 fHistFull[k].fHistFullHar[j].fHistV0sbv2D->SetYTitle("Pt (GeV/c)");
3298 fHistFull[k].fHistFullHar[j].fHistV0sbv2D->SetZTitle("v (%)");
3299 delete histTitle;
3300 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0sbv2D);
3301
3302 // Create the 1D v histograms (V0s sidebands) - v,Eta
3303 histTitle = new TString("FlowV0sb_vEta_Sel");
3304 *histTitle += k+1;
3305 histTitle->Append("_Har");
3306 *histTitle += j+1;
3307 fHistFull[k].fHistFullHar[j].fHistV0sbvEta = fHistFull[k].fHistFullHar[j].fHistV0sbvObsEta->ProjectionX(histTitle->Data());
3308 fHistFull[k].fHistFullHar[j].fHistV0sbvEta->SetTitle(histTitle->Data());
3309 fHistFull[k].fHistFullHar[j].fHistV0sbvEta->SetXTitle((char*)fLabel.Data());
3310 fHistFull[k].fHistFullHar[j].fHistV0sbvEta->SetYTitle("v (%)");
3311 delete histTitle;
3312 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0sbvEta);
3313
3314 // (V0s sidebands) v,Pt
3315 histTitle = new TString("FlowV0sb_vPt_Sel");
3316 *histTitle += k+1;
3317 histTitle->Append("_Har");
3318 *histTitle += j+1;
3319 fHistFull[k].fHistFullHar[j].fHistV0sbvPt = fHistFull[k].fHistFullHar[j].fHistV0sbvObsPt->ProjectionX(histTitle->Data());
3320 fHistFull[k].fHistFullHar[j].fHistV0sbvPt->SetTitle(histTitle->Data());
3321 fHistFull[k].fHistFullHar[j].fHistV0sbvPt->SetXTitle("Pt (GeV/c)");
3322 fHistFull[k].fHistFullHar[j].fHistV0sbvPt->SetYTitle("v (%)");
3323 delete histTitle;
3324 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0sbvPt);
3325
3326 // Calulate v = vObs / Resolution
3327 if(fRes[k][j])
3328 {
3329 cout << "***** Resolution of the " << j+1 << "th harmonic = " << fRes[k][j] << " +/- " << fResErr[k][j] << endl;
3330 // selected tracks
3331 fHistFull[k].fHistFullHar[j].fHistv2D->Scale(1. / fRes[k][j]);
3332 fHistFull[k].fHistFullHar[j].fHistvEta->Scale(1. / fRes[k][j]);
3333 fHistFull[k].fHistFullHar[j].fHistvPt->Scale(1. / fRes[k][j]);
3334 content = fHistFull[k].fHistv->GetBinContent(j+1);
3335 content /= fRes[k][j];
3336 fHistFull[k].fHistv->SetBinContent(j+1, content);
3337 // selected V0s
3338 fHistFull[k].fHistFullHar[j].fHistV0v2D->Scale(1. / fRes[k][j]);
3339 fHistFull[k].fHistFullHar[j].fHistV0vEta->Scale(1. / fRes[k][j]);
3340 fHistFull[k].fHistFullHar[j].fHistV0vPt->Scale(1. / fRes[k][j]);
3341 contentV0 = fHistFull[k].fHistV0v->GetBinContent(j+1);
3342 contentV0 /= fRes[k][j];
3343 fHistFull[k].fHistV0v->SetBinContent(j+1, contentV0);
3344 // V0s sidebands
3345 fHistFull[k].fHistFullHar[j].fHistV0sbv2D->Scale(1. / fRes[k][j]);
3346 fHistFull[k].fHistFullHar[j].fHistV0sbvEta->Scale(1. / fRes[k][j]);
3347 fHistFull[k].fHistFullHar[j].fHistV0sbvPt->Scale(1. / fRes[k][j]);
3348 // The systematic error of the resolution is folded in.
3349 // tracks
3350 error = fHistFull[k].fHistv->GetBinError(j+1) ;
3351 error /= fRes[k][j];
3352 totalError = fabs(content) * TMath::Sqrt((error/content)*(error/content) + (fResErr[k][j]/fRes[k][j])*(fResErr[k][j]/fRes[k][j]));
3353 fHistFull[k].fHistv->SetBinError(j+1, totalError);
3354 cout << "***** v" << j+1 << "= (" << content << " +/- " << error << " +/- " << totalError << "(with syst.)) %" << endl;
3355 // V0s
3356 errorV0 = fHistFull[k].fHistV0v->GetBinError(j+1) ;
3357 errorV0 /= fRes[k][j];
3358 if (contentV0>0)
3359 totalError = fabs(contentV0) * TMath::Sqrt((errorV0/contentV0)*(errorV0/contentV0) + (fResErr[k][j]/fRes[k][j])*(fResErr[k][j]/fRes[k][j]));
3360 else
3361 totalError = 0;
3362 fHistFull[k].fHistV0v->SetBinError(j+1, totalError);
3363 }
3364 else
3365 {
3366 cout << "##### Resolution of the " << j+1 << "th harmonic was zero." << endl;
3367 // tracks hist
3368 fHistFull[k].fHistFullHar[j].fHistv2D->Reset();
3369 fHistFull[k].fHistFullHar[j].fHistvEta->Reset();
3370 fHistFull[k].fHistFullHar[j].fHistvPt->Reset();
3371 fHistFull[k].fHistv->SetBinContent(j+1, 0.);
3372 fHistFull[k].fHistv->SetBinError(j+1, 0.);
3373 // v0s hist
3374 fHistFull[k].fHistFullHar[j].fHistV0v2D->Reset();
3375 fHistFull[k].fHistFullHar[j].fHistV0vEta->Reset();
3376 fHistFull[k].fHistFullHar[j].fHistV0vPt->Reset();
3377 fHistFull[k].fHistFullHar[j].fHistV0sbv2D->Reset();
3378 fHistFull[k].fHistFullHar[j].fHistV0sbvEta->Reset();
3379 fHistFull[k].fHistFullHar[j].fHistV0sbvPt->Reset();
3380 fHistFull[k].fHistV0v->SetBinContent(j+1, 0.);
3381 fHistFull[k].fHistV0v->SetBinError(j+1, 0.);
3382 }
3383 }
3384 }
3385
3386 return kTRUE ;
3387}
3388//-----------------------------------------------------------------------
3389Double_t AliFlowAnalyser::Chi(Double_t res)
3390{
3391 // chi from the event plane resolution
3392
3393 Double_t chi = 2.0;
3394 Double_t delta = 1.0;
3395 for(int i = 0; i < 15; i++)
3396 {
3397 if(ResEventPlane(chi) < res) { chi = chi + delta ; }
3398 else { chi = chi - delta ; }
3399 delta = delta / 2.;
3400 }
3401
3402 return chi ;
3403}
3404//-----------------------------------------------------------------------
3405Double_t AliFlowAnalyser::ResEventPlane(Double_t chi)
3406{
3407 // plane resolution as function of chi
3408
3409 Double_t con = TMath::Sqrt(TMath::Pi()/2)/2 ; // ~ 0.626657
3410 Double_t arg = chi * chi / 4.;
3411 Double_t res = con * chi * exp(-arg) * (TMath::BesselI0(arg) + TMath::BesselI1(arg));
3412
3413 return res ;
3414}
3415//-----------------------------------------------------------------------
3416Double_t AliFlowAnalyser::ResEventPlaneK2(Double_t chi)
3417{
3418 // ... for the case k=2.
3419
3420 Double_t con = TMath::Sqrt(TMath::Pi()/2)/2 ; // ~ 0.626657
3421 Double_t arg = chi * chi / 4.;
3422
3423 Double_t besselOneHalf = TMath::Sqrt(arg/(TMath::Pi()/2)) * sinh(arg)/arg;
3424 Double_t besselThreeHalfs = TMath::Sqrt(arg/(TMath::Pi()/2)) * (cosh(arg)/arg - sinh(arg)/(arg*arg));
3425 Double_t res = con * chi * exp(-arg) * (besselOneHalf + besselThreeHalfs);
3426
3427 return res ;
3428}
3429//-----------------------------------------------------------------------
3430Double_t AliFlowAnalyser::ResEventPlaneK3(Double_t chi)
3431{
3432 // ...for the case k=3.
3433
3434 Double_t con = TMath::Sqrt(TMath::Pi()/2)/2 ; // ~ 0.626657
3435 Double_t arg = chi * chi / 4.;
3436 Double_t res = con * chi * exp(-arg) * (TMath::BesselI1(arg) + TMath::BesselI(2, arg));
3437
3438 return res ;
3439}
3440//-----------------------------------------------------------------------
3441Double_t AliFlowAnalyser::ResEventPlaneK4(Double_t chi)
3442{
3443 // ...for the case k=4.
3444
3445 Double_t con = TMath::Sqrt(TMath::Pi()/2)/2 ; // ~ 0.626657
3446 Double_t arg = chi * chi / 4.;
3447
3448 Double_t besselOneHalf = TMath::Sqrt(arg/(TMath::Pi()/2)) * sinh(arg)/arg;
3449 Double_t besselThreeHalfs = TMath::Sqrt(arg/(TMath::Pi()/2)) * (cosh(arg)/arg - sinh(arg)/(arg*arg));
3450 Double_t besselFiveHalfs = besselOneHalf - 3*besselThreeHalfs/arg;
3451 Double_t res = con * chi * exp(-arg) * (besselThreeHalfs + besselFiveHalfs);
3452
3453 return res ;
3454}
3455//-------------------------------------------------------------
3456// ###
3457//-----------------------------------------------------------------------
3458
3459
3460#endif