1 //////////////////////////////////////////////////////////////////////
5 // Author: Emanuele Simili
7 //////////////////////////////////////////////////////////////////////
8 //_____________________________________________________________
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.
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).
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
44 //////////////////////////////////////////////////////////////////////
46 #ifndef ALIFLOWANALYSER_CXX
47 #define ALIFLOWANALYSER_CXX
58 #include <TObjArray.h>
59 #include <TOrdCollection.h>
60 #include <TParticle.h>
61 #include <TParticlePDG.h>
62 #include <TDatabasePDG.h>
67 #include <TProfile2D.h>
71 #include "AliFlowEvent.h"
72 #include "AliFlowTrack.h"
73 #include "AliFlowV0.h"
74 #include "AliFlowConstants.h"
75 #include "AliFlowSelection.h"
76 #include "AliFlowAnalyser.h"
83 using namespace std; //required for resolving the 'cout' symbol
85 ClassImp(AliFlowAnalyser) ;
86 //-----------------------------------------------------------------------
87 AliFlowAnalyser::AliFlowAnalyser(const AliFlowSelection* flowSelect)
89 // default constructor (selection given or default selection)
91 if(flowSelect) { fFlowSelect = new AliFlowSelection(*flowSelect) ; }
92 else { fFlowSelect = new AliFlowSelection() ; }
94 // output file (histograms)
95 fHistFileName = "flowAnalysPlot.root" ;
100 fPhiBins = Flow::nPhiBins ;
102 fPhiMax = 2*TMath::Pi() ;
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
118 fPhiWgtHistList = 0 ;
121 //-----------------------------------------------------------------------
122 AliFlowAnalyser::~AliFlowAnalyser()
124 // default distructor (no actions)
126 //-----------------------------------------------------------------------
127 Bool_t AliFlowAnalyser::Init()
129 // sets some defaults for the analysis
131 cout << "* FlowAnalysis * - Init()" << endl ; cout << endl ;
133 // Open output files (->plots)
134 fHistFile = new TFile(fHistFileName.Data(), "RECREATE");
137 // counters and pointers to 0
142 //fNumberOfEvents = 0 ;
143 fNumberOfTracks = 0 ;
150 for(Int_t ii=0;ii<3;ii++) { fVertex[ii] = 0 ; }
152 // Check for Reading weights: if weight file is not there, wgt arrays are not used (phiwgt & bayesian)
155 fReadPhiWgt = kFALSE ;
157 cout << " No wgt file. Phi & Bayesian weights will not be used ! " << endl ;
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 ;
171 SetPtRangevEta(1.,0.); // (Flow::fPtMin , Flow::fPtMax) ;
172 SetEtaRangevPt(1.,0.); // (Flow::fEtaMin , Flow::fEtaMax) ;
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. ;
234 enum { nTriggerBins = 12,
238 nFitPtsBinsTpc = 161,
240 nFitPtsBinsTrd = 131,
242 nFitOverMaxBins = 55,
243 nMultOverOrigBins =100,
245 nVertexXYBins = 1000,
263 // Histograms Booking ...
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");
270 fHistCharge = new TH1F("Flow_Charge", "Flow_Charge", nChargeBins, chargeMin, chargeMax);
271 fHistCharge->SetXTitle("Charge");
272 fHistCharge->SetYTitle("Counts");
274 // Reconstructed Momentum (constrained, if so)
275 fHistPtot = new TH1F("Flow_totalP","Flow_totalP", fPtBins, fPtMin, fPtMax);
277 fHistPtot->SetXTitle("P (GeV/c)");
278 fHistPtot->SetYTitle("Counts");
280 // Reconstructed pT (constrained, if so)
281 fHistPt = new TH1F("Flow_Pt","Flow_Pt", nMomenBins, pMin, pMax);
283 fHistPt->SetXTitle("Pt (GeV/c)");
284 fHistPt->SetYTitle("Counts");
286 // Reconstructed P.id vs pT
287 fHistPidPt = new TH2F("Flow_PidPt","Flow_PidPt", 12, 0., 12., fPtBins, fPtMin, fPtMax);
289 fHistPidPt->SetXTitle("e-, e+, mu-, mu+, pi-, pi+, K-, K+, pr-, pr+, d-, d+");
290 fHistPidPt->SetYTitle("Pt (GeV/c)");
292 // Distance of closest approach - Transverse, Unsigned
293 fHistDca = new TH1F("Flow_TransDCA", "Flow_TransverseDCA", nDcaBins, dcaMin, dcaMax);
295 fHistDca->SetXTitle("|Track's Signed dca (cm)|");
296 fHistDca->SetYTitle("Counts");
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");
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");
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");
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");
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)");
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)");
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");
340 fHistChi2TPC = new TH1F("Flow_Chi2_TPC", "Flow_Chi2_TPC", nChi2Bins, chi2Min, chi2Max);
341 fHistChi2TPC->SetXTitle("Chi square for TPC");
342 fHistChi2TPC->SetYTitle("Counts");
344 fHistChi2ITS = new TH1F("Flow_Chi2_ITS", "Flow_Chi2_ITS", nChi2Bins, chi2Min, chi2MaxI);
345 fHistChi2ITS->SetXTitle("Chi square for ITS");
346 fHistChi2ITS->SetYTitle("Counts");
348 fHistChi2TRD = new TH1F("Flow_Chi2_TRD", "Flow_Chi2_TRD", nChi2Bins, chi2Min, chi2Max);
349 fHistChi2TRD->SetXTitle("Chi square for TRD");
350 fHistChi2TRD->SetYTitle("Counts");
352 fHistChi2TOF = new TH1F("Flow_Chi2_TOF", "Flow_Chi2_TOF", nChi2Bins, chi2Min, chi2Max);
353 fHistChi2TOF->SetXTitle("Chi square for TOF");
354 fHistChi2TOF->SetYTitle("Counts");
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");
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");
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");
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");
376 fHistFitPtsTPC = new TH1F("Flow_FitPts_TPC", "Flow_FitPts_TPC", nFitPtsBinsTpc, fitPtsMin, fitPtsMaxTpc);
377 fHistFitPtsTPC->SetXTitle("Fit Points");
378 fHistFitPtsTPC->SetYTitle("Counts");
380 fHistFitPtsITS = new TH1F("Flow_HitPts_ITS", "Flow_HitPts_ITS", nFitPtsBinsIts, fitPtsMin, fitPtsMaxIts);
381 fHistFitPtsITS->SetXTitle("Fit Points");
382 fHistFitPtsITS->SetYTitle("Counts");
384 fHistFitPtsTRD = new TH1F("Flow_FitPts_TRD", "Flow_FitPts_TRD", nFitPtsBinsTrd, fitPtsMin, fitPtsMaxTrd);
385 fHistFitPtsTRD->SetXTitle("Fit Points");
386 fHistFitPtsTRD->SetYTitle("Counts");
388 fHistFitPtsTOF = new TH1F("Flow_HitPts_TOF", "Flow_HitPts_TOF", nFitPtsBinsTof, fitPtsMin, fitPtsMaxTof);
389 fHistFitPtsTOF->SetXTitle("Fit Points");
390 fHistFitPtsTOF->SetYTitle("Counts");
394 fHistMaxPtsTPC = new TH1F("Flow_MaxPts_TPC", "Flow_MaxPts_TPC", nFitPtsBinsTpc, fitPtsMin, fitPtsMaxTpc);
395 fHistMaxPtsTPC->SetXTitle("Max Points");
396 fHistMaxPtsTPC->SetYTitle("Counts");
398 fHistMaxPtsITS = new TH1F("Flow_MaxPts_ITS", "Flow_MaxPts_ITS", nFitPtsBinsIts, fitPtsMin, fitPtsMaxIts);
399 fHistMaxPtsITS->SetXTitle("Max Points");
400 fHistMaxPtsITS->SetYTitle("Counts");
402 fHistMaxPtsTRD = new TH1F("Flow_MaxPts_TRD", "Flow_MaxPts_TRD", nFitPtsBinsTrd, fitPtsMin, fitPtsMaxTrd);
403 fHistMaxPtsTRD->SetXTitle("Max Points");
404 fHistMaxPtsTRD->SetYTitle("Counts");
406 fHistMaxPtsTOF = new TH1F("Flow_MaxPts_TOF", "Flow_MaxPts_TOF", nFitPtsBinsTof, fitPtsMin, fitPtsMaxTof);
407 fHistMaxPtsTOF->SetXTitle("Max Points");
408 fHistMaxPtsTOF->SetYTitle("Counts");
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");
416 fHistFitOverMax = new TH1F("Flow_FitOverMax", "Flow_FitOverMax", nFitOverMaxBins, fitOverMaxMin, fitOverMaxMax);
417 fHistFitOverMax->SetXTitle("(Fit Points - 1) / Max Points");
418 fHistFitOverMax->SetYTitle("Counts");
421 fHistLenght = new TH1F("Flow_TrackLenght", "Flow_TrackLenght", nLgBins, lgMin, lgMax);
422 fHistLenght->SetXTitle("Lenght of the Track (cm)");
423 fHistLenght->SetYTitle("Counts");
426 fHistOrigMult = new TH1F("Flow_OrigMult", "Flow_OrigMult", nMultBins, multMin, multMax);
427 fHistOrigMult->SetXTitle("Original Mult");
428 fHistOrigMult->SetYTitle("Counts");
431 fHistMultEta = new TH1F("Flow_MultEta", "Flow_MultEta", nMultBins, multMin, multMax);
432 fHistMultEta->SetXTitle("Mult for Centrality");
433 fHistMultEta->SetYTitle("Counts");
436 fHistMult = new TH1F("Flow_Mult", "Flow_Mult", nMultBins, multMin, multMax);
437 fHistMult->SetXTitle("Mult");
438 fHistMult->SetYTitle("Counts");
441 fHistV0Mult = new TH1F("FlowV0_Mult","FlowV0_Mult", nMultBins, multMin, multV0);
442 fHistV0Mult->SetXTitle("V0s Multiplicity");
443 fHistV0Mult->SetYTitle("Counts");
446 fHistMultOverOrig = new TH1F("Flow_MultOverOrig", "Flow_MultOverOrig", nMultOverOrigBins, multOverOrigMin, multOverOrigMax);
447 fHistMultOverOrig->SetXTitle("Mult / Orig. Mult");
448 fHistMultOverOrig->SetYTitle("Counts");
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");
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");
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");
466 fHistVertexZ = new TH1F("Flow_VertexZ", "Flow_VertexZ", nVertexZBins, vertexZMin, vertexZMax);
467 fHistVertexZ->SetXTitle("Vertex Z (cm)");
468 fHistVertexZ->SetYTitle("Counts");
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)");
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");
481 fHistEtaSym = new TH1F("Flow_EtaSym_TPC", "Flow_EtaSym_TPC", nEtaSymBins, etaSymMin, etaSymMax);
482 fHistEtaSym->SetXTitle("Eta Symmetry Ratio TPC");
483 fHistEtaSym->SetYTitle("Counts");
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");
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");
496 fHistPhi = new TH1F("Flow_xPhi", "Flow_xPhi (all)", fPhiBins, fPhiMin, fPhiMax);
497 fHistPhi->SetXTitle("Phi (rad)");
498 fHistPhi->SetYTitle("Counts");
501 fHistPhiCons = new TH1F("Flow_cPhi", "Flow_cPhi", fPhiBins, fPhiMin, fPhiMax);
502 fHistPhiCons->SetXTitle("cPhi (rad)");
503 fHistPhiCons->SetYTitle("Counts");
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)");
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)");
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)");
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)");
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)");
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)");
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)");
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)");
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)");
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)");
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)");
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)");
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)");
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");
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");
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");
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());
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)");
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)> (%)");
614 fHistPidPiPlus = new TH1F("Flow_PidPiPlus", "Flow_PidPiPlus", nPidBins, pidMin, pidMax);
615 fHistPidPiPlus->SetXTitle("ALICE P.Id.");
616 fHistPidPiPlus->SetYTitle("Counts");
619 fHistPidPiMinus = new TH1F("Flow_PidPiMinus", "Flow_PidPiMinus", nPidBins, pidMin, pidMax);
620 fHistPidPiMinus->SetXTitle("ALICE P.Id.");
621 fHistPidPiMinus->SetYTitle("Counts");
624 fHistPidProton = new TH1F("Flow_PidProton", "Flow_PidProton", nPidBins, pidMin, pidMax);
625 fHistPidProton->SetXTitle("ALICE P.Id.");
626 fHistPidProton->SetYTitle("Counts");
629 fHistPidAntiProton = new TH1F("Flow_PidAntiProton", "Flow_PidAntiProton", nPidBins, pidMin, pidMax);
630 fHistPidAntiProton->SetXTitle("ALICE P.Id.");
631 fHistPidAntiProton->SetYTitle("Counts");
634 fHistPidKplus = new TH1F("Flow_PidKplus", "Flow_PidKplus", nPidBins, pidMin, pidMax);
635 fHistPidKplus->SetXTitle("ALICE P.Id.");
636 fHistPidKplus->SetYTitle("Counts");
639 fHistPidKminus = new TH1F("Flow_PidKminus", "Flow_PidKminus", nPidBins, pidMin, pidMax);
640 fHistPidKminus->SetXTitle("ALICE P.Id.");
641 fHistPidKminus->SetYTitle("Counts");
644 fHistPidDeuteron = new TH1F("Flow_PidDeuteron", "Flow_PidDeuteron", nPidBins, pidMin, pidMax);
645 fHistPidDeuteron->SetXTitle("ALICE P.Id.");
646 fHistPidDeuteron->SetYTitle("Counts");
649 fHistPidAntiDeuteron = new TH1F("Flow_PidAntiDeuteron", "Flow_PidAntiDeuteron", nPidBins, pidMin, pidMax);
650 fHistPidAntiDeuteron->SetXTitle("ALICE P.Id.");
651 fHistPidAntiDeuteron->SetYTitle("Counts");
654 fHistPidElectron = new TH1F("Flow_PidElectron", "Flow_PidElectron", nPidBins, pidMin, pidMax);
655 fHistPidElectron->SetXTitle("ALICE P.Id.");
656 fHistPidElectron->SetYTitle("Counts");
659 fHistPidPositron = new TH1F("Flow_PidPositron", "Flow_PidPositron", nPidBins, pidMin, pidMax);
660 fHistPidPositron->SetXTitle("ALICE P.Id.");
661 fHistPidPositron->SetYTitle("Counts");
664 fHistPidMuonPlus = new TH1F("Flow_PidMuonPlus", "Flow_PidMuonPlus", nPidBins, pidMin, pidMax);
665 fHistPidMuonPlus->SetXTitle("ALICE P.Id.");
666 fHistPidMuonPlus->SetYTitle("Counts");
669 fHistPidMuonMinus = new TH1F("Flow_PidMuonMinus", "Flow_PidMuonMinus", nPidBins, pidMin, pidMax);
670 fHistPidMuonMinus->SetXTitle("ALICE P.Id.");
671 fHistPidMuonMinus->SetYTitle("Counts");
674 fHistPidPiPlusPart = new TH1F("Flow_PidPiPlusPart", "Flow_PidPiPlusPart", nPidBins, pidMin, pidMax);
675 fHistPidPiPlusPart->SetXTitle("ALICE P.Id. (pid = pi+)");
676 fHistPidPiPlusPart->SetYTitle("Counts");
679 fHistPidPiMinusPart = new TH1F("Flow_PidPiMinusPart", "Flow_PidPiMinusPart", nPidBins, pidMin, pidMax);
680 fHistPidPiMinusPart->SetXTitle("ALICE P.Id. (pid = pi-)");
681 fHistPidPiMinusPart->SetYTitle("Counts");
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");
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");
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");
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");
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");
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");
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");
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");
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");
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");
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");
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");
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");
751 fHistCent = new TH1F("Flow_Cent", "Flow_Cent", nCentBins, centMin, centMax);
752 fHistCent->SetXTitle("Centrality Bin");
753 fHistCent->SetYTitle("Counts");
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");
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");
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 (+)");
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 (-)");
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 (+)");
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 (-)");
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");
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");
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");
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");
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 (+)");
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 (-)");
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 (+)");
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 (-)");
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+)");
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-)");
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+)");
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-)");
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+)");
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-)");
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-)");
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-)");
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+)");
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-)");
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+)");
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-)");
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+)");
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-)");
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+)");
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-)");
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-)");
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-)");
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+)");
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-)");
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+)");
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-)");
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+)");
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-)");
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+)");
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-)");
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-)");
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-)");
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+)");
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+)");
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-)");
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+)");
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-)");
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+)");
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-)");
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-)");
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-)");
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+)");
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-)");
1028 for (int n = 0; n < Flow::nSubs; n++) // for sub-events
1030 for (int k = 0; k < Flow::nSels; k++)
1032 for (int j = 0; j < Flow::nHars; j++)
1034 float order = (float)(j + 1);
1035 int i = Flow::nSubs * k + n ;
1038 histTitle = new TString("Flow_Psi_Sub");
1040 histTitle->Append("_Sel");
1042 histTitle->Append("_Har");
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");
1052 if(fV0loop) // All V0s (if there, if flag on)
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");
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");
1071 fHistV0Chi2 = new TH1F("FlowV0_Chi2", "FlowV0_Chi2", nChi2Bins, chi2Min, chi2MaxC);
1072 fHistV0Chi2->SetXTitle("Chi square at Main Vertex");
1073 fHistV0Chi2->SetYTitle("Counts");
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)");
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)");
1102 fHistV0MassWin = new TH1F("FlowV0_MassWinPart", "FlowV0_MassWinPart", nMassBins, massMin, massMax);
1103 fHistV0MassWin->SetXTitle("Invariant Mass (GeV)");
1104 fHistV0MassWin->SetYTitle("Counts");
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");
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");
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)");
1156 for (int k = 0; k < Flow::nSels; k++) // for each selection
1159 histTitle = new TString("Flow_Cos_Sel");
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)>");
1167 histTitle = new TString("Flow_Res_Sel");
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");
1175 histTitle = new TString("Flow_vObs_Sel");
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 (%)");
1183 histTitle = new TString("FlowV0_vObs_Sel");
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 (%)");
1190 // vObs V0 sideband SX
1191 histTitle = new TString("FlowV0sb_vObs_sx_Sel");
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 (%)");
1198 // vObs V0 sideband DX
1199 histTitle = new TString("FlowV0sb_vObs_dx_Sel");
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 (%)");
1206 // PID for tracks used in R.P.
1207 histTitle = new TString("Flow_BayPidMult_Sel");
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");
1215 for (int j = 0; j < Flow::nHars; j++) // for each harmonic
1217 float order = (float)(j+1);
1220 histTitle = new TString("Flow_Mul_Sel");
1222 histTitle->Append("_Har");
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");
1230 histTitle = new TString("Flow_Psi_Sel");
1232 histTitle->Append("_Har");
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");
1239 // event plane difference of two selections
1240 histTitle = new TString("Flow_Psi_Diff_Sel");
1242 histTitle->Append("_Har");
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.);
1252 fHistFull[k].fHistFullHar[j].fHistPsiDiff = new TH1F(histTitle->Data(), histTitle->Data(), nPsiBins, -psiMax/2., psiMax/2.);
1258 fHistFull[k].fHistFullHar[j].fHistPsiDiff->SetXTitle("#Psi_{1,Sel1} - #Psi_{1,Sel2}(rad)");
1262 fHistFull[k].fHistFullHar[j].fHistPsiDiff->SetXTitle("#Psi_{2,Sel1} - #Psi_{2,Sel2}(rad)");
1269 fHistFull[k].fHistFullHar[j].fHistPsiDiff->SetXTitle("#Psi_{1,Sel1} - #Psi_{2,Sel2}(rad)");
1273 fHistFull[k].fHistFullHar[j].fHistPsiDiff->SetXTitle("#Psi_{1,Sel1} - #Psi_{2,Sel1}(rad)");
1276 fHistFull[k].fHistFullHar[j].fHistPsiDiff->SetYTitle("Counts");
1279 // correlation of sub-event planes
1280 histTitle = new TString("Flow_Psi_Sub_Corr_Sel");
1282 histTitle->Append("_Har");
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");
1290 // correlation of sub-event planes of different order
1291 histTitle = new TString("Flow_Psi_Sub_Corr_Diff_Sel");
1293 histTitle->Append("_Har");
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");
1302 histTitle = new TString("Flow_NormQ_Sel");
1304 histTitle->Append("_Har");
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");
1312 // particle-plane azimuthal correlation
1313 histTitle = new TString("Flow_Phi_Corr_Sel");
1315 histTitle->Append("_Har");
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");
1323 // neutral particle-plane azimuthal correlation
1324 histTitle = new TString("FlowV0_Phi_Corr_Sel");
1326 histTitle->Append("_Har");
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");
1334 // neutral sidebands-plane azimuthal correlation
1335 histTitle = new TString("FlowV0sb_Phi_Corr_Sel");
1337 histTitle->Append("_Har");
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");
1346 histTitle = new TString("Flow_YieldPt_Sel");
1348 histTitle->Append("_Har");
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");
1357 histTitle = new TString("Flow_EtaPtPhi3D_Sel");
1359 histTitle->Append("_Har");
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)");
1368 histTitle = new TString("Flow_Yield2D_Sel");
1370 histTitle->Append("_Har");
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)");
1379 histTitle = new TString("Flow_3dDca_Sel") ;
1381 histTitle->Append("_Har");
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)|");
1388 // Yield(pt) - excluded from R.P.
1389 histTitle = new TString("Flow_YieldPt_outSel");
1391 histTitle->Append("_Har");
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");
1399 // EtaPtPhi - excluded from R.P.
1400 histTitle = new TString("Flow_EtaPtPhi3D_outSel");
1402 histTitle->Append("_Har");
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)");
1410 // Yield(eta,pt) - excluded from R.P.
1411 histTitle = new TString("Flow_Yield2D_outSel");
1413 histTitle->Append("_Har");
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)");
1421 // Dca - 3D - excluded from R.P.
1422 histTitle = new TString("Flow_3dDca_outSel") ;
1424 histTitle->Append("_Har");
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)|");
1431 // Flow observed - v_obs,Pt,Eta
1432 histTitle = new TString("Flow_vObs2D_Sel");
1434 histTitle->Append("_Har");
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)");
1442 histTitle = new TString("Flow_vObsEta_Sel");
1444 histTitle->Append("_Har");
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 (%)");
1452 histTitle = new TString("Flow_vObsPt_Sel");
1454 histTitle->Append("_Har");
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 (%)");
1461 // neutral Flow observed - Pt,Eta
1462 histTitle = new TString("FlowV0_vObs2D_Sel");
1464 histTitle->Append("_Har");
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)");
1471 // neutral Flow observed - Eta
1472 histTitle = new TString("FlowV0_vObsEta_Sel");
1474 histTitle->Append("_Har");
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 (%)");
1481 // neutral Flow observed - Pt
1482 histTitle = new TString("FlowV0_vObsPt_Sel");
1484 histTitle->Append("_Har");
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 (%)");
1491 // neutral sidebands Flow observed - Pt,Eta
1492 histTitle = new TString("FlowV0sb_vObs2D_Sel");
1494 histTitle->Append("_Har");
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)");
1501 // neutral sidebands Flow observed - Eta
1502 histTitle = new TString("FlowV0sb_vObsEta_Sel");
1504 histTitle->Append("_Har");
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 (%)");
1511 // neutral sidebands Flow observed - Pt
1512 histTitle = new TString("FlowV0sb_vObsPt_Sel");
1514 histTitle->Append("_Har");
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 (%)");
1521 // SX neutral sidebands Flow observed - Eta
1522 histTitle = new TString("FlowV0sb_vObsEta_sx_Sel");
1524 histTitle->Append("_Har");
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 (%)");
1531 // SX neutral sidebands Flow observed - Pt
1532 histTitle = new TString("FlowV0sb_vObsPt_sx_Sel");
1534 histTitle->Append("_Har");
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 (%)");
1541 // DX neutral sidebands Flow observed - Eta
1542 histTitle = new TString("FlowV0sb_vObsEta_dx_Sel");
1544 histTitle->Append("_Har");
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 (%)");
1551 // DX neutral sidebands Flow observed - Pt
1552 histTitle = new TString("FlowV0sb_vObsPt_dx_Sel");
1554 histTitle->Append("_Har");
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 (%)");
1563 histTitle = new TString("Flow_Phi_TPCplus_Sel");
1565 histTitle->Append("_Har");
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");
1573 histTitle = new TString("Flow_Phi_TPCminus_Sel");
1575 histTitle->Append("_Har");
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");
1583 histTitle = new TString("Flow_Phi_TPCcross_Sel");
1585 histTitle->Append("_Har");
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");
1593 histTitle = new TString("Flow_Phi_TPC_Sel");
1595 histTitle->Append("_Har");
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");
1602 // Phi lab flattened
1604 histTitle = new TString("Flow_Phi_Flat_TPCplus_Sel");
1606 histTitle->Append("_Har");
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");
1614 histTitle = new TString("Flow_Phi_Flat_TPCminus_Sel");
1616 histTitle->Append("_Har");
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");
1624 histTitle = new TString("Flow_Phi_Flat_TPCcross_Sel");
1626 histTitle->Append("_Har");
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");
1634 histTitle = new TString("Flow_Phi_Flat_TPC_Sel");
1636 histTitle->Append("_Har");
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");
1644 cout << " Init() - Histograms booked" << endl ;
1648 //-----------------------------------------------------------------------
1650 //-----------------------------------------------------------------------
1651 Bool_t AliFlowAnalyser::Finish()
1653 // Close the analysis and saves the histograms on the histFile .
1655 cout << "* FlowAnalysis * - Finish()" << endl ; cout << endl ;
1657 // Write all histograms
1659 fHistFile->Write() ;
1661 // Write Resolution corrected histograms
1664 fVnResHistList->Write();
1665 delete fVnResHistList ;
1667 else { cout << " E.P. resolution has not been calculated. No v_n histograms!" << endl ; }
1669 // Write PhiWgt histograms
1672 fPhiWgtHistList->Write();
1673 delete fPhiWgtHistList ;
1676 fFlowSelect->Write();
1677 // delete fFlowSelect ;
1679 fHistFile->Close() ;
1681 cout << " Finish() - Histograms saved : " << fHistFileName.Data() << endl ; cout << endl ;
1685 //-----------------------------------------------------------------------
1687 //----------------------------------------------------------------------
1688 Float_t AliFlowAnalyser::GetRunBayesian(Int_t nPid, Int_t selN)
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.
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. ; }
1699 //-----------------------------------------------------------------------
1700 void AliFlowAnalyser::PrintRunBayesian(Int_t selN)
1702 // Prints the normalized particle abundance of all the analysed events
1703 // (in selection selN).
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++)
1711 bayes = GetRunBayesian(i, selN) ;
1712 cout << bayes << "_" << names[i] << " ; " ;
1718 //-----------------------------------------------------------------------
1719 void AliFlowAnalyser::FillWgtArrays(TFile* wgtFile)
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).
1725 fPhiWgtFile = wgtFile ;
1727 TString* histTitle ;
1728 TH1D* TPC_all ; TH1D* TPC_plus ; TH1D* TPC_minus ; TH1D* TPC_cross ;
1730 for(int k=0;k<Flow::nSels;k++)
1732 for(int j=0;j<Flow::nHars;j++)
1735 histTitle = new TString("Flow_Phi_Weight_TPCplus_Sel");
1737 histTitle->Append("_Har");
1739 TPC_plus = (TH1D*)fPhiWgtFile->Get(histTitle->Data());
1742 histTitle = new TString("Flow_Phi_Weight_TPCminus_Sel");
1744 histTitle->Append("_Har");
1746 TPC_minus = (TH1D*)fPhiWgtFile->Get(histTitle->Data());
1749 histTitle = new TString("Flow_Phi_Weight_TPCcross_Sel");
1751 histTitle->Append("_Har");
1753 TPC_cross = (TH1D*)fPhiWgtFile->Get(histTitle->Data());
1757 histTitle = new TString("Flow_Phi_Weight_TPC_Sel");
1759 histTitle->Append("_Har");
1761 TPC_all = (TH1D*)fPhiWgtFile->Get(histTitle->Data());
1764 for(int n=0;n<fPhiBins;n++)
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 ;
1775 histTitle = new TString("Flow_BayPidMult_Sel");
1777 PID_bay = (TH1D*)fPhiWgtFile->Get(histTitle->Data());
1779 Double_t totCount = PID_bay->GetSumOfWeights() ;
1780 for (int n=0;n<Flow::nPid;n++)
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 ;
1788 delete TPC_all ; delete TPC_plus ; delete TPC_minus ; delete TPC_cross ;
1793 //-----------------------------------------------------------------------
1794 void AliFlowAnalyser::FillEvtPhiWgt(AliFlowEvent* fFlowEvent)
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(...).
1800 fFlowEvent->SetPhiWeight(fPhiWgt);
1801 fFlowEvent->SetPhiWeightPlus(fPhiWgtPlus);
1802 fFlowEvent->SetPhiWeightMinus(fPhiWgtMinus);
1803 fFlowEvent->SetPhiWeightCross(fPhiWgtCross);
1807 //-----------------------------------------------------------------------
1808 void AliFlowAnalyser::FillBayesianWgt(AliFlowEvent* fFlowEvent)
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!).
1816 Double_t bayes[Flow::nPid] ;
1817 Double_t bayCheck = 0. ;
1818 for (int n=0;n<Flow::nPid;n++)
1820 bayes[n] = fBayesianWgt[0][n] ;
1821 bayCheck += bayes[n] ;
1822 // cout << "Bayesian V[" << n << "] = " << fBayesianWgt[0][n] << endl ;
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 ; }
1829 //-----------------------------------------------------------------------
1830 void AliFlowAnalyser::Weightening()
1832 // Calculates weights, and fills PhiWgt histograms .
1833 // This is called at the end of the event analysis.
1835 cout << " AliFlowAnalyser::Weightening() " << endl ; cout << endl ;
1837 // PhiWgt histogram collection
1838 fPhiWgtHistList = new TOrdCollection(4*Flow::nSels*Flow::nHars) ;
1840 // Creates PhiWgt Histograms
1841 TString* histTitle ;
1842 for(int k = 0; k < Flow::nSels; k++)
1844 for(int j = 0; j < Flow::nHars; j++)
1847 histTitle = new TString("Flow_Phi_Weight_TPCplus_Sel");
1849 histTitle->Append("_Har");
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");
1857 histTitle = new TString("Flow_Phi_Weight_TPCminus_Sel");
1859 histTitle->Append("_Har");
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");
1867 histTitle = new TString("Flow_Phi_Weight_TPCcross_Sel");
1869 histTitle->Append("_Har");
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");
1877 histTitle = new TString("Flow_Phi_Weight_TPC_Sel");
1879 histTitle->Append("_Har");
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");
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 ;
1894 for (int i=0;i<fPhiBins;i++)
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.);
1906 if(meanTPC==0) { cout << " Sel." << k << " , Har." << j << " : empty phi histogram ! " << endl ; }
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);
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);
1924 //-----------------------------------------------------------------------
1926 //-----------------------------------------------------------------------
1927 Bool_t AliFlowAnalyser::Analyse(AliFlowEvent* flowEvent)
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.
1934 cout << " AliFlowAnalyser::Analyze(" << fFlowEvent << " ) - " << fEventNumber << endl ;
1935 if(!flowEvent) { return kFALSE ; }
1936 else { fFlowEvent = flowEvent ; }
1938 if(fFlowSelect->Select(fFlowEvent)) // event selected - here below the ANALYSIS FLAGS are setted -
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 ;
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
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
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
1958 if(fShuffle) { fFlowEvent->RandomShuffle() ; } // tracks re-shuffling
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
1964 cout << " * 2 . Calculating event quantities all in one shoot . " << endl ;
1965 fFlowEvent->MakeAll() ;
1967 if(FillFromFlowEvent(fFlowEvent)) // calculates event quantities
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)
1977 cout << " * 3 . Event psi = 0 - Skipping! " << endl ;
1983 cout << " * 0 . Event " << fEventNumber << " (event ID = " << fFlowEvent->EventID() << ") discarded . " << endl ;
1984 delete fFlowEvent ; fFlowEvent = 0 ;
1991 //-----------------------------------------------------------------------
1992 Bool_t AliFlowAnalyser::FillFromFlowEvent(AliFlowEvent* fFlowEvent)
1994 // gets event quantities, returns kFALSE if Q vector is always 0
1996 Int_t selCheck = 0 ;
1997 for(Int_t k = 0; k < Flow::nSels; k++)
1999 fFlowSelect->SetSelection(k) ;
2000 for(Int_t j = 0; j < Flow::nHars; j++)
2002 fFlowSelect->SetHarmonic(j) ;
2003 for(Int_t n = 0; n < Flow::nSubs; n++)
2005 fFlowSelect->SetSubevent(n) ;
2006 fPsiSub[n][k][j] = fFlowEvent->Psi(fFlowSelect) ; // sub-event quantities
2007 fMultSub[n][k][j] = fFlowEvent->Mult(fFlowSelect) ;
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] ;
2017 if(!selCheck) { return kFALSE ; } // if there are no particles in the selection -> skip the event
2021 //-----------------------------------------------------------------------
2022 void AliFlowAnalyser::FillEventHistograms(AliFlowEvent* fFlowEvent)
2024 // fill event histograms
2026 cout << " Fill Event Histograms ... " << endl ;
2028 float trigger = (float)fFlowEvent->L0TriggerWord() ;
2029 fHistTrigger->Fill(trigger);
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());
2036 int cent = fFlowEvent->Centrality();
2037 fHistCent->Fill((float)cent);
2039 fHistMult->Fill((float)fNumberOfTracks) ;
2040 fHistV0Mult->Fill((float)fNumberOfV0s) ;
2041 if(origMult) { fHistMultOverOrig->Fill((float)fNumberOfTracks/(float)origMult) ; }
2043 fFlowEvent->VertexPos(fVertex) ;
2044 fHistVertexZ->Fill(fVertex[2]) ;
2045 fHistVertexXY2D->Fill(fVertex[0],fVertex[1]) ;
2048 fHistPartZDC->Fill(fFlowEvent->ZDCpart()) ;
2049 for(int ii=0;ii<3;ii++) { fHistEnergyZDC->Fill(ii,fFlowEvent->ZDCenergy(ii)) ; }
2051 // sub-event Psi_Subs
2052 for(int k = 0; k < Flow::nSels; k++)
2054 for(int j = 0; j < Flow::nHars; j++)
2056 for(int n = 0; n < Flow::nSubs; n++)
2058 int iii = Flow::nSubs * k + n ; //cout << " " << k << j << n << " , " << iii << endl ;
2059 fHistSub[iii].fHistSubHar[j].fHistPsiSubs->Fill(fPsiSub[n][k][j]) ;
2064 // full event Psi, PsiSubCorr, PsiSubCorrDiff, cos, mult, q
2065 for(int k = 0; k < Flow::nSels; k++)
2067 for(int j = 0; j < Flow::nHars; j++)
2069 float order = (float)(j+1);
2070 fHistFull[k].fHistFullHar[j].fHistPsi->Fill(fPsi[k][j]);
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
2084 float psi1 ; float psi2 ;
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
2093 if(fPsiSub[0][k][j] != 0. && fPsiSub[1][k][j] != 0.)
2095 float psiSubCorr; // this is: delta_Psi
2096 if(fV1Ep1Ep2 == kFALSE || order != 1)
2098 psiSubCorr = fPsiSub[0][k][j] - fPsiSub[1][k][j];
2100 else // i.e. (fV1Ep1Ep2 == kTRUE && order == 1)
2102 psiSubCorr = fPsiSub[0][k][0] + fPsiSub[1][k][0] - 2*fPsi[k][1];
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);
2110 if(j < Flow::nHars - 1) // subevents of different harmonics
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) ;
2125 fHistFull[k].fHistFullHar[j].fHistMult->Fill((float)fMult[k][j]) ;
2126 fHistFull[k].fHistFullHar[j].fHistQnorm->Fill(fQnorm[k][j]) ;
2132 //-----------------------------------------------------------------------
2133 void AliFlowAnalyser::FillParticleHistograms(TObjArray* fFlowTracks)
2135 // fills tracks histograms
2137 cout << " Tracks Loop . " << endl ;
2139 float corrMultUnit = 0. ;
2140 float corrMultN = 0. ;
2141 float etaSymPosTpcN = 0. ;
2142 float etaSymNegTpcN = 0. ;
2143 float etaSymPosTpcNpart = 0. ;
2144 float etaSymNegTpcNpart = 0. ;
2146 float hMinusN = 0. ;
2147 float piPlusN = 0. ;
2148 float piMinusN = 0. ;
2149 float protonN = 0. ;
2151 float kMinusN = 0. ;
2153 float deuteronN = 0. ;
2155 float electronN = 0. ;
2156 float positronN = 0. ;
2157 float muonMinusN = 0. ;
2158 float muonPlusN = 0. ;
2160 for(fTrackNumber=0;fTrackNumber<fNumberOfTracks;fTrackNumber++)
2162 fFlowTrack = (AliFlowTrack*)fFlowTracks->At(fTrackNumber) ;
2163 //cout << "Track n. " << fTrackNumber << endl ; fFlowTrack->Dump() ;
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
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);
2214 // - here ITS (chi2 & nHits)
2215 fHistChi2ITS->Fill(chi2ITS);
2217 fHistChi2normITS->Fill(chi2ITS/((float)fitPtsITS));
2218 fHistFitPtsITS->Fill((float)fitPtsITS);
2219 fHistMaxPtsITS->Fill((float)maxPtsITS);
2221 // - here TPC (chi2 & nHits)
2222 fHistChi2TPC->Fill(chi2TPC);
2224 fHistChi2normTPC->Fill(chi2TPC/((float)fitPtsTPC));
2225 fHistFitPtsTPC->Fill((float)fitPtsTPC);
2226 fHistMaxPtsTPC->Fill((float)maxPtsTPC);
2228 fHistFitOverMaxTPC->Fill((float)(fitPtsTPC)/(float)maxPtsTPC);
2230 // - here TRD (chi2 & nHits)
2231 fHistChi2TRD->Fill(chi2TRD);
2233 fHistChi2normTRD->Fill(chi2TRD/((float)fitPtsTRD));
2234 fHistFitPtsTRD->Fill((float)fitPtsTRD);
2235 fHistMaxPtsTRD->Fill((float)maxPtsTRD);
2237 // - here TOF (chi2 & nHits)
2238 fHistChi2TOF->Fill(chi2TOF);
2240 fHistChi2normTOF->Fill(chi2TOF/((float)fitPtsTOF));
2241 fHistFitPtsTOF->Fill((float)fitPtsTOF);
2242 fHistMaxPtsTOF->Fill((float)maxPtsTOF);
2244 // fit over max (all)
2245 int maxPts = maxPtsITS + maxPtsTPC + maxPtsTRD + maxPtsTOF ;
2246 int fitPts = fitPtsITS + fitPtsTPC + fitPtsTRD + fitPtsTOF ;
2248 fHistFitOverMax->Fill((float)(fitPts)/(float)maxPts) ;
2251 fHistLenght->Fill(lenght) ;
2252 fHistInvMass->Fill(invMass) ;
2254 // PID histograms & multiplicity count (for bayesian histogram)
2258 fHistMeanDedxPos2D->Fill(lpTPC, dEdx) ;
2259 fHistMeanDedxPos2DITS->Fill(lpITS, its) ;
2260 fHistMeanDedxPos2DTRD->Fill(lpTRD, trd) ;
2261 fHistMeanDedxPos2DTOF->Fill(lpTOF, tof) ;
2263 float positron = fFlowTrack->ElectronPositronProb() ;
2264 fHistPidPositron->Fill(positron) ;
2265 if(strcmp(pid, "e+") == 0)
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) ;
2274 float muonPlus = fFlowTrack->MuonPlusMinusProb() ;
2275 fHistPidMuonPlus->Fill(muonPlus) ;
2276 if(strcmp(pid, "mu+") == 0)
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) ;
2285 float piPlus = fFlowTrack->PionPlusMinusProb() ;
2286 fHistPidPiPlus->Fill(piPlus) ;
2287 if(strcmp(pid, "pi+") == 0)
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) ;
2296 float kplus = fFlowTrack->KaonPlusMinusProb() ;
2297 fHistPidKplus->Fill(kplus) ;
2298 if(strcmp(pid, "k+") == 0)
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) ;
2307 float proton = fFlowTrack->ProtonPbarProb() ;
2308 fHistPidProton->Fill(proton) ;
2309 if(strcmp(pid, "pr+") == 0)
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) ;
2318 float deuteron = fFlowTrack->DeuteriumAntiDeuteriumProb() ;
2319 fHistPidDeuteron->Fill(deuteron) ;
2320 if(strcmp(pid, "d+") == 0)
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) ;
2330 else if(charge == -1)
2333 fHistMeanDedxNeg2D->Fill(lpTPC, dEdx) ;
2334 fHistMeanDedxNeg2DITS->Fill(lpITS, its) ;
2335 fHistMeanDedxNeg2DTRD->Fill(lpTRD, trd) ;
2336 fHistMeanDedxNeg2DTOF->Fill(lpTOF, tof) ;
2338 float electron = fFlowTrack->ElectronPositronProb() ;
2339 fHistPidElectron->Fill(electron);
2340 if(strcmp(pid, "e-") == 0)
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);
2349 float muonMinus = fFlowTrack->MuonPlusMinusProb() ;
2350 fHistPidMuonMinus->Fill(muonMinus) ;
2351 if(strcmp(pid, "mu-") == 0)
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) ;
2360 float piMinus = fFlowTrack->PionPlusMinusProb() ;
2361 fHistPidPiMinus->Fill(piMinus) ;
2362 if(strcmp(pid, "pi-") == 0)
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);
2371 float kminus = fFlowTrack->KaonPlusMinusProb() ;
2372 fHistPidKminus->Fill(kminus);
2373 if(strcmp(pid, "k-") == 0)
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);
2382 float antiproton = fFlowTrack->ProtonPbarProb() ;
2383 fHistPidAntiProton->Fill(antiproton);
2384 if(strcmp(pid, "pr-") == 0)
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);
2393 float antideuteron = fFlowTrack->DeuteriumAntiDeuteriumProb() ;
2394 fHistPidAntiDeuteron->Fill(antideuteron);
2395 if(strcmp(pid, "d-") == 0)
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);
2406 // Yield3D, Yield2D, Eta, Pt, Phi, bayP.Id.
2407 fHistPtot->Fill(totalp) ;
2409 fHistPhi->Fill(phi);
2410 fHistAllEtaPtPhi3D->Fill(eta, pt, phi) ;
2411 fHistYieldAll2D->Fill(eta, pt) ;
2412 fHistBayPidMult->Fill(fPidId) ;
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) ;
2423 fHistYieldUnc2D->Fill(etaGlob, ptGlob) ;
2424 fHistUncEtaPtPhi3D->Fill(etaGlob, ptGlob, phiGlob) ;
2425 fHistPhiPtUnc->Fill(phiGlob, ptGlob) ;
2427 if(fFlowTrack->Charge()>0) { fHistPtPhiPos->Fill(phi, pt); }
2428 else if(fFlowTrack->Charge()<0) { fHistPtPhiNeg->Fill(phi, pt); }
2430 // fills selected part histograms
2431 if(fFlowSelect->SelectPart(fFlowTrack))
2433 if(strlen(fFlowSelect->PidPart()) != 0)
2435 float rapidity = fFlowTrack->Y();
2436 fHistBinEta->Fill(rapidity, rapidity);
2437 fHistYieldPart2D->Fill(rapidity, pt);
2441 fHistBinEta->Fill(eta, eta) ;
2442 fHistYieldPart2D->Fill(eta, pt) ;
2444 fHistBayPidMultPart->Fill(fPidId) ;
2445 fHistBinPt->Fill(pt, pt) ;
2446 fHistEtaPtPhi3DPart->Fill(eta,pt,phi) ;
2447 fHistDcaGlobalPart->Fill(dcaGlobal) ;
2448 fHistInvMassPart->Fill(invMass) ;
2451 fHistMeanDedxPos3DPart->Fill(lpITS, dEdx, fPidId) ;
2452 fHistMeanDedxPos3DPartITS->Fill(lpITS, its, fPidId) ;
2454 else if(charge == -1)
2456 fHistMeanDedxNeg3DPart->Fill(lpITS, dEdx, fPidId) ;
2457 fHistMeanDedxNeg3DPartITS->Fill(lpITS, its, fPidId) ;
2459 if(eta > 0.) { etaSymPosTpcNpart++ ; } // for fHistEtaSymPart
2460 else { etaSymNegTpcNpart++ ; }
2464 fHistEtaPtPhi3DOut->Fill(eta,pt,phi) ;
2465 fHistYieldOut2D->Fill(eta, pt) ;
2466 fHistDcaGlobalOut->Fill(dcaGlobal) ;
2467 fHistInvMassOut->Fill(invMass) ;
2471 for(int j = 0; j < Flow::nHars; j++)
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);
2480 //For Eta symmetry TPC
2481 if(fFlowTrack->FitPtsTPC())
2483 if(eta > 0.) { etaSymPosTpcN++ ; } // for fHistEtaSym
2484 else { etaSymNegTpcN++ ; }
2487 // not to call it twice ...
2488 int nEtaS = HarmonicsLoop(fFlowTrack) ;
2490 //For Correlated Multiplicity
2491 corrMultN += (float)nEtaS ;
2493 //For Correlated Multiplicity in 1 unit rapidity
2494 if(TMath::Abs(eta) <= 0.5) { corrMultUnit += (float)nEtaS ; }
2498 } // end of tracks loop
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] ;
2507 fHistEtaSym->Fill(etaSymTpc);
2508 fHistEtaSymVerZ2D->Fill(vertexZ , etaSymTpc);
2510 fHistEtaSymPart->Fill(etaSymTpc);
2511 fHistEtaSymVerZ2DPart->Fill(vertexZ , etaSymTpcPart);
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);
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) ;
2540 //-----------------------------------------------------------------------
2541 Int_t AliFlowAnalyser::HarmonicsLoop(AliFlowTrack* fFlowTrack)
2543 // HarmonicsLoop, does the correlation analysis, fills histograms
2544 // (internally called by ::FillParticleHistograms(..) loop) .
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() ;
2552 // Looping over Selections and Harmonics
2554 for (int k = 0; k < Flow::nSels; k++)
2556 fFlowSelect->SetSelection(k) ;
2557 for (int j = 0; j < Flow::nHars; j++)
2559 bool oddHar = (j+1) % 2;
2560 fFlowSelect->SetHarmonic(j);
2561 double order = (double)(j+1);
2563 if(fFlowEvent->EtaSubs()) // particles with the opposite subevent
2565 if(eta > 0) { psi_i = fPsiSub[1][k][j] ; } //check
2566 else { psi_i = fPsiSub[0][k][j] ; }
2568 else if(order > 3. && !oddHar)
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 ; }
2574 else // random subevents
2576 psi_i = fPsi[k][j] ;
2579 if(fFlowSelect->Select(fFlowTrack)) // Get detID
2581 Bool_t kTpcPlus = kFALSE ;
2582 Bool_t kTpcMinus = kFALSE ;
2583 Bool_t kTpcAll = kFALSE ;
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()));
2590 // Set Tpc (+ and -)
2591 if(fFlowTrack->FitPtsTPC()) //OR*: AliESDtrack:: "TBits& GetTPCClusterMap()" or "Int_t GetTPCclusters(Int_t* idx)" ...
2593 if(zFirstPoint >= 0. && eta > 0.) { kTpcPlus = kTRUE ; }
2594 else if(zFirstPoint <= 0. && eta < 0.) { kTpcMinus = kTRUE ; }
2595 else { kTpcAll = kTRUE ; }
2598 // PID Multiplicities (particle for R.P.) - done just one time for each selection
2599 if(j==0) { fHistFull[k].fHistBayPidMult->Fill(fPidId) ; }
2601 // Calculate weights for filling histograms
2602 float wt = 1. ; // TMath::Abs(fFlowEvent->Weight(k, j, fFlowTrack)) ;
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) ;
2610 // Get phiWgt from file
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
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) ;
2623 if(oddHar && eta<0.) { phiWgt *= -1. ; } // restore value
2625 // Remove autocorrelations
2627 if(!fFlowEvent->EtaSubs()) // random subevents
2629 if(order > 3. && !oddHar) // 2nd harmonic event plane
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() ; }
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 ; }
2645 // Remove autocorrelations of the second order 'particles' which are used for v1{EP1,EP2}.
2646 if (fV1Ep1Ep2 == kTRUE && order == 1)
2648 AliFlowSelection usedForPsi2 = *fFlowSelect ;
2649 usedForPsi2.SetHarmonic(1);
2650 if(usedForPsi2.Select(fFlowTrack)) // particle was used for Psi2
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() ; }
2657 else // particle was not used for Psi2
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()));
2671 // Caculate v for all particles selected for correlation analysis
2672 if(fFlowSelect->SelectPart(fFlowTrack))
2677 if (fV1Ep1Ep2 == kFALSE || order != 1)
2679 v = 100 * cos(order * (phi - psi_i)) ;
2681 else // i.e. (fV1Ep1Ep2 == kTRUE && order == 1)
2683 v = 100 * cos(phi + psi_i - 2*psi_2) ;
2687 if(eta < 0 && oddHar) { vFlip *= -1 ; }
2688 if(strlen(fFlowSelect->PidPart()) != 0) // pid, fill rapidity
2690 float rapidity = fFlowTrack->Y();
2691 fHistFull[k].fHistFullHar[j].fHistvObs2D->Fill(rapidity, pt, v);
2693 if(fPtRangevEta[1] > fPtRangevEta[0]) // cut is used
2695 if(pt < fPtRangevEta[1] && pt >= fPtRangevEta[0]) // check cut range, fill if in range
2697 fHistFull[k].fHistFullHar[j].fHistvObsEta->Fill(rapidity, v);
2700 else // cut is not used, fill in any case
2702 fHistFull[k].fHistFullHar[j].fHistvObsEta->Fill(rapidity, v);
2705 else // no pid, fill eta
2707 fHistFull[k].fHistFullHar[j].fHistvObs2D->Fill(eta, pt, v);
2709 if(fPtRangevEta[1] > fPtRangevEta[0]) // cut is used
2711 if(pt < fPtRangevEta[1] && pt >= fPtRangevEta[0]) // check cut range, fill if in range
2713 fHistFull[k].fHistFullHar[j].fHistvObsEta->Fill(eta, v);
2716 else // cut is not used, fill in any case
2718 fHistFull[k].fHistFullHar[j].fHistvObsEta->Fill(eta, v);
2722 if(fEtaRangevPt[1] > fEtaRangevPt[0]) // cut is used
2724 if(TMath::Abs(eta) < fEtaRangevPt[1] && TMath::Abs(eta) >= fEtaRangevPt[0]) // check cut range, fill if in range
2726 fHistFull[k].fHistFullHar[j].fHistvObsPt->Fill(pt, vFlip); // for odd harmonis /-/
2729 else // cut is not used, fill in any case
2731 fHistFull[k].fHistFullHar[j].fHistvObsPt->Fill(pt, vFlip);
2735 Bool_t etaPtNoCut = kTRUE;
2736 if(fPtRangevEta[1] > fPtRangevEta[0] && (pt < fPtRangevEta[0] || pt >= fPtRangevEta[1]))
2738 etaPtNoCut = kFALSE;
2740 if(fEtaRangevPt[1] > fEtaRangevPt[0] && (TMath::Abs(eta) < fEtaRangevPt[0] || TMath::Abs(eta) >= fEtaRangevPt[1]))
2742 etaPtNoCut = kFALSE;
2744 if(etaPtNoCut) { fHistFull[k].fHistvObs->Fill(order, vFlip) ; }
2746 // Correlation of Phi of selected particles with Psi
2748 if(eta < 0 && oddHar)
2750 phi_i += TMath::Pi() ; // backward particle and odd harmonic
2751 if(phi_i > 2*TMath::Pi()) { phi_i -= 2*TMath::Pi() ; }
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));
2762 //-----------------------------------------------------------------------
2763 void AliFlowAnalyser::FillV0Histograms(TObjArray* fFlowV0s)
2767 cout << " V0s Loop . " << endl ;
2769 int corrMultV0 = 0 ;
2770 for(fV0Number=0;fV0Number<fNumberOfV0s;fV0Number++)
2772 fFlowV0 = (AliFlowV0*)fFlowV0s->At(fV0Number) ;
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() ; }
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) ;
2803 if(fFlowSelect->SelectPart(fFlowV0))
2805 bool inWin = fFlowSelect->SelectV0Part(fFlowV0) ;
2806 bool sx = fFlowSelect->SelectV0sxSide(fFlowV0) ;
2807 bool dx = fFlowSelect->SelectV0dxSide(fFlowV0) ;
2810 fHistV0YieldPart2D->Fill(eta, pt) ;
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);
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);
2830 for(int k = 0; k < Flow::nSels; k++) // sort of HarmonicsLoop - selection number used
2832 fFlowSelect->SetSelection(k) ;
2833 for(Int_t j=0;j<Flow::nHars;j++)
2835 Bool_t oddHar = (j+1) % 2 ;
2836 Float_t order = (Float_t)(j+1) ;
2837 fFlowSelect->SetHarmonic(j);
2839 // Remove autocorrelations
2840 Float_t psi_i, psi_2 ;
2842 Float_t phiDaughter1 = 0. ; Float_t phiDaughter2 = 0. ;
2843 Double_t phiWgt1 = 0. ; Double_t phiWgt2 = 0. ;
2847 fFlowTrack = daughterPlus ;
2848 phiDaughter1 = fFlowTrack->Phi() ;
2849 if(fFlowSelect->Select(fFlowTrack)) // Get phiWgt from file
2851 if(order > 3. && !oddHar) { phiWgt1 = fFlowEvent->PhiWeight(k, 1, fFlowTrack) ; }
2852 else { phiWgt1 = fFlowEvent->PhiWeight(k, j, fFlowTrack) ; }
2857 fFlowTrack = daughterMinus ;
2858 phiDaughter2 = fFlowTrack->Phi() ;
2859 if(fFlowSelect->Select(fFlowTrack)) // Get phiWgt from file
2861 if(order > 3. && !oddHar) { phiWgt2 = fFlowEvent->PhiWeight(k, 1, fFlowTrack) ; }
2862 else { phiWgt2 = fFlowEvent->PhiWeight(k, j, fFlowTrack) ; }
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() ; }
2875 if(order > 3. && !oddHar) { psi_i = psi_2 ; }
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 ; }
2886 // Caculate v for all V0s selected for correlation analysis
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 ; }
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) ; }
2896 if(fPtRangevEta[1] > fPtRangevEta[0]) // cut is used
2898 if(pt < fPtRangevEta[1] && pt >= fPtRangevEta[0]) // check cut range, fill if in range
2900 if(inWin) { fHistFull[k].fHistFullHar[j].fHistV0vObsEta->Fill(eta, v) ; }
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) ; }
2909 else // cut is not used, fill in any case
2911 if(inWin) { fHistFull[k].fHistFullHar[j].fHistV0vObsEta->Fill(eta, v) ; }
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) ; }
2919 if(fEtaRangevPt[1] > fEtaRangevPt[0]) // cut is used
2921 if(TMath::Abs(eta) < fEtaRangevPt[1] && TMath::Abs(eta) >= fEtaRangevPt[0]) // check cut range, fill if in range
2923 if(inWin) { fHistFull[k].fHistFullHar[j].fHistV0vObsPt->Fill(pt, vFlip) ; } // for odd harmonis /-/
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) ; }
2932 else // cut is not used, fill in any case
2934 if(inWin) { fHistFull[k].fHistFullHar[j].fHistV0vObsPt->Fill(pt, vFlip) ; }
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) ; }
2943 Bool_t etaPtNoCut = kTRUE;
2944 if(fPtRangevEta[1] > fPtRangevEta[0] && (pt < fPtRangevEta[0] || pt >= fPtRangevEta[1]))
2946 etaPtNoCut = kFALSE;
2948 if(fEtaRangevPt[1] > fEtaRangevPt[0] && (TMath::Abs(eta) < fEtaRangevPt[0] || TMath::Abs(eta) >= fEtaRangevPt[1]))
2950 etaPtNoCut = kFALSE;
2954 if(inWin) { fHistFull[k].fHistV0vObs->Fill(order, vFlip) ; }
2957 if(sx) { fHistFull[k].fHistV0sbvObsSx->Fill(order, vFlip) ; }
2958 else if(dx) { fHistFull[k].fHistV0sbvObsDx->Fill(order, vFlip) ; }
2962 // Correlation of Phi of selected v0s with Psi
2964 if(eta < 0 && oddHar)
2966 phi_i += TMath::Pi() ; // backward particle and odd harmonic
2967 if(phi_i > 2*TMath::Pi()) { phi_i -= 2*TMath::Pi() ; }
2969 float dPhi = phi_i - psi_i;
2970 if(dPhi < 0.) { dPhi += 2*TMath::Pi() ; }
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)) ; }
2977 //delete fFlowV0 ; fFlowV0 = 0 ;
2978 //delete daughterPlus ; daughterPlus = 0 ;
2979 //delete daughterMinus ; daughterMinus = 0 ;
2981 fHistV0MultPart->Fill(corrMultV0) ;
2985 //-----------------------------------------------------------------------
2986 // void AliFlowAnalyser::FillLabels()
2988 // // Reads the tracks label (i.e. the link from ESD tracking to KineTree)
2991 // cout << " Fill Labels . " << endl ;
2995 //-----------------------------------------------------------------------
2996 void AliFlowAnalyser::PrintEventQuantities()
2998 // prints event by event calculated quantities
3000 cout << endl ; cout << " # " << fEventNumber << " - Event quantities : " << endl ; cout << endl ;
3002 cout << "fQ[k][j] = " ;
3003 for(int k=0;k<Flow::nSels;k++)
3005 for(int j=0;j<Flow::nHars;j++)
3007 cout << (Float_t)fQ[k][j].X() << "," << (Float_t)fQ[k][j].Y() << " ; " ;
3012 cout << endl ; cout << "fPsi[k][j] = " ;
3013 for(int k=0;k<Flow::nSels;k++)
3015 for(int j=0;j<Flow::nHars;j++) { Float_t aaa = (Float_t)fPsi[k][j] ; cout << aaa << " , " ; }
3019 cout << endl ; cout << "fQnorm[k][j] = " ;
3020 for(int k=0;k<Flow::nSels;k++)
3022 for(int j=0;j<Flow::nHars;j++) { Float_t aaa = (Float_t)fQnorm[k][j] ; cout << aaa << " , " ; }
3026 cout << endl ; cout << "fMult[k][j] = " ;
3027 for(int k=0;k<Flow::nSels;k++)
3029 for(int j=0;j<Flow::nHars;j++) { Float_t aaa = (Float_t)fMult[k][j] ; cout << aaa << " , " ; }
3033 cout << endl ; cout << "fMultSub[n][k][j] = " ;
3034 for(int n=0;n<Flow::nSubs;n++)
3036 for(int k=0;k<Flow::nSels;k++)
3038 for(int j=0;j<Flow::nHars;j++) { Float_t aaa = fMultSub[n][k][j] ; cout << aaa << " , " ; }
3043 cout << endl ; cout << "fPsiSub[n][k][j] = " ;
3044 for(int n=0;n<Flow::nSubs;n++)
3046 for(int k=0;k<Flow::nSels;k++)
3048 for(int j=0;j<Flow::nHars;j++) { Float_t aaa = fPsiSub[n][k][j] ; cout << aaa << " , " ; }
3053 cout << endl ; cout << "Delta_PsiSub[k][j] = " ;
3054 for(int k=0;k<Flow::nSels;k++)
3056 for(int j=0;j<Flow::nHars;j++) { Float_t aaa = fPsiSub[0][k][j]-fPsiSub[1][k][j] ; cout << aaa << " , " ; }
3062 //----------------------------------------------------------------------
3064 //----------------------------------------------------------------------
3065 Bool_t AliFlowAnalyser::Resolution()
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.
3071 cout << " AliFlowAnalyser::Resolution()" << endl ;
3073 // VnRes histogram collection
3074 fVnResHistList = new TOrdCollection(Flow::nSels*Flow::nHars);
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;
3084 for(int k = 0; k < Flow::nSels; k++)
3086 // v for tracks (corrected for resolution)
3087 histTitle = new TString("Flow_v_Sel");
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 (%)");
3094 fVnResHistList->AddLast(fHistFull[k].fHistv);
3096 // v for v0s (corrected for resolution)
3097 histTitle = new TString("FlowV0_v_Sel");
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 (%)");
3104 fVnResHistList->AddLast(fHistFull[k].fHistV0v);
3106 for(int j = 0; j < Flow::nHars; j++)
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.)
3113 double resSub = 0. ;
3114 double resSubErr = 0. ;
3115 if(fV1Ep1Ep2 == kTRUE && order == 1) // calculate resolution of second order event plane first
3118 double res2error = 0. ;
3119 if(fHistFull[k].fHistCos->GetBinContent(2) > 0.)
3121 if(fEtaSub) // sub res only
3123 res2 = TMath::Sqrt(fHistFull[k].fHistCos->GetBinContent(2));
3124 res2error = fHistFull[k].fHistCos->GetBinError(2) / (2. * res2);
3128 if (fHistFull[k].fHistCos->GetBinContent(2) > 0.92) // resolution saturates
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;
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
3155 else if(fEtaSub) // sub res only
3157 resSub = TMath::Sqrt(cosPair[k][j]);
3158 resSubErr = cosPairErr[k][j] / (2. * resSub);
3159 fRes[k][j] = resSub;
3160 fResErr[k][j] = resSubErr;
3162 else if(order==4. || order==6.|| order==8.) // 2nd harmonic event plane
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);
3172 fRes[k][j] = ResEventPlaneK2(TMath::Sqrt(2.) * chiSub); // full event plane res.
3173 mResDelta = ResEventPlaneK2(TMath::Sqrt(2.) * chiSubDelta);
3177 fRes[k][j] = ResEventPlaneK3(TMath::Sqrt(2.) * chiSub); // full event plane res.
3178 mResDelta = ResEventPlaneK3(TMath::Sqrt(2.) * chiSubDelta);
3182 fRes[k][j] = ResEventPlaneK4(TMath::Sqrt(2.) * chiSub); // full event plane res.
3183 mResDelta = ResEventPlaneK4(TMath::Sqrt(2.) * chiSubDelta);
3185 fResErr[k][j] = resSubErr * fabs((double)fRes[k][j] - mResDelta) / deltaResSub;
3189 if(cosPair[k][j] > 0.92) // resolution saturates
3192 fResErr[k][j] = 0.007;
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;
3207 else // subevent correlation must be positive
3212 fHistFull[k].fHistRes->SetBinContent(j+1, fRes[k][j]);
3213 fHistFull[k].fHistRes->SetBinError(j+1, fResErr[k][j]);
3215 // Create the v 2D histogram (Flow corrected for res.) - v,Pt,Eta
3216 histTitle = new TString("Flow_v2D_Sel");
3218 histTitle->Append("_Har");
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 (%)");
3226 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistv2D);
3228 // Create the 1D v histograms - v,Eta
3229 histTitle = new TString("Flow_vEta_Sel");
3231 histTitle->Append("_Har");
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 (%)");
3238 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistvEta);
3241 histTitle = new TString("Flow_vPt_Sel");
3243 histTitle->Append("_Har");
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 (%)");
3250 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistvPt);
3252 // Create the v 2D histogram (V0s Flow corrected for res.) - v,Pt,Eta
3253 histTitle = new TString("FlowV0_v2D_Sel");
3255 histTitle->Append("_Har");
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 (%)");
3263 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0v2D);
3265 // Create the 1D v histograms (V0s) - v,Eta
3266 histTitle = new TString("FlowV0_vEta_Sel");
3268 histTitle->Append("_Har");
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 (%)");
3275 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0vEta);
3278 histTitle = new TString("FlowV0_vPt_Sel");
3280 histTitle->Append("_Har");
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 (%)");
3287 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0vPt);
3289 // Create the v 2D histogram (V0s sidebands Flow corrected for res.) - v,Pt,Eta
3290 histTitle = new TString("FlowV0sb_v2D_Sel");
3292 histTitle->Append("_Har");
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 (%)");
3300 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0sbv2D);
3302 // Create the 1D v histograms (V0s sidebands) - v,Eta
3303 histTitle = new TString("FlowV0sb_vEta_Sel");
3305 histTitle->Append("_Har");
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 (%)");
3312 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0sbvEta);
3314 // (V0s sidebands) v,Pt
3315 histTitle = new TString("FlowV0sb_vPt_Sel");
3317 histTitle->Append("_Har");
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 (%)");
3324 fVnResHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistV0sbvPt);
3326 // Calulate v = vObs / Resolution
3329 cout << "***** Resolution of the " << j+1 << "th harmonic = " << fRes[k][j] << " +/- " << fResErr[k][j] << endl;
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);
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);
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.
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;
3356 errorV0 = fHistFull[k].fHistV0v->GetBinError(j+1) ;
3357 errorV0 /= fRes[k][j];
3359 totalError = fabs(contentV0) * TMath::Sqrt((errorV0/contentV0)*(errorV0/contentV0) + (fResErr[k][j]/fRes[k][j])*(fResErr[k][j]/fRes[k][j]));
3362 fHistFull[k].fHistV0v->SetBinError(j+1, totalError);
3366 cout << "##### Resolution of the " << j+1 << "th harmonic was zero." << endl;
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.);
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.);
3388 //-----------------------------------------------------------------------
3389 Double_t AliFlowAnalyser::Chi(Double_t res)
3391 // chi from the event plane resolution
3394 Double_t delta = 1.0;
3395 for(int i = 0; i < 15; i++)
3397 if(ResEventPlane(chi) < res) { chi = chi + delta ; }
3398 else { chi = chi - delta ; }
3404 //-----------------------------------------------------------------------
3405 Double_t AliFlowAnalyser::ResEventPlane(Double_t chi)
3407 // plane resolution as function of chi
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));
3415 //-----------------------------------------------------------------------
3416 Double_t AliFlowAnalyser::ResEventPlaneK2(Double_t chi)
3418 // ... for the case k=2.
3420 Double_t con = TMath::Sqrt(TMath::Pi()/2)/2 ; // ~ 0.626657
3421 Double_t arg = chi * chi / 4.;
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);
3429 //-----------------------------------------------------------------------
3430 Double_t AliFlowAnalyser::ResEventPlaneK3(Double_t chi)
3432 // ...for the case k=3.
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));
3440 //-----------------------------------------------------------------------
3441 Double_t AliFlowAnalyser::ResEventPlaneK4(Double_t chi)
3443 // ...for the case k=4.
3445 Double_t con = TMath::Sqrt(TMath::Pi()/2)/2 ; // ~ 0.626657
3446 Double_t arg = chi * chi / 4.;
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);
3455 //-------------------------------------------------------------
3457 //-----------------------------------------------------------------------